/*
 * $Id: mod_pmin.cpp 494 2012-06-13 20:41:16Z johnoel $
 *
 * Author: David Fournier
 * Copyright (c) 2008-2012 Regents of the University of California 
 */
#include <admodel.h>

void get_confidence_interval(const dvector& left_bd, const dvector& right_bd,
  dmatrix& ms, const dvector& xs, const dvector& siglevel, const int& level_index,
  dvector& xdist,int index);
void get_onesided_intervals(const dvector& left_bd, const dvector& right_bd,
  dmatrix& ms, const dvector& xs, const dvector& siglevel, const int& level_index,
  dvector& xdist,int index);
void report_confidence_limits(const ofstream& ofs3,int numsig_levels,
  dvector& siglevel, const dvector& left_bd, const dvector& right_bd);
void report_onesided_confidence_limits(const ofstream& ofs3,int numsig_levels,
  dvector& siglevel, const dvector& left_bd, const dvector& right_bd,int ip);

void get_ee(const dmatrix& hh, const ofstream& of5); 

ofstream of5("eigv.rpt");

dmatrix trans(const dvector& x)
{
  int mmin=x.indexmin();
  int mmax=x.indexmax();
  dmatrix tmp(mmin,mmax,1,1);
  for (int i=mmin;i<=mmax;i++)
  {
    tmp(i,1)=x(i);
  }
  return tmp;
}
  double mult_factor(int j)
  {
    switch(j)
    {
      case 0:
        return .0;
      case 1:
        return .5;
      case 2:
        return 1.0;
      case 3:
        return 1.5;
      default:
        return 2.0;
    }
  }

  double trimax(double x,double y,double z);

#if defined(USE_ADPVM)
  void function_minimizer::pvm_slave_likeprof_routine(void)
  {
    do
    {
      int prof_switch=get_int_from_master();
      if (!prof_switch) break;
      if (prof_switch !=3)
      {
        cerr << "Error in prof_switch " << prof_switch << endl;
        ad_exit(1);
      }
      int underflow_flag=get_int_from_master();
      pvm_slave_prof_minimize(underflow_flag);
    }
    while(1);
  }
#endif
  void function_minimizer::likeprof_routine(double global_min)
  {
   // dvector siglevel(sshhiitt);
    dvector siglevel("{.90,.95,.975}");
    int num_pp=likeprof_params::likeprofptr[0]->get_stepnumber();
    {
      for (int ip=1;ip<likeprof_params::num_likeprof_params;ip++)
      {
        int sno=likeprof_params::likeprofptr[ip]->get_stepnumber();
        if (sno)
        {
	  if (sno>num_pp) num_pp=sno;
        }
      }
    }
    double relsig=.5;
    initial_params::current_phase = initial_params::max_number_phases;
    // DF NOV 28 11
    if (random_effects_flag)
    {
     initial_params::set_inactive_only_random_effects();
    }
    int nvar=initial_params::nvarcalc();
    dvector xsave(1,nvar);
    int ii=1;
    initial_params::copy_all_values(xsave,ii);
    double old_value;  // this is where we were
    double new_value;  // this is where we want to go
    double fprof;//
    double current_value;
    dmatrix lprof(0,likeprof_params::num_likeprof_params-1,-num_pp,num_pp);
    dmatrix ldet(0,likeprof_params::num_likeprof_params-1,-num_pp,num_pp);
    dmatrix ln_det_proj_jac(0,likeprof_params::num_likeprof_params-1,
      -num_pp,num_pp);
    dmatrix pdet(0,likeprof_params::num_likeprof_params-1,-num_pp,num_pp);
    dmatrix actual_value(0,likeprof_params::num_likeprof_params-1,
      -num_pp,num_pp);
    dvector all_values(-num_pp,num_pp);
    dvector all_num_sigs(-num_pp,num_pp);
    dvector xxxtmp(-num_pp,num_pp);
    //d3_array hesses(-num_pp,num_pp,1,nvar,1,nvar);
    dmatrix lg_jacob(0,likeprof_params::num_likeprof_params-1,-num_pp,num_pp);
    dmatrix lg_prjacob(0,likeprof_params::num_likeprof_params-1,-num_pp,num_pp);
    dmatrix xxtmp(-num_pp,num_pp,1,nvar);
    dvector xvector(1,nvar);
    dmatrix xmax(-num_pp,num_pp,1,nvar);  // this holds the conditional max
    dmatrix gprof(-num_pp,num_pp,1,nvar);  // this holds the conditional max
    dmatrix fgrads(-num_pp,num_pp,1,nvar);  // this holds the conditional max
    dmatrix xdist(0,likeprof_params::num_likeprof_params-1,-num_pp,num_pp);
    int sign=0;
    double sigma;
    ofstream ofs2((char*) (ad_comm::adprogram_name + adstring(".prf")) );
    double udet =unrestricted_hess_determinant();
    const int offset=0;
    dvector xscale(1,nvar);   // need to get scale from somewhere
    dvector likeprof_save(0,likeprof_params::num_likeprof_params-1);
    dmatrix penalties(0,likeprof_params::num_likeprof_params-1,-num_pp,num_pp);
    penalties.initialize();
#ifdef CURVE_CORRECT
    int nlp=likeprof_params::num_likeprof_params;
    d3_array eigenvals(0,nlp-1,-num_pp,num_pp,1,nvar-1);
    d3_array curvcor(0,nlp-1,-num_pp,num_pp,1,nvar-1);
#endif

    int ip;
    for (ip=0;ip<likeprof_params::num_likeprof_params;ip++)
    {
      likeprof_save(ip)=likeprof_params::likeprofptr[ip]->get_value();
    }
    double final_weight;
    for (ip=0;ip<likeprof_params::num_likeprof_params;ip++)
    {
      int sno=likeprof_params::likeprofptr[ip]->get_stepnumber();
      double snz=likeprof_params::likeprofptr[ip]->get_stepsize();
      if (sno)
      {
	num_pp=sno;
      }
      if (snz)
      {
	relsig=snz;
      }
      if (ip>0)
      {
	int ii=1;
	initial_params::restore_all_values(xsave,ii);
      }
      sigma=likeprof_params::likeprofptr[ip]->get_sigma(); // this is the
      if (sigma==0.0)
        cerr << "error standard dev of likeporf parameter is 0" << endl;
					// estimated sd
      old_value=likeprof_save(ip);
      old_value=old_value+offset*relsig*sigma;  // this is where we
      int bigint_flag=0;
      int bigint_flag1=0;
      int lastj=0;
      int lastj1=0;
      double ldiff=0.0;
      double num_sigs;
      for (int i=1;i<=2;i++)  // go in positive and negative directions
      {
        num_sigs=0.0;
	bigint_flag=0;
	lastj=0;
	bigint_flag1=0;
	lastj1=0;
	int underflow_flag=0;
	if (i>1) // get the parameter values at the global minimum
	{
	 int ii=1;
	 initial_params::restore_all_values(xsave,ii);
	}
	if (i==1)
	{
	  sign=1;
	}
	else
	{
	  sign=-1;
	}
	current_value=old_value;  // initialize at the minimum
	for (int j=0;j<=num_pp;j++)  // go in positive and negative directions
	{
	 if (j!=0 || sign > 0)
	 {
	  if (bigint_flag==0)
	  {
	    num_sigs+=mult_factor(j)*relsig*sign;
	    current_value+=mult_factor(j)*relsig*sign*sigma;
	    new_value=current_value;
	    // new_value=current_value+j*relsig*sign*sigma;
	  }
	  else
	  {
	    if (bigint_flag1==0)
	    {
	      num_sigs+=1.5*relsig*sign;
	      current_value+=1.5*relsig*sign*sigma;
	      new_value=current_value;
	    }
	    else
	    {
	      num_sigs+=2.5*relsig*sign;
	      current_value+=2.5*relsig*sign*sigma;
	      new_value=current_value;
	    }
	  }
          if (!ad_comm::pvm_manager)
          {
            if (random_effects_flag==0)
            {
	      prof_minimize(ip,sigma,new_value,fprof,underflow_flag,
	        global_min,penalties(ip,j*sign),final_weight); // get the 
                                                        // conditional max
            }
            else
            {
	      prof_minimize_re(ip,sigma,new_value,fprof,underflow_flag,
	        global_min,penalties(ip,j*sign),final_weight); // get the 
                                                        // conditional max
            }
          }
          else
          {
#if defined(USE_ADPVM)
            if (ad_comm::pvm_manager->mode==1) // master
            {
              send_int_to_slaves(3);
	      pvm_master_prof_minimize(ip,sigma,new_value,fprof,underflow_flag,
	        global_min,penalties(ip,j*sign),final_weight); // get the 
            }
#else
            {
              cerr << "PVM not included with this distribution" << endl;
              ad_exit(1);
            }
#endif
          }
          all_num_sigs(j*sign)=num_sigs;
	  initial_params::xinit(xvector); // save the
	  int ic=1;
	  initial_params::copy_all_values(xmax(sign*j),ic); // save the
							// conditional maximum
	  /*int check=*/initial_params::stddev_scale(xscale,xvector);
        //#if defined(DO_PROFILE)
	  //dvector curvscale(1,nvar);   // need to get scale from somewhere
        //#endif

       // #if defined(DO_PROFILE)
	 // check=initial_params::stddev_curvscale(curvscale,xvector);
       // #endif
	  //cout << "xscale = " << endl << xscale << endl;
          //{
           // ofstream ofs("xscale");
           // ofs << xscale << endl;
         // }
#if defined(USE_DDOUBLE)
	  lg_jacob(ip,sign*j)=sum(log(fabs(xscale)+double(1.e-60)));
#else
	  lg_jacob(ip,sign*j)=sum(log(fabs(xscale)+1.e-60));
#endif
	  if (!underflow_flag)
	  {
	    lprof(ip,sign*j)=fprof;
	  }
	  else
	  {
	    lprof(ip,sign*j)=lprof(ip,sign*(j-1))+2;
	  }
	  double xx=likeprof_params::likeprofptr[ip]->get_value();
	  if (!underflow_flag)
	  {
	    actual_value(ip,sign*j)=xx;
	  }
	  else
	  {
	    actual_value(ip,sign*j)=new_value;  // this is where we
	  }
	
	  ldiff=fprof-lprof(ip,0);
	  if ( ldiff > 40.0)
	  {
	    underflow_flag=1;
	  }

	  //if ( (fprof-lprof(ip,0)) > 3.0 && bigint_flag ==0)
	  if ( abs(j) >4 )
	  {
	    bigint_flag=1;
	    lastj=5;
	  }
	  if ( abs(j) >5 )
	  {
	    bigint_flag1=1;
	    lastj1=6;
	  }
	  // get the gradient for the profile likelihood variable at
	  // the conditional maximum
	  initial_params::current_phase = initial_params::max_number_phases;
	  int nvar=initial_params::nvarcalc();
	  dvector g(1,nvar);
	  dvector fg(1,nvar);
	  if (!underflow_flag)
	  {
	    // g is the grad. wrt the prof lik var
	    // fg is the grad. wrt the log-likelihood
	    get_particular_grad(ip,nvar,fg,g);
	    gprof(sign*j)=g;

	    fgrads(sign*j)=fg;
            xdist(ip,sign*j)=norm(elem_div(gprof(sign*j),xscale));
	  }
	  else
	  {
	    gprof(sign*j)=gprof(sign*(j-1));
            xdist(ip,sign*j)=xdist(ip,sign*(j-1));
	  }
         
        //#if defined(DO_PROFILE)
	/*
	  if (!underflow_flag)
	  {
            hess_routine_and_constraint(ip,g,fg);  // calculate the hessian 
                                            // at the conditional max
          }

	  if (!underflow_flag)
	  {
            ldet(ip,sign*j)=projected_hess_determinant(g,underflow_flag,
              xscale,ln_det_proj_jac(ip,sign*j));
          }
          else
	  {
            ldet(ip,sign*j)=ldet(ip,sign*(j-1));
            ln_det_proj_jac(ip,sign*j)=ln_det_proj_jac(ip,sign*(j-1));
          }
	 */
        //#endif
	 }   // end of the j loop
	}
      }
      //for (int ix=-num_pp;ix<=num_pp;ix++)
      //{
       // xdist(ip,ix)=norm(elem_div(gprof(ix),xscale));
     // }
      if ( (option_match(ad_comm::argc,ad_comm::argv,"-prsave"))>-1)
      {
        adstring profrep_name=(likeprof_params::likeprofptr[ip]->label());
        ofstream ofs3((char*) (profrep_name+adstring(".pvl")));
        for (int ix=-num_pp;ix<=num_pp;ix++)
	{
	  ofs3 << "#Step " << ix << endl;
	  ofs3 << "#num sigmas " << all_num_sigs(ix) << endl;
	  ofs3 << xmax(ix) << endl;
	}
      }
    }
  #if defined(DO_PROFILE)
    for (ip=0;ip<likeprof_params::num_likeprof_params;ip++)
    {
      ldet(ip)=ldet(ip)-ln_det_proj_jac(ip);
    }
    {
    
      ofstream ofs("det.tmp");
      for (ip=0;ip<likeprof_params::num_likeprof_params;ip++)
      {
        ofs << "the log dets" << endl;
        ofs << "ldet" << endl;
        ofs << ldet(ip) << endl << endl;
        ofs << "lndet_proj_jac" << endl;
        ofs << ln_det_proj_jac(ip) << endl << endl;
        ofs << "ldet-lndet_proj_jac" << endl;
        ofs << ldet(ip)-ln_det_proj_jac(ip) +ln_det_proj_jac(ip,0) 
            << endl << endl;
      }
    }
    //for (ip=0;ip<likeprof_params::num_likeprof_params;ip++)
    //{
     // ldet(ip)=ldet(ip)-(ln_det_proj_jac(ip)-ln_det_proj_jac(ip));
    //}
  #endif
      sigma=likeprof_params::likeprofptr[0]->get_sigma(); // this is the
      if (sigma==0.0)
        cerr << "Error standard dev of likeprof parameter is 0" << endl;
#ifndef CURVE_CORRECT
    normalize_posterior_distribution(udet,siglevel,ofs2,num_pp,
      all_values,actual_value,global_min,offset,lprof,ldet,xdist,
      penalties);
#else
    normalize_posterior_distribution(udet,siglevel,ofs2,num_pp,
      all_values,actual_value,global_min,offset,lprof,ldet,xdist,
      eigenvals,curvcor);
#endif
  }

void get_ee(const dmatrix& hh, const ofstream& _of5) 
{
  ofstream& of5= (ofstream&) _of5; 
  int mmin=hh.rowmin();
  int mmax=hh.rowmax();
  dvector l(mmin,mmax);
  dmatrix tmp(mmin,mmax,1,2);
  dvector ll(mmin,mmax);
  dmatrix e=eigenvectors(hh,l);
  for (int i=mmin;i<=mmax;i++)
  {
    ll(i)=e(i)*(hh*e(i)); 
    of5 << l(i) << "  " << ll(i) << endl;
  }
}
