ADMB Documentation  11.1x.2738
 All Classes Files Functions Variables Typedefs Friends Defines
df1b2impf.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: df1b2impf.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_funnel(const dvector& x,const dvector& u0,
00020   const dmatrix& Hess,const dvector& _xadjoint,const dvector& _uadjoint,
00021   const dmatrix& _Hessadjoint,function_minimizer * pmin)
00022 {
00023   ADUNCONST(dvector,xadjoint)
00024   ADUNCONST(dvector,uadjoint)
00025   ADUNCONST(dmatrix,Hessadjoint)
00026   const int xs=x.size();
00027   const int us=u0.size();
00028   gradient_structure::set_YES_DERIVATIVES();
00029   int nvar=x.size()+u0.size()+u0.size()*u0.size();
00030   independent_variables y(1,nvar);
00031 
00032   // need to set random effects active together with whatever
00033   // init parameters should be active in this phase
00034   initial_params::set_inactive_only_random_effects();
00035   initial_params::set_active_random_effects();
00036   /*int onvar=*/initial_params::nvarcalc();
00037   initial_params::xinit(y);    // get the initial values into the
00038   // do we need this next line?
00039   y(1,xs)=x;
00040 
00041   int i,j;
00042   dvar_vector d(1,xs+us);
00043 
00044   // contribution for quadratic prior
00045   if (quadratic_prior::get_num_quadratic_prior()>0)
00046   {
00047     //Hess+=quadratic_prior::get_cHessian_contribution();
00048     int & vxs = (int&)(xs);
00049     quadratic_prior::get_cHessian_contribution(Hess,vxs);
00050   }
00051  // Here need hooks for sparse matrix structures
00052   int ii=xs+us+1;
00053   for (i=1;i<=us;i++)
00054     for (j=1;j<=us;j++)
00055     y(ii++)=Hess(i,j);
00056 
00057   dvar_vector vy=dvar_vector(y);
00058   initial_params::stddev_vscale(d,vy);
00059   dvar_matrix vHess(1,us,1,us);
00060   ii=xs+us+1;
00061   if (initial_df1b2params::have_bounded_random_effects)
00062   {
00063     cerr << "can't do importance sampling with bounded random effects"
00064      " at present" << endl;
00065     ad_exit(1);
00066   }
00067   else
00068   {
00069     for (i=1;i<=us;i++)
00070     {
00071       for (j=1;j<=us;j++)
00072       {
00073         vHess(i,j)=vy(ii++);
00074       }
00075     }
00076   }
00077 
00078    int nsamp=pmin->lapprox->num_importance_samples;
00079    dvar_vector sample_value(1,nsamp);
00080    sample_value.initialize();
00081    dvar_matrix ch=choleski_decomp(inv(vHess));
00082 
00083    dvariable vf=0.0;
00084 
00085    // **************************************************************
00086    // **************************************************************
00087    // **************************************************************
00088    int& nfunnelblocks=pmin->lapprox->nfunnelblocks;
00089    int lbound=1;
00090    int samplesize=nsamp/nfunnelblocks;
00091    int ubound=samplesize;
00092    int mean_count=0;
00093    dvar_vector fvalues(1,nfunnelblocks);
00094    ivector blocksizes(1,nfunnelblocks);
00095    for (int iblock=1;iblock<=nfunnelblocks;iblock++)
00096    {
00097      funnel_dvariable fdv;
00098      if (lbound>nsamp) break;
00099      ad_begin_funnel();
00100      int icount=0;
00101 
00102      for (int is=lbound;is<=ubound;is++)
00103      {
00104        if (is>nsamp) break;
00105        icount++;
00106        dvar_vector tau=ch*pmin->lapprox->epsilon(is);
00107 
00108        vy(xs+1,xs+us).shift(1)+=tau;
00109        initial_params::reset(vy);    // get the values into the model
00110        vy(xs+1,xs+us).shift(1)-=tau;
00111 
00112        *objective_function_value::pobjfun=0.0;
00113        pmin->AD_uf_outer();
00114 
00115        if (pmin->lapprox->use_outliers==0)
00116        {
00117          sample_value(icount)=*objective_function_value::pobjfun
00118            -0.5*norm2(pmin->lapprox->epsilon(is));
00119        }
00120        else
00121        {
00122          dvector& e=pmin->lapprox->epsilon(is);
00123 
00124          sample_value(icount)=*objective_function_value::pobjfun
00125            +sum(log(.95*exp(-0.5*square(e))+.05/3.0*exp(-square(e)/18.0)));
00126        }
00127      }
00128 
00129      if (icount>0)
00130      {
00131        mean_count+=1;
00132        dvar_vector tsp=sample_value(1,icount);
00133        double min_vf=min(value(tsp));
00134        fdv=log(mean(exp(min_vf-tsp)))-min_vf;
00135        dvariable tmp;
00136        tmp=fdv;
00137        fvalues(mean_count)=tmp;
00138        blocksizes(mean_count)=icount;
00139      }
00140      lbound+=samplesize;
00141      ubound+=samplesize;
00142    }
00143 
00144    double fm=mean(value(fvalues(1,mean_count)));
00145    dvar_vector nfval=exp(fvalues(1,mean_count)-fm);
00146    vf=-fm-log(nfval*blocksizes(1,mean_count)/sum(blocksizes(1,mean_count)));
00147    //vf-=us*.91893853320467241;
00148 
00149    int sgn=0;
00150    dvariable ld;
00151    if (ad_comm::no_ln_det_choleski_flag)
00152      ld=0.5*ln_det(vHess,sgn);
00153    else
00154      ld=0.5*ln_det_choleski(vHess);
00155 
00156    vf+=ld;
00157    double f=value(vf);
00158    dvector g(1,nvar);
00159    gradcalc(nvar,g);
00160 
00161    // put uhat back into the model
00162    gradient_structure::set_NO_DERIVATIVES();
00163    vy(xs+1,xs+us).shift(1)=u0;
00164    initial_params::reset(vy);    // get the values into the model
00165    gradient_structure::set_YES_DERIVATIVES();
00166 
00167   ii=1;
00168   for (i=1;i<=xs;i++)
00169     xadjoint(i)=g(ii++);
00170   for (i=1;i<=us;i++)
00171     uadjoint(i)=g(ii++);
00172   for (i=1;i<=us;i++)
00173     for (j=1;j<=us;j++)
00174       Hessadjoint(i,j)=g(ii++);
00175   return f;
00176 }