ADMB Documentation  11.1.2503
 All Classes Files Functions Variables Typedefs Friends Defines
df1b2ghmult.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: df1b2ghmult.cpp 1935 2014-04-26 02:02:58Z 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   int i,j;
00046 
00047   // contribution for quadratic prior
00048   if (quadratic_prior::get_num_quadratic_prior()>0)
00049   {
00050     //Hess+=quadratic_prior::get_cHessian_contribution();
00051     int & vxs = (int&)(xs);
00052     quadratic_prior::get_cHessian_contribution(Hess,vxs);
00053   }
00054  // Here need hooks for sparse matrix structures
00055 
00056   dvar3_array & block_diagonal_vhessian=
00057     *pmin->lapprox->block_diagonal_vhessian;
00058   block_diagonal_vhessian.initialize();
00059   dvar3_array& block_diagonal_ch=
00060     *pmin->lapprox->block_diagonal_vch;
00061     //dvar3_array(*pmin->lapprox->block_diagonal_ch);
00062   int ii=xs+us+1;
00063   d3_array& bdH=(*pmin->lapprox->block_diagonal_hessian);
00064   int ic;
00065   for (ic=1;ic<=nsc;ic++)
00066   {
00067     int lus=lrea(ic);
00068     for (i=1;i<=lus;i++)
00069       for (j=1;j<=lus;j++)
00070         y(ii++)=bdH(ic)(i,j);
00071   }
00072 
00073   dvector g(1,nvar);
00074   gradcalc(0,g);
00075   gradient_structure::set_YES_DERIVATIVES();
00076   dvar_vector vy=dvar_vector(y);
00077   //initial_params::stddev_vscale(d,vy);
00078   ii=xs+us+1;
00079   if (initial_df1b2params::have_bounded_random_effects)
00080   {
00081     cerr << "can't do importance sampling with bounded random effects"
00082      " at present" << endl;
00083     ad_exit(1);
00084   }
00085   else
00086   {
00087     for (int ic=1;ic<=nsc;ic++)
00088     {
00089       int lus=lrea(ic);
00090       if (lus>0)
00091       {
00092         for (i=1;i<=lus;i++)
00093         {
00094           for (j=1;j<=lus;j++)
00095           {
00096             block_diagonal_vhessian(ic,i,j)=vy(ii++);
00097           }
00098         }
00099         block_diagonal_ch(ic)=
00100           choleski_decomp(inv(block_diagonal_vhessian(ic)));
00101       }
00102     }
00103   }
00104 
00105    int nsamp=pmin->lapprox->use_gauss_hermite;
00106    pmin->lapprox->in_gauss_hermite_phase=1;
00107    dvar_vector sample_value(1,nsamp);
00108    sample_value.initialize();
00109 
00110    dvar_vector tau(1,us);;
00111    // !!! This only works for one random efect in each separable call
00112    // at present.
00113 
00114    if (pmin->lapprox->gh->mi)
00115    {
00116      delete pmin->lapprox->gh->mi;
00117      pmin->lapprox->gh->mi=0;
00118    }
00119 
00120    pmin->lapprox->gh->mi=new multi_index(1,nsamp,
00121     pmin->lapprox->multi_random_effects);
00122 
00123    multi_index & mi = *(pmin->lapprox->gh->mi);
00124 
00125    //for (int is=1;is<=nsamp;is++)
00126    dvector& xx=pmin->lapprox->gh->x;
00127    do
00128    {
00129      int offset=0;
00130      pmin->lapprox->num_separable_calls=0;
00131      //pmin->lapprox->gh->is=is;
00132      for (ic=1;ic<=nsc;ic++)
00133      {
00134        int lus=lrea(ic);
00135        // will need vector stuff here when more than one random effect
00136        if (lus>0)
00137        {
00138          //tau(offset+1,offset+lus).shift(1)=block_diagonal_ch(ic)(1,1)*
00139          //  pmin->lapprox->gh->x(is);
00140          dvector xv(1,lus);
00141          for (int iu=1;iu<=lus;iu++)
00142          {
00143            xv(iu)= xx(mi()(iu));
00144          }
00145          tau(offset+1,offset+lus).shift(1)=block_diagonal_ch(ic)*xv;
00146 
00147          offset+=lus;
00148        }
00149      }
00150 
00151      // have to reorder the terms to match the block diagonal hessian
00152      imatrix & ls=*(pmin->lapprox->block_diagonal_re_list);
00153      int mmin=ls.indexmin();
00154      int mmax=ls.indexmax();
00155 
00156      int ii=1;
00157      int i;
00158      for (i=mmin;i<=mmax;i++)
00159      {
00160        int cmin=ls(i).indexmin();
00161        int cmax=ls(i).indexmax();
00162        for (int j=cmin;j<=cmax;j++)
00163        {
00164          vy(ls(i,j))+=tau(ii++);
00165        }
00166      }
00167      if (ii-1 != us)
00168      {
00169        cerr << "error in interface" << endl;
00170        ad_exit(1);
00171      }
00172      initial_params::reset(vy);    // get the values into the model
00173      ii=1;
00174      for (i=mmin;i<=mmax;i++)
00175      {
00176        int cmin=ls(i).indexmin();
00177        int cmax=ls(i).indexmax();
00178        for (int j=cmin;j<=cmax;j++)
00179        {
00180          vy(ls(i,j))-=tau(ii++);
00181        }
00182      }
00183 
00184      *objective_function_value::pobjfun=0.0;
00185      pmin->AD_uf_outer();
00186      ++mi;
00187    }
00188    while(mi.get_depth()<=pmin->lapprox->multi_random_effects);
00189 
00190    nsc=pmin->lapprox->num_separable_calls;
00191 
00192    dvariable vf=pmin->do_gauss_hermite_integration();
00193 
00194    int sgn=0;
00195    dvariable ld=0.0;
00196    if (ad_comm::no_ln_det_choleski_flag)
00197    {
00198      for (int ic=1;ic<=nsc;ic++)
00199      {
00200        if (allocated(block_diagonal_vhessian(ic)))
00201        {
00202          ld+=ln_det(block_diagonal_vhessian(ic),sgn);
00203        }
00204      }
00205      ld*=0.5;
00206    }
00207    else
00208    {
00209      for (int ic=1;ic<=nsc;ic++)
00210      {
00211        if (allocated(block_diagonal_vhessian(ic)))
00212        {
00213          ld+=ln_det_choleski(block_diagonal_vhessian(ic));
00214        }
00215      }
00216      ld*=0.5;
00217    }
00218 
00219    vf+=ld;
00220    //vf+=us*0.91893853320467241;
00221 
00222    double f=value(vf);
00223    gradcalc(nvar,g);
00224 
00225    // put uhat back into the model
00226    gradient_structure::set_NO_DERIVATIVES();
00227    vy(xs+1,xs+us).shift(1)=u0;
00228    initial_params::reset(vy);    // get the values into the model
00229    gradient_structure::set_YES_DERIVATIVES();
00230 
00231    pmin->lapprox->in_gauss_hermite_phase=0;
00232 
00233   ii=1;
00234   for (i=1;i<=xs;i++)
00235     xadjoint(i)=g(ii++);
00236   for (i=1;i<=us;i++)
00237     uadjoint(i)=g(ii++);
00238   for (ic=1;ic<=nsc;ic++)
00239   {
00240     int lus=lrea(ic);
00241     for (i=1;i<=lus;i++)
00242     {
00243       for (j=1;j<=lus;j++)
00244       {
00245         (*pmin->lapprox->block_diagonal_vhessianadjoint)(ic)(i,j)=g(ii++);
00246       }
00247     }
00248   }
00249   return f;
00250 }