ADMB Documentation  11.1.1916
 All Classes Files Functions Variables Typedefs Friends Defines
df1b2im3.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: df1b2im3.cpp 1667 2014-02-24 23:53:59Z 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 
00016 static void xxx(void){;}
00017 
00022 double calculate_importance_sample_block_diagonal(const dvector& x,
00023   const dvector& u0,const dmatrix& Hess,const dvector& _xadjoint,
00024   const dvector& _uadjoint,const dmatrix& _Hessadjoint,
00025   function_minimizer * pmin)
00026 {
00027   ADUNCONST(dvector,xadjoint)
00028   ADUNCONST(dvector,uadjoint)
00029   //ADUNCONST(dmatrix,Hessadjoint)
00030   const int xs=x.size();
00031   const int us=u0.size();
00032   gradient_structure::set_NO_DERIVATIVES();
00033   int nsc=pmin->lapprox->num_separable_calls;
00034   const ivector lrea = (*pmin->lapprox->num_local_re_array)(1,nsc);
00035   int hroom =  int(sum(square(lrea)));
00036   int nvar=x.size()+u0.size()+hroom;
00037   independent_variables y(1,nvar);
00038 
00039   // need to set random effects active together with whatever
00040   // init parameters should be active in this phase
00041   initial_params::set_inactive_only_random_effects();
00042   initial_params::set_active_random_effects();
00043   /*int onvar=*/initial_params::nvarcalc();
00044   initial_params::xinit(y);    // get the initial values into the
00045   // do we need this next line?
00046   y(1,xs)=x;
00047 
00048   int i,j;
00049 
00050   // contribution for quadratic prior
00051   if (quadratic_prior::get_num_quadratic_prior()>0)
00052   {
00053     //Hess+=quadratic_prior::get_cHessian_contribution();
00054     int & vxs = (int&)(xs);
00055     quadratic_prior::get_cHessian_contribution(Hess,vxs);
00056   }
00057  // Here need hooks for sparse matrix structures
00058 
00059   dvar3_array & block_diagonal_vhessian=
00060     *pmin->lapprox->block_diagonal_vhessian;
00061   block_diagonal_vhessian.initialize();
00062   dvar3_array& block_diagonal_ch=
00063     *pmin->lapprox->block_diagonal_vch;
00064     //dvar3_array(*pmin->lapprox->block_diagonal_ch);
00065   int ii=xs+us+1;
00066   d3_array& bdH=(*pmin->lapprox->block_diagonal_hessian);
00067   int ic;
00068   for (ic=1;ic<=nsc;ic++)
00069   {
00070     int lus=lrea(ic);
00071     for (i=1;i<=lus;i++)
00072       for (j=1;j<=lus;j++)
00073         y(ii++)=bdH(ic)(i,j);
00074   }
00075 
00076   dvector g(1,nvar);
00077   gradcalc(0,g);
00078   gradient_structure::set_YES_DERIVATIVES();
00079   dvar_vector vy=dvar_vector(y);
00080   //initial_params::stddev_vscale(d,vy);
00081   ii=xs+us+1;
00082   if (initial_df1b2params::have_bounded_random_effects)
00083   {
00084     cerr << "can't do importance sampling with bounded random effects"
00085      " at present" << endl;
00086     ad_exit(1);
00087   }
00088   else
00089   {
00090     for (int ic=1;ic<=nsc;ic++)
00091     {
00092       int lus=lrea(ic);
00093       for (i=1;i<=lus;i++)
00094       {
00095         for (j=1;j<=lus;j++)
00096         {
00097           block_diagonal_vhessian(ic,i,j)=vy(ii++);
00098         }
00099       }
00100       block_diagonal_ch(ic)=
00101         choleski_decomp(inv(block_diagonal_vhessian(ic)));
00102     }
00103   }
00104 
00105    int nsamp=pmin->lapprox->num_importance_samples;
00106 
00107    dvariable vf=0.0;
00108 
00109    dvar_vector sample_value(1,nsamp);
00110    sample_value.initialize();
00111 
00112    dvar_vector tau(1,us);;
00113    int is;
00114    for (is=1;is<=nsamp;is++)
00115    {
00116      int offset=0;
00117      pmin->lapprox->importance_sampling_counter=is;
00118      for (ic=1;ic<=nsc;ic++)
00119      {
00120        int lus=lrea(ic);
00121        tau(offset+1,offset+lus).shift(1)=block_diagonal_ch(ic)*
00122          pmin->lapprox->epsilon(is)(offset+1,offset+lus).shift(1);
00123        offset+=lus;
00124      }
00125 
00126      // have to reorder the terms to match the block diagonal hessian
00127      imatrix & ls=*(pmin->lapprox->block_diagonal_re_list);
00128      int mmin=ls.indexmin();
00129      int mmax=ls.indexmax();
00130 
00131      int ii=1;
00132      int i;
00133      for (i=mmin;i<=mmax;i++)
00134      {
00135        int cmin=ls(i).indexmin();
00136        int cmax=ls(i).indexmax();
00137        for (int j=cmin;j<=cmax;j++)
00138        {
00139          vy(ls(i,j))+=tau(ii++);
00140        }
00141      }
00142      if (ii-1 != us)
00143      {
00144        cerr << "error in interface" << endl;
00145        ad_exit(1);
00146      }
00147      initial_params::reset(vy);    // get the values into the model
00148      ii=1;
00149      for (i=mmin;i<=mmax;i++)
00150      {
00151        int cmin=ls(i).indexmin();
00152        int cmax=ls(i).indexmax();
00153        for (int j=cmin;j<=cmax;j++)
00154        {
00155          vy(ls(i,j))-=tau(ii++);
00156        }
00157      }
00158 
00159      *objective_function_value::pobjfun=0.0;
00160      //int istop=0;
00161      if (is==65)
00162      {
00163         xxx();
00164      }
00165      pmin->AD_uf_outer();
00166 
00167      if (pmin->lapprox->use_outliers==0)
00168      {
00169        double neps=0.5*norm2(pmin->lapprox->epsilon(is));
00170 
00171        (*pmin->lapprox->importance_sampling_values)(is)=
00172            neps-value(*objective_function_value::pobjfun);
00173 
00174        (*pmin->lapprox->importance_sampling_weights)(is)=neps;
00175 
00176         sample_value(is)=*objective_function_value::pobjfun
00177           -neps;
00178      }
00179      else
00180      {
00181        dvector& e=pmin->lapprox->epsilon(is);
00182        double neps=-sum(log(.95*exp(-0.5*square(e))+
00183           0.0166666667*exp(-square(e)/18.0)));
00184 
00185        (*pmin->lapprox->importance_sampling_values)(is)=
00186          neps-value(*objective_function_value::pobjfun);
00187 
00188        sample_value(is)=*objective_function_value::pobjfun
00189          -neps;
00190      }
00191    }
00192 
00193    nsc=pmin->lapprox->num_separable_calls;
00194    dmatrix weights(1,nsc,1,nsamp);
00195    for (is=1;is<=nsamp;is++)
00196    {
00197      int offset=0;
00198      for (ic=1;ic<=nsc;ic++)
00199      {
00200        int lus=lrea(ic);
00201        dvector e= pmin->lapprox->epsilon(is)(offset+1,offset+lus).shift(1);
00202        offset+=lus;
00203        if (pmin->lapprox->use_outliers==0)
00204        {
00205          weights(ic,is)=0.5*norm2(e);
00206        }
00207        else
00208        {
00209          weights(ic,is)=-sum(log(.95*exp(-0.5*square(e))+
00210           0.0166666667*exp(-square(e)/18.0)));
00211        }
00212      }
00213    }
00214    dvar_vector lcomp(1,nsc);
00215    for (int isc=1;isc<=nsc;isc++)
00216    {
00217      dvar_vector & comps=
00218        (*pmin->lapprox->importance_sampling_components)(isc);
00219 
00220      dvar_vector diff=comps-weights(isc);
00221      double dmin=min(value(diff));
00222      lcomp(isc)=-dmin+log(mean(exp(-diff+dmin)));
00223 
00224      for (int is=1;is<=nsamp;is++)
00225      {
00226      }
00227    }
00228 
00229    if (laplace_approximation_calculator::
00230      print_importance_sampling_weights_flag==1)
00231    {
00232      double min_vf=min(value(sample_value));
00233      dvector tmp=exp(value(sample_value)-min_vf);
00234      cout << "The unsorted normalized importance samplng weights are " << endl
00235           << tmp << endl;
00236      cout << "The sorted normalized importance samplng weights are " << endl
00237           << sort(tmp) << endl;
00238      ad_exit(1);
00239    }
00240 
00241 
00242    double min_vf=min(value(sample_value));
00243    vf=min_vf-log(mean(exp(min_vf-sample_value)));
00244    vf-=us*0.91893853320467241;
00245 
00246    int sgn=0;
00247    dvariable ld=0.0;
00248    if (ad_comm::no_ln_det_choleski_flag)
00249    {
00250      for (int ic=1;ic<=nsc;ic++)
00251      {
00252        ld+=ln_det(block_diagonal_vhessian(ic),sgn);
00253      }
00254      ld*=0.5;
00255    }
00256    else
00257    {
00258      for (int ic=1;ic<=nsc;ic++)
00259      {
00260        ld+=ln_det_choleski(block_diagonal_vhessian(ic));
00261      }
00262      ld*=0.5;
00263    }
00264 
00265    vf+=ld;
00266 
00267    double f=value(vf);
00268    gradcalc(nvar,g);
00269 
00270    // put uhat back into the model
00271    gradient_structure::set_NO_DERIVATIVES();
00272    vy(xs+1,xs+us).shift(1)=u0;
00273    initial_params::reset(vy);    // get the values into the model
00274    gradient_structure::set_YES_DERIVATIVES();
00275 
00276   ii=1;
00277   for (i=1;i<=xs;i++)
00278     xadjoint(i)=g(ii++);
00279   for (i=1;i<=us;i++)
00280     uadjoint(i)=g(ii++);
00281   for (ic=1;ic<=nsc;ic++)
00282   {
00283     int lus=lrea(ic);
00284     for (i=1;i<=lus;i++)
00285     {
00286       for (j=1;j<=lus;j++)
00287       {
00288         (*pmin->lapprox->block_diagonal_vhessianadjoint)(ic)(i,j)=g(ii++);
00289       }
00290     }
00291   }
00292   return f;
00293 }
00294 #endif