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