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