ADMB Documentation  11.1.1932
 All Classes Files Functions Variables Typedefs Friends Defines
df1b2impf.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: df1b2impf.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_funnel(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 
00086    // **************************************************************
00087    // **************************************************************
00088    // **************************************************************
00089    int& nfunnelblocks=pmin->lapprox->nfunnelblocks;
00090    int lbound=1;
00091    int samplesize=nsamp/nfunnelblocks;
00092    int ubound=samplesize;
00093    int mean_count=0;
00094    dvar_vector fvalues(1,nfunnelblocks);
00095    ivector blocksizes(1,nfunnelblocks);
00096    for (int iblock=1;iblock<=nfunnelblocks;iblock++)
00097    {
00098      funnel_dvariable fdv;
00099      if (lbound>nsamp) break;
00100      ad_begin_funnel();
00101      int icount=0;
00102 
00103      for (int is=lbound;is<=ubound;is++)
00104      {
00105        if (is>nsamp) break;
00106        icount++;
00107        dvar_vector tau=ch*pmin->lapprox->epsilon(is);
00108 
00109        vy(xs+1,xs+us).shift(1)+=tau;
00110        initial_params::reset(vy);    // get the values into the model
00111        vy(xs+1,xs+us).shift(1)-=tau;
00112 
00113        *objective_function_value::pobjfun=0.0;
00114        pmin->AD_uf_outer();
00115 
00116        if (pmin->lapprox->use_outliers==0)
00117        {
00118          sample_value(icount)=*objective_function_value::pobjfun
00119            -0.5*norm2(pmin->lapprox->epsilon(is));
00120        }
00121        else
00122        {
00123          dvector& e=pmin->lapprox->epsilon(is);
00124 
00125          sample_value(icount)=*objective_function_value::pobjfun
00126            +sum(log(.95*exp(-0.5*square(e))+.05/3.0*exp(-square(e)/18.0)));
00127        }
00128      }
00129 
00130      if (icount>0)
00131      {
00132        mean_count+=1;
00133        dvar_vector tsp=sample_value(1,icount);
00134        double min_vf=min(value(tsp));
00135        fdv=log(mean(exp(min_vf-tsp)))-min_vf;
00136        dvariable tmp;
00137        tmp=fdv;
00138        fvalues(mean_count)=tmp;
00139        blocksizes(mean_count)=icount;
00140      }
00141      lbound+=samplesize;
00142      ubound+=samplesize;
00143    }
00144 
00145    double fm=mean(value(fvalues(1,mean_count)));
00146    dvar_vector nfval=exp(fvalues(1,mean_count)-fm);
00147    vf=-fm-log(nfval*blocksizes(1,mean_count)/sum(blocksizes(1,mean_count)));
00148    //vf-=us*.91893853320467241;
00149 
00150    int sgn=0;
00151    dvariable ld;
00152    if (ad_comm::no_ln_det_choleski_flag)
00153      ld=0.5*ln_det(vHess,sgn);
00154    else
00155      ld=0.5*ln_det_choleski(vHess);
00156 
00157    vf+=ld;
00158    double f=value(vf);
00159    dvector g(1,nvar);
00160    gradcalc(nvar,g);
00161 
00162    // put uhat back into the model
00163    gradient_structure::set_NO_DERIVATIVES();
00164    vy(xs+1,xs+us).shift(1)=u0;
00165    initial_params::reset(vy);    // get the values into the model
00166    gradient_structure::set_YES_DERIVATIVES();
00167 
00168   ii=1;
00169   for (i=1;i<=xs;i++)
00170     xadjoint(i)=g(ii++);
00171   for (i=1;i<=us;i++)
00172     uadjoint(i)=g(ii++);
00173   for (i=1;i<=us;i++)
00174     for (j=1;j<=us;j++)
00175       Hessadjoint(i,j)=g(ii++);
00176   return f;
00177 }
00178 #endif