ADMB Documentation  11.1.2497
 All Classes Files Functions Variables Typedefs Friends Defines
df1b2lp2.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: df1b2lp2.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 void get_second_ders(int xs,int us,const init_df1b2vector y,dmatrix& Hess,
00015   dmatrix& Dux, df1b2_gradlist * f1b2gradlist,function_minimizer * pfmin,
00016   laplace_approximation_calculator* lap);
00017 double calculate_laplace_approximation(const dvector& x,const dvector& u0,
00018   const dmatrix& Hess,const dvector& _xadjoint,const dvector& _uadjoint,
00019   const dmatrix& _Hessadjoint,function_minimizer * pmin);
00020 
00021 static void xxx(ivector re_list,ivector fe_list){}
00022 
00027 dvector laplace_approximation_calculator::block_diagonal_calculations
00028   (const dvector& _x,const double& _f,function_minimizer * pfmin)
00029 {
00030   // for use when there is no separability
00031   ADUNCONST(dvector,x)
00032   ADUNCONST(double,f)
00033   //int i,j;
00034   int i;
00035 
00036   initial_params::set_inactive_only_random_effects();
00037   gradient_structure::set_NO_DERIVATIVES();
00038   initial_params::reset(x);    // get current x values into the model
00039   gradient_structure::set_YES_DERIVATIVES();
00040 
00041   initial_params::set_active_only_random_effects();
00042   //int lmn_flag=0;
00043   if (!inner_lmnflag)
00044   {
00045     if (!ADqd_flag)
00046       uhat=get_uhat_quasi_newton_block_diagonal(x,pfmin);
00047     else
00048       uhat=get_uhat_quasi_newton_qd(x,pfmin);
00049   }
00050   else
00051   {
00052     uhat=get_uhat_lm_newton(x,pfmin);
00053   }
00054   if (!allocated(scale))
00055   {
00056     scale.allocate(1,uhat.indexmax());
00057   }
00058   else
00059   {
00060     if (scale.indexmax() != uhat.indexmax())
00061     {
00062       scale.deallocate();
00063       scale.allocate(1,uhat.indexmax());
00064     }
00065   }
00066 
00067   if (!allocated(curv))
00068   {
00069     curv.allocate(1,uhat.indexmax());
00070   }
00071   else
00072   {
00073     if (curv.indexmax() != uhat.indexmax())
00074     {
00075       curv.deallocate();
00076       curv.allocate(1,uhat.indexmax());
00077     }
00078   }
00079 
00080   if (sparse_hessian_flag==0)
00081   {
00082     for (i=1;i<=xsize;i++)
00083     {
00084       y(i)=x(i);
00085     }
00086     for (int i=1;i<=usize;i++)
00087     {
00088       y(i+xsize)=uhat(i);
00089     }
00090   }
00091   else
00092   {
00093     for (i=1;i<=xsize;i++)
00094     {
00095       value(y(i))=x(i);
00096     }
00097     for (int i=1;i<=usize;i++)
00098     {
00099       value(y(i+xsize))=uhat(i);
00100     }
00101   }
00102 
00103   for(int ii=1;ii<=num_nr_iters;ii++)
00104   {
00105     {
00106       // test newton raphson
00107       //Hess.initialize();
00108       /*int check=*/initial_params::stddev_scale(scale,uhat);
00109       /*check=*/initial_params::stddev_curvscale(curv,uhat);
00110       max_separable_g=0.0;
00111       pmin->inner_opt_flag=1;
00112       step=get_newton_raphson_info_block_diagonal(pfmin);
00113       cout << "max separable g " << max_separable_g << endl;
00114       cout << "Newton raphson " << ii << endl;
00115       uhat+=step;
00116 
00117       evaluate_function(uhat,pfmin);
00118       pmin->inner_opt_flag=0;
00119     }
00120 
00121     if (sparse_hessian_flag==0)
00122     {
00123       for (int i=1;i<=usize;i++)
00124       {
00125         y(i+xsize)=uhat(i);
00126       }
00127     }
00128     else
00129     {
00130       for (int i=1;i<=usize;i++)
00131       {
00132         value(y(i+xsize))=uhat(i);
00133       }
00134     }
00135   }
00136 
00137   cout << initial_df1b2params::cobjfun << endl;
00138   xadjoint.initialize();
00139   uadjoint.initialize();
00140   block_diagonal_flag=2;
00141   used_flags.initialize();
00142   funnel_init_var::lapprox=this;
00143   if (use_gauss_hermite>0)
00144   {
00145     df1b2variable pen=0.0;
00146     initial_df1b2params::reset(y,pen);
00147     initial_params::straight_through_flag=0;
00148     block_diagonal_flag=6;
00149     num_separable_calls=0;
00150     // get the block diagonal hessians to use in importance sampling
00151     pfmin->user_function();
00152     //cout << (*pfmin->lapprox->block_diagonal_hessian) << endl;
00153     block_diagonal_flag=2;
00154     initial_params::straight_through_flag=0;
00155 
00156     // do importance sampling and get ders bakc to Hessian adjoint
00157     // new stuff for more than one random effect in each separable call
00158     //  Apr 17 07
00159     if (multi_random_effects==0)
00160     {
00161       f=do_gauss_hermite_block_diagonal(x,uhat,Hess,xadjoint,
00162         uadjoint,Hessadjoint,pfmin);
00163     }
00164     else
00165     {
00166       f=do_gauss_hermite_block_diagonal_multi(x,uhat,Hess,xadjoint,
00167         uadjoint,Hessadjoint,pfmin);
00168     }
00169     int xmax=xadjoint.indexmax();
00170     dvector x_con(1,xmax);
00171     x_con.initialize();
00172     dvector dscale(1,nvar);   // need to get scale from somewhere
00173     dscale=1.0;
00174     /*int check=*/initial_params::stddev_scale(dscale,x);
00175     dvector sscale=dscale(1,xsize);
00176     // *******************************************************
00177     // *******************************************************
00178     // *******************************************************
00179     // derivatives from hessian adjoint back
00180     {
00181       x_con.initialize();
00182 
00183       for (int i=1;i<=num_separable_calls;i++)
00184       {
00185         ivector& re_list=(*block_diagonal_re_list)(i);
00186         ivector& fe_list=(*block_diagonal_fe_list)(i);
00187         dmatrix& Dux=(*block_diagonal_Dux)(i);
00188         dmatrix& H=(*block_diagonal_hessian)(i);
00189         xxx(re_list,fe_list);
00190         int mmax=re_list.indexmax();
00191         dvector tmp(1,mmax);
00192 
00193         int j;
00194         for (j=1;j<=re_list.indexmax();j++)
00195         {
00196           tmp(j)=uadjoint(re_list(j)-xmax);
00197         }
00198 
00199         if (allocated(fe_list))
00200         {
00201           if (allocated(H))
00202           {
00203             dvector tmp1=solve(H,tmp);
00204             dvector xtmp=tmp1*Dux;
00205             for (j=1;j<=fe_list.indexmax();j++)
00206             {
00207               x_con(fe_list(j))+=xtmp(j);
00208             }
00209           }
00210         }
00211       }
00212       if (initial_df1b2params::separable_flag)
00213       {
00214         x_con=elem_prod(x_con,sscale);
00215       }
00216     }
00217     xadjoint-=x_con;
00218     // *******************************************************
00219     // *******************************************************
00220     // *******************************************************
00221 
00222     block_diagonal_flag=3;
00223     //pfmin->lapprox->xadjoint.initialize();
00224     //pfmin->lapprox->uadjoint.initialize();
00225     pfmin->lapprox->num_separable_calls=0;
00226     pfmin->lapprox->check_local_xadjoint.initialize();
00227     pfmin->lapprox->check_local_xadjoint2.initialize();
00228     pfmin->lapprox->check_local_uadjoint.initialize();
00229     pfmin->lapprox->check_local_uadjoint2.initialize();
00230     //df1b2_gradlist::set_yes_derivatives();
00231     //initial_df1b2params::reset(y,pen);
00232     pfmin->user_function();
00233     dvector lcx=elem_prod(check_local_xadjoint2,sscale);
00234     xadjoint+=lcx;
00235     //df1b2_gradlist::set_no_derivatives();
00236     funnel_init_var::lapprox=0;
00237     block_diagonal_flag=0;
00238     initial_params::set_inactive_only_random_effects();
00239   }
00240   else if (num_importance_samples>0)
00241   {
00242     df1b2variable pen=0.0;
00243     initial_df1b2params::reset(y,pen);
00244     initial_params::straight_through_flag=0;
00245     block_diagonal_flag=6;
00246     num_separable_calls=0;
00247     // get the block diagonal hessians to use in importance sampling
00248     pfmin->user_function();
00249     //cout << (*pfmin->lapprox->block_diagonal_hessian) << endl;
00250     block_diagonal_flag=2;
00251     initial_params::straight_through_flag=0;
00252     // do importance sampling and get ders bakc to Hessian adjoint
00253     if (isfunnel_flag==0)
00254     {
00255       if (antiflag==0)
00256       {
00257         f=calculate_importance_sample_block_diagonal_option2(x,uhat,Hess,
00258           xadjoint,uadjoint,Hessadjoint,pfmin);
00259       }
00260       else
00261       {
00262         f=calculate_importance_sample_block_diagonal_option_antithetical(x,
00263           uhat,Hess,xadjoint,uadjoint,Hessadjoint,pfmin);
00264       }
00265     }
00266     else
00267     {
00268       f=calculate_importance_sample_block_diagonal_funnel(x,uhat,Hess,xadjoint,
00269         uadjoint,Hessadjoint,pfmin);
00270     }
00271 
00272     int xmax=xadjoint.indexmax();
00273     dvector x_con(1,xmax);
00274     x_con.initialize();
00275     dvector dscale(1,nvar);   // need to get scale from somewhere
00276     dscale=1.0;
00277     /*int check=*/initial_params::stddev_scale(dscale,x);
00278     dvector sscale=dscale(1,xsize);
00279     // *******************************************************
00280     // *******************************************************
00281     // *******************************************************
00282     // derivatives from hessian adjoint back
00283     {
00284       x_con.initialize();
00285 
00286       for (int i=1;i<=num_separable_calls;i++)
00287       {
00288         dmatrix& H=(*block_diagonal_hessian)(i);
00289         if (allocated(H))
00290         {
00291           ivector& re_list=(*block_diagonal_re_list)(i);
00292           ivector& fe_list=(*block_diagonal_fe_list)(i);
00293           dmatrix& Dux=(*block_diagonal_Dux)(i);
00294           xxx(re_list,fe_list);
00295           int mmax=re_list.indexmax();
00296           dvector tmp(1,mmax);
00297 
00298           int j;
00299           for (j=1;j<=re_list.indexmax();j++)
00300           {
00301             tmp(j)=uadjoint(re_list(j)-xmax);
00302           }
00303 
00304           if (allocated(fe_list))
00305           {
00306             dvector tmp1=solve(H,tmp);
00307             dvector xtmp=tmp1*Dux;
00308             for (j=1;j<=fe_list.indexmax();j++)
00309             {
00310               x_con(fe_list(j))+=xtmp(j);
00311             }
00312           }
00313         }
00314       }
00315       if (initial_df1b2params::separable_flag)
00316       {
00317         //for (i=1;i<=usize;i++)
00318         //{
00319         //  Dux(i)=elem_prod(Dux(i),sscale);
00320         //}
00321         x_con=elem_prod(x_con,sscale);
00322       }
00323     }
00324     xadjoint-=x_con;
00325     // *******************************************************
00326     // *******************************************************
00327     // *******************************************************
00328 
00329     block_diagonal_flag=3;
00330     //pfmin->lapprox->xadjoint.initialize();
00331     //pfmin->lapprox->uadjoint.initialize();
00332     pfmin->lapprox->num_separable_calls=0;
00333     pfmin->lapprox->check_local_xadjoint.initialize();
00334     pfmin->lapprox->check_local_xadjoint2.initialize();
00335     pfmin->lapprox->check_local_uadjoint.initialize();
00336     pfmin->lapprox->check_local_uadjoint2.initialize();
00337     //df1b2_gradlist::set_yes_derivatives();
00338     //initial_df1b2params::reset(y,pen);
00339     pfmin->user_function();
00340     dvector lcx=elem_prod(check_local_xadjoint2,sscale);
00341     xadjoint+=lcx;
00342     //df1b2_gradlist::set_no_derivatives();
00343     funnel_init_var::lapprox=0;
00344     block_diagonal_flag=0;
00345     initial_params::set_inactive_only_random_effects();
00346   }
00347   else
00348   {
00349     if (function_minimizer::first_hessian_flag)
00350     {
00351       // need hessin of random effects for stddeV report
00352       df1b2variable pen=0.0;
00353       initial_df1b2params::reset(y,pen);
00354       initial_params::straight_through_flag=0;
00355       block_diagonal_flag=6;
00356       allocate_block_diagonal_stuff();
00357       num_separable_calls=0;
00358       // get the block diagonal hessians to use in importance sampling
00359       pfmin->user_function();
00360       //cout << (*pfmin->lapprox->block_diagonal_hessian) << endl;
00361       block_diagonal_flag=2;
00362     }
00363     initial_params::straight_through_flag=1;
00364     df1b2variable pen=0.0;
00365     if (saddlepointflag==2)
00366     {
00367       pmin->inner_opt_flag=0;
00368       f=get_fx_fu(pfmin);
00369     }
00370     initial_df1b2params::reset(y,pen);
00371     pmin->inner_opt_flag=1;
00372     pfmin->pre_user_function();
00373     pmin->inner_opt_flag=0;
00374     initial_params::straight_through_flag=0;
00375     if (saddlepointflag!=2)
00376     {
00377       f=initial_df1b2params::cobjfun;
00378     }
00379     else
00380     {
00381       xadjoint=(*grad_x_u)(1,xsize)-(*grad_x);
00382     }
00383 
00384     if (saddlepointflag!=2 && pmin->multinomial_weights==0)
00385     {
00386       f-=usize*.91893853320467241;
00387     }
00388 
00389     funnel_init_var::lapprox=0;
00390     block_diagonal_flag=0;
00391     dvector scale1(1,xsize);   // need to get scale from somewhere
00392     initial_params::set_inactive_only_random_effects();
00393     /*int check=*/initial_params::stddev_scale(scale1,x);
00394     for (i=1;i<=xadjoint.indexmax();i++)
00395       xadjoint(i)*=scale1(i);
00396   }
00397   //cout << initial_df1b2params::cobjfun << endl;
00398   //f=initial_df1b2params::cobjfun;
00399   return xadjoint;
00400 }
00401 
00406 dvector laplace_approximation_calculator::get_newton_raphson_info_block_diagonal
00407   (function_minimizer * pfmin)
00408 {
00409   int nv=initial_df1b2params::set_index();
00410   if (allocated(used_flags))
00411   {
00412     if (used_flags.indexmax() != nv)
00413     {
00414       used_flags.safe_deallocate();
00415     }
00416   }
00417   if (!allocated(used_flags))
00418   {
00419     used_flags.safe_allocate(1,nv);
00420   }
00421 
00422   //for (ip=1;ip<=num_der_blocks;ip++)
00423   {
00424     used_flags.initialize();
00425 
00426     // do we need to reallocate memory for df1b2variables?
00427     int ip = 0;
00428     check_for_need_to_reallocate(ip);
00429 
00430     df1b2_gradlist::set_no_derivatives();
00431     //cout << re_objective_function_value::pobjfun << endl;
00432     //cout << re_objective_function_value::pobjfun->ptr << endl;
00433     (*re_objective_function_value::pobjfun)=0;
00434     df1b2variable pen=0.0;
00435     df1b2variable zz=0.0;
00436 
00437     initial_df1b2params::reset(y,pen);
00438 
00439     // call function to do block diagonal newton-raphson
00440     // the step vector from the newton-raphson is in the vector step
00441     df1b2_gradlist::set_yes_derivatives();
00442 
00443     funnel_init_var::lapprox=this;
00444     //cout << funnel_init_var::lapprox << endl;
00445     block_diagonal_flag=1;
00446     pfmin->pre_user_function();
00447     //pfmin->user_function();
00448     funnel_init_var::lapprox=0;
00449     block_diagonal_flag=0;
00450   }
00451   return step;
00452 }