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