ADMB Documentation  11.1.2394
 All Classes Files Functions Variables Typedefs Friends Defines
df1b2im5.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: df1b2im5.cpp 2294 2014-09-04 23:29:42Z 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 calculate_importance_sample_block_diagonal_option_antithetical
00020   (const dvector& x,const dvector& u0,const dmatrix& Hess,
00021   const dvector& _xadjoint,const dvector& _uadjoint,
00022   const dmatrix& _Hessadjoint,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     int ic;
00088     for (ic=1;ic<=nsc;ic++)
00089     {
00090       int lus=lrea(ic);
00091       for (i=1;i<=lus;i++)
00092       {
00093         for (j=1;j<=lus;j++)
00094         {
00095           block_diagonal_vhessian(ic,i,j)=vy(ii++);
00096         }
00097       }
00098       block_diagonal_ch(ic)=
00099         choleski_decomp(inv(block_diagonal_vhessian(ic)));
00100     }
00101   }
00102 
00103    int nsamp=pmin->lapprox->num_importance_samples;
00104 
00105    dvariable vf=0.0;
00106 
00107    dvar_vector sample_value(1,nsamp);
00108    sample_value.initialize();
00109 
00110    dvar_vector tau(1,us);;
00111    int is;
00112    for (is=1;is<=nsamp;is++)
00113    {
00114      int offset=0;
00115      pmin->lapprox->importance_sampling_counter=is;
00116      int ic;
00117      for (ic=1;ic<=nsc;ic++)
00118      {
00119        int lus=lrea(ic);
00120        tau(offset+1,offset+lus).shift(1)=block_diagonal_ch(ic)*
00121          (*pmin->lapprox->antiepsilon)(is);
00122        offset+=lus;
00123      }
00124 
00125      // have to reorder the terms to match the block diagonal hessian
00126      imatrix & ls=*(pmin->lapprox->block_diagonal_re_list);
00127      int mmin=ls.indexmin();
00128      int mmax=ls.indexmax();
00129 
00130      int ii=1;
00131      int i;
00132      for (i=mmin;i<=mmax;i++)
00133      {
00134        int cmin=ls(i).indexmin();
00135        int cmax=ls(i).indexmax();
00136        for (int j=cmin;j<=cmax;j++)
00137        {
00138          vy(ls(i,j))+=tau(ii++);
00139        }
00140      }
00141      if (ii-1 != us)
00142      {
00143        cerr << "error in interface" << endl;
00144        ad_exit(1);
00145      }
00146      initial_params::reset(vy);    // get the values into the model
00147      ii=1;
00148      for (i=mmin;i<=mmax;i++)
00149      {
00150        int cmin=ls(i).indexmin();
00151        int cmax=ls(i).indexmax();
00152        for (int j=cmin;j<=cmax;j++)
00153        {
00154          vy(ls(i,j))-=tau(ii++);
00155        }
00156      }
00157 
00158      *objective_function_value::pobjfun=0.0;
00159      pmin->AD_uf_outer();
00160 
00161      if (pmin->lapprox->use_outliers==0)
00162      {
00163        // assumes that all separable calls have the same number
00164        // of random effects
00165        double neps=0.5*nsc*norm2((*pmin->lapprox->antiepsilon)(is));
00166 
00167        (*pmin->lapprox->importance_sampling_values)(is)=
00168            neps-value(*objective_function_value::pobjfun);
00169 
00170        (*pmin->lapprox->importance_sampling_weights)(is)=neps;
00171 
00172         sample_value(is)=*objective_function_value::pobjfun
00173           -neps;
00174      }
00175      else
00176      {
00177        dvector& e=pmin->lapprox->epsilon(is);
00178        double neps=-sum(log(.95*exp(-0.5*square(e))+
00179           0.0166666667*exp(-square(e)/18.0)));
00180 
00181        (*pmin->lapprox->importance_sampling_values)(is)=
00182          neps-value(*objective_function_value::pobjfun);
00183 
00184        sample_value(is)=*objective_function_value::pobjfun
00185          -neps;
00186      }
00187    }
00188 
00189    nsc=pmin->lapprox->num_separable_calls;
00190    dmatrix weights(1,nsc,1,nsamp);
00191    for (is=1;is<=nsamp;is++)
00192    {
00193      int offset=0;
00194      int ic;
00195      for (ic=1;ic<=nsc;ic++)
00196      {
00197        int lus=lrea(ic);
00198        // assumes that all spearable calls have the same number of
00199        // random effects
00200        dvector e= (*pmin->lapprox->antiepsilon)(is);
00201        offset+=lus;
00202        if (pmin->lapprox->use_outliers==0)
00203        {
00204          weights(ic,is)=0.5*norm2(e);
00205        }
00206        else
00207        {
00208          weights(ic,is)=-sum(log(.95*exp(-0.5*square(e))+
00209           0.0166666667*exp(-square(e)/18.0)));
00210        }
00211      }
00212    }
00213    dvar_vector lcomp(1,nsc);
00214    for (int isc=1;isc<=nsc;isc++)
00215    {
00216      dvar_vector & comps=
00217        (*pmin->lapprox->importance_sampling_components)(isc);
00218 
00219      dvar_vector diff=comps-weights(isc);
00220      double dmin=min(value(diff));
00221      lcomp(isc)=dmin-log(mean(exp(dmin-diff)));
00222    }
00223 
00224    //double ns=lcomp.indexmax()-lcomp.indexmin()+1;
00225    //double min_vf=min(value(lcomp));
00226    vf= sum(lcomp);
00227    vf-=us*0.91893853320467241;
00228 
00229    int sgn=0;
00230    dvariable ld=0.0;
00231    if (ad_comm::no_ln_det_choleski_flag)
00232    {
00233      for (int ic=1;ic<=nsc;ic++)
00234      {
00235        ld+=ln_det(block_diagonal_vhessian(ic),sgn);
00236      }
00237      ld*=0.5;
00238    }
00239    else
00240    {
00241      for (int ic=1;ic<=nsc;ic++)
00242      {
00243        ld+=ln_det_choleski(block_diagonal_vhessian(ic));
00244      }
00245      ld*=0.5;
00246    }
00247 
00248    vf+=ld;
00249 
00250    double f=value(vf);
00251    gradcalc(nvar,g);
00252 
00253    // put uhat back into the model
00254    gradient_structure::set_NO_DERIVATIVES();
00255    vy(xs+1,xs+us).shift(1)=u0;
00256    initial_params::reset(vy);    // get the values into the model
00257    gradient_structure::set_YES_DERIVATIVES();
00258 
00259   ii=1;
00260   for (i=1;i<=xs;i++)
00261     xadjoint(i)=g(ii++);
00262   for (i=1;i<=us;i++)
00263     uadjoint(i)=g(ii++);
00264   for (ic=1;ic<=nsc;ic++)
00265   {
00266     int lus=lrea(ic);
00267     for (i=1;i<=lus;i++)
00268     {
00269       for (j=1;j<=lus;j++)
00270       {
00271         (*pmin->lapprox->block_diagonal_vhessianadjoint)(ic)(i,j)=g(ii++);
00272       }
00273     }
00274   }
00275   return f;
00276 }