ADMB Documentation  11.6rc.3396
 All Classes Files Functions Variables Typedefs Friends Defines
df1b2lap.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id$
00003  *
00004  * Author: David Fournier
00005  * Copyright (c) 2008-2012 Regents of the University of California
00006  */
00011 #include <sstream>
00012 using std::istringstream;
00013 #include <admodel.h>
00014 #include <df1b2fun.h>
00015 #include <adrndeff.h>
00016 
00017 #ifndef OPT_LIB
00018   #include <cassert>
00019 #endif
00020         int fcount =0;
00021 static int no_stuff=0;
00022 static int write_sparse_flag=0;
00023     static void trapper(void)
00024     {
00025       //int x=5;
00026     }
00027 int noboundepen_flag=1;
00028 unsigned int global_nvar=0;
00029 
00030 double CHECK_HESSIAN_PENALTY=0.0;
00031 
00032 double evaluate_function(const dvector& x,function_minimizer * pfmin);
00033 void get_newton_raphson_info(int xs,int us,const init_df1b2vector _y,
00034   dmatrix& Hess, dvector& grad,
00035   df1b2_gradlist* f1b2gradlist,function_minimizer * pfmin);
00036 
00037 //dvariable AD_uf_inner(const dvector& x,const dvar_vector& u);
00038 void get_second_ders(int xs,int us,const init_df1b2vector y,dmatrix& Hess,
00039   dmatrix& Dux, df1b2_gradlist * f1b2gradlist,function_minimizer * pfmin,
00040   laplace_approximation_calculator* lap);
00041 
00042   double re_objective_function_value::fun_without_pen=0;
00043 
00044 int laplace_approximation_calculator::antiflag=0;
00045 int laplace_approximation_calculator::alternative_user_function_flag=0;
00046 int laplace_approximation_calculator::saddlepointflag=0;
00047 int laplace_approximation_calculator::print_importance_sampling_weights_flag=0;
00048 int laplace_approximation_calculator::sparse_hessian_flag=0;
00049 
00050 int laplace_approximation_calculator::where_are_we_flag=0;
00051 dvar_vector *
00052   laplace_approximation_calculator::variance_components_vector=0;
00053 
00058 dvector laplace_approximation_calculator::get_uhat_quasi_newton
00059   (const dvector& x,function_minimizer * pfmin)
00060 {
00061   //int on,nopt;
00062   pfmin->inner_opt_flag=1;
00063   double f=0.0;
00064   double fb=1.e+100;
00065   dvector g(1,usize);
00066   dvector ub(1,usize);
00067   independent_variables u(1,usize);
00068   gradcalc(0,g);
00069   fmc1.itn=0;
00070   fmc1.ifn=0;
00071   fmc1.ireturn=0;
00072   initial_params::xinit(u);    // get the initial values into the
00073   initial_params::xinit(ubest);    // get the initial values into the
00074   fmc1.ialph=0;
00075   fmc1.ihang=0;
00076   fmc1.ihflag=0;
00077   fmc1.use_control_c=0;
00078 
00079   if (init_switch)
00080   {
00081     u.initialize();
00082   }
00083   else
00084   {
00085     u=ubest;
00086   }
00087 
00088   fmc1.dfn=1.e-2;
00089   dvariable pen=0.0;
00090   //cout << "starting  norm(u) = " << norm(u) << endl;
00091   while (fmc1.ireturn>=0)
00092   {
00093    /*
00094     double nu=norm(value(u));
00095     if (nu>400)
00096     {
00097       cout << "U norm(u) = " << nu  << endl;
00098     }
00099     cout << "V norm(u) = " << nu
00100          << " f = " << f  << endl;
00101     */
00102     fmc1.fmin(f,u,g);
00103     //cout << "W norm(u) = " << norm(value(u)) << endl;
00104     if (fmc1.ireturn>0)
00105     {
00106       dvariable vf=0.0;
00107       pen=initial_params::reset(dvar_vector(u));
00108       *objective_function_value::pobjfun=0.0;
00109       pfmin->AD_uf_inner();
00110       if (saddlepointflag)
00111       {
00112         *objective_function_value::pobjfun*=-1.0;
00113       }
00114       if ( no_stuff==0 && quadratic_prior::get_num_quadratic_prior()>0)
00115       {
00116         quadratic_prior::get_M_calculations();
00117       }
00118       vf+=*objective_function_value::pobjfun;
00119 
00120      /*  this is now done in the operator = function
00121       if (quadratic_prior::get_num_quadratic_prior()>0)
00122       {
00123         vf+= quadratic_prior::get_quadratic_priors();
00124       }
00125       */
00126 
00127       objective_function_value::fun_without_pen=value(vf);
00128 
00129       //cout << " pen = " << pen << endl;
00130       if (noboundepen_flag==0)
00131       {
00132         vf+=pen;
00133       }
00134       f=value(vf);
00135       if (f<fb)
00136       {
00137         fb=f;
00138         ub=u;
00139       }
00140       gradcalc(usize,g);
00141       //cout << " f = " << setprecision(17) << f << " " << norm(g)
00142        // << " " << norm(u) << endl;
00143     }
00144     u=ub;
00145   }
00146   cout <<  " inner maxg = " <<  fmc1.gmax;
00147   if (fabs(fmc1.gmax)>1.e+3)
00148     trapper();
00149 
00150   if (fabs(fmc1.gmax)>1.e-4)
00151   {
00152     fmc1.itn=0;
00153     //fmc1.crit=1.e-9;
00154     fmc1.ifn=0;
00155     fmc1.ireturn=0;
00156     fmc1.ihang=0;
00157     fmc1.ihflag=0;
00158     fmc1.ialph=0;
00159     initial_params::xinit(u);    // get the initial values into the
00160     //u.initialize();
00161     while (fmc1.ireturn>=0)
00162     {
00163       fmc1.fmin(f,u,g);
00164       if (fmc1.ireturn>0)
00165       {
00166         dvariable vf=0.0;
00167         pen=initial_params::reset(dvar_vector(u));
00168         *objective_function_value::pobjfun=0.0;
00169         pfmin->AD_uf_inner();
00170         if ( no_stuff==0 && quadratic_prior::get_num_quadratic_prior()>0)
00171         {
00172           quadratic_prior::get_M_calculations();
00173         }
00174         vf+=*objective_function_value::pobjfun;
00175         objective_function_value::fun_without_pen=value(vf);
00176 
00177         vf+=pen;
00178         f=value(vf);
00179         if (f<fb)
00180         {
00181           fb=f;
00182           ub=u;
00183         }
00184         gradcalc(usize,g);
00185         //cout << " f = " << setprecision(15) << f << " " << norm(g) << endl;
00186       }
00187     }
00188     u=ub;
00189     cout <<  "  Inner second time = " << fmc1.gmax;
00190   }
00191   cout << "  Inner f = " << fb << endl;
00192   fmc1.ireturn=0;
00193   fmc1.fbest=fb;
00194   gradient_structure::set_NO_DERIVATIVES();
00195   *objective_function_value::pobjfun=0.0;
00196   pfmin->AD_uf_inner();
00197   if ( no_stuff==0 && quadratic_prior::get_num_quadratic_prior()>0)
00198   {
00199     quadratic_prior::get_M_calculations();
00200   }
00201   gradient_structure::set_YES_DERIVATIVES();
00202   pfmin->inner_opt_flag=0;
00203   return u;
00204 }
00205 
00210 dvector laplace_approximation_calculator::get_uhat_lm_newton
00211   (const dvector& x,function_minimizer * pfmin)
00212 {
00213   double f=0.0;
00214   //dvector g(1,usize);
00215   //g.initialize();
00216   independent_variables u(1,usize);
00217   if (init_switch)
00218   {
00219     initial_params::xinit(u);    // get the initial values into the
00220   }
00221   else
00222   {
00223     u=ubest;
00224   }
00225 
00226   int iprint_save=pfmin->iprint;
00227   pfmin->iprint=inner_iprint;
00228   pfmin->limited_memory_quasi_newton(f,u,5,1,inner_maxfn,inner_crit);
00229   pfmin->iprint=iprint_save;
00230   if (initial_params::mc_phase==0)
00231   {
00232     cout << "  Inner f = " << f << endl;
00233   }
00234   //gradient_structure::set_NO_DERIVATIVES();
00235   //initial_params::set_inactive_only_random_effects();
00236   //initial_params::reset(dvar_vector(x));
00237   return u;
00238 }
00239 
00244 void set_partition_sizes(int & num_der_blocks,ivector& minder,
00245   ivector& maxder,int nvariables)
00246 {
00247   if (num_der_blocks==0)
00248   {
00249     num_der_blocks=1;
00250   }
00251 #if defined(USE_ADPVM)
00252   if (ad_comm::pvm_manager)
00253   {
00254     switch (ad_comm::pvm_manager->mode)
00255     {
00256     case 1: //master
00257       {
00258         int i;
00259         int ndb=ad_comm::pvm_manager->num_slave_processes+1;
00260 
00261         minder.allocate(1,ndb);
00262         maxder.allocate(1,ndb);
00263 
00264         int nd=nvariables/ndb;
00265         int r= nvariables - nd * ndb;
00266         ivector partition(1,ndb);
00267         partition=nd;
00268         partition(1,r)+=1;
00269         minder(1)=1;
00270         maxder(1)=partition(1);
00271 
00272         for (i=2;i<=ndb;i++)
00273         {
00274           minder(i)=maxder(i-1)+1;
00275           maxder(i)=minder(i)+partition(i)-1;
00276         }
00277         send_int_to_slaves(minder(2,ndb));
00278         send_int_to_slaves(maxder(2,ndb));
00279         break;
00280       }
00281 
00282     case 2: //slave
00283       minder.allocate(1,1);
00284       maxder.allocate(1,1);
00285       minder(1)=get_int_from_master();
00286       maxder(1)=get_int_from_master();
00287       break;
00288     }
00289   }
00290   else
00291 #endif // #USE_ADPVM
00292   {
00293     minder.allocate(1,num_der_blocks);
00294     maxder.allocate(1,num_der_blocks);
00295 
00296     int nd=nvariables/num_der_blocks;
00297     int r= nvariables - nd * num_der_blocks;
00298     ivector partition(1,num_der_blocks);
00299     partition=nd;
00300     partition(1,r)+=1;
00301 
00302     minder(1)=1;
00303     maxder(1)=partition(1);
00304     for (int i=2;i<=num_der_blocks;i++)
00305     {
00306       minder(i)=maxder(i-1)+1;
00307       maxder(i)=minder(i)+partition(i)-1;
00308     }
00309   }
00310 }
00311 
00316 laplace_approximation_calculator::laplace_approximation_calculator
00317 (
00318   int _xsize,
00319   int _usize,
00320   int _minder,
00321   int _maxder,
00322   function_minimizer* _pmin
00323 ):
00324   init_switch(1),
00325   separable_call_level(0),
00326   triplet_information(0),
00327   compressed_triplet_information(0),
00328   nested_separable_calls_counter(1,10),
00329   nested_tree_position(1,5),
00330   local_dtemp(1,_xsize),
00331   nr_debug(0),
00332   pmin(_pmin),
00333   block_diagonal_flag(0),
00334   xsize(_xsize),
00335   usize(_usize),
00336   bHess_pd_flag(0),
00337   sparse_triplet(0),
00338   sparse_iterator(0),
00339   sparse_count(0),
00340   sparse_count_adjoint(0),
00341   sparse_triplet2(0),
00342   vsparse_triplet(0),
00343   vsparse_triplet_adjoint(0),
00344   sparse_symbolic(0),
00345   sparse_symbolic2(0),
00346   fmc1(usize,1),
00347   fmc(_xsize),
00348   xadjoint(1,_xsize),
00349   check_local_uadjoint(1,_usize),
00350   check_local_uadjoint2(1,_usize),
00351   check_local_xadjoint(1,_xsize),
00352   check_local_xadjoint2(1,_xsize),
00353   uadjoint(1,_usize),
00354   uhat(1,_usize),
00355   ubest(1,_usize)
00356 {
00357   ubest.initialize();
00358   nested_shape.allocate(100,100,10);
00359   nested_shape.initialize();
00360   nested_tree_position.initialize();
00361   nested_separable_calls_counter.initialize();
00362   calling_set=0;
00363   antiepsilon=0;
00364   dd_nr_flag=0;
00365   no_re_ders_flag =0;
00366   importance_sampling_components=0;
00367   is_diagnostics_flag=0;
00368   importance_sampling_values = 0;
00369   importance_sampling_weights = 0;
00370   no_function_component_flag=0;
00371   uhat.initialize();
00372   hesstype=1;
00373   use_gauss_hermite=0;
00374   in_gauss_hermite_phase=0;
00375   multi_random_effects=0;
00376   //var_flag=1;
00377   var_flag=0;
00378   num_separable_calls=0;
00379   separable_calls_counter=0;
00380   importance_sampling_counter=0;
00381   num_local_re_array=0;
00382   num_local_fixed_array=0;
00383   isfunnel_flag=0;
00384   antiflag=0;
00385   rseed=3457;
00386   nfunnelblocks=0;
00387   separable_function_difference=0;
00388   gh=0;
00389   block_diagonal_hessian=0;
00390   block_diagonal_Dux=0;
00391   block_diagonal_re_list=0;
00392   block_diagonal_fe_list=0;
00393   block_diagonal_vhessian=0;
00394   block_diagonal_vhessianadjoint=0;
00395   block_diagonal_ch=0;
00396   block_diagonal_vch=0;
00397   have_users_hesstype=0;
00398   pHess_non_quadprior_part=0;
00399   bw=0;
00400   bHess=0;
00401   grad_x_u=0;
00402   grad_x=0;
00403   int ndi=20000;
00404   int nopt=0;
00405   int on=-1;
00406   if ( (inner_lmnflag=option_match(ad_comm::argc,ad_comm::argv,"-ndi",nopt))
00407     >-1)
00408   {
00409     if (!nopt)
00410     {
00411       cerr << "Usage -ndi option needs integer  -- set to default 20000.\n";
00412     }
00413     else
00414     {
00415       int jj=atoi(ad_comm::argv[inner_lmnflag+1]);
00416       if (jj<=0)
00417       {
00418         cerr << "Usage -ndi option needs positive integer"
00419         "  -- set to defalt 20000" << endl;
00420       }
00421       else
00422       {
00423         ndi=jj;
00424       }
00425     }
00426   }
00427   if ( (on=option_match(ad_comm::argc,ad_comm::argv,"-noreders",nopt))>-1)
00428   {
00429     no_re_ders_flag=1;
00430   }
00431   derindex=new imatrix(1,ndi);
00432   bHessadjoint=0;
00433   if (initial_df1b2params::separable_flag)
00434   {
00435     if (initial_df1b2params::pointer_table)
00436     {
00437       delete initial_df1b2params::pointer_table;
00438     }
00439     initial_df1b2params::pointer_table=0;
00440   }
00441   inner_lmnflag=0;
00442   inner_lmnsteps=10;
00443   if ( (on=option_match(ad_comm::argc,ad_comm::argv,"-noinit",nopt))>-1)
00444   {
00445     init_switch=0;
00446   }
00447   if ( (inner_lmnflag=option_match(ad_comm::argc,ad_comm::argv,"-ilmn",nopt))
00448     >-1)
00449   {
00450     if (!nopt)
00451     {
00452       cerr << "Usage -ilmn option needs integer  -- set to default 10" << endl;
00453     }
00454     else
00455     {
00456       int jj=atoi(ad_comm::argv[inner_lmnflag+1]);
00457       if (jj<=0)
00458       {
00459         cerr << "Usage -ilmn option needs positive integer"
00460           "  -- set to defalt 10" << endl;
00461       }
00462       else
00463       {
00464         inner_lmnsteps=jj;
00465       }
00466     }
00467   }
00468   else
00469   {
00470     inner_lmnflag=0;
00471   }
00472   inner_noprintx=0;
00473 
00474   num_der_blocks=1; // default value
00475   if ( (on=option_match(ad_comm::argc,ad_comm::argv,"-ndb",nopt))>-1)
00476   {
00477     if (!nopt)
00478     {
00479       cerr << "Usage -ndb option needs non-negative integer  -- ignored.\n";
00480     }
00481     else
00482     {
00483       int _num_der_blocks=atoi(ad_comm::argv[on+1]);
00484       if (_num_der_blocks<=0)
00485       {
00486         cerr << "Usage -ndb option needs positive integer  -- ignored" << endl;
00487       }
00488       else
00489       {
00490         num_der_blocks=_num_der_blocks;
00491       }
00492     }
00493   }
00494 
00495   if ( (on=option_match(ad_comm::argc,ad_comm::argv,"-isf",nopt))>-1)
00496   {
00497     if (!nopt)
00498     {
00499       cerr << "Usage -isf option needs non-negative integer  -- ignored.\n";
00500     }
00501     else
00502     {
00503       int _nfunnelblocks=atoi(ad_comm::argv[on+1]);
00504       if (_nfunnelblocks<=0)
00505       {
00506         cerr << "Usage -isf option needs positive integer  -- ignored" << endl;
00507       }
00508       else
00509       {
00510         nfunnelblocks=_nfunnelblocks;
00511         isfunnel_flag=1;
00512       }
00513     }
00514   }
00515 
00516   antiflag=0;
00517   if ( (on=option_match(ad_comm::argc,ad_comm::argv,"-anti",nopt))>-1)
00518   {
00519     antiflag=1;
00520   }
00521 
00522   if ( (on=option_match(ad_comm::argc,ad_comm::argv,"-nrdbg",nopt))>-1)
00523   {
00524     nr_debug=1;
00525   }
00526 
00527   if ( (on=option_match(ad_comm::argc,ad_comm::argv,"-ddnr",nopt))>-1)
00528   {
00529     dd_nr_flag=1;
00530   }
00531 
00532   nvariables=xsize+usize;
00533   /*
00534   int rem=0;
00535   if (nvariables%xsize!=0)
00536     rem=1;
00537   */
00538 
00539   if (nvariables/num_der_blocks<xsize)
00540   {
00541     cerr << " Error -- the number of der blocks (-ndb) can not be larger than "
00542          << nvariables/xsize << endl;
00543     ad_exit(1);
00544   }
00545 
00546   set_partition_sizes(num_der_blocks,minder,maxder,nvariables);
00547 
00548  // !!! remove
00549   //maxder(1)=5;
00550 
00551   fmc1.crit=1.e-3;
00552 
00553   //cout << "Need to fix Hess" << endl;
00554   for (int i=1;i<=num_der_blocks;i++)
00555   {
00556     if (minder(i)<1 || maxder(i) > nvariables || maxder(i) < minder(i))
00557     {
00558       cerr << " minder or maxder value out of bounds in"
00559         "laplace_approximation_calculator constructor "
00560        << endl << " values are minder = " << minder
00561        << " maxder = " << maxder << endl;
00562       ad_exit(1);
00563     }
00564   }
00565 
00566   num_nr_iters=2;
00567   if ( (on=option_match(ad_comm::argc,ad_comm::argv,"-nr",nopt))>-1)
00568   {
00569     if (!nopt)
00570     {
00571       cerr << "Usage -nr option needs non-negative integer  -- ignored" << endl;
00572     }
00573     else
00574     {
00575       int _num_nr_iters=atoi(ad_comm::argv[on+1]);
00576       if (_num_nr_iters<0)
00577       {
00578         cerr << "Usage -nr option needs non-negative integer  -- ignored.\n";
00579       }
00580       else
00581       {
00582         num_nr_iters=_num_nr_iters;
00583       }
00584     }
00585   }
00586 
00587   if ( (on=option_match(ad_comm::argc,ad_comm::argv,"-nochol",nopt))>-1)
00588   {
00589     ad_comm::no_ln_det_choleski_flag=1;
00590   }
00591 
00592   ad_comm::print_hess_and_exit_flag=0;
00593   if ( (on=option_match(ad_comm::argc,ad_comm::argv,"-phe",nopt))>-1)
00594   {
00595     ad_comm::print_hess_and_exit_flag=1;
00596   }
00597 
00598   double _nr_crit=1.e-11;
00599   nr_crit=1.e-11;
00600   if ( (on=option_match(ad_comm::argc,ad_comm::argv,"-nrcrit",nopt))>-1)
00601   {
00602     if (!nopt)
00603     {
00604       cerr << "Usage -nrcrit option needs number  -- ignored" << endl;
00605     }
00606     else
00607     {
00608       istringstream ist(ad_comm::argv[on+1]);
00609       ist >> _nr_crit;
00610 
00611       if (_nr_crit<=0)
00612       {
00613         cerr << "Usage -nrcrit option needs positive number  -- ignored.\n";
00614         _nr_crit=0.0;
00615       }
00616     }
00617     if (_nr_crit>0)
00618     {
00619       nr_crit=_nr_crit;
00620     }
00621   }
00622   if ( (on=option_match(ad_comm::argc,ad_comm::argv,"-gh",nopt))>-1)
00623   {
00624     if (!nopt)
00625     {
00626       cerr << "Usage -gh option needs positive integer  -- ignored" << endl;
00627     }
00628     else
00629     {
00630       int _inner_gh=atoi(ad_comm::argv[on+1]);
00631       if (_inner_gh<=0)
00632       {
00633         cerr << "Usage -gh option needs positive integer  -- ignored" << endl;
00634       }
00635       else
00636       {
00637         use_gauss_hermite=_inner_gh;
00638       }
00639     }
00640   }
00641 
00642   inner_crit=1.e-3;
00643   if ( (on=option_match(ad_comm::argc,ad_comm::argv,"-icrit",nopt))>-1)
00644   {
00645     double _inner_crit=0.0;;
00646     if (!nopt)
00647     {
00648       cerr << "Usage -icrit option needs number  -- ignored" << endl;
00649     }
00650     else
00651     {
00652       istringstream ist(ad_comm::argv[on+1]);
00653       ist >> _inner_crit;
00654 
00655       if (_inner_crit<=0)
00656       {
00657         cerr << "Usage -icrit option needs positive number  -- ignored" << endl;
00658         _inner_crit=0.0;
00659       }
00660     }
00661     if (_inner_crit>0)
00662     {
00663       inner_crit=_inner_crit;
00664     }
00665   }
00666   fmc1.crit=inner_crit;
00667 
00668   fmc.dfn=.01;
00669 
00670   inner_iprint=0;
00671   if ( (on=option_match(ad_comm::argc,ad_comm::argv,"-iiprint",nopt))>-1)
00672   {
00673     if (!nopt)
00674     {
00675       cerr << "Usage -iprint option needs non-negative integer  -- ignored.\n";
00676     }
00677     else
00678     {
00679       int _inner_iprint=atoi(ad_comm::argv[on+1]);
00680       if (_inner_iprint<=0)
00681       {
00682         cerr << "Usage -iip option needs non-negative integer  -- ignored.\n";
00683       }
00684       else
00685       {
00686         inner_iprint=_inner_iprint;
00687       }
00688     }
00689   }
00690   fmc1.iprint=inner_iprint;
00691 
00692   inner_maxfn=1000;
00693   if ( (on=option_match(ad_comm::argc,ad_comm::argv,"-imaxfn",nopt))>-1)
00694   {
00695     if (!nopt)
00696     {
00697       cerr << "Usage -maxfn option needs non-negative integer  -- ignored.\n";
00698     }
00699     else
00700     {
00701       int _inner_maxfn=atoi(ad_comm::argv[on+1]);
00702       if (_inner_maxfn<0)
00703       {
00704         cerr << "Usage -iip option needs non-negative integer  -- ignored.\n";
00705       }
00706       else
00707       {
00708         inner_maxfn=_inner_maxfn;
00709       }
00710     }
00711   }
00712   fmc1.maxfn=inner_maxfn;
00713   // what sort of structure on the Hessian do we have
00714   nopt=0;
00715 
00716   rseed=3456;
00717   num_importance_samples=0;
00718   if ( (on=option_match(ad_comm::argc,ad_comm::argv,"-is",nopt))>-1)
00719   {
00720     if (!nopt)
00721     {
00722       cerr << "Usage -is option needs positive integer  -- ignored" << endl;
00723     }
00724     else
00725     {
00726       int tht=atoi(ad_comm::argv[on+1]);
00727       if (tht<=0)
00728       {
00729         cerr << "Usage -is option needs non-negative integer  -- ignored.\n";
00730       }
00731       else
00732       {
00733         num_importance_samples=tht;
00734       }
00735       if (nopt==2)
00736       {
00737         int rseed1=atoi(ad_comm::argv[on+2]);
00738         if (rseed1<=0)
00739         {
00740           cerr << "Usage -is option needs non-negative integer  -- ignored.\n";
00741         }
00742         else
00743         {
00744           rseed=rseed1;
00745         }
00746       }
00747     }
00748   }
00749   int use_balanced=0;
00750   if ( (on=option_match(ad_comm::argc,ad_comm::argv,"-isb",nopt))>-1)
00751   {
00752     use_balanced=1;
00753     if (!nopt)
00754     {
00755       cerr << "Usage -isb option needs positive integer  -- ignored" << endl;
00756     }
00757     else
00758     {
00759       int tht=atoi(ad_comm::argv[on+1]);
00760       if (tht<=0)
00761       {
00762         cerr << "Usage -isb option needs non-negative integer  -- ignored.\n";
00763       }
00764       else
00765       {
00766         num_importance_samples=2*tht;
00767       }
00768       if (nopt==2)
00769       {
00770         int rseed1=atoi(ad_comm::argv[on+2]);
00771         if (rseed1<=0)
00772         {
00773           cerr << "Usage -isb option needs non-negative integer  -- ignored.\n";
00774         }
00775         else
00776         {
00777           rseed=rseed1;
00778         }
00779       }
00780     }
00781   }
00782   use_outliers=0;
00783   if ( (on=option_match(ad_comm::argc,ad_comm::argv,"-iso",nopt))>-1)
00784   {
00785     use_outliers=1;
00786   }
00787   if (num_importance_samples)
00788   {
00789     if ( (on=option_match(ad_comm::argc,ad_comm::argv,"-isdiag",nopt))>-1)
00790     {
00791       is_diagnostics_flag=1;
00792     }
00793     if (importance_sampling_values)
00794     {
00795       delete importance_sampling_values;
00796       importance_sampling_values=0;
00797     }
00798     importance_sampling_values =
00799       new dvector(1,num_importance_samples);
00800 
00801     if (importance_sampling_weights)
00802     {
00803       delete importance_sampling_weights;
00804       importance_sampling_weights=0;
00805     }
00806     importance_sampling_weights =
00807       new dvector(1,num_importance_samples);
00808 
00809     random_number_generator rng(rseed);
00810     if (allocated(epsilon)) epsilon.deallocate();
00811     if (use_balanced)
00812     {
00813       // check for eveb num_importance samples
00814       if (num_importance_samples%2)
00815         num_importance_samples+=1;
00816     }
00817     epsilon.allocate(1,num_importance_samples,1,usize);
00818     if (use_balanced)
00819     {
00820       int n2=num_importance_samples/2;
00821       epsilon.sub(1,n2).fill_randn(rng);
00822       if (use_outliers)
00823       {
00824         dmatrix os(1,num_importance_samples,1,usize);
00825         os.fill_randu(rng);
00826         for (int i=1;i<=num_importance_samples;i++)
00827         {
00828           for (int j=1;j<=usize;j++)
00829           {
00830             if (os(i,j)<0.05) epsilon(i,j)*=3.0;
00831           }
00832         }
00833       }
00834       for (int i=1;i<=n2;i++)
00835       {
00836         epsilon(i+n2)=-epsilon(i);
00837       }
00838     }
00839     else
00840     {
00841       epsilon.fill_randn(rng);
00842       if (use_outliers)
00843       {
00844         dmatrix os(1,num_importance_samples,1,usize);
00845         os.fill_randu(rng);
00846         for (int i=1;i<=num_importance_samples;i++)
00847         {
00848           for (int j=1;j<=usize;j++)
00849           {
00850             if (os(i,j)<0.05) epsilon(i,j)*=3.0;
00851           }
00852         }
00853       }
00854     }
00855 
00856     double eps_mult=1.0;
00857     if ( (on=option_match(ad_comm::argc,ad_comm::argv,"-epsmult",nopt))>-1)
00858     {
00859       if (!nopt)
00860       {
00861         cerr << "Usage -epsmult option needs number  -- ignored" << endl;
00862       }
00863       else
00864       {
00865         istringstream ist(ad_comm::argv[on+1]);
00866         ist >> eps_mult;
00867 
00868         if (eps_mult<=0.0 || eps_mult>1.0)
00869         {
00870           cerr << "Usage -epsmult option needs positive number between 0 and 1 "
00871               "-- ignored" << endl;
00872           eps_mult=1.0;
00873         }
00874       }
00875       for (int i=1;i<=num_importance_samples;i++)
00876       {
00877         epsilon(i)*=eps_mult;
00878       }
00879     }
00880   }
00881 
00882   nopt=0;
00883   if ( (on=option_match(ad_comm::argc,ad_comm::argv,"-ht",nopt))>-1)
00884   {
00885     if (!nopt)
00886     {
00887       cerr << "Usage -ht option needs positive integer  -- ignored" << endl;
00888       set_default_hessian_type();
00889     }
00890     else
00891     {
00892       int tht=atoi(ad_comm::argv[on+1]);
00893       if (tht<=0)
00894       {
00895         cerr << "Usage -ht option needs non-negative integer  -- ignored.\n";
00896         set_default_hessian_type();
00897       }
00898       else
00899       {
00900         have_users_hesstype=1;
00901         hesstype=tht;
00902       }
00903       if (nopt>=2)
00904       {
00905         int tbw=atoi(ad_comm::argv[on+2]);
00906         if (tbw<=0)
00907         {
00908           cerr << "Usage -ht option needs non-negative bandwidth"
00909              "  -- ignored" << endl;
00910         }
00911         else
00912         {
00913           ad_comm::bandwidth=tbw;
00914         }
00915       }
00916     }
00917   }
00918   else
00919   {
00920     set_default_hessian_type();
00921   }
00922 
00923   // !! need to check nvar calculation
00924 #ifndef OPT_LIB
00925   assert(maxder(1) >= minder(1));
00926 #endif
00927   nvar = (unsigned int)(maxder(1) - minder(1) + 1);
00928 
00929   switch (hesstype)
00930   {
00931   case 0:
00932   case 1:
00933   case 4:
00934     grad.allocate(1,usize);
00935     Hess.allocate(1,usize,1,usize);
00936     Hessadjoint.allocate(1,usize,1,usize);
00937     Dux.allocate(1,usize,1,xsize);
00938     break;
00939   case 3:
00940     {
00941       int bw=2;
00942       if (ad_comm::bandwidth) bw=ad_comm::bandwidth;
00943       bHess=new banded_symmetric_dmatrix(1,usize,bw);
00944       bHessadjoint=new banded_symmetric_dmatrix(1,usize,bw);
00945       grad.allocate(1,usize);
00946       Dux.allocate(1,usize,1,xsize);
00947     }
00948     break;
00949   default:
00950     break;
00951   }
00952 
00953   step.allocate(1,usize);
00954   // !!! nov 12
00955   f1b2gradlist = new df1b2_gradlist(4000000U,200000U,8000000U,400000U,
00956     2000000U,100000U,adstring("f1b2list1"));
00957 
00958   if (hesstype==2)
00959   {
00960     ad_dstar::allocate(nvar);
00961     global_nvar=nvar;
00962     df1b2variable::set_nvar(nvar);
00963     df1b2variable::set_minder(minder(1));
00964     df1b2variable::set_maxder(maxder(1));
00965     df1b2variable::set_blocksize();
00966   }
00967   if (hesstype!=2)
00968   {
00969     if (sparse_hessian_flag==0)
00970     {
00971       ad_dstar::allocate(nvar);
00972       global_nvar=nvar;
00973       df1b2variable::set_nvar(nvar);
00974       df1b2variable::set_minder(minder(1));
00975       df1b2variable::set_maxder(maxder(1));
00976       df1b2variable::set_blocksize();
00977       y.allocate(1,nvariables);
00978     }
00979     else
00980     {
00981       unsigned int nsave=nvar;
00982       nvar=1;
00983       ad_dstar::allocate(nvar);
00984       global_nvar=nvar;
00985       df1b2variable::set_nvar(1);
00986       df1b2variable::set_minder(minder(1));
00987       df1b2variable::set_maxder(maxder(1));
00988       df1b2variable::set_blocksize();
00989       y.allocate(1,nvariables);
00990       nvar=nsave;
00991       global_nvar=nvar;
00992     }
00993   }
00994   if (df1b2variable::adpool_counter !=0)
00995   {
00996     cout << "this can't happen" << endl;
00997     ad_exit(1);
00998   }
00999   df1b2variable::adpool_vector[df1b2variable::adpool_counter]=
01000     df1b2variable::pool;
01001 #ifndef OPT_LIB
01002   assert(nvariables >= 0);
01003 #endif
01004   df1b2variable::nvar_vector[df1b2variable::adpool_counter]=
01005     (unsigned int)nvariables;
01006   //df1b2variable::adpool_counter++;
01007   df1b2variable::increment_adpool_counter();
01008 }
01009 
01014 laplace_approximation_calculator::laplace_approximation_calculator(
01015   int _xsize,
01016   int _usize,
01017   ivector _minder,
01018   ivector _maxder,
01019   function_minimizer* _pmin
01020 ):
01021   separable_call_level(1),
01022   triplet_information(0),
01023   compressed_triplet_information(0),
01024   nested_separable_calls_counter(1,10),
01025   nested_tree_position(1,5),
01026   nr_debug(0),
01027   pmin(_pmin),
01028   xsize(_xsize),
01029   usize(_usize),
01030   bHess_pd_flag(0),
01031   sparse_triplet(0),
01032   sparse_iterator(0),
01033   sparse_count(0),
01034   sparse_count_adjoint(0),
01035   sparse_triplet2(0),
01036   vsparse_triplet(0),
01037   vsparse_triplet_adjoint(0),
01038   sparse_symbolic(0),
01039   sparse_symbolic2(0),
01040   fmc1(0),
01041   fmc(_xsize),
01042   xadjoint(1,_xsize),
01043   check_local_uadjoint(1,_usize),
01044   check_local_uadjoint2(1,_usize),
01045   check_local_xadjoint(1,_xsize),
01046   check_local_xadjoint2(1,_xsize),
01047   uadjoint(1,_usize),
01048   uhat(1,_usize)
01049 {
01050   nested_tree_position.initialize();
01051   nested_separable_calls_counter.initialize();
01052   //hesstype=1;
01053   uhat.initialize();
01054   //var_flag=1;
01055   var_flag=0;
01056   calling_set=0;
01057   antiepsilon=0;
01058   num_separable_calls=0;
01059   num_local_re_array=0;
01060   num_local_fixed_array=0;
01061   isfunnel_flag=0;
01062   antiflag=0;
01063   rseed=3457;
01064   nfunnelblocks=0;
01065   block_diagonal_hessian=0;
01066   block_diagonal_Dux=0;
01067   block_diagonal_re_list=0;
01068   block_diagonal_fe_list=0;
01069   block_diagonal_vhessian=0;
01070   block_diagonal_vhessianadjoint=0;
01071   pHess_non_quadprior_part=0;
01072   bw=0;
01073   bHess=0;
01074   grad_x_u=0;
01075   grad_x=0;
01076   have_users_hesstype=0;
01077   int mmin=_minder.indexmin();
01078   int mmax=_minder.indexmax();
01079   num_der_blocks= mmax-mmin+1;
01080   minder.allocate(mmin,mmax);
01081   maxder.allocate(mmin,mmax);
01082   minder=_minder;
01083   maxder=_maxder;
01084   inner_iprint = 0;
01085   fmc1.iprint = inner_iprint;
01086   fmc1.crit=1.e-3;
01087   nvariables=xsize+usize;
01088   Hess.allocate(1,usize,1,usize);
01089   for (int i=1;i<=num_der_blocks;i++)
01090   {
01091     if (minder(i)<1 || maxder(i) > nvariables || maxder(i) < minder(i))
01092     {
01093       cerr << " minder or maxder value out of bounds in"
01094         "laplace_approximation_calculator constructor "
01095        << endl << " values are minder = " << minder
01096        << " maxder = " << maxder << endl;
01097       ad_exit(1);
01098     }
01099   }
01100 #ifndef OPT_LIB
01101   assert(maxder(1) >= minder(1));
01102 #endif
01103   nvar = (unsigned int)(maxder(1) - minder(1) + 1);
01104   Hessadjoint.allocate(1,usize,1,usize);
01105   Dux.allocate(1,usize,1,xsize);
01106   // !!! nov 12
01107   f1b2gradlist = new df1b2_gradlist(4000000U,200000U,8000000U,400000U,
01108     2000000U,100000U,adstring("f1b2list1"));
01109   ad_dstar::allocate(nvar);
01110   global_nvar=nvar;
01111   df1b2variable::set_nvar(nvar);
01112   df1b2variable::set_minder(minder(1));
01113   df1b2variable::set_maxder(maxder(1));
01114   df1b2variable::set_blocksize();
01115   y.allocate(1,nvariables);
01116 }
01117 
01122 laplace_approximation_calculator::~laplace_approximation_calculator()
01123 {
01124   if(importance_sampling_weights)
01125   {
01126     delete importance_sampling_weights;
01127     importance_sampling_weights = 0;
01128   }
01129   if(importance_sampling_components)
01130   {
01131     delete importance_sampling_components;
01132     importance_sampling_components = 0;
01133   }
01134   if(importance_sampling_values)
01135   {
01136     delete importance_sampling_values;
01137     importance_sampling_values = 0;
01138   }
01139   if (block_diagonal_vch)
01140   {
01141     delete block_diagonal_vch;
01142     block_diagonal_vch=0;
01143   }
01144   if (block_diagonal_ch)
01145   {
01146     delete block_diagonal_ch;
01147     block_diagonal_ch=0;
01148   }
01149   if (block_diagonal_hessian)
01150   {
01151     delete block_diagonal_hessian;
01152     block_diagonal_hessian=0;
01153   }
01154   if (block_diagonal_Dux )
01155   {
01156     delete block_diagonal_Dux;
01157     block_diagonal_Dux =0;
01158   }
01159   if (block_diagonal_re_list )
01160   {
01161     delete block_diagonal_re_list;
01162     block_diagonal_re_list =0;
01163   }
01164   if (block_diagonal_fe_list )
01165   {
01166     delete block_diagonal_fe_list;
01167     block_diagonal_fe_list =0;
01168   }
01169   if (block_diagonal_vhessian )
01170   {
01171     delete block_diagonal_vhessian;
01172     block_diagonal_vhessian =0;
01173   }
01174   if (block_diagonal_vhessianadjoint )
01175   {
01176     delete block_diagonal_vhessianadjoint;
01177     block_diagonal_vhessianadjoint =0;
01178   }
01179   if (separable_function_difference)
01180   {
01181     delete separable_function_difference;
01182     separable_function_difference=0;
01183   }
01184 
01185   if (derindex)
01186   {
01187     delete derindex;
01188     derindex=0;
01189   }
01190 
01191   if (pHess_non_quadprior_part)
01192   {
01193     delete pHess_non_quadprior_part;
01194     pHess_non_quadprior_part=0;
01195   }
01196 
01197   if (triplet_information)
01198   {
01199     delete triplet_information;
01200     triplet_information=0;
01201   }
01202 
01203   if (bHessadjoint)
01204   {
01205     delete bHessadjoint;
01206     bHessadjoint=0;
01207   }
01208   if (grad_x)
01209   {
01210     delete grad_x;
01211     grad_x=0;
01212   }
01213   if (grad_x_u)
01214   {
01215     delete grad_x_u;
01216     grad_x_u=0;
01217   }
01218   if (bHess)
01219   {
01220     delete bHess;
01221     bHess=0;
01222   }
01223   if (f1b2gradlist)
01224   {
01225     df1b2_gradlist::set_no_derivatives();
01226     delete f1b2gradlist;
01227     f1b2gradlist=0;
01228   }
01229   if (df1b2variable::pool)
01230   {
01231     //df1b2variable::pool->set_size(-1);
01232   }
01233   if (vsparse_triplet_adjoint)
01234   {
01235     delete vsparse_triplet_adjoint;
01236     vsparse_triplet_adjoint=0;
01237   }
01238   if (vsparse_triplet)
01239   {
01240     delete vsparse_triplet;
01241     vsparse_triplet=0;
01242   }
01243   if (sparse_triplet2)
01244   {
01245     delete sparse_triplet2;
01246     sparse_triplet2=0;
01247   }
01248   if (sparse_triplet)
01249   {
01250     delete sparse_triplet;
01251     sparse_triplet=0;
01252   }
01253   if (sparse_symbolic)
01254   {
01255     delete sparse_symbolic;
01256     sparse_symbolic=0;
01257   }
01258   if (sparse_symbolic2)
01259   {
01260     delete sparse_symbolic2;
01261     sparse_symbolic2=0;
01262   }
01263   if (num_local_re_array)
01264   {
01265     delete num_local_re_array;
01266     num_local_re_array = NULL;
01267   }
01268   if (num_local_fixed_array)
01269   {
01270     delete num_local_fixed_array;
01271     num_local_fixed_array = NULL;
01272   }
01273 }
01274 
01279 dvector laplace_approximation_calculator::operator()
01280   (const dvector& _x, const double& _f, function_minimizer * pfmin)
01281 {
01282   dvector g;
01283 
01284 #if defined(USE_ADPVM)
01285   if (pfmin->test_trust_flag)
01286   {
01287     return test_trust_region_method(_x,_f,pfmin);
01288   }
01289   if (ad_comm::pvm_manager)
01290   {
01291     switch (ad_comm::pvm_manager->mode)
01292     {
01293     case 1:
01294       return default_calculations_parallel_master(_x,_f,pfmin);
01295       break;
01296     case 2:
01297     {
01298       dvector g(1,1);
01299       default_calculations_parallel_slave(_x,_f,pfmin);
01300       return g;
01301       break;
01302     }
01303     default:
01304       cerr << "illegal value for mode " << endl;
01305       ad_exit(1);
01306     }
01307   }
01308   else
01309 #endif  //# if defined(USE_ADPVM)
01310 
01311   {
01312     //hesstype=1;
01313     //cout << hesstype << endl;
01314     switch (hesstype)
01315     {
01316       case 1:
01317       {
01318         int check_der_flag=0;
01319         int on=-1;
01320         int nopt=0;
01321         if ((on=option_match(ad_comm::argc,ad_comm::argv,"-chkder",nopt))>-1)
01322         {
01323           check_der_flag=1;
01324         }
01325         if (check_der_flag==1)
01326         {
01327           g = default_calculations_check_derivatives(_x,pfmin,_f);
01328         }
01329         else
01330         {
01331           g = default_calculations(_x,_f,pfmin);
01332         }
01333         break;
01334       }
01335       case 2:
01336       {
01337         // Hessian for the random effects is block diagonal
01338         g = block_diagonal_calculations(_x,_f,pfmin);
01339         break;
01340       }
01341       case 3:
01342       case 4:  // not banded but full hessian
01343       {
01344         // Hessian for the random effects is block diagonal
01345         if (laplace_approximation_calculator::variance_components_vector)
01346         {
01347           g = banded_calculations_lme(_x,_f,pfmin);
01348         }
01349         else
01350         {
01351           g = banded_calculations(_x,_f,pfmin);
01352         }
01353         break;
01354       }
01355       default:
01356       {
01357         cerr << "illegal value for hesstype " << endl;
01358         ad_exit(1);
01359       }
01360     }
01361   }
01362 
01363   return g;
01364 }
01365 
01366 void random_effects_userfunction(double f,const dvector& x, const dvector& g);
01367 
01372 void get_second_ders(int xs,int us,const init_df1b2vector _y,dmatrix& Hess,
01373   dmatrix& Dux, df1b2_gradlist * f1b2gradlist,function_minimizer * pfmin,
01374   laplace_approximation_calculator * lpc)
01375 {
01376   // Note that xs is the number of active non random effects
01377   // parameters
01378   // us is the number of random effects parameters
01379   int j;
01380   int i;
01381   ADUNCONST(init_df1b2vector,y)
01382   int num_der_blocks=lpc->num_der_blocks;
01383   int xsize=lpc->xsize;
01384   int usize=lpc->usize;
01385 
01386   for (int ip=1;ip<=num_der_blocks;ip++)
01387   {
01388     df1b2variable::minder=lpc->minder(ip);
01389     df1b2variable::maxder=lpc->maxder(ip);
01390     lpc->set_u_dot(ip);
01391     df1b2_gradlist::set_yes_derivatives();
01392     (*re_objective_function_value::pobjfun)=0;
01393     df1b2variable pen=0.0;
01394     df1b2variable zz=0.0;
01395     initial_params::straight_through_flag=0;
01396     //initial_params::straight_through_flag=1;
01397     initial_df1b2params::reset(y,pen);
01398     initial_params::straight_through_flag=0;
01399     if (initial_df1b2params::separable_flag)
01400     {
01401       initial_df1b2params::separable_calculation_type=2;
01402       Hess.initialize();
01403       Dux.initialize();
01404     }
01405 
01406     //cout << "2D" << endl;
01407     pfmin->user_function();
01408 
01409     //pfmin->user_function(y,zz);
01410     (*re_objective_function_value::pobjfun)+=pen;
01411     (*re_objective_function_value::pobjfun)+=zz;
01412     if (!initial_df1b2params::separable_flag)
01413     {
01414       set_dependent_variable(*re_objective_function_value::pobjfun);
01415       df1b2_gradlist::set_no_derivatives();
01416       df1b2variable::passnumber=1;
01417       df1b2_gradcalc1();
01418       int mind=y(1).minder;
01419       int jmin=max(mind,xsize+1);
01420       int jmax=min(y(1).maxder,xsize+usize);
01421       for (i=1;i<=usize;i++)
01422         for (j=jmin;j<=jmax;j++){
01423           //cout<<i<<" "<<j<<" "<<mind<<" "<<xsize<<" "<<usize<<endl; 
01424           Hess(i,j-xsize)=y(i+xsize).u_bar[j-mind];
01425         }
01426       jmin=max(mind,1);
01427       jmax=min(y(1).maxder,xsize);
01428       for (i=1;i<=usize;i++)
01429         for (j=jmin;j<=jmax;j++)
01430           Dux(i,j)=y(i+xsize).u_bar[j-1];
01431     }
01432     if (ip<num_der_blocks)
01433       f1b2gradlist->reset();
01434   }
01435 
01436   if  (ad_comm::print_hess_and_exit_flag)
01437   {
01438     for (i=1;i<=usize;i++)Hess(i,i) -= 2*CHECK_HESSIAN_PENALTY;
01439     cout<< "--------------------------------------------------------------------"<<endl;
01440     cout << "param:" << endl;
01441     for (i=1;i<=usize;i++)cout << y(i+xsize)<<" ";
01442     cout << endl;
01443     cout<< "--------------------------------------------------------------------"<<endl;
01444     cout << "Hess:" << endl;
01445     cout <<  Hess << endl;
01446     cout<< "--------------------------------------------------------------------"<<endl;
01447     ad_exit(0);
01448   }
01449 
01450   //cout << "Dux" << endl;
01451   //cout << setprecision(16) << Dux << endl;
01452 }
01453 
01458 double calculate_laplace_approximation(const dvector& x,const dvector& u0,
01459   const dmatrix& _Hess,const dvector& _xadjoint,const dvector& _uadjoint,
01460   const dmatrix& _Hessadjoint,function_minimizer * pmin)
01461 {
01462   ADUNCONST(dvector,xadjoint)
01463   ADUNCONST(dvector,uadjoint)
01464   ADUNCONST(dmatrix,Hessadjoint)
01465   ADUNCONST(dmatrix,Hess)
01466   const int xs=x.size();
01467   const int us=u0.size();
01468   gradient_structure::set_YES_DERIVATIVES();
01469   int nvar;
01470   //dcompressed_triplet & lst2 = *(pmin->lapprox->sparse_triplet);
01471   //hs_symbolic & ssymb=*(pmin->lapprox->sparse_symbolic);
01472   //dcompressed_triplet & xxxt = *(pmin->lapprox->sparse_triplet);
01473   dcompressed_triplet & lst = *(pmin->lapprox->sparse_triplet2);
01474   hs_symbolic & ssymb=*(pmin->lapprox->sparse_symbolic2);
01475   {
01476   /*
01477     cout << norm2(make_dmatrix(lst)-make_dmatrix(xxxt)) << endl;
01478     ofstream ofs1("m1");
01479     ofs1 << setfixed() << setprecision(3) << setw(10)
01480          << make_dmatrix(lst) << endl;
01481     ofstream ofs2("m2");
01482     ofs2 << setfixed() << setprecision(3) << setw(10)
01483          << make_dmatrix(xxxt) << endl;
01484    */
01485   }
01486 
01487   if (laplace_approximation_calculator::alternative_user_function_flag==0)
01488   {
01489     if (pmin->lapprox->sparse_hessian_flag==0)
01490     {
01491       nvar=x.size()+u0.size()+u0.size()*u0.size();
01492     }
01493     else
01494     {
01495       int sz= lst.indexmax()-lst.indexmin()+1;
01496       nvar=x.size()+u0.size()+sz;
01497     }
01498   }
01499   else
01500   {
01501     nvar=x.size()+u0.size();
01502   }
01503   independent_variables y(1,nvar);
01504 
01505   // need to set random effects active together with whatever
01506   // init parameters should be active in this phase
01507   initial_params::set_inactive_only_random_effects();
01508   initial_params::set_active_random_effects();
01509   /*int onvar=*/initial_params::nvarcalc();
01510   initial_params::xinit(y);    // get the initial values into the
01511   y(1,xs)=x;
01512 
01513   int i,j;
01514   dvar_vector d(1,xs+us);
01515 
01516   dmatrix Hess_save;
01517   // contribution for quadratic prior
01518   if (quadratic_prior::get_num_quadratic_prior()>0)
01519   {
01520     if (allocated(Hess_save)) Hess_save.deallocate();
01521     int mmin=Hess.indexmin();
01522     int mmax=Hess.indexmax();
01523     Hess_save.allocate(mmin,mmax,mmin,mmax);
01524     Hess_save=Hess;
01525     int & vxs = (int&)(xs);
01526     quadratic_prior::get_cHessian_contribution(Hess,vxs);
01527   }
01528 
01529   dvar_matrix vHess;
01530   int ii=xs+us+1;
01531   if (laplace_approximation_calculator::alternative_user_function_flag!=0)
01532   {
01533     dvar_vector vy=dvar_vector(y);
01534     initial_params::reset(vy);    // get the initial values into the
01535   }
01536   else
01537   {
01538     // Here need hooks for sparse matrix structures
01539     if (pmin->lapprox->sparse_hessian_flag==0)
01540     {
01541       for (i=1;i<=us;i++)
01542         for (j=1;j<=us;j++)
01543           y(ii++)=Hess(i,j);
01544     }
01545     else
01546     {
01547       int smin=lst.indexmin();
01548       int smax=lst.indexmax();
01549       for (i=smin;i<=smax;i++)
01550         y(ii++)=lst(i);
01551     }
01552 
01553     if (quadratic_prior::get_num_quadratic_prior()>0)
01554     {
01555       Hess=Hess_save;
01556     }
01557 
01558     dvar_vector vy=dvar_vector(y);
01559     initial_params::stddev_vscale(d,vy);
01560     if (pmin->lapprox->sparse_hessian_flag==0)
01561     {
01562       vHess.allocate(1,us,1,us);
01563     }
01564     initial_params::reset(vy);    // get the initial values into the
01565     ii=xs+us+1;
01566     if (initial_df1b2params::have_bounded_random_effects)
01567     {
01568       for (i=1;i<=us;i++)
01569       {
01570         if (d(i+xs)>0.0)
01571         {
01572           for (j=1;j<=us;j++)
01573           {
01574             if (d(j+xs)>0.0)
01575               vHess(i,j)=vy(ii++)/(d(i+xs)*d(j+xs));
01576             else
01577               vHess(i,j)=vy(ii++)/d(i+xs);
01578           }
01579         }
01580         else
01581         {
01582           for (j=1;j<=us;j++)
01583           {
01584             if (d(j+xs)>0.0)
01585               vHess(i,j)=vy(ii++)/d(j+xs);
01586             else
01587               vHess(i,j)=vy(ii++);
01588           }
01589         }
01590       }
01591     }
01592     else
01593     {
01594       if (laplace_approximation_calculator::sparse_hessian_flag==0)
01595       {
01596         for (i=1;i<=us;i++)
01597         {
01598           for (j=1;j<=us;j++)
01599           {
01600             vHess(i,j)=vy(ii++);
01601           }
01602         }
01603       }
01604       else
01605       {
01606         int mmin=lst.indexmin();
01607         int mmax=lst.indexmax();
01608         dvar_compressed_triplet * vsparse_triplet =
01609         pmin->lapprox->vsparse_triplet;
01610 
01611         if (vsparse_triplet==0)
01612         {
01613           pmin->lapprox->vsparse_triplet=
01614             new dvar_compressed_triplet(mmin,mmax,us,us);
01615           vsparse_triplet = pmin->lapprox->vsparse_triplet;
01616           for (i=mmin;i<=mmax;i++)
01617           {
01618             (*vsparse_triplet)(1,i)=lst(1,i);
01619             (*vsparse_triplet)(2,i)=lst(2,i);
01620           }
01621         }
01622         else
01623         {
01624           if (!allocated(*vsparse_triplet))
01625           {
01626             (*vsparse_triplet).allocate(mmin,mmax,us,us);
01627             for (i=mmin;i<=mmax;i++)
01628             {
01629               (*vsparse_triplet)(1,i)=lst(1,i);
01630               (*vsparse_triplet)(2,i)=lst(2,i);
01631             }
01632           }
01633         }
01634         dcompressed_triplet * vsparse_triplet_adjoint =
01635         pmin->lapprox->vsparse_triplet_adjoint;
01636 
01637         if (vsparse_triplet_adjoint==0)
01638         {
01639           pmin->lapprox->vsparse_triplet_adjoint=
01640           new dcompressed_triplet(mmin,mmax,us,us);
01641           vsparse_triplet_adjoint = pmin->lapprox->vsparse_triplet_adjoint;
01642           for (i=mmin;i<=mmax;i++)
01643           {
01644             (*vsparse_triplet_adjoint)(1,i)=lst(1,i);
01645             (*vsparse_triplet_adjoint)(2,i)=lst(2,i);
01646           }
01647         }
01648         else
01649         {
01650           if (!allocated(*vsparse_triplet_adjoint))
01651           {
01652             (*vsparse_triplet_adjoint).allocate(mmin,mmax,us,us);
01653             for (i=mmin;i<=mmax;i++)
01654             {
01655               (*vsparse_triplet_adjoint)(1,i)=lst(1,i);
01656               (*vsparse_triplet_adjoint)(2,i)=lst(2,i);
01657             }
01658           }
01659         }
01660         vsparse_triplet->get_x()=vy(ii,ii+mmax-mmin).shift(1);
01661       }
01662     }
01663   }
01664 
01665   dvariable vf=0.0;
01666 
01667   *objective_function_value::pobjfun=0.0;
01668   if (laplace_approximation_calculator::alternative_user_function_flag==1)
01669   {
01670     laplace_approximation_calculator::alternative_user_function_flag=2;
01671   }
01672   pmin->AD_uf_outer();
01673   if (laplace_approximation_calculator::alternative_user_function_flag==2)
01674   {
01675     laplace_approximation_calculator::alternative_user_function_flag=1;
01676   }
01677   if ( no_stuff==0 && quadratic_prior::get_num_quadratic_prior()>0)
01678   {
01679     quadratic_prior::get_M_calculations();
01680   }
01681   vf+=*objective_function_value::pobjfun;
01682   //cout << setprecision(15) << vf << endl;
01683   // *********************************************
01684   // *********************************************
01685   // *********************************************
01686   // dvector tmpg(1,nvar);
01687   // tmpg.initialize();
01688   // gradcalc(nvar,tmpg);
01689   // *********************************************
01690   // *********************************************
01691   // *********************************************
01692 
01693    int sgn=0;
01694    dvariable ld = 0;
01695    double f=0.0;
01696    if (laplace_approximation_calculator::alternative_user_function_flag==0)
01697    {
01698      if (ad_comm::no_ln_det_choleski_flag)
01699      {
01700        if(laplace_approximation_calculator::saddlepointflag==0)
01701        {
01702          ld=0.5*ln_det(vHess,sgn);
01703        }
01704        else
01705        {
01706          ld=0.5*ln_det(-vHess,sgn);
01707        }
01708      }
01709      else
01710      {
01711        if(laplace_approximation_calculator::saddlepointflag==0)
01712        {
01713          int ierr=0;
01714          if (pmin->lapprox->sparse_hessian_flag==0)
01715          {
01716            ld=0.5*ln_det_choleski_error(vHess,ierr);
01717            if (ierr==1)
01718            {
01719              ofstream ofs("hessian.diag");
01720              ofs << vHess << endl;
01721              ofs <<  x << endl;
01722              ofs <<  u0 << endl;
01723              ofs << "Matrix not positive definite in Ln_det_choleski"
01724                  << endl;
01725              cerr << "Matrix not positive definite in Ln_det_choleski\n"
01726                   << "see file hessian.diag for details"
01727                  << endl;
01728              ad_exit(1);
01729            }
01730          }
01731          else
01732          {
01733            //double ld1=0.5*ln_det(*(pmin->lapprox->sparse_triplet),
01734            //  *(pmin->lapprox->sparse_symbolic));
01735 
01736            if (write_sparse_flag)
01737            {
01738              //ofstream ofs("sparse");
01739              //ofs << *(pmin->lapprox->vsparse_triplet) << endl;
01740            }
01741            ld=0.5*ln_det(*(pmin->lapprox->vsparse_triplet),
01742              ssymb,*(pmin->lapprox->sparse_triplet2));
01743              //*(pmin->lapprox->sparse_symbolic2),pmin->lapprox);
01744            //cout << ld-ld1 << endl;
01745          }
01746        }
01747        else
01748        {
01749          if (pmin->lapprox->sparse_hessian_flag==0)
01750          {
01751            ld=0.5*ln_det_choleski(-vHess);
01752          }
01753          else
01754          {
01755            cerr << "need to fix this" << endl;
01756            ad_exit(1);
01757            //ld+=ln_det(-*(pmin->lapprox->vsparse_triplet));
01758          }
01759        }
01760      }
01761 
01762 #ifdef DIAG
01763      int ps1=0;
01764      if (ps1)
01765      {
01766        dmatrix cHess=value(vHess);
01767        cout << " ln_det = " << ld << " ";
01768        dvector eig=eigenvalues(cHess);
01769        dmatrix r(1,2,1,eig.indexmax());
01770        dvector num(1,eig.indexmax());
01771        num.fill_seqadd(1,1);
01772        r(1)=eig;
01773        r(2)=num;
01774        dmatrix tt=trans(r);
01775        dmatrix ss=trans(sort(tt,1));
01776        cout << ss(1,1)  << " " << ss(1,eig.indexmax()) << " " ;
01777        int nnn=(int)ss(2,1);
01778        dmatrix d=eigenvectors(cHess);
01779        //cout << " d " << d(nnn) << " " << d(nnn)*cHess*d(nnn) << endl;
01780        dmatrix t=trans(d);
01781        r(1)=t(nnn);
01782        r(2)=num;
01783        dmatrix tt2=trans(r);
01784        dmatrix ss2=trans(sort(tt2,1));
01785 
01786        cout << " t " << setprecision(3) << ss2(1)(1,5) << " --- "
01787             << t(nnn)*cHess*t(nnn) << endl;
01788        cout << "   " << setprecision(3) << ss2(2)(1,5) << endl;
01789        //cout << " t " << t(1) << " " << t(1)*cHess*t(2) << endl;
01790      }
01791 #endif
01792 
01793      int nx=0;
01794      if (nx==0)
01795      {
01796        vf+=ld;
01797      }
01798      f=value(vf);
01799      f-=us*0.5*log(2.0*PI);
01800    }
01801    else
01802    {
01803      f=value(vf);
01804    }
01805 
01806    dvector g(1,nvar);
01807    gradcalc(nvar,g);
01808 
01809   ii=1;
01810   for (i=1;i<=xs;i++)
01811     xadjoint(i)=g(ii++);
01812   for (i=1;i<=us;i++)
01813     uadjoint(i)=g(ii++);
01814 
01815   if (laplace_approximation_calculator::alternative_user_function_flag==0)
01816   {
01817     if (pmin->lapprox->sparse_hessian_flag==0)
01818     {
01819       if (allocated(Hessadjoint))
01820       {
01821         for (i=1;i<=us;i++)
01822           for (j=1;j<=us;j++)
01823             Hessadjoint(i,j)=g(ii++);
01824       }
01825     }
01826     else
01827     {
01828       dcompressed_triplet * vsparse_triplet_adjoint =
01829         pmin->lapprox->vsparse_triplet_adjoint;
01830 
01831       int smin=lst.indexmin();
01832       int smax=lst.indexmax();
01833       for (i=smin;i<=smax;i++)
01834       {
01835         (*vsparse_triplet_adjoint)(i)=g(ii);
01836         ii++;
01837       }
01838     }
01839   }
01840 
01841    /*
01842    if (quadratic_prior::get_num_quadratic_prior()>0)
01843    {
01844      // *******************************************************
01845      // *******************************************************
01846      // for quadratic prior option
01847      // temporarily get der wrt x of x ---> ln_det(F_uu + inv(A(x)))
01848      int nvar=x.size()+u0.size();
01849      independent_variables y(1,nvar);
01850      initial_params::xinit(y);    // get the initial values into the
01851      y(1,xs)=x;
01852      dvar_vector vy=dvar_vector(y);
01853      initial_params::reset(vy);    // get the initial values into the
01854 
01855      pmin->AD_uf_outer();
01856 
01857      dvar_matrix tmpH=quadratic_prior::get_vHessian_contribution();
01858 
01859      tmpH+=Hess;
01860      dvariable ld;
01861      ld=0.5*ln_det(tmpH,sgn);
01862      dvector g(1,nvar);
01863      gradcalc(nvar,g);
01864      int ii=1;
01865      double checksum=0.0;
01866      for (i=1;i<=xs;i++)
01867        xadjoint(i)+=g(ii++);
01868      for (i=1;i<=us;i++)
01869        checksum+=g(ii++);
01870 
01871      if (fabs(checksum)> 1.e-10)
01872      {
01873        cerr << "checksum too big " << checksum << endl;
01874        ad_exit(1);
01875      }
01876 
01877    }
01878   */
01879 
01880 
01881    // *******************************************************
01882    // *******************************************************
01883 
01884   return f;
01885 }
01886 
01891 void get_newton_raphson_info(int xs,int us,const init_df1b2vector _y,
01892   dmatrix& Hess, dvector& grad, df1b2_gradlist * f1b2gradlist,
01893   function_minimizer * pfmin)
01894 {
01895   // Note that xs is the number of active non random effects
01896   // parameters
01897   // us is the number of random effects parameters
01898   int j;
01899   int i;
01900   ADUNCONST(init_df1b2vector,y)
01901   {
01902     df1b2_gradlist::set_yes_derivatives();
01903     //cout << re_objective_function_value::pobjfun << endl;
01904     //cout << re_objective_function_value::pobjfun->ptr << endl;
01905     (*re_objective_function_value::pobjfun)=0;
01906     df1b2variable pen=0.0;
01907     df1b2variable zz=0.0;
01908     initial_df1b2params::reset(y,pen);
01909     pfmin->user_function();
01910     //pfmin->user_function(y,zz);
01911     (*re_objective_function_value::pobjfun)+=pen;
01912     (*re_objective_function_value::pobjfun)+=zz;
01913     set_dependent_variable(*re_objective_function_value::pobjfun);
01914     df1b2_gradlist::set_no_derivatives();
01915     df1b2variable::passnumber=1;
01916   }
01917 
01918   int mind=y(1).minder;
01919   int jmin=max(mind,xs+1);
01920   int jmax=min(y(1).maxder,xs+us);
01921   if (!initial_df1b2params::separable_flag)
01922   {
01923     df1b2_gradcalc1();
01924     for (i=1;i<=us;i++)
01925       for (j=jmin;j<=jmax;j++)
01926         Hess(i,j-xs)=y(i+xs).u_bar[j-mind];
01927     for (j=jmin;j<=jmax;j++)
01928       grad(j-xs)= re_objective_function_value::pobjfun->u_dot[j-mind];
01929   }
01930 }
01931 
01935 void laplace_approximation_calculator::check_for_need_to_reallocate(int ip)
01936 {
01937 }
01938   //void laplace_approximation_calculator::get_newton_raphson_info
01939   //  (function_minimizer * pfmin)
01940   //{
01941   //  int i,j,ip;
01942   //
01943   //  for (ip=1;ip<=num_der_blocks;ip++)
01944   //  {
01945   //    // do we need to reallocate memory for df1b2variables?
01946   //    check_for_need_to_reallocate(ip);
01947   //    df1b2_gradlist::set_yes_derivatives();
01948   //    //cout << re_objective_function_value::pobjfun << endl;
01949   //    //cout << re_objective_function_value::pobjfun->ptr << endl;
01950   //    (*re_objective_function_value::pobjfun)=0;
01951   //    df1b2variable pen=0.0;
01952   //    df1b2variable zz=0.0;
01953   //    initial_df1b2params::reset(y,pen);
01954   //    if (initial_df1b2params::separable_flag)
01955   //    {
01956   //      Hess.initialize();
01957   //      grad.initialize();
01958   //    }
01959   //    pfmin->user_function();
01960   //    //pfmin->user_function(y,zz);
01961   //    (*re_objective_function_value::pobjfun)+=pen;
01962   //    (*re_objective_function_value::pobjfun)+=zz;
01963   //
01964   //    if (!initial_df1b2params::separable_flag)
01965   //    {
01966   //      set_dependent_variable(*re_objective_function_value::pobjfun);
01967   //      df1b2_gradlist::set_no_derivatives();
01968   //      df1b2variable::passnumber=1;
01969   //      df1b2_gradcalc1();
01970   //      int mind=y(1).minder;
01971   //      int jmin=max(mind,xsize+1);
01972   //      int jmax=min(y(1).maxder,xsize+usize);
01973   //      for (i=1;i<=usize;i++)
01974   //        for (j=jmin;j<=jmax;j++)
01975   //          Hess(i,j-xsize)=y(i+xsize).u_bar[j-mind];
01976   //
01977   //      for (j=jmin;j<=jmax;j++)
01978   //        grad(j-xsize)= re_objective_function_value::pobjfun->u_dot[j-mind];
01979   //    }
01980   //  }
01981   //  {
01982   //    //ofstream ofs("hess.tmp");
01983   //    //ofs << Hess << endl;
01984   //    //ofs << grad << endl;
01985   //  }
01986   //}
01987   //
01988 
01993 double evaluate_function(const dvector& x,function_minimizer * pfmin)
01994 {
01995   int usize=initial_params::nvarcalc();
01996   //double f=0.0;
01997   dvector g(1,usize);
01998   independent_variables u(1,usize);
01999   u=x;
02000   dvariable vf=0.0;
02001   vf=initial_params::reset(dvar_vector(u));
02002   //vf=0.0;
02003   *objective_function_value::pobjfun=0.0;
02004   laplace_approximation_calculator::where_are_we_flag=2;
02005   pfmin->AD_uf_inner();
02006   if ( no_stuff==0 && quadratic_prior::get_num_quadratic_prior()>0)
02007   {
02008     quadratic_prior::get_M_calculations();
02009   }
02010   vf+=*objective_function_value::pobjfun;
02011   laplace_approximation_calculator::where_are_we_flag=0;
02012   initial_df1b2params::cobjfun=value(vf);
02013   gradcalc(usize,g);
02014   double maxg=max(fabs(g));
02015   if (!initial_params::mc_phase)
02016   {
02017     cout << setprecision(16) << " f = " << vf
02018          << " max g = " << maxg << endl;
02019   }
02020   return maxg;
02021 }
02022 
02027 double evaluate_function(double& fval,const dvector& x,
02028   function_minimizer* pfmin)
02029 {
02030   int usize=initial_params::nvarcalc();
02031   //double f=0.0;
02032   dvector g(1,usize);
02033   independent_variables u(1,usize);
02034   u=x;
02035   dvariable vf=0.0;
02036   vf=initial_params::reset(dvar_vector(u));
02037   //vf=0.0;
02038   *objective_function_value::pobjfun=0.0;
02039   laplace_approximation_calculator::where_are_we_flag=2;
02040   pfmin->AD_uf_inner();
02041   if ( no_stuff==0 && quadratic_prior::get_num_quadratic_prior()>0)
02042   {
02043     quadratic_prior::get_M_calculations();
02044   }
02045   vf+=*objective_function_value::pobjfun;
02046   laplace_approximation_calculator::where_are_we_flag=0;
02047   initial_df1b2params::cobjfun=value(vf);
02048   gradcalc(usize,g);
02049   double maxg=max(fabs(g));
02050   fval=value(vf);
02051   if (!initial_params::mc_phase)
02052   {
02053     cout << setprecision(10) << " f = " << vf
02054          << " max g = " << maxg << endl;
02055   }
02056   return maxg;
02057 }
02058 
02063 double evaluate_function(double& fval,const dvector& x,const dvector& g,
02064   function_minimizer * pfmin)
02065 {
02066   int usize=initial_params::nvarcalc();
02067   //double f=0.0;
02068   //dvector g(1,usize);
02069   independent_variables u(1,usize);
02070   u=x;
02071   dvariable vf=0.0;
02072   vf=initial_params::reset(dvar_vector(u));
02073   //vf=0.0;
02074   *objective_function_value::pobjfun=0.0;
02075   laplace_approximation_calculator::where_are_we_flag=2;
02076   pfmin->AD_uf_inner();
02077   if ( no_stuff==0 && quadratic_prior::get_num_quadratic_prior()>0)
02078   {
02079     quadratic_prior::get_M_calculations();
02080   }
02081   vf+=*objective_function_value::pobjfun;
02082   laplace_approximation_calculator::where_are_we_flag=0;
02083   initial_df1b2params::cobjfun=value(vf);
02084   gradcalc(usize,g);
02085   double maxg=max(fabs(g));
02086   fval=value(vf);
02087   if (!initial_params::mc_phase)
02088   {
02089     cout << setprecision(15) << " f = " << vf
02090          << " max g = " << maxg << endl;
02091   }
02092   return maxg;
02093 }
02094 
02099 double evaluate_function_quiet(const dvector& x,function_minimizer * pfmin)
02100 {
02101   int usize=initial_params::nvarcalc();
02102   //double f=0.0;
02103   dvector g(1,usize);
02104   independent_variables u(1,usize);
02105   u=x;
02106   dvariable vf=0.0;
02107   vf=initial_params::reset(dvar_vector(u));
02108   //vf=0.0;
02109   *objective_function_value::pobjfun=0.0;
02110   laplace_approximation_calculator::where_are_we_flag=2;
02111   pfmin->AD_uf_inner();
02112   if ( no_stuff==0 && quadratic_prior::get_num_quadratic_prior()>0)
02113   {
02114     quadratic_prior::get_M_calculations();
02115   }
02116   vf+=*objective_function_value::pobjfun;
02117   laplace_approximation_calculator::where_are_we_flag=0;
02118   initial_df1b2params::cobjfun=value(vf);
02119   gradcalc(usize,g);
02120   double maxg=max(fabs(g));
02121   return maxg;
02122 }
02123 
02128 void evaluate_function_gradient(double& f,const dvector& x,
02129   function_minimizer * pfmin,dvector& xadjoint,dvector& uadjoint)
02130 {
02131   int usize=initial_params::nvarcalc();
02132   dvector g(1,usize);
02133   independent_variables u(1,usize);
02134   u=x;
02135   dvariable vf=0.0;
02136   vf=initial_params::reset(dvar_vector(u));
02137   //vf=0.0;
02138   *objective_function_value::pobjfun=0.0;
02139   pfmin->AD_uf_inner();
02140   if ( no_stuff==0 && quadratic_prior::get_num_quadratic_prior()>0)
02141   {
02142     quadratic_prior::get_M_calculations();
02143   }
02144   vf+=*objective_function_value::pobjfun;
02145   initial_df1b2params::cobjfun=value(vf);
02146   f=value(vf);
02147   gradcalc(usize,g);
02148   int xsize=xadjoint.indexmax();
02149   int us=uadjoint.indexmax();
02150   xadjoint=g(1,xsize);
02151   uadjoint=g(xsize+1,xsize+us).shift(1);
02152 }
02153 
02158 double evaluate_function_no_derivatives(const dvector& x,
02159   function_minimizer* pfmin)
02160 {
02161   double fval;
02162   gradient_structure::set_NO_DERIVATIVES();
02163   int usize=initial_params::nvarcalc();
02164   //double f=0.0;
02165   dvector g(1,usize);
02166   independent_variables u(1,usize);
02167   u=x;
02168   dvariable vf=0.0;
02169   vf=initial_params::reset(dvar_vector(u));
02170   //vf=0.0;
02171   *objective_function_value::pobjfun=0.0;
02172   pfmin->AD_uf_inner();
02173   if ( no_stuff==0 && quadratic_prior::get_num_quadratic_prior()>0)
02174   {
02175     quadratic_prior::get_M_calculations();
02176   }
02177   vf+=*objective_function_value::pobjfun;
02178   initial_df1b2params::cobjfun=value(vf);
02179   fval=value(vf);
02180   gradient_structure::set_YES_DERIVATIVES();
02181   return fval;
02182 }
02183 
02188 void cleanup_laplace_stuff(laplace_approximation_calculator * l)
02189 {
02190   if (l)
02191   {
02192     delete l;
02193     l=0;
02194     df1b2variable::pool->deallocate();
02195 
02196     for (int i=0;i<df1b2variable::adpool_counter;i++)
02197     {
02198       delete df1b2variable::adpool_vector[i];
02199       df1b2variable::adpool_vector[i]=0;
02200       df1b2variable::nvar_vector[i]=0;
02201       df1b2variable::adpool_counter=0;
02202     }
02203   }
02204 }
02205 
02210 gauss_hermite_stuff::gauss_hermite_stuff
02211   (laplace_approximation_calculator * lapprox,
02212   int use_gauss_hermite,int num_separable_calls ,const ivector& itmp) :
02213   x(1,use_gauss_hermite), w(1,use_gauss_hermite), mi(0)
02214 {
02215   // for now we must have the same number of random effects in each
02216   // separable function call"
02217   int i;
02218   for (i=2;i<=num_separable_calls;i++)
02219   {
02220     if (itmp(i)!=itmp(i-1))
02221     {
02222       cerr << " At present for the adaptive gauss_hermite must have the same"
02223            << endl
02224            << " number of random effects in each separable function call"
02225            << endl;
02226       ad_exit(1);
02227     }
02228   }
02229   for (i=1;i<=num_separable_calls;i++)
02230   {
02231     if (itmp(i)>1)
02232     {
02233       lapprox->multi_random_effects=
02234         max(lapprox->multi_random_effects,itmp(i));
02235      /*
02236       cerr << "At present gauss-hermite is only implemented for"
02237         " one random effect per separable function call "
02238        << endl;
02239       ad_exit(1);
02240      */
02241     }
02242   }
02243   if (allocated(gauss_hermite_values))
02244     gauss_hermite_values.deallocate();
02245   if (lapprox->multi_random_effects==0)
02246   {
02247     gauss_hermite_values.allocate(1,num_separable_calls,1,use_gauss_hermite);
02248   }
02249   else
02250   {
02251    ivector indx=pow(use_gauss_hermite,itmp);
02252     gauss_hermite_values.allocate(1,num_separable_calls,1,indx);
02253   }
02254 
02255   normalized_gauss_hermite(x,w);
02256   is=0;
02257 }
02258 
02263 void laplace_approximation_calculator::set_default_hessian_type(void )
02264 {
02265   //if (function_minimizer::hesstype==0)
02266   {
02267     if (initial_df1b2params::separable_flag)
02268     {
02269      // have SEPARABLE_FUNCTION but no object of type quadratic_penalty or
02270      // normal prior -- can't tell if we should use hesstype 2 or 3
02271      // until we run
02272       hesstype=3;
02273     }
02274     else
02275     {
02276       hesstype=1;  // assume block diagonal
02277     }
02278   }
02279  /*
02280   else
02281   {
02282     hesstype=function_minimizer::hesstype;
02283   }
02284  */
02285 }
02286 
02291 double laplace_approximation_calculator::get_fx_fu(function_minimizer * pfmin)
02292 {
02293   initial_params::set_inactive_only_random_effects();
02294   initial_params::set_active_random_effects();
02295 
02296   if (grad_x==0)
02297   {
02298     grad_x=new dvector(1,xsize);
02299   }
02300   grad_x->initialize();
02301 
02302   if (grad_x_u==0)
02303   {
02304     grad_x_u=new dvector(1,xsize+usize);
02305   }
02306   pfmin->inner_opt_flag=0;
02307   independent_variables u(1,xsize+usize);
02308   //gradcalc(0,*grad_x_u);
02309   initial_params::xinit(u);    // get the initial values into the
02310 
02311   dvariable pen=0.0;
02312   dvariable vf=0.0;
02313   pen=initial_params::reset(dvar_vector(u));
02314   *objective_function_value::pobjfun=0.0;
02315   pfmin->AD_uf_inner();
02316   vf+=*objective_function_value::pobjfun;
02317   objective_function_value::fun_without_pen=value(vf);
02318   vf+=pen;
02319   gradcalc(xsize+usize,*grad_x_u);
02320   double f=value(vf);
02321   return f;
02322 }
02323 
02328   void laplace_approximation_calculator::begin_separable_call_stuff(void)
02329   {
02330     separable_call_level++;
02331     //build_up_nested_shape();
02332     //clean(nested_tree_position,separable_call_level);
02333     //nested_separable_calls_counter(separable_call_level)++;
02334     //nested_tree_position(separable_call_level)++;
02335   }
02336 
02341   void laplace_approximation_calculator::end_separable_call_stuff(void)
02342   {
02343     //build_up_nested_shape();
02344     //clean(nested_tree_position,separable_call_level);
02345     separable_call_level--;
02346   }
02347 
02352 void laplace_approximation_calculator::build_up_nested_shape(void)
02353 {
02354   int ll;
02355   ivector& a =nested_separable_calls_counter;
02356   int clean_level=0;
02357   switch(separable_call_level)
02358   {
02359   case 1:
02360     ll=1;
02361     if (nested_separable_calls_counter(ll)>0)
02362     {
02363       if (clean_level==0) clean_level=ll+1;
02364       if (nested_separable_calls_counter(ll+1)>0)
02365       {
02366         nested_shape(a(ll))=a(ll+1);
02367       }
02368       else
02369       {
02370         break;
02371       }
02372     }
02373     else
02374     {
02375       break;
02376     }
02377    case 2:
02378     ll=2;
02379     if (nested_separable_calls_counter(ll)>0)
02380     {
02381       if (clean_level==0) clean_level=ll+1;
02382       if (nested_separable_calls_counter(ll+1)>0)
02383       {
02384         nested_shape(a(1),a(2))=a(ll+1);
02385       }
02386       else
02387       {
02388         break;
02389       }
02390     }
02391     else
02392     {
02393       break;
02394     }
02395    case 3:
02396     ll=3;
02397     if (nested_separable_calls_counter(ll)>0)
02398     {
02399       if (clean_level==0) clean_level=ll+1;
02400       if (nested_separable_calls_counter(ll+1)>0)
02401       {
02402         nested_shape(a(1),a(2),a(3))=a(ll+1);
02403       }
02404       else
02405       {
02406         break;
02407       }
02408     }
02409     else
02410     {
02411       break;
02412     }
02413    case 4:
02414     ll=4;
02415     if (nested_separable_calls_counter(ll)>0)
02416     {
02417       if (clean_level==0) clean_level=ll+1;
02418       if (nested_separable_calls_counter(ll+1)>0)
02419       {
02420         nested_shape(a(1),a(2),a(3),a(4))=a(ll+1);
02421       }
02422       else
02423       {
02424         break;
02425       }
02426     }
02427     else
02428     {
02429       break;
02430     }
02431    default:
02432      cerr << "illegal value in " <<
02433        "laplace_approximation_calculator::build_up_nested_shape"
02434        << endl;
02435   }
02436   if (clean_level>0)
02437   {
02438     int mmax=a.indexmax();
02439     for (int i=clean_level;i<=mmax;i++)
02440     {
02441       a(i)=0;
02442     }
02443   }
02444 }
02445 
02450 ivector * nested_calls_shape::get_ptr1(void)
02451 {
02452   return ptr1;
02453 }
02454 
02459 imatrix * nested_calls_shape::get_ptr2(void)
02460 {
02461   return ptr2;
02462 }
02463 
02468 i3_array * nested_calls_shape::get_ptr3(void)
02469 {
02470   return ptr3;
02471 }
02472 
02477 i4_array * nested_calls_shape::get_ptr4(void)
02478 {
02479   return ptr4;
02480 }
02481 
02486 ostream &  operator << (const ostream& _s,const nested_calls_shape& _m)
02487 {
02488   ADUNCONST(nested_calls_shape,m)
02489   ADUNCONST(ofstream,s)
02490   if (m.get_ptr1())
02491     s<< *(m.get_ptr1()) << endl << endl;
02492 
02493   if (m.get_ptr2())
02494     s<< *(m.get_ptr2()) << endl << endl;
02495 
02496   if (m.get_ptr3())
02497     s<< *(m.get_ptr3()) << endl << endl;
02498 
02499   if (m.get_ptr4())
02500     s<< *(m.get_ptr4()) << endl << endl;
02501 
02502   return s;
02503 }
02504 
02509 void nested_calls_shape::trim(void)
02510 {
02511   int mmax1=0;
02512   if (ptr1)
02513   {
02514     int imin=ptr1->indexmin();
02515     int imax=ptr1->indexmax();
02516     for (int i=imin;i<=imax;i++)
02517     {
02518       if ( (*ptr1)(i)==0) break;
02519       mmax1++;
02520     }
02521   }
02522   if (mmax1==0)
02523   {
02524     delete ptr1;
02525     ptr1=0;
02526   }
02527   else
02528   {
02529     ivector * tmp=new ivector(1,mmax1);
02530     (*tmp)=(*ptr1)(1,mmax1);
02531     delete ptr1;
02532     ptr1=tmp;
02533   }
02534   if (ptr2)
02535   {
02536     if (!ptr1)
02537     {
02538       delete ptr2;
02539       ptr2=0;
02540     }
02541     else
02542     {
02543       imatrix * tmp=new imatrix(1,mmax1);
02544       int imin=tmp->indexmin();
02545       int imax=tmp->indexmax();
02546       for (int i=imin;i<=imax;i++)
02547       {
02548         int jmin=(*ptr2)(i).indexmin();
02549         int jmax=(*ptr2)(i).indexmax();
02550         int mmax2=0;
02551         for (int j=jmin;j<=jmax;j++)
02552         {
02553           if ((*ptr2)(i,j)==0) break;
02554           mmax2++;
02555         }
02556         if (mmax2>0)
02557         {
02558           (*tmp)(i).allocate(1,mmax2);
02559           (*tmp)(i)=(*ptr2)(i)(1,mmax2);
02560         }
02561       }
02562       delete ptr2;
02563       ptr2=tmp;
02564     }
02565   }
02566   if (ptr3)
02567   {
02568     cerr << "warning not dealitn with prt3" << endl;
02569     delete ptr3;
02570     ptr3=0;
02571   }
02572 }
02573 
02578 nested_calls_shape::~nested_calls_shape()
02579 {
02580   if (ptr1)
02581   {
02582     delete ptr1;
02583     ptr1=0;
02584   }
02585   if (ptr2)
02586   {
02587     delete ptr2;
02588     ptr2=0;
02589   }
02590   if (ptr3)
02591   {
02592     delete ptr3;
02593     ptr3=0;
02594   }
02595   if (ptr4)
02596   {
02597     delete ptr4;
02598     ptr4=0;
02599   }
02600 }
02601 
02606 void nested_calls_shape::allocate(int n,int m,int p )
02607 {
02608   if (ptr1)
02609   {
02610     delete ptr1;
02611     ptr1=0;
02612   }
02613   ptr1=new ivector(1,n);
02614 
02615   if (ptr2)
02616   {
02617     delete ptr2;
02618     ptr2=0;
02619   }
02620   ptr2=new imatrix(1,n,1,m);
02621   if (ptr3)
02622   {
02623     delete ptr3;
02624     ptr3=0;
02625   }
02626   ptr3=new i3_array(1,n,1,m,1,p);
02627   /*
02628   if (ptr4)
02629   {
02630     delete ptr4;
02631     ptr4=0;
02632   }
02633   */
02634 }
02635 
02640 void nested_calls_shape::allocate(int n,int m)
02641 {
02642   if (ptr1)
02643   {
02644     delete ptr1;
02645     ptr1=0;
02646   }
02647   ptr1=new ivector(1,n);
02648 
02649   if (ptr2)
02650   {
02651     delete ptr2;
02652     ptr2=0;
02653   }
02654   ptr2=new imatrix(1,n,1,m);
02655   /*
02656   if (ptr3)
02657   {
02658     delete ptr3;
02659     ptr3=0;
02660   }
02661   ptr=new i3_array(1,n,1,m,1,p);
02662   if (ptr4)
02663   {
02664     delete ptr4;
02665     ptr4=0;
02666   }
02667   */
02668 }
02669 
02674 void nested_calls_shape::initialize(void)
02675 {
02676   if (ptr1)
02677   {
02678     ptr1->initialize();
02679   }
02680 
02681   if (ptr2)
02682   {
02683     ptr2->initialize();
02684   }
02685 
02686   if (ptr3)
02687   {
02688     ptr3->initialize();
02689   }
02690 
02691   if (ptr4)
02692   {
02693     ptr4->initialize();
02694   }
02695 }
02696 
02701 void nested_calls_shape::allocate(int n)
02702 {
02703   if (ptr1)
02704   {
02705     delete ptr1;
02706     ptr1=0;
02707   }
02708   ptr1=new ivector(1,n);
02709 
02710   /*
02711   if (ptr2)
02712   {
02713     delete ptr2;
02714     ptr2=0;
02715   }
02716   ptr=new imatrix(1,n,1,m);
02717   if (ptr3)
02718   {
02719     delete ptr3;
02720     ptr3=0;
02721   }
02722   ptr=new i3_array(1,n,1,m,1,p);
02723   if (ptr4)
02724   {
02725     delete ptr4;
02726     ptr4=0;
02727   }
02728   */
02729 }
02730 
02735 void nested_calls_indices::allocate(const nested_calls_shape& _nsc)
02736 {
02737   int mmax=0;
02738   ADUNCONST(nested_calls_shape,nsc)
02739   mmax=nsc.get_ptr1()->indexmax();
02740   if (ptr1)
02741   {
02742     delete ptr1;
02743     ptr1=0;
02744   }
02745   if (nsc.get_ptr1())
02746   {
02747     ptr1=new imatrix(1,mmax,1,*nsc.get_ptr1());
02748   }
02749   if (ptr2)
02750   {
02751     delete ptr2;
02752     ptr2=0;
02753   }
02754   if (nsc.get_ptr2())
02755   {
02756     ptr2=new i3_array(1,mmax,1,*nsc.get_ptr1(),1,*nsc.get_ptr2());
02757   }
02758   if (ptr3)
02759   {
02760     delete ptr3;
02761     ptr3=0;
02762   }
02763   if (nsc.get_ptr3())
02764   {
02765     ptr3=new i4_array(1,mmax,1,*nsc.get_ptr1(),1,*nsc.get_ptr2(),
02766       1,*nsc.get_ptr3());
02767   }
02768 }
02769 
02774 dvector laplace_approximation_calculator::get_uhat_lm_newton2
02775   (const dvector& x,function_minimizer * pfmin)
02776 {
02777   //int on,nopt;
02778   pfmin->inner_opt_flag=1;
02779   double f=0.0;
02780   double fb=1.e+100;
02781   dvector g(1,usize);
02782   dvector ub(1,usize);
02783   independent_variables u(1,usize);
02784   gradcalc(0,g);
02785   fmc1.itn=0;
02786   fmc1.ifn=0;
02787   fmc1.ireturn=0;
02788   initial_params::xinit(u);    // get the initial values into the
02789   initial_params::xinit(ubest);    // get the initial values into the
02790   fmc1.ialph=0;
02791   fmc1.ihang=0;
02792   fmc1.ihflag=0;
02793 
02794   if (init_switch)
02795   {
02796     u.initialize();
02797   }
02798   else
02799   {
02800     u=ubest;
02801   }
02802 
02803   fmc1.dfn=1.e-2;
02804   dvariable pen=0.0;
02805   //cout << "starting  norm(u) = " << norm(u) << endl;
02806   while (fmc1.ireturn>=0)
02807   {
02808    /*
02809     double nu=norm(value(u));
02810     if (nu>400)
02811     {
02812       cout << "U norm(u) = " << nu  << endl;
02813     }
02814     cout << "V norm(u) = " << nu
02815          << " f = " << f  << endl;
02816     */
02817     fmc1.fmin(f,u,g);
02818     //cout << "W norm(u) = " << norm(value(u)) << endl;
02819     if (fmc1.ireturn>0)
02820     {
02821       dvariable vf=0.0;
02822       pen=initial_params::reset(dvar_vector(u));
02823       *objective_function_value::pobjfun=0.0;
02824       pfmin->AD_uf_inner();
02825       if (saddlepointflag)
02826       {
02827         *objective_function_value::pobjfun*=-1.0;
02828       }
02829       if ( no_stuff==0 && quadratic_prior::get_num_quadratic_prior()>0)
02830       {
02831         quadratic_prior::get_M_calculations();
02832       }
02833       vf+=*objective_function_value::pobjfun;
02834 
02835      /*  this is now done in the operator = function
02836       if (quadratic_prior::get_num_quadratic_prior()>0)
02837       {
02838         vf+= quadratic_prior::get_quadratic_priors();
02839       }
02840       */
02841 
02842       objective_function_value::fun_without_pen=value(vf);
02843 
02844       //cout << " pen = " << pen << endl;
02845       if (noboundepen_flag==0)
02846       {
02847         vf+=pen;
02848       }
02849       f=value(vf);
02850       if (f<fb)
02851       {
02852         fb=f;
02853         ub=u;
02854       }
02855       gradcalc(usize,g);
02856       //cout << " f = " << setprecision(17) << f << " " << norm(g)
02857        // << " " << norm(u) << endl;
02858     }
02859     u=ub;
02860   }
02861   cout <<  " inner maxg = " <<  fmc1.gmax;
02862   if (fabs(fmc1.gmax)>1.e+3)
02863     trapper();
02864 
02865   if (fabs(fmc1.gmax)>1.e-4)
02866   {
02867     fmc1.itn=0;
02868     //fmc1.crit=1.e-9;
02869     fmc1.ifn=0;
02870     fmc1.ireturn=0;
02871     fmc1.ihang=0;
02872     fmc1.ihflag=0;
02873     fmc1.ialph=0;
02874     initial_params::xinit(u);    // get the initial values into the
02875     //u.initialize();
02876     while (fmc1.ireturn>=0)
02877     {
02878       fmc1.fmin(f,u,g);
02879       if (fmc1.ireturn>0)
02880       {
02881         dvariable vf=0.0;
02882         pen=initial_params::reset(dvar_vector(u));
02883         *objective_function_value::pobjfun=0.0;
02884         pfmin->AD_uf_inner();
02885         if ( no_stuff==0 && quadratic_prior::get_num_quadratic_prior()>0)
02886         {
02887           quadratic_prior::get_M_calculations();
02888         }
02889         vf+=*objective_function_value::pobjfun;
02890         objective_function_value::fun_without_pen=value(vf);
02891 
02892         vf+=pen;
02893         f=value(vf);
02894         if (f<fb)
02895         {
02896           fb=f;
02897           ub=u;
02898         }
02899         gradcalc(usize,g);
02900         //cout << " f = " << setprecision(15) << f << " " << norm(g) << endl;
02901       }
02902     }
02903     u=ub;
02904     cout <<  "  Inner second time = " << fmc1.gmax;
02905   }
02906   cout << "  Inner f = " << fb << endl;
02907   fmc1.ireturn=0;
02908   fmc1.fbest=fb;
02909   gradient_structure::set_NO_DERIVATIVES();
02910   *objective_function_value::pobjfun=0.0;
02911   pfmin->AD_uf_inner();
02912   if ( no_stuff==0 && quadratic_prior::get_num_quadratic_prior()>0)
02913   {
02914     quadratic_prior::get_M_calculations();
02915   }
02916   gradient_structure::set_YES_DERIVATIVES();
02917   pfmin->inner_opt_flag=0;
02918   return u;
02919 }