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