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