ADMB Documentation  11.1x.2711
 All Classes Files Functions Variables Typedefs Friends Defines
df1b2ghmult.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: df1b2ghmult.cpp 2684 2014-11-18 19:50:21Z johnoel $
00003  *
00004  * Author: David Fournier
00005  * Copyright (c) 2008-2012 Regents of the University of California
00006  */
00011 #  include <admodel.h>
00012 #  include <df1b2fun.h>
00013 #  include <adrndeff.h>
00014 
00019 double do_gauss_hermite_block_diagonal_multi(const dvector& x,
00020   const dvector& u0,const dmatrix& Hess,const dvector& _xadjoint,
00021   const dvector& _uadjoint,const dmatrix& _Hessadjoint,
00022   function_minimizer * pmin)
00023 {
00024   ADUNCONST(dvector,xadjoint)
00025   ADUNCONST(dvector,uadjoint)
00026   //ADUNCONST(dmatrix,Hessadjoint)
00027   const int xs=x.size();
00028   const int us=u0.size();
00029   gradient_structure::set_NO_DERIVATIVES();
00030   int nsc=pmin->lapprox->num_separable_calls;
00031   const ivector lrea = (*pmin->lapprox->num_local_re_array)(1,nsc);
00032   int hroom =  sum(square(lrea));
00033   int nvar=x.size()+u0.size()+hroom;
00034   independent_variables y(1,nvar);
00035 
00036   // need to set random effects active together with whatever
00037   // init parameters should be active in this phase
00038   initial_params::set_inactive_only_random_effects();
00039   initial_params::set_active_random_effects();
00040   /*int onvar=*/initial_params::nvarcalc();
00041   initial_params::xinit(y);    // get the initial values into the
00042   // do we need this next line?
00043   y(1,xs)=x;
00044 
00045   // contribution for quadratic prior
00046   if (quadratic_prior::get_num_quadratic_prior()>0)
00047   {
00048     //Hess+=quadratic_prior::get_cHessian_contribution();
00049     int & vxs = (int&)(xs);
00050     quadratic_prior::get_cHessian_contribution(Hess,vxs);
00051   }
00052  // Here need hooks for sparse matrix structures
00053 
00054   dvar3_array & block_diagonal_vhessian=
00055     *pmin->lapprox->block_diagonal_vhessian;
00056   block_diagonal_vhessian.initialize();
00057   dvar3_array& block_diagonal_ch=
00058     *pmin->lapprox->block_diagonal_vch;
00059     //dvar3_array(*pmin->lapprox->block_diagonal_ch);
00060   int ii=xs+us+1;
00061   d3_array& bdH=(*pmin->lapprox->block_diagonal_hessian);
00062   for (int ic=1;ic<=nsc;ic++)
00063   {
00064     int lus=lrea(ic);
00065     for (int i=1;i<=lus;i++)
00066       for (int j=1;j<=lus;j++)
00067         y(ii++)=bdH(ic)(i,j);
00068   }
00069 
00070   dvector g(1,nvar);
00071   gradcalc(0,g);
00072   gradient_structure::set_YES_DERIVATIVES();
00073   dvar_vector vy=dvar_vector(y);
00074   //initial_params::stddev_vscale(d,vy);
00075   ii=xs+us+1;
00076   if (initial_df1b2params::have_bounded_random_effects)
00077   {
00078     cerr << "can't do importance sampling with bounded random effects"
00079      " at present" << endl;
00080     ad_exit(1);
00081   }
00082   else
00083   {
00084     for (int ic=1;ic<=nsc;ic++)
00085     {
00086       int lus=lrea(ic);
00087       if (lus>0)
00088       {
00089         for (int i=1;i<=lus;i++)
00090         {
00091           for (int j=1;j<=lus;j++)
00092           {
00093             block_diagonal_vhessian(ic,i,j)=vy(ii++);
00094           }
00095         }
00096         block_diagonal_ch(ic)=
00097           choleski_decomp(inv(block_diagonal_vhessian(ic)));
00098       }
00099     }
00100   }
00101 
00102    int nsamp=pmin->lapprox->use_gauss_hermite;
00103    pmin->lapprox->in_gauss_hermite_phase=1;
00104    dvar_vector sample_value(1,nsamp);
00105    sample_value.initialize();
00106 
00107    dvar_vector tau(1,us);;
00108    // !!! This only works for one random efect in each separable call
00109    // at present.
00110 
00111    if (pmin->lapprox->gh->mi)
00112    {
00113      delete pmin->lapprox->gh->mi;
00114      pmin->lapprox->gh->mi=0;
00115    }
00116 
00117    pmin->lapprox->gh->mi=new multi_index(1,nsamp,
00118     pmin->lapprox->multi_random_effects);
00119 
00120    multi_index & mi = *(pmin->lapprox->gh->mi);
00121 
00122    //for (int is=1;is<=nsamp;is++)
00123    dvector& xx=pmin->lapprox->gh->x;
00124    do
00125    {
00126      int offset=0;
00127      pmin->lapprox->num_separable_calls=0;
00128      //pmin->lapprox->gh->is=is;
00129      for (int ic=1;ic<=nsc;ic++)
00130      {
00131        int lus=lrea(ic);
00132        // will need vector stuff here when more than one random effect
00133        if (lus>0)
00134        {
00135          //tau(offset+1,offset+lus).shift(1)=block_diagonal_ch(ic)(1,1)*
00136          //  pmin->lapprox->gh->x(is);
00137          dvector xv(1,lus);
00138          for (int iu=1;iu<=lus;iu++)
00139          {
00140            xv(iu)= xx(mi()(iu));
00141          }
00142          tau(offset+1,offset+lus).shift(1)=block_diagonal_ch(ic)*xv;
00143 
00144          offset+=lus;
00145        }
00146      }
00147 
00148      // have to reorder the terms to match the block diagonal hessian
00149      imatrix & ls=*(pmin->lapprox->block_diagonal_re_list);
00150      int mmin=ls.indexmin();
00151      int mmax=ls.indexmax();
00152 
00153      ii=1;
00154      for (int i=mmin;i<=mmax;i++)
00155      {
00156        int cmin=ls(i).indexmin();
00157        int cmax=ls(i).indexmax();
00158        for (int j=cmin;j<=cmax;j++)
00159        {
00160          vy(ls(i,j))+=tau(ii++);
00161        }
00162      }
00163      if (ii-1 != us)
00164      {
00165        cerr << "error in interface" << endl;
00166        ad_exit(1);
00167      }
00168      initial_params::reset(vy);    // get the values into the model
00169      ii=1;
00170      for (int i=mmin;i<=mmax;i++)
00171      {
00172        int cmin=ls(i).indexmin();
00173        int cmax=ls(i).indexmax();
00174        for (int j=cmin;j<=cmax;j++)
00175        {
00176          vy(ls(i,j))-=tau(ii++);
00177        }
00178      }
00179 
00180      *objective_function_value::pobjfun=0.0;
00181      pmin->AD_uf_outer();
00182      ++mi;
00183    }
00184    while(mi.get_depth()<=pmin->lapprox->multi_random_effects);
00185 
00186    nsc=pmin->lapprox->num_separable_calls;
00187 
00188    dvariable vf=pmin->do_gauss_hermite_integration();
00189 
00190    int sgn=0;
00191    dvariable ld=0.0;
00192    if (ad_comm::no_ln_det_choleski_flag)
00193    {
00194      for (int ic=1;ic<=nsc;ic++)
00195      {
00196        if (allocated(block_diagonal_vhessian(ic)))
00197        {
00198          ld+=ln_det(block_diagonal_vhessian(ic),sgn);
00199        }
00200      }
00201      ld*=0.5;
00202    }
00203    else
00204    {
00205      for (int ic=1;ic<=nsc;ic++)
00206      {
00207        if (allocated(block_diagonal_vhessian(ic)))
00208        {
00209          ld+=ln_det_choleski(block_diagonal_vhessian(ic));
00210        }
00211      }
00212      ld*=0.5;
00213    }
00214 
00215    vf+=ld;
00216    //vf+=us*0.91893853320467241;
00217 
00218    double f=value(vf);
00219    gradcalc(nvar,g);
00220 
00221    // put uhat back into the model
00222    gradient_structure::set_NO_DERIVATIVES();
00223    vy(xs+1,xs+us).shift(1)=u0;
00224    initial_params::reset(vy);    // get the values into the model
00225    gradient_structure::set_YES_DERIVATIVES();
00226 
00227    pmin->lapprox->in_gauss_hermite_phase=0;
00228 
00229   ii=1;
00230   for (int i=1;i<=xs;i++)
00231     xadjoint(i)=g(ii++);
00232   for (int i=1;i<=us;i++)
00233     uadjoint(i)=g(ii++);
00234   for (int ic=1;ic<=nsc;ic++)
00235   {
00236     int lus=lrea(ic);
00237     for (int i=1;i<=lus;i++)
00238     {
00239       for (int j=1;j<=lus;j++)
00240       {
00241         (*pmin->lapprox->block_diagonal_vhessianadjoint)(ic)(i,j)=g(ii++);
00242       }
00243     }
00244   }
00245   return f;
00246 }