ADMB Documentation  11.1.1650
 All Classes Files Functions Variables Typedefs Friends Defines
df1b2lp5.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: df1b2lp5.cpp 1109 2013-07-12 00:22:26Z johnoel $
00003  *
00004  * Author: David Fournier
00005  * Copyright (c) 2008-2012 Regents of the University of California
00006  */
00012 #if defined(USE_LAPLACE)
00013 #  include <admodel.h>
00014 #  include <df1b2fun.h>
00015 #  include <adrndeff.h>
00016 void get_second_ders(int xs,int us,const init_df1b2vector y,dmatrix& Hess,
00017   dmatrix& Dux, df1b2_gradlist * f1b2gradlist,function_minimizer * pfmin,
00018   laplace_approximation_calculator* lap);
00019 double calculate_laplace_approximation(const dvector& x,const dvector& u0,
00020   const dmatrix& Hess,const dvector& _xadjoint,const dvector& _uadjoint,
00021   const dmatrix& _Hessadjoint,function_minimizer * pmin);
00022 
00023 #if defined(USE_ADPVM)
00024 
00028 void get_second_ders_slave(int xs,int us,const init_df1b2vector _y,
00029   dmatrix& Hess,
00030   dmatrix& Dux, df1b2_gradlist * f1b2gradlist,function_minimizer * pfmin,
00031   laplace_approximation_calculator * lpc)
00032 {
00033   // Note that xs is the number of active non random effects
00034   // parameters
00035   // us is the number of random effects parameters
00036   int j;
00037   int i;
00038   ADUNCONST(init_df1b2vector,y)
00039   int num_der_blocks=lpc->num_der_blocks;
00040   int xsize=lpc->xsize;
00041   int usize=lpc->usize;
00042 
00043   for (int ip=1;ip<=num_der_blocks;ip++)
00044   {
00045     df1b2variable::minder=lpc->minder(ip);
00046     df1b2variable::maxder=lpc->maxder(ip);
00047     lpc->set_u_dot(ip);
00048     df1b2_gradlist::set_yes_derivatives();
00049     (*re_objective_function_value::pobjfun)=0;
00050     df1b2variable pen=0.0;
00051     df1b2variable zz=0.0;
00052     initial_df1b2params::reset(y,pen);
00053     if (initial_df1b2params::separable_flag)
00054     {
00055       initial_df1b2params::separable_calculation_type=2;
00056       Hess.initialize();
00057       Dux.initialize();
00058     }
00059 
00060     //cout << "2D" << endl;
00061     pfmin->user_function();
00062 
00063     //pfmin->user_function(y,zz);
00064     (*re_objective_function_value::pobjfun)+=pen;
00065     (*re_objective_function_value::pobjfun)+=zz;
00066 
00067     if (!initial_df1b2params::separable_flag)
00068     {
00069       set_dependent_variable(*re_objective_function_value::pobjfun);
00070       df1b2_gradlist::set_no_derivatives();
00071       df1b2variable::passnumber=1;
00072       df1b2_gradcalc1();
00073 
00074       int mind=y(1).minder;
00075       int jmin=max(mind,xsize+1);
00076       int jmax=min(y(1).maxder,xsize+usize);
00077       dmatrix tHess(1,usize,jmin,jmax);
00078       for (i=1;i<=usize;i++)
00079         for (j=jmin;j<=jmax;j++)
00080         {
00081           //Hess(i,j-xsize)=y(i+xsize).u_bar[j-mind];
00082           tHess(i,j)=y(i+xsize).u_bar[j-mind];
00083         }
00084       send_dmatrix_to_master(tHess);
00085 
00086       jmin=max(mind,1);
00087       jmax=min(y(1).maxder,xsize);
00088       if (jmax>=jmin)
00089       {
00090         dmatrix tDux(1,usize,jmin,jmax);
00091         for (i=1;i<=usize;i++)
00092           for (j=jmin;j<=jmax;j++)
00093           {
00094             //Dux(i,j)=y(i+xsize).u_bar[j-1];
00095             tDux(i,j)=y(i+xsize).u_bar[j-1];
00096           }
00097 
00098         send_int_to_master(1);
00099         send_dmatrix_to_master(tDux);
00100       }
00101       else
00102       {
00103         send_int_to_master(0);
00104       }
00105     }
00106     if (ip<num_der_blocks)
00107       f1b2gradlist->reset();
00108   }
00109   //cout << "Hess" << endl;
00110   //cout << Hess << endl;
00111   //cout << "Dux" << endl;
00112   //cout << setprecision(16) << Dux << endl;
00113 }
00114 
00119 void laplace_approximation_calculator::default_calculations_parallel_slave
00120   (const dvector& _x,const double& _f,function_minimizer * pfmin)
00121 {
00122   // for use when there is no separability
00123   ADUNCONST(dvector,x)
00124   ADUNCONST(double,f)
00125   int i,j;
00126   x=get_dvector_from_master();
00127 
00128   initial_params::set_inactive_only_random_effects();
00129   gradient_structure::set_NO_DERIVATIVES();
00130   initial_params::reset(x);    // get current x values into the model
00131   gradient_structure::set_YES_DERIVATIVES();
00132 
00133   initial_params::set_active_only_random_effects();
00134   int lmn_flag=0;
00135   if (!inner_lmnflag)
00136   {
00137     if (!ADqd_flag)
00138       uhat=get_uhat_quasi_newton(x,pfmin);
00139     else
00140       uhat=get_uhat_quasi_newton_qd(x,pfmin);
00141   }
00142   else
00143   {
00144     uhat=get_uhat_lm_newton(x,pfmin);
00145   }
00146 
00147   for (i=1;i<=xsize;i++)
00148   {
00149     y(i)=x(i);
00150   }
00151   for (i=1;i<=usize;i++)
00152   {
00153     y(i+xsize)=uhat(i);
00154   }
00155   //cout << y << endl;
00156 
00157   for(int ii=1;ii<=num_nr_iters;ii++)
00158   {
00159     {
00160       // test newton raphson
00161       Hess.initialize();
00162      cout << "Newton raphson " << ii << endl;
00163       get_newton_raphson_info_slave(pfmin);
00164 
00165       step=get_dvector_from_master();
00166 
00167       f1b2gradlist->reset();
00168       f1b2gradlist->list.initialize();
00169       f1b2gradlist->list2.initialize();
00170       f1b2gradlist->list3.initialize();
00171       f1b2gradlist->nlist.initialize();
00172       f1b2gradlist->nlist2.initialize();
00173       f1b2gradlist->nlist3.initialize();
00174 
00175       uhat+=step;
00176 
00177       evaluate_function(uhat,pfmin);
00178     }
00179 
00180     for (i=1;i<=usize;i++)
00181     {
00182       y(i+xsize)=uhat(i);
00183     }
00184   }
00185 
00186   get_second_ders_slave(xsize,usize,y,Hess,Dux,f1b2gradlist,pfmin,this);
00187   int sgn=0;
00188 
00189   //f=calculate_laplace_approximation(x,uhat,Hess,xadjoint,uadjoint,
00190    // Hessadjoint,pfmin);
00191 
00192 
00193   for (int ip=num_der_blocks;ip>=1;ip--)
00194   {
00195     df1b2variable::minder=minder(ip);
00196     df1b2variable::maxder=maxder(ip);
00197     int mind=y(1).minder;
00198     int jmin=max(mind,xsize+1);
00199     int jmax=min(y(1).maxder,xsize+usize);
00200     dmatrix Hessadjoint=get_dmatrix_from_master();
00201     //for (i=1;i<=usize;i++)
00202 
00203     if (initial_df1b2params::separable_flag)
00204     {
00205       for (j=1;j<=xsize+usize;j++)
00206       {
00207         *y(j).get_u_tilde()=0;
00208       }
00209       Hess.initialize();
00210       initial_df1b2params::separable_calculation_type=3;
00211       pfmin->user_function();
00212     }
00213     else
00214     {
00215       if (ip<num_der_blocks)
00216       {
00217         f1b2gradlist->reset();
00218         set_u_dot(ip);
00219         df1b2_gradlist::set_yes_derivatives();
00220         (*re_objective_function_value::pobjfun)=0;
00221         df1b2variable pen=0.0;
00222         df1b2variable zz=0.0;
00223 
00224         initial_df1b2params::reset(y,pen);
00225         pfmin->user_function();
00226         (*re_objective_function_value::pobjfun)+=pen;
00227         (*re_objective_function_value::pobjfun)+=zz;
00228 
00229         set_dependent_variable(*re_objective_function_value::pobjfun);
00230         df1b2_gradlist::set_no_derivatives();
00231         df1b2variable::passnumber=1;
00232         df1b2_gradcalc1();
00233       }
00234 
00235       for (i=1;i<=usize;i++)
00236       {
00237         for (j=1;j<=usize;j++)
00238         {
00239           y(i+xsize).get_u_bar_tilde()[j-mind]=0;
00240         }
00241       }
00242 
00243       for (i=1;i<=usize;i++)
00244       {
00245         for (j=jmin;j<=jmax;j++)
00246         {
00247           //Hess(i,j-xsize)=y(i+xsize).u_bar[j-mind];
00248           y(i+xsize).get_u_bar_tilde()[j-mind]=Hessadjoint(i,j-xsize);
00249         }
00250       }
00251 
00252       int mind=y(1).minder;
00253       df1b2variable::passnumber=2;
00254       df1b2_gradcalc1();
00255 
00256       df1b2variable::passnumber=3;
00257       df1b2_gradcalc1();
00258 
00259       f1b2gradlist->reset();
00260       f1b2gradlist->list.initialize();
00261       f1b2gradlist->list2.initialize();
00262       f1b2gradlist->list3.initialize();
00263       f1b2gradlist->nlist.initialize();
00264       f1b2gradlist->nlist2.initialize();
00265       f1b2gradlist->nlist3.initialize();
00266     }
00267 
00268     dvector dtmp(1,xsize);
00269     for (i=1;i<=xsize;i++)
00270     {
00271       dtmp(i)=*y(i).get_u_tilde();
00272     }
00273     if (initial_df1b2params::separable_flag)
00274     {
00275       dvector scale(1,nvar);   // need to get scale from somewhere
00276       int check=initial_params::stddev_scale(scale,x);
00277       dvector sscale=scale(1,Dux(1).indexmax());
00278       for (i=1;i<=usize;i++)
00279       {
00280         Dux(i)=elem_prod(Dux(i),sscale);
00281       }
00282       dtmp=elem_prod(dtmp,sscale);
00283     }
00284 
00285     send_dvector_to_master(dtmp);
00286     dvector utmp(1,usize);
00287     for (i=1;i<=usize;i++)
00288     {
00289       utmp(i)=*y(i+xsize).get_u_tilde();
00290     }
00291     send_dvector_to_master(utmp);
00292   }
00293 }
00294 
00299 void laplace_approximation_calculator::get_newton_raphson_info_slave
00300   (function_minimizer * pfmin)
00301 {
00302   int i,j,ip;
00303 
00304   for (ip=1;ip<=num_der_blocks;ip++)
00305   {
00306     df1b2variable::minder=minder(ip);
00307     df1b2variable::maxder=maxder(ip);
00308 
00309     set_u_dot(ip);
00310 
00311     // do we need to reallocate memory for df1b2variables?
00312     check_for_need_to_reallocate(ip);
00313     df1b2_gradlist::set_yes_derivatives();
00314     //cout << re_objective_function_value::pobjfun << endl;
00315     //cout << re_objective_function_value::pobjfun->ptr << endl;
00316     (*re_objective_function_value::pobjfun)=0;
00317     df1b2variable pen=0.0;
00318     df1b2variable zz=0.0;
00319     initial_df1b2params::reset(y,pen);
00320     if (initial_df1b2params::separable_flag)
00321     {
00322       Hess.initialize();
00323       grad.initialize();
00324     }
00325     pfmin->user_function();
00326     (*re_objective_function_value::pobjfun)+=pen;
00327     (*re_objective_function_value::pobjfun)+=zz;
00328     //for (j=0;j<xsize+usize;j++)
00329      // cout << re_objective_function_value::pobjfun->u_dot[j] << " ";
00330 
00331     if (!initial_df1b2params::separable_flag)
00332     {
00333       set_dependent_variable(*re_objective_function_value::pobjfun);
00334       df1b2_gradlist::set_no_derivatives();
00335       df1b2variable::passnumber=1;
00336       df1b2_gradcalc1();
00337       int mind=y(1).minder;
00338       int jmin=max(mind,xsize+1);
00339       int jmax=min(y(1).maxder,xsize+usize);
00340       dmatrix tHess(1,usize,jmin,jmax);
00341       dvector tgrad(jmin,jmax);
00342       for (i=1;i<=usize;i++)
00343         for (j=jmin;j<=jmax;j++)
00344         {
00345           //Hess(i,j-xsize)=y(i+xsize).u_bar[j-mind];
00346           tHess(i,j)=y(i+xsize).u_bar[j-mind];
00347         }
00348 
00349       for (j=jmin;j<=jmax;j++)
00350       {
00351         //grad(j-xsize)= re_objective_function_value::pobjfun->u_dot[j-mind];
00352         tgrad(j)= re_objective_function_value::pobjfun->u_dot[j-mind];
00353       }
00354       send_dmatrix_to_master(tHess);
00355       send_dvector_to_master(tgrad);
00356     }
00357     if (ip<num_der_blocks)
00358       f1b2gradlist->reset();
00359   }
00360   {
00361     //for (j=0;j<xsize+usize;j++)
00362     //  cout << re_objective_function_value::pobjfun->u_dot[j] << " ";
00363     //cout << endl;
00364     //ofstream ofs("hess.tmp");
00365     //cout << Hess << endl;
00366     //cout << grad << endl;
00367   }
00368 }
00369 #else
00370 void laplace_approximation_calculator::default_calculations_parallel_slave
00371   (const dvector& _x,const double& _f,function_minimizer * pfmin){}
00372 void laplace_approximation_calculator::get_newton_raphson_info_slave
00373   (function_minimizer * pfmin){}
00374 #endif  // #if defined(USE_ADPVM)
00375 
00376 #endif //#if defined(USE_LAPLACE)