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