ADMB Documentation  11.1.1897
 All Classes Files Functions Variables Typedefs Friends Defines
df1b2impspf.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: df1b2impspf.cpp 1755 2014-03-05 22:21:42Z 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 
00016 class dvar_hs_smatrix;
00017 
00022 double calculate_importance_sample_shess(const dvector& x,const dvector& u0,
00023   const dmatrix& Hess,const dvector& _xadjoint,const dvector& _uadjoint,
00024   const dmatrix& _Hessadjoint,function_minimizer * pmin)
00025 {
00026   ADUNCONST(dvector,xadjoint)
00027   ADUNCONST(dvector,uadjoint)
00028   const int xs=x.size();
00029   const int us=u0.size();
00030   gradient_structure::set_YES_DERIVATIVES();
00031 
00032   int nvar=0;
00033   if (pmin->lapprox->sparse_hessian_flag==0)
00034   {
00035     cerr << "Error we should not be here" << endl;
00036     ad_exit(1);
00037   }
00038 
00039   dcompressed_triplet & lst = *(pmin->lapprox->sparse_triplet2);
00040   int smin=lst.indexmin();
00041   int smax=lst.indexmax();
00042   int sz= smax-smin+1;
00043   nvar=x.size()+u0.size()+sz;
00044   independent_variables y(1,nvar);
00045 
00046   // need to set random effects active together with whatever
00047   // init parameters should be active in this phase
00048   initial_params::set_inactive_only_random_effects();
00049   initial_params::set_active_random_effects();
00050   /*int onvar=*/initial_params::nvarcalc();
00051   initial_params::xinit(y);    // get the initial values into the
00052   // do we need this next line?
00053   y(1,xs)=x;
00054 
00055   dvar_vector d(1,xs+us);
00056 
00057   // contribution for quadratic prior
00058   if (quadratic_prior::get_num_quadratic_prior()>0)
00059   {
00060     //Hess+=quadratic_prior::get_cHessian_contribution();
00061     int & vxs = (int&)(xs);
00062     quadratic_prior::get_cHessian_contribution(Hess,vxs);
00063   }
00064  // Here need hooks for sparse matrix structures
00065 
00066   int ii=xs+us+1;
00067   for (int i=smin;i<=smax;i++)
00068     y(ii++)=lst(i);
00069 
00070   dvar_vector vy=dvar_vector(y);
00071   initial_params::stddev_vscale(d,vy);
00072 
00073 
00074   ii=xs+us+1;
00075   if (initial_df1b2params::have_bounded_random_effects)
00076   {
00077     cerr << "can't do importance sampling with bounded random effects"
00078      " at present" << endl;
00079     ad_exit(1);
00080   }
00081   else
00082   {
00083     dcompressed_triplet & lst = *(pmin->lapprox->sparse_triplet2);
00084     int mmin=lst.indexmin();
00085     int mmax=lst.indexmax();
00086     dvar_compressed_triplet * vsparse_triplet =
00087       pmin->lapprox->vsparse_triplet;
00088 
00089     if (vsparse_triplet==0)
00090     {
00091       pmin->lapprox->vsparse_triplet=
00092         new dvar_compressed_triplet(mmin,mmax,us,us);
00093       vsparse_triplet = pmin->lapprox->vsparse_triplet;
00094       for (int i=mmin;i<=mmax;i++)
00095       {
00096         (*vsparse_triplet)(1,i)=lst(1,i);
00097         (*vsparse_triplet)(2,i)=lst(2,i);
00098       }
00099     }
00100     else
00101     {
00102       if (!allocated(*vsparse_triplet))
00103       {
00104         (*vsparse_triplet).allocate(mmin,mmax,us,us);
00105         for (int i=mmin;i<=mmax;i++)
00106         {
00107           (*vsparse_triplet)(1,i)=lst(1,i);
00108           (*vsparse_triplet)(2,i)=lst(2,i);
00109         }
00110       }
00111     }
00112     dcompressed_triplet * vsparse_triplet_adjoint =
00113       pmin->lapprox->vsparse_triplet_adjoint;
00114 
00115     if (vsparse_triplet_adjoint==0)
00116     {
00117       pmin->lapprox->vsparse_triplet_adjoint=
00118         new dcompressed_triplet(mmin,mmax,us,us);
00119       vsparse_triplet_adjoint = pmin->lapprox->vsparse_triplet_adjoint;
00120       for (int i=mmin;i<=mmax;i++)
00121       {
00122         (*vsparse_triplet_adjoint)(1,i)=lst(1,i);
00123         (*vsparse_triplet_adjoint)(2,i)=lst(2,i);
00124       }
00125     }
00126     else
00127     {
00128       if (!allocated(*vsparse_triplet_adjoint))
00129       {
00130         (*vsparse_triplet_adjoint).allocate(mmin,mmax,us,us);
00131         for (int i=mmin;i<=mmax;i++)
00132         {
00133           (*vsparse_triplet_adjoint)(1,i)=lst(1,i);
00134           (*vsparse_triplet_adjoint)(2,i)=lst(2,i);
00135         }
00136       }
00137     }
00138     vsparse_triplet->get_x()=vy(ii,ii+mmax-mmin).shift(1);
00139   }
00140 
00141    int nsamp=pmin->lapprox->num_importance_samples;
00142    dvar_vector sample_value(1,nsamp);
00143    sample_value.initialize();
00144    //dvar_matrix ch=choleski_decomp(inv(vHess));
00145 
00146    dvar_compressed_triplet * vsparse_triplet
00147      = pmin->lapprox->vsparse_triplet;
00148 
00149    dvar_hs_smatrix * dhs= return_choleski_decomp(*vsparse_triplet);
00150 
00151    dvariable vf=0.0;
00152 
00153    // get fhat for ocputing weighted sums.
00154    initial_params::reset(vy);    // get the values into the model
00155    *objective_function_value::pobjfun=0.0;
00156     pmin->AD_uf_outer();
00157     double fhat=value(*objective_function_value::pobjfun);
00158    // check if we do this in a group of funnels
00159    if (pmin->lapprox->isfunnel_flag==0)
00160    {
00161      for (int is=1;is<=nsamp;is++)
00162      {
00163        //dvar_vector tau=ch*pmin->lapprox->epsilon(is);
00164        dvar_vector tau=return_choleski_factor_solve(dhs,
00165          pmin->lapprox->epsilon(is));
00166 
00167        vy(xs+1,xs+us).shift(1)+=tau;
00168        initial_params::reset(vy);    // get the values into the model
00169        vy(xs+1,xs+us).shift(1)-=tau;
00170 
00171        *objective_function_value::pobjfun=0.0;
00172        pmin->AD_uf_outer();
00173 
00174       /*
00175        if (pmin->lapprox->istudent_flag==0)
00176       */
00177        {
00178          sample_value(is)=*objective_function_value::pobjfun
00179            -0.5*norm2(pmin->lapprox->epsilon(is));
00180        }
00181       /*
00182        else
00183        {
00184          sample_value(is)=*objective_function_value::pobjfun
00185            +sum(log_student_density(pmin->lapprox->istudent_flag,
00186            pmin->lapprox->epsilon(is)));
00187        }
00188       */
00189        if (pmin->lapprox->importance_sampling_values)
00190          (*(pmin->lapprox->importance_sampling_values))(is)
00191           =value(sample_value(is))-fhat;
00192      }
00193      dvariable min_vf=min(sample_value);
00194      vf=min_vf-log(mean(exp(min_vf-sample_value)));
00195    }
00196    else
00197    {
00198      int& nfunnelblocks=pmin->lapprox->nfunnelblocks;
00199      int lbound=1;
00200      int samplesize=nsamp/nfunnelblocks;
00201      int ubound=samplesize;
00202      int mean_count=0;
00203      dvar_vector fvalues(1,nfunnelblocks);
00204      ivector blocksizes(1,nfunnelblocks);
00205      for (int iblock=1;iblock<=nfunnelblocks;iblock++)
00206      {
00207        funnel_dvariable fdv;
00208        if (lbound>nsamp) break;
00209        ad_begin_funnel();
00210        int icount=0;
00211 
00212 
00213        for (int is=lbound;is<=ubound;is++)
00214        {
00215          if (is>nsamp) break;
00216          icount++;
00217          dvar_vector tau=return_choleski_factor_solve(dhs,
00218            pmin->lapprox->epsilon(is));
00219          vy(xs+1,xs+us).shift(1)+=tau;
00220          initial_params::reset(vy);    // get the values into the model
00221          vy(xs+1,xs+us).shift(1)-=tau;
00222 
00223          *objective_function_value::pobjfun=0.0;
00224          pmin->AD_uf_outer();
00225        /*
00226          if (pmin->lapprox->istudent_flag==0)
00227        */
00228          {
00229            sample_value(icount)=*objective_function_value::pobjfun
00230              -0.5*norm2(pmin->lapprox->epsilon(is));
00231          }
00232         /*
00233          else
00234          {
00235            sample_value(icount)=*objective_function_value::pobjfun
00236              +sum(log_student_density(pmin->lapprox->istudent_flag,
00237               pmin->lapprox->epsilon(is)));
00238          }
00239        */
00240          if (pmin->lapprox->importance_sampling_values)
00241            (*(pmin->lapprox->importance_sampling_values))(is)
00242              =value(sample_value(icount))-fhat;
00243        }
00244 
00245        if (icount>0)
00246        {
00247          mean_count+=1;
00248          dvar_vector tsp=sample_value(1,icount);
00249          double min_vf=min(value(tsp));
00250          fdv=log(mean(exp(min_vf-tsp)))-min_vf;
00251          dvariable tmp;
00252          tmp=fdv;
00253          fvalues(mean_count)=tmp;
00254          blocksizes(mean_count)=icount;
00255        }
00256        lbound+=samplesize;
00257        ubound+=samplesize;
00258      }
00259 
00260      double fm=mean(value(fvalues(1,mean_count)));
00261      dvar_vector nfval=exp(fvalues(1,mean_count)-fm);
00262      vf=-fm-log(nfval*blocksizes(1,mean_count)/sum(blocksizes(1,mean_count)));
00263    }
00264   /*
00265    if (pmin->lapprox->istudent_flag==0)
00266   */
00267    {
00268      vf-=us*0.91893853320467241;
00269    }
00270   /*
00271    else
00272    {
00273      (*(pmin->lapprox->importance_sampling_values))
00274         +=us*0.91893853320467241;
00275    }
00276   */
00277 
00278    dvariable ld;
00279 
00280    //dvariable sld=0.5*ln_det_choleski(make_dvar_matrix(*vsparse_triplet));
00281    dvariable sld=0.5*ln_det(*vsparse_triplet);
00282    vf+=sld;
00283 
00284    double f=value(vf);
00285    dvector g(1,nvar);
00286    gradcalc(nvar,g);
00287 
00288    // put uhat back into the model
00289    gradient_structure::set_NO_DERIVATIVES();
00290    vy(xs+1,xs+us).shift(1)=u0;
00291    initial_params::reset(vy);    // get the values into the model
00292    gradient_structure::set_YES_DERIVATIVES();
00293 
00294   ii=1;
00295   for (int i=1;i<=xs;i++)
00296     xadjoint(i)=g(ii++);
00297   for (int i=1;i<=us;i++)
00298     uadjoint(i)=g(ii++);
00299 
00300   dcompressed_triplet * vsparse_triplet_adjoint =
00301     pmin->lapprox->vsparse_triplet_adjoint;
00302   int mmin=vsparse_triplet_adjoint->indexmin();
00303   int mmax=vsparse_triplet_adjoint->indexmax();
00304   vsparse_triplet_adjoint->get_x()
00305     =g(ii,ii+mmax-mmin).shift(1);
00306 
00307   return f;
00308 }
00309 #endif