ADMB Documentation  11.1.1927
 All Classes Files Functions Variables Typedefs Friends Defines
df1b2im2.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: df1b2im2.cpp 1919 2014-04-22 22:02:01Z 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 banded_symmetric_dmatrix& bHess,const dvector& _xadjoint,
00022   const dvector& _uadjoint,
00023   const banded_symmetric_dmatrix& _bHessadjoint,function_minimizer * pmin)
00024 {
00025   ADUNCONST(dvector,xadjoint)
00026   ADUNCONST(dvector,uadjoint)
00027   ADUNCONST(banded_symmetric_dmatrix,bHessadjoint)
00028   int bw=bHess.bandwidth();
00029   const int xs=x.size();
00030   const int us=u0.size();
00031   gradient_structure::set_YES_DERIVATIVES();
00032   int nvar=x.size()+u0.size()+((bw+1)*(2*u0.size()-bw))/2;
00033   independent_variables y(1,nvar);
00034 
00035   // need to set random effects active together with whatever
00036   // init parameters should be active in this phase
00037   initial_params::set_inactive_only_random_effects();
00038   initial_params::set_active_random_effects();
00039   /*int onvar=*/initial_params::nvarcalc();
00040   initial_params::xinit(y);    // get the initial values into the
00041   y(1,xs)=x;
00042 
00043   int i,j;
00044 
00045  // Here need hooks for sparse matrix structures
00046   int ii=xs+us+1;
00047   for (i=1;i<=us;i++)
00048   {
00049     int jmin=admax(1,i-bw+1);
00050     for (j=jmin;j<=i;j++)
00051     y(ii++)=bHess(i,j);
00052   }
00053 
00054   dvar_vector vy=dvar_vector(y);
00055   banded_symmetric_dvar_matrix vHess(1,us,bw);
00056   initial_params::reset(vy);    // get the initial values into the
00057   ii=xs+us+1;
00058   for (i=1;i<=us;i++)
00059   {
00060     int jmin=admax(1,i-bw+1);
00061     for (j=jmin;j<=i;j++)
00062       vHess(i,j)=vy(ii++);
00063   }
00064 
00065    int nsamp=pmin->lapprox->num_importance_samples;
00066    dvar_vector sample_value(1,nsamp);
00067    sample_value.initialize();
00068 
00069    int ierr = 0;
00070    banded_lower_triangular_dvar_matrix ch=choleski_decomp(vHess,ierr);
00071    if (ierr)
00072    {
00073      cerr << "error in choleski decomp" << endl;
00074      ad_exit(1);
00075    }
00076 
00077    dvariable vf=0.0;
00078    if (laplace_approximation_calculator::
00079      print_importance_sampling_weights_flag==0)
00080    {
00081      for (int is=1;is<=nsamp;is++)
00082      {
00083        dvar_vector tau=solve_trans(ch,pmin->lapprox->epsilon(is));
00084 
00085        vy(xs+1,xs+us).shift(1)+=tau;
00086        initial_params::reset(vy);    // get the values into the model
00087        vy(xs+1,xs+us).shift(1)-=tau;
00088 
00089        *objective_function_value::pobjfun=0.0;
00090        pmin->AD_uf_outer();
00091 
00092        sample_value(is)=*objective_function_value::pobjfun
00093          -0.5*norm2(pmin->lapprox->epsilon(is));
00094      }
00095    }
00096    else
00097    {
00098      dvector normal_weight(1,nsamp);
00099      int is;
00100      for (is=1;is<=nsamp;is++)
00101      {
00102        dvar_vector tau=solve_trans(ch,pmin->lapprox->epsilon(is));
00103 
00104        vy(xs+1,xs+us).shift(1)+=tau;
00105        initial_params::reset(vy);    // get the values into the model
00106        vy(xs+1,xs+us).shift(1)-=tau;
00107 
00108        *objective_function_value::pobjfun=0.0;
00109        pmin->AD_uf_outer();
00110 
00111        sample_value(is)=*objective_function_value::pobjfun;
00112        normal_weight(is)=0.5*norm2(pmin->lapprox->epsilon(is));
00113      }
00114      dvector tmp1(value(sample_value)-normal_weight);
00115      double min_vf=min(tmp1);
00116      dvector tmp=exp(tmp1-min_vf);
00117      cout << "The unsorted normalized importance samplng weights are " << endl
00118           << tmp << endl;
00119      cout << "The sorted normalized importance samplng weights are " << endl
00120           << sort(tmp) << endl;
00121      cout << "The sample value normal weight pairs are " << endl;
00122      for (is=1;is<=nsamp;is++)
00123      {
00124        cout << sample_value(is) << "  " << normal_weight(is) << endl;
00125      }
00126      cout << "The normalized sample value normal weight pairs are " << endl;
00127      for (is=1;is<=nsamp;is++)
00128      {
00129        cout << normal_weight(is) << "  "
00130             << sample_value(is)-normal_weight(is) << endl;
00131      }
00132      ad_exit(1);
00133    }
00134 
00135    dvariable min_vf=min(sample_value);
00136    vf=min_vf-log(mean(exp(min_vf-sample_value)));
00137    vf-=us*0.91893853320467241;
00138 
00139    int sgn=0;
00140    dvariable ld;
00141 
00142    ld=0.5*ln_det_choleski(vHess,sgn);
00143 
00144    vf+=ld;
00145    double f=value(vf);
00146    dvector g(1,nvar);
00147    gradcalc(nvar,g);
00148 
00149   ii=1;
00150   for (i=1;i<=xs;i++)
00151     xadjoint(i)=g(ii++);
00152   for (i=1;i<=us;i++)
00153     uadjoint(i)=g(ii++);
00154   for (i=1;i<=us;i++)
00155   {
00156     int jmin=admax(1,i-bw+1);
00157     for (j=jmin;j<=i;j++)
00158       bHessadjoint(i,j)=g(ii++);
00159   }
00160   return f;
00161 }
00162 #endif //#if defined(USE_LAPLACE)