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