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