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