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