ADMB Documentation  11.1.2192
 All Classes Files Functions Variables Typedefs Friends Defines
df1b2imp.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: df1b2imp.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 calculate_importance_sample(const dvector& x,const dvector& u0,
00020   const dmatrix& Hess,const dvector& _xadjoint,const dvector& _uadjoint,
00021   const dmatrix& _Hessadjoint,function_minimizer * pmin)
00022 {
00023   ADUNCONST(dvector,xadjoint)
00024   ADUNCONST(dvector,uadjoint)
00025   ADUNCONST(dmatrix,Hessadjoint)
00026   const int xs=x.size();
00027   const int us=u0.size();
00028   gradient_structure::set_YES_DERIVATIVES();
00029   int nvar=x.size()+u0.size()+u0.size()*u0.size();
00030   independent_variables y(1,nvar);
00031 
00032   // need to set random effects active together with whatever
00033   // init parameters should be active in this phase
00034   initial_params::set_inactive_only_random_effects();
00035   initial_params::set_active_random_effects();
00036   /*int onvar=*/initial_params::nvarcalc();
00037   initial_params::xinit(y);    // get the initial values into the
00038   // do we need this next line?
00039   y(1,xs)=x;
00040 
00041   int i,j;
00042   dvar_vector d(1,xs+us);
00043 
00044   // contribution for quadratic prior
00045   if (quadratic_prior::get_num_quadratic_prior()>0)
00046   {
00047     //Hess+=quadratic_prior::get_cHessian_contribution();
00048     int & vxs = (int&)(xs);
00049     quadratic_prior::get_cHessian_contribution(Hess,vxs);
00050   }
00051  // Here need hooks for sparse matrix structures
00052   int ii=xs+us+1;
00053   for (i=1;i<=us;i++)
00054     for (j=1;j<=us;j++)
00055     y(ii++)=Hess(i,j);
00056 
00057   dvar_vector vy=dvar_vector(y);
00058   initial_params::stddev_vscale(d,vy);
00059   dvar_matrix vHess(1,us,1,us);
00060   ii=xs+us+1;
00061   if (initial_df1b2params::have_bounded_random_effects)
00062   {
00063     cerr << "can't do importance sampling with bounded random effects"
00064      " at present" << endl;
00065     ad_exit(1);
00066   }
00067   else
00068   {
00069     for (i=1;i<=us;i++)
00070     {
00071       for (j=1;j<=us;j++)
00072       {
00073         vHess(i,j)=vy(ii++);
00074       }
00075     }
00076   }
00077 
00078    int nsamp=pmin->lapprox->num_importance_samples;
00079    dvar_vector sample_value(1,nsamp);
00080    sample_value.initialize();
00081    dvar_matrix ch=choleski_decomp(inv(vHess));
00082 
00083    dvariable vf=0.0;
00084    for (int is=1;is<=nsamp;is++)
00085    {
00086      dvar_vector tau=ch*pmin->lapprox->epsilon(is);
00087 
00088      vy(xs+1,xs+us).shift(1)+=tau;
00089      initial_params::reset(vy);    // get the values into the model
00090      vy(xs+1,xs+us).shift(1)-=tau;
00091 
00092      *objective_function_value::pobjfun=0.0;
00093      pmin->AD_uf_outer();
00094 
00095      sample_value(is)=*objective_function_value::pobjfun
00096        -0.5*norm2(pmin->lapprox->epsilon(is));
00097    }
00098 
00099 
00100    if (laplace_approximation_calculator::
00101      print_importance_sampling_weights_flag==1)
00102    {
00103      double min_vf=min(value(sample_value));
00104      dvector tmp=exp(value(sample_value)-min_vf);
00105      cout << "The unsorted normalized importance samplng weights are " << endl
00106           << tmp << endl;
00107      cout << "The sorted normalized importance samplng weights are " << endl
00108           << sort(tmp) << endl;
00109      ad_exit(1);
00110    }
00111 
00112 
00113    dvariable min_vf=min(sample_value);
00114    vf=min_vf-log(mean(exp(min_vf-sample_value)));
00115    vf-=us*0.91893853320467241;
00116 
00117    int sgn=0;
00118    dvariable ld;
00119    if (ad_comm::no_ln_det_choleski_flag)
00120      ld=0.5*ln_det(vHess,sgn);
00121    else
00122      ld=0.5*ln_det_choleski(vHess);
00123 
00124    vf+=ld;
00125    double f=value(vf);
00126    dvector g(1,nvar);
00127    gradcalc(nvar,g);
00128 
00129    // put uhat back into the model
00130    gradient_structure::set_NO_DERIVATIVES();
00131    vy(xs+1,xs+us).shift(1)=u0;
00132    initial_params::reset(vy);    // get the values into the model
00133    gradient_structure::set_YES_DERIVATIVES();
00134 
00135   ii=1;
00136   for (i=1;i<=xs;i++)
00137     xadjoint(i)=g(ii++);
00138   for (i=1;i<=us;i++)
00139     uadjoint(i)=g(ii++);
00140   for (i=1;i<=us;i++)
00141     for (j=1;j<=us;j++)
00142       Hessadjoint(i,j)=g(ii++);
00143   return f;
00144 }