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