ADMB Documentation  11.1.2534
 All Classes Files Functions Variables Typedefs Friends Defines
df1b2lap.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: df1b2lap.cpp 2397 2014-09-20 08:10:20Z 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   nopt=0;
00705 
00706   rseed=3456;
00707   num_importance_samples=0;
00708   if ( (on=option_match(ad_comm::argc,ad_comm::argv,"-is",nopt))>-1)
00709   {
00710     if (!nopt)
00711     {
00712       cerr << "Usage -is option needs positive integer  -- ignored" << endl;
00713     }
00714     else
00715     {
00716       int tht=atoi(ad_comm::argv[on+1]);
00717       if (tht<=0)
00718       {
00719         cerr << "Usage -is option needs non-negative integer  -- ignored.\n";
00720       }
00721       else
00722       {
00723         num_importance_samples=tht;
00724       }
00725       if (nopt==2)
00726       {
00727         int rseed1=atoi(ad_comm::argv[on+2]);
00728         if (rseed1<=0)
00729         {
00730           cerr << "Usage -is option needs non-negative integer  -- ignored.\n";
00731         }
00732         else
00733         {
00734           rseed=rseed1;
00735         }
00736       }
00737     }
00738   }
00739   int use_balanced=0;
00740   if ( (on=option_match(ad_comm::argc,ad_comm::argv,"-isb",nopt))>-1)
00741   {
00742     use_balanced=1;
00743     if (!nopt)
00744     {
00745       cerr << "Usage -isb option needs positive integer  -- ignored" << endl;
00746     }
00747     else
00748     {
00749       int tht=atoi(ad_comm::argv[on+1]);
00750       if (tht<=0)
00751       {
00752         cerr << "Usage -isb option needs non-negative integer  -- ignored.\n";
00753       }
00754       else
00755       {
00756         num_importance_samples=2*tht;
00757       }
00758       if (nopt==2)
00759       {
00760         int rseed1=atoi(ad_comm::argv[on+2]);
00761         if (rseed1<=0)
00762         {
00763           cerr << "Usage -isb option needs non-negative integer  -- ignored.\n";
00764         }
00765         else
00766         {
00767           rseed=rseed1;
00768         }
00769       }
00770     }
00771   }
00772   use_outliers=0;
00773   if ( (on=option_match(ad_comm::argc,ad_comm::argv,"-iso",nopt))>-1)
00774   {
00775     use_outliers=1;
00776   }
00777   if (num_importance_samples)
00778   {
00779     if ( (on=option_match(ad_comm::argc,ad_comm::argv,"-isdiag",nopt))>-1)
00780     {
00781       is_diagnostics_flag=1;
00782     }
00783     if (importance_sampling_values)
00784     {
00785       delete importance_sampling_values;
00786       importance_sampling_values=0;
00787     }
00788     importance_sampling_values =
00789       new dvector(1,num_importance_samples);
00790 
00791     if (importance_sampling_weights)
00792     {
00793       delete importance_sampling_weights;
00794       importance_sampling_weights=0;
00795     }
00796     importance_sampling_weights =
00797       new dvector(1,num_importance_samples);
00798 
00799     random_number_generator rng(rseed);
00800     if (allocated(epsilon)) epsilon.deallocate();
00801     if (use_balanced)
00802     {
00803       // check for eveb num_importance samples
00804       if (num_importance_samples%2)
00805         num_importance_samples+=1;
00806     }
00807     epsilon.allocate(1,num_importance_samples,1,usize);
00808     if (use_balanced)
00809     {
00810       int n2=num_importance_samples/2;
00811       epsilon.sub(1,n2).fill_randn(rng);
00812       if (use_outliers)
00813       {
00814         dmatrix os(1,num_importance_samples,1,usize);
00815         os.fill_randu(rng);
00816         for (i=1;i<=num_importance_samples;i++)
00817         {
00818           for (int j=1;j<=usize;j++)
00819           {
00820             if (os(i,j)<0.05) epsilon(i,j)*=3.0;
00821           }
00822         }
00823       }
00824       for (int i=1;i<=n2;i++)
00825       {
00826         epsilon(i+n2)=-epsilon(i);
00827       }
00828     }
00829     else
00830     {
00831       epsilon.fill_randn(rng);
00832       if (use_outliers)
00833       {
00834         dmatrix os(1,num_importance_samples,1,usize);
00835         os.fill_randu(rng);
00836         for (i=1;i<=num_importance_samples;i++)
00837         {
00838           for (int j=1;j<=usize;j++)
00839           {
00840             if (os(i,j)<0.05) epsilon(i,j)*=3.0;
00841           }
00842         }
00843       }
00844     }
00845 
00846     double eps_mult=1.0;
00847     if ( (on=option_match(ad_comm::argc,ad_comm::argv,"-epsmult",nopt))>-1)
00848     {
00849       if (!nopt)
00850       {
00851         cerr << "Usage -epsmult option needs number  -- ignored" << endl;
00852       }
00853       else
00854       {
00855         istringstream ist(ad_comm::argv[on+1]);
00856         ist >> eps_mult;
00857 
00858         if (eps_mult<=0.0 || eps_mult>1.0)
00859         {
00860           cerr << "Usage -epsmult option needs positive number between 0 and 1 "
00861               "-- ignored" << endl;
00862           eps_mult=1.0;
00863         }
00864       }
00865       for (i=1;i<=num_importance_samples;i++)
00866       {
00867         epsilon(i)*=eps_mult;
00868       }
00869     }
00870   }
00871 
00872   nopt=0;
00873   if ( (on=option_match(ad_comm::argc,ad_comm::argv,"-ht",nopt))>-1)
00874   {
00875     if (!nopt)
00876     {
00877       cerr << "Usage -ht option needs positive integer  -- ignored" << endl;
00878       set_default_hessian_type();
00879     }
00880     else
00881     {
00882       int tht=atoi(ad_comm::argv[on+1]);
00883       if (tht<=0)
00884       {
00885         cerr << "Usage -ht option needs non-negative integer  -- ignored.\n";
00886         set_default_hessian_type();
00887       }
00888       else
00889       {
00890         have_users_hesstype=1;
00891         hesstype=tht;
00892       }
00893       if (nopt>=2)
00894       {
00895         int tbw=atoi(ad_comm::argv[on+2]);
00896         if (tbw<=0)
00897         {
00898           cerr << "Usage -ht option needs non-negative bandwidth"
00899              "  -- ignored" << endl;
00900         }
00901         else
00902         {
00903           ad_comm::bandwidth=tbw;
00904         }
00905       }
00906     }
00907   }
00908   else
00909   {
00910     set_default_hessian_type();
00911   }
00912 
00913   // !! need to check nvar calculation
00914   nvar=maxder(1)-minder(1)+1;
00915 
00916   switch (hesstype)
00917   {
00918   case 0:
00919   case 1:
00920   case 4:
00921     grad.allocate(1,usize);
00922     Hess.allocate(1,usize,1,usize);
00923     Hessadjoint.allocate(1,usize,1,usize);
00924     Dux.allocate(1,usize,1,xsize);
00925     break;
00926   case 3:
00927     {
00928       int bw=2;
00929       if (ad_comm::bandwidth) bw=ad_comm::bandwidth;
00930       bHess=new banded_symmetric_dmatrix(1,usize,bw);
00931       bHessadjoint=new banded_symmetric_dmatrix(1,usize,bw);
00932       grad.allocate(1,usize);
00933       Dux.allocate(1,usize,1,xsize);
00934     }
00935     break;
00936   default:
00937     break;
00938   }
00939 
00940   step.allocate(1,usize);
00941   // !!! nov 12
00942   f1b2gradlist = new df1b2_gradlist(4000000U,200000U,8000000U,400000U,
00943     2000000U,100000U,adstring("f1b2list1"));
00944 
00945   if (hesstype==2)
00946   {
00947     ad_dstar::allocate(nvar);
00948     global_nvar=nvar;
00949     df1b2variable::set_nvar(nvar);
00950     df1b2variable::set_minder(minder(1));
00951     df1b2variable::set_maxder(maxder(1));
00952     df1b2variable::set_blocksize();
00953   }
00954   if (hesstype!=2)
00955   {
00956     if (sparse_hessian_flag==0)
00957     {
00958       ad_dstar::allocate(nvar);
00959       global_nvar=nvar;
00960       df1b2variable::set_nvar(nvar);
00961       df1b2variable::set_minder(minder(1));
00962       df1b2variable::set_maxder(maxder(1));
00963       df1b2variable::set_blocksize();
00964       y.allocate(1,nvariables);
00965     }
00966     else
00967     {
00968       int nsave=nvar;
00969       nvar=1;
00970       ad_dstar::allocate(nvar);
00971       global_nvar=nvar;
00972       df1b2variable::set_nvar(1);
00973       df1b2variable::set_minder(minder(1));
00974       df1b2variable::set_maxder(maxder(1));
00975       df1b2variable::set_blocksize();
00976       y.allocate(1,nvariables);
00977       nvar=nsave;
00978       global_nvar=nvar;
00979     }
00980   }
00981   if (df1b2variable::adpool_counter !=0)
00982   {
00983     cout << "this can't happen" << endl;
00984     ad_exit(1);
00985   }
00986   df1b2variable::adpool_vector[df1b2variable::adpool_counter]=
00987     df1b2variable::pool;
00988   df1b2variable::nvar_vector[df1b2variable::adpool_counter]=
00989       nvariables;
00990   //df1b2variable::adpool_counter++;
00991   df1b2variable::increment_adpool_counter();
00992 }
00993 
00998 laplace_approximation_calculator::laplace_approximation_calculator(
00999   int _xsize,
01000   int _usize,
01001   ivector _minder,
01002   ivector _maxder,
01003   function_minimizer* _pmin
01004 ):
01005   separable_call_level(1),
01006   triplet_information(0),
01007   compressed_triplet_information(0),
01008   nested_separable_calls_counter(1,10),
01009   nested_tree_position(1,5),
01010   nr_debug(0),
01011   pmin(_pmin),
01012   xsize(_xsize),
01013   usize(_usize),
01014   bHess_pd_flag(0),
01015   sparse_triplet(0),
01016   sparse_iterator(0),
01017   sparse_count(0),
01018   sparse_count_adjoint(0),
01019   sparse_triplet2(0),
01020   vsparse_triplet(0),
01021   vsparse_triplet_adjoint(0),
01022   sparse_symbolic(0),
01023   sparse_symbolic2(0),
01024   fmc1(0),
01025   fmc(_xsize),
01026   xadjoint(1,_xsize),
01027   check_local_uadjoint(1,_usize),
01028   check_local_uadjoint2(1,_usize),
01029   check_local_xadjoint(1,_xsize),
01030   check_local_xadjoint2(1,_xsize),
01031   uadjoint(1,_usize),
01032   uhat(1,_usize)
01033 {
01034   nested_tree_position.initialize();
01035   nested_separable_calls_counter.initialize();
01036   //hesstype=1;
01037   uhat.initialize();
01038   //var_flag=1;
01039   var_flag=0;
01040   calling_set=0;
01041   antiepsilon=0;
01042   num_separable_calls=0;
01043   num_local_re_array=0;
01044   num_local_fixed_array=0;
01045   isfunnel_flag=0;
01046   antiflag=0;
01047   rseed=3457;
01048   nfunnelblocks=0;
01049   block_diagonal_hessian=0;
01050   block_diagonal_Dux=0;
01051   block_diagonal_re_list=0;
01052   block_diagonal_fe_list=0;
01053   block_diagonal_vhessian=0;
01054   block_diagonal_vhessianadjoint=0;
01055   pHess_non_quadprior_part=0;
01056   bw=0;
01057   bHess=0;
01058   grad_x_u=0;
01059   grad_x=0;
01060   have_users_hesstype=0;
01061   int mmin=_minder.indexmin();
01062   int mmax=_minder.indexmax();
01063   num_der_blocks= mmax-mmin+1;
01064   minder.allocate(mmin,mmax);
01065   maxder.allocate(mmin,mmax);
01066   minder=_minder;
01067   maxder=_maxder;
01068   fmc1.iprint=inner_iprint;
01069   fmc1.crit=1.e-3;
01070   nvariables=xsize+usize;
01071   Hess.allocate(1,usize,1,usize);
01072   for (int i=1;i<=num_der_blocks;i++)
01073   {
01074     if (minder(i)<1 || maxder(i) > nvariables || maxder(i) < minder(i))
01075     {
01076       cerr << " minder or maxder value out of bounds in"
01077         "laplace_approximation_calculator constructor "
01078        << endl << " values are minder = " << minder
01079        << " maxder = " << maxder << endl;
01080       ad_exit(1);
01081     }
01082   }
01083   nvar=maxder(1)-minder(1)+1;
01084   Hessadjoint.allocate(1,usize,1,usize);
01085   Dux.allocate(1,usize,1,xsize);
01086   // !!! nov 12
01087   f1b2gradlist = new df1b2_gradlist(4000000U,200000U,8000000U,400000U,
01088     2000000U,100000U,adstring("f1b2list1"));
01089   ad_dstar::allocate(nvar);
01090   global_nvar=nvar;
01091   df1b2variable::set_nvar(nvar);
01092   df1b2variable::set_minder(minder(1));
01093   df1b2variable::set_maxder(maxder(1));
01094   df1b2variable::set_blocksize();
01095   y.allocate(1,nvariables);
01096 }
01097 
01102 laplace_approximation_calculator::~laplace_approximation_calculator()
01103 {
01104   if(importance_sampling_weights)
01105   {
01106     delete importance_sampling_weights;
01107     importance_sampling_weights = 0;
01108   }
01109   if(importance_sampling_components)
01110   {
01111     delete importance_sampling_components;
01112     importance_sampling_components = 0;
01113   }
01114   if(importance_sampling_values)
01115   {
01116     delete importance_sampling_values;
01117     importance_sampling_values = 0;
01118   }
01119   if (block_diagonal_vch)
01120   {
01121     delete block_diagonal_vch;
01122     block_diagonal_vch=0;
01123   }
01124   if (block_diagonal_ch)
01125   {
01126     delete block_diagonal_ch;
01127     block_diagonal_ch=0;
01128   }
01129   if (block_diagonal_hessian)
01130   {
01131     delete block_diagonal_hessian;
01132     block_diagonal_hessian=0;
01133   }
01134   if (block_diagonal_Dux )
01135   {
01136     delete block_diagonal_Dux;
01137     block_diagonal_Dux =0;
01138   }
01139   if (block_diagonal_re_list )
01140   {
01141     delete block_diagonal_re_list;
01142     block_diagonal_re_list =0;
01143   }
01144   if (block_diagonal_fe_list )
01145   {
01146     delete block_diagonal_fe_list;
01147     block_diagonal_fe_list =0;
01148   }
01149   if (block_diagonal_vhessian )
01150   {
01151     delete block_diagonal_vhessian;
01152     block_diagonal_vhessian =0;
01153   }
01154   if (block_diagonal_vhessianadjoint )
01155   {
01156     delete block_diagonal_vhessianadjoint;
01157     block_diagonal_vhessianadjoint =0;
01158   }
01159   if (separable_function_difference)
01160   {
01161     delete separable_function_difference;
01162     separable_function_difference=0;
01163   }
01164 
01165   if (derindex)
01166   {
01167     delete derindex;
01168     derindex=0;
01169   }
01170 
01171   if (pHess_non_quadprior_part)
01172   {
01173     delete pHess_non_quadprior_part;
01174     pHess_non_quadprior_part=0;
01175   }
01176 
01177   if (triplet_information)
01178   {
01179     delete triplet_information;
01180     triplet_information=0;
01181   }
01182 
01183   if (bHessadjoint)
01184   {
01185     delete bHessadjoint;
01186     bHessadjoint=0;
01187   }
01188   if (grad_x)
01189   {
01190     delete grad_x;
01191     grad_x=0;
01192   }
01193   if (grad_x_u)
01194   {
01195     delete grad_x_u;
01196     grad_x_u=0;
01197   }
01198   if (bHess)
01199   {
01200     delete bHess;
01201     bHess=0;
01202   }
01203   if (f1b2gradlist)
01204   {
01205     df1b2_gradlist::set_no_derivatives();
01206     delete f1b2gradlist;
01207     f1b2gradlist=0;
01208   }
01209   if (df1b2variable::pool)
01210   {
01211     //df1b2variable::pool->set_size(-1);
01212   }
01213   if (vsparse_triplet_adjoint)
01214   {
01215     delete vsparse_triplet_adjoint;
01216     vsparse_triplet_adjoint=0;
01217   }
01218   if (vsparse_triplet)
01219   {
01220     delete vsparse_triplet;
01221     vsparse_triplet=0;
01222   }
01223   if (sparse_triplet2)
01224   {
01225     delete sparse_triplet2;
01226     sparse_triplet2=0;
01227   }
01228   if (sparse_triplet)
01229   {
01230     delete sparse_triplet;
01231     sparse_triplet=0;
01232   }
01233   if (sparse_symbolic)
01234   {
01235     delete sparse_symbolic;
01236     sparse_symbolic=0;
01237   }
01238   if (sparse_symbolic2)
01239   {
01240     delete sparse_symbolic2;
01241     sparse_symbolic2=0;
01242   }
01243 }
01244 
01249 dvector laplace_approximation_calculator::operator()
01250   (const dvector& _x, const double& _f, function_minimizer * pfmin)
01251 {
01252   dvector g;
01253 
01254 #if defined(USE_ADPVM)
01255   if (pfmin->test_trust_flag)
01256   {
01257     return test_trust_region_method(_x,_f,pfmin);
01258   }
01259   if (ad_comm::pvm_manager)
01260   {
01261     switch (ad_comm::pvm_manager->mode)
01262     {
01263     case 1:
01264       return default_calculations_parallel_master(_x,_f,pfmin);
01265       break;
01266     case 2:
01267     {
01268       dvector g(1,1);
01269       default_calculations_parallel_slave(_x,_f,pfmin);
01270       return g;
01271       break;
01272     }
01273     default:
01274       cerr << "illegal value for mode " << endl;
01275       ad_exit(1);
01276     }
01277   }
01278   else
01279 #endif  //# if defined(USE_ADPVM)
01280 
01281   {
01282     //hesstype=1;
01283     //cout << hesstype << endl;
01284     switch (hesstype)
01285     {
01286       case 1:
01287       {
01288         int check_der_flag=0;
01289         int on=-1;
01290         int nopt=0;
01291         if ((on=option_match(ad_comm::argc,ad_comm::argv,"-chkder",nopt))>-1)
01292         {
01293           check_der_flag=1;
01294         }
01295         if (check_der_flag==1)
01296         {
01297           g = default_calculations_check_derivatives(_x,pfmin,_f);
01298         }
01299         else
01300         {
01301           g = default_calculations(_x,_f,pfmin);
01302         }
01303         break;
01304       }
01305       case 2:
01306       {
01307         // Hessian for the random effects is block diagonal
01308         g = block_diagonal_calculations(_x,_f,pfmin);
01309         break;
01310       }
01311       case 3:
01312       case 4:  // not banded but full hessian
01313       {
01314         // Hessian for the random effects is block diagonal
01315         if (laplace_approximation_calculator::variance_components_vector)
01316         {
01317           g = banded_calculations_lme(_x,_f,pfmin);
01318         }
01319         else
01320         {
01321           g = banded_calculations(_x,_f,pfmin);
01322         }
01323         break;
01324       }
01325       default:
01326       {
01327         cerr << "illegal value for hesstype " << endl;
01328         ad_exit(1);
01329       }
01330     }
01331   }
01332 
01333   return g;
01334 }
01335 
01336 void random_effects_userfunction(double f,const dvector& x, const dvector& g);
01337 
01342 void get_second_ders(int xs,int us,const init_df1b2vector _y,dmatrix& Hess,
01343   dmatrix& Dux, df1b2_gradlist * f1b2gradlist,function_minimizer * pfmin,
01344   laplace_approximation_calculator * lpc)
01345 {
01346   // Note that xs is the number of active non random effects
01347   // parameters
01348   // us is the number of random effects parameters
01349   int j;
01350   int i;
01351   ADUNCONST(init_df1b2vector,y)
01352   int num_der_blocks=lpc->num_der_blocks;
01353   int xsize=lpc->xsize;
01354   int usize=lpc->usize;
01355 
01356   for (int ip=1;ip<=num_der_blocks;ip++)
01357   {
01358     df1b2variable::minder=lpc->minder(ip);
01359     df1b2variable::maxder=lpc->maxder(ip);
01360     lpc->set_u_dot(ip);
01361     df1b2_gradlist::set_yes_derivatives();
01362     (*re_objective_function_value::pobjfun)=0;
01363     df1b2variable pen=0.0;
01364     df1b2variable zz=0.0;
01365     initial_params::straight_through_flag=0;
01366     //initial_params::straight_through_flag=1;
01367     initial_df1b2params::reset(y,pen);
01368     initial_params::straight_through_flag=0;
01369     if (initial_df1b2params::separable_flag)
01370     {
01371       initial_df1b2params::separable_calculation_type=2;
01372       Hess.initialize();
01373       Dux.initialize();
01374     }
01375 
01376     //cout << "2D" << endl;
01377     pfmin->user_function();
01378 
01379     //pfmin->user_function(y,zz);
01380     (*re_objective_function_value::pobjfun)+=pen;
01381     (*re_objective_function_value::pobjfun)+=zz;
01382 
01383     if (!initial_df1b2params::separable_flag)
01384     {
01385       set_dependent_variable(*re_objective_function_value::pobjfun);
01386       df1b2_gradlist::set_no_derivatives();
01387       df1b2variable::passnumber=1;
01388       df1b2_gradcalc1();
01389 
01390       int mind=y(1).minder;
01391       int jmin=max(mind,xsize+1);
01392       int jmax=min(y(1).maxder,xsize+usize);
01393       for (i=1;i<=usize;i++)
01394         for (j=jmin;j<=jmax;j++)
01395           Hess(i,j-xsize)=y(i+xsize).u_bar[j-mind];
01396 
01397       jmin=max(mind,1);
01398       jmax=min(y(1).maxder,xsize);
01399       for (i=1;i<=usize;i++)
01400         for (j=jmin;j<=jmax;j++)
01401           Dux(i,j)=y(i+xsize).u_bar[j-1];
01402     }
01403     if (ip<num_der_blocks)
01404       f1b2gradlist->reset();
01405   }
01406 
01407   if  (ad_comm::print_hess_and_exit_flag)
01408   {
01409     cout << "Hess" << endl;
01410     cout << Hess << endl;
01411     ad_exit(0);
01412   }
01413   //cout << "Dux" << endl;
01414   //cout << setprecision(16) << Dux << endl;
01415 }
01416 
01421 double calculate_laplace_approximation(const dvector& x,const dvector& u0,
01422   const dmatrix& _Hess,const dvector& _xadjoint,const dvector& _uadjoint,
01423   const dmatrix& _Hessadjoint,function_minimizer * pmin)
01424 {
01425   ADUNCONST(dvector,xadjoint)
01426   ADUNCONST(dvector,uadjoint)
01427   ADUNCONST(dmatrix,Hessadjoint)
01428   ADUNCONST(dmatrix,Hess)
01429   const int xs=x.size();
01430   const int us=u0.size();
01431   gradient_structure::set_YES_DERIVATIVES();
01432   int nvar;
01433   //dcompressed_triplet & lst2 = *(pmin->lapprox->sparse_triplet);
01434   //hs_symbolic & ssymb=*(pmin->lapprox->sparse_symbolic);
01435   //dcompressed_triplet & xxxt = *(pmin->lapprox->sparse_triplet);
01436   dcompressed_triplet & lst = *(pmin->lapprox->sparse_triplet2);
01437   hs_symbolic & ssymb=*(pmin->lapprox->sparse_symbolic2);
01438   {
01439   /*
01440     cout << norm2(make_dmatrix(lst)-make_dmatrix(xxxt)) << endl;
01441     ofstream ofs1("m1");
01442     ofs1 << setfixed() << setprecision(3) << setw(10)
01443          << make_dmatrix(lst) << endl;
01444     ofstream ofs2("m2");
01445     ofs2 << setfixed() << setprecision(3) << setw(10)
01446          << make_dmatrix(xxxt) << endl;
01447    */
01448   }
01449 
01450   if (pmin->lapprox->sparse_hessian_flag==0)
01451   {
01452     nvar=x.size()+u0.size()+u0.size()*u0.size();
01453   }
01454   else
01455   {
01456     int sz= lst.indexmax()-lst.indexmin()+1;
01457     nvar=x.size()+u0.size()+sz;
01458   }
01459   independent_variables y(1,nvar);
01460 
01461   // need to set random effects active together with whatever
01462   // init parameters should be active in this phase
01463   initial_params::set_inactive_only_random_effects();
01464   initial_params::set_active_random_effects();
01465   /*int onvar=*/initial_params::nvarcalc();
01466   initial_params::xinit(y);    // get the initial values into the
01467   y(1,xs)=x;
01468 
01469   int i,j;
01470   dvar_vector d(1,xs+us);
01471 
01472   dmatrix Hess_save;
01473   // contribution for quadratic prior
01474   if (quadratic_prior::get_num_quadratic_prior()>0)
01475   {
01476     if (allocated(Hess_save)) Hess_save.deallocate();
01477     int mmin=Hess.indexmin();
01478     int mmax=Hess.indexmax();
01479     Hess_save.allocate(mmin,mmax,mmin,mmax);
01480     Hess_save=Hess;
01481     int & vxs = (int&)(xs);
01482     quadratic_prior::get_cHessian_contribution(Hess,vxs);
01483   }
01484  // Here need hooks for sparse matrix structures
01485   int ii=xs+us+1;
01486   if (pmin->lapprox->sparse_hessian_flag==0)
01487   {
01488     for (i=1;i<=us;i++)
01489       for (j=1;j<=us;j++)
01490       y(ii++)=Hess(i,j);
01491   }
01492   else
01493   {
01494     int smin=lst.indexmin();
01495     int smax=lst.indexmax();
01496     for (i=smin;i<=smax;i++)
01497       y(ii++)=lst(i);
01498   }
01499 
01500   if (quadratic_prior::get_num_quadratic_prior()>0)
01501   {
01502     Hess=Hess_save;
01503   }
01504 
01505   dvar_vector vy=dvar_vector(y);
01506   initial_params::stddev_vscale(d,vy);
01507   dvar_matrix vHess;
01508   if (pmin->lapprox->sparse_hessian_flag==0)
01509   {
01510     vHess.allocate(1,us,1,us);
01511   }
01512   initial_params::reset(vy);    // get the initial values into the
01513   ii=xs+us+1;
01514   if (initial_df1b2params::have_bounded_random_effects)
01515   {
01516     for (i=1;i<=us;i++)
01517     {
01518       if (d(i+xs)>0.0)
01519       {
01520         for (j=1;j<=us;j++)
01521         {
01522           if (d(j+xs)>0.0)
01523             vHess(i,j)=vy(ii++)/(d(i+xs)*d(j+xs));
01524           else
01525             vHess(i,j)=vy(ii++)/d(i+xs);
01526         }
01527       }
01528       else
01529       {
01530         for (j=1;j<=us;j++)
01531         {
01532           if (d(j+xs)>0.0)
01533             vHess(i,j)=vy(ii++)/d(j+xs);
01534           else
01535             vHess(i,j)=vy(ii++);
01536         }
01537       }
01538     }
01539   }
01540   else
01541   {
01542     if (laplace_approximation_calculator::sparse_hessian_flag==0)
01543     {
01544       for (i=1;i<=us;i++)
01545       {
01546         for (j=1;j<=us;j++)
01547         {
01548           vHess(i,j)=vy(ii++);
01549         }
01550       }
01551     }
01552     else
01553     {
01554       int mmin=lst.indexmin();
01555       int mmax=lst.indexmax();
01556       dvar_compressed_triplet * vsparse_triplet =
01557         pmin->lapprox->vsparse_triplet;
01558 
01559       if (vsparse_triplet==0)
01560       {
01561         pmin->lapprox->vsparse_triplet=
01562           new dvar_compressed_triplet(mmin,mmax,us,us);
01563         vsparse_triplet = pmin->lapprox->vsparse_triplet;
01564         for (i=mmin;i<=mmax;i++)
01565         {
01566           (*vsparse_triplet)(1,i)=lst(1,i);
01567           (*vsparse_triplet)(2,i)=lst(2,i);
01568         }
01569       }
01570       else
01571       {
01572         if (!allocated(*vsparse_triplet))
01573         {
01574           (*vsparse_triplet).allocate(mmin,mmax,us,us);
01575           for (i=mmin;i<=mmax;i++)
01576           {
01577             (*vsparse_triplet)(1,i)=lst(1,i);
01578             (*vsparse_triplet)(2,i)=lst(2,i);
01579           }
01580         }
01581       }
01582       dcompressed_triplet * vsparse_triplet_adjoint =
01583         pmin->lapprox->vsparse_triplet_adjoint;
01584 
01585       if (vsparse_triplet_adjoint==0)
01586       {
01587         pmin->lapprox->vsparse_triplet_adjoint=
01588           new dcompressed_triplet(mmin,mmax,us,us);
01589         vsparse_triplet_adjoint = pmin->lapprox->vsparse_triplet_adjoint;
01590         for (i=mmin;i<=mmax;i++)
01591         {
01592           (*vsparse_triplet_adjoint)(1,i)=lst(1,i);
01593           (*vsparse_triplet_adjoint)(2,i)=lst(2,i);
01594         }
01595       }
01596       else
01597       {
01598         if (!allocated(*vsparse_triplet_adjoint))
01599         {
01600           (*vsparse_triplet_adjoint).allocate(mmin,mmax,us,us);
01601           for (i=mmin;i<=mmax;i++)
01602           {
01603             (*vsparse_triplet_adjoint)(1,i)=lst(1,i);
01604             (*vsparse_triplet_adjoint)(2,i)=lst(2,i);
01605           }
01606         }
01607       }
01608       vsparse_triplet->get_x()=vy(ii,ii+mmax-mmin).shift(1);
01609     }
01610   }
01611    dvariable vf=0.0;
01612 
01613    *objective_function_value::pobjfun=0.0;
01614    pmin->AD_uf_outer();
01615    if ( no_stuff==0 && quadratic_prior::get_num_quadratic_prior()>0)
01616    {
01617      quadratic_prior::get_M_calculations();
01618    }
01619    vf+=*objective_function_value::pobjfun;
01620    //cout << setprecision(15) << vf << endl;
01621   // *********************************************
01622   // *********************************************
01623   // *********************************************
01624   // dvector tmpg(1,nvar);
01625   // tmpg.initialize();
01626   // gradcalc(nvar,tmpg);
01627   // *********************************************
01628   // *********************************************
01629   // *********************************************
01630 
01631    int sgn=0;
01632    dvariable ld = 0;
01633    if (ad_comm::no_ln_det_choleski_flag)
01634    {
01635      if(laplace_approximation_calculator::saddlepointflag==0)
01636      {
01637        ld=0.5*ln_det(vHess,sgn);
01638      }
01639      else
01640      {
01641        ld=0.5*ln_det(-vHess,sgn);
01642      }
01643    }
01644    else
01645    {
01646      if(laplace_approximation_calculator::saddlepointflag==0)
01647      {
01648        int ierr=0;
01649        if (pmin->lapprox->sparse_hessian_flag==0)
01650        {
01651          ld=0.5*ln_det_choleski_error(vHess,ierr);
01652          if (ierr==1)
01653          {
01654            ofstream ofs("hessian.diag");
01655            ofs << vHess << endl;
01656            ofs <<  x << endl;
01657            ofs <<  u0 << endl;
01658            ofs << "Matrix not positive definite in Ln_det_choleski"
01659                << endl;
01660            ad_exit(1);
01661          }
01662        }
01663        else
01664        {
01665          //double ld1=0.5*ln_det(*(pmin->lapprox->sparse_triplet),
01666          //  *(pmin->lapprox->sparse_symbolic));
01667 
01668          if (write_sparse_flag)
01669          {
01670            //ofstream ofs("sparse");
01671            //ofs << *(pmin->lapprox->vsparse_triplet) << endl;
01672          }
01673          ld=0.5*ln_det(*(pmin->lapprox->vsparse_triplet),
01674            ssymb,*(pmin->lapprox->sparse_triplet2));
01675            //*(pmin->lapprox->sparse_symbolic2),pmin->lapprox);
01676          //cout << ld-ld1 << endl;
01677        }
01678      }
01679      else
01680      {
01681        if (pmin->lapprox->sparse_hessian_flag==0)
01682        {
01683          ld=0.5*ln_det_choleski(-vHess);
01684        }
01685        else
01686        {
01687          cerr << "need to fix this" << endl;
01688          ad_exit(1);
01689          //ld+=ln_det(-*(pmin->lapprox->vsparse_triplet));
01690        }
01691      }
01692    }
01693 
01694    int ps1=0;
01695    if (ps1)
01696    {
01697      dmatrix cHess=value(vHess);
01698      cout << " ln_det = " << ld << " ";
01699      dvector eig=eigenvalues(cHess);
01700      dmatrix r(1,2,1,eig.indexmax());
01701      dvector num(1,eig.indexmax());
01702      num.fill_seqadd(1,1);
01703      r(1)=eig;
01704      r(2)=num;
01705      dmatrix tt=trans(r);
01706      dmatrix ss=trans(sort(tt,1));
01707      cout << ss(1,1)  << " " << ss(1,eig.indexmax()) << " " ;
01708      int nnn=(int)ss(2,1);
01709      dmatrix d=eigenvectors(cHess);
01710      //cout << " d " << d(nnn) << " " << d(nnn)*cHess*d(nnn) << endl;
01711      dmatrix t=trans(d);
01712      r(1)=t(nnn);
01713      r(2)=num;
01714      dmatrix tt2=trans(r);
01715      dmatrix ss2=trans(sort(tt2,1));
01716 
01717      cout << " t " << setprecision(3) << ss2(1)(1,5) << " --- "
01718           << t(nnn)*cHess*t(nnn) << endl;
01719      cout << "   " << setprecision(3) << ss2(2)(1,5) << endl;
01720      //cout << " t " << t(1) << " " << t(1)*cHess*t(2) << endl;
01721    }
01722 
01723    int nx=0;
01724    if (nx==0)
01725    {
01726      vf+=ld;
01727    }
01728    double f=value(vf);
01729    f-=us*0.5*log(2.0*PI);
01730    dvector g(1,nvar);
01731    gradcalc(nvar,g);
01732 
01733   ii=1;
01734   for (i=1;i<=xs;i++)
01735     xadjoint(i)=g(ii++);
01736   for (i=1;i<=us;i++)
01737     uadjoint(i)=g(ii++);
01738   if (pmin->lapprox->sparse_hessian_flag==0)
01739   {
01740     if (allocated(Hessadjoint))
01741     {
01742       for (i=1;i<=us;i++)
01743         for (j=1;j<=us;j++)
01744           Hessadjoint(i,j)=g(ii++);
01745     }
01746   }
01747   else
01748   {
01749     dcompressed_triplet * vsparse_triplet_adjoint =
01750       pmin->lapprox->vsparse_triplet_adjoint;
01751 
01752     int smin=lst.indexmin();
01753     int smax=lst.indexmax();
01754     for (i=smin;i<=smax;i++)
01755     {
01756       (*vsparse_triplet_adjoint)(i)=g(ii);
01757       ii++;
01758     }
01759   }
01760 
01761    /*
01762    if (quadratic_prior::get_num_quadratic_prior()>0)
01763    {
01764      // *******************************************************
01765      // *******************************************************
01766      // for quadratic prior option
01767      // temporarily get der wrt x of x ---> ln_det(F_uu + inv(A(x)))
01768      int nvar=x.size()+u0.size();
01769      independent_variables y(1,nvar);
01770      initial_params::xinit(y);    // get the initial values into the
01771      y(1,xs)=x;
01772      dvar_vector vy=dvar_vector(y);
01773      initial_params::reset(vy);    // get the initial values into the
01774 
01775      pmin->AD_uf_outer();
01776 
01777      dvar_matrix tmpH=quadratic_prior::get_vHessian_contribution();
01778 
01779      tmpH+=Hess;
01780      dvariable ld;
01781      ld=0.5*ln_det(tmpH,sgn);
01782      dvector g(1,nvar);
01783      gradcalc(nvar,g);
01784      int ii=1;
01785      double checksum=0.0;
01786      for (i=1;i<=xs;i++)
01787        xadjoint(i)+=g(ii++);
01788      for (i=1;i<=us;i++)
01789        checksum+=g(ii++);
01790 
01791      if (fabs(checksum)> 1.e-10)
01792      {
01793        cerr << "checksum too big " << checksum << endl;
01794        ad_exit(1);
01795      }
01796 
01797    }
01798   */
01799 
01800 
01801    // *******************************************************
01802    // *******************************************************
01803 
01804 
01805   return f;
01806 }
01807 
01812 void get_newton_raphson_info(int xs,int us,const init_df1b2vector _y,
01813   dmatrix& Hess, dvector& grad, df1b2_gradlist * f1b2gradlist,
01814   function_minimizer * pfmin)
01815 {
01816   // Note that xs is the number of active non random effects
01817   // parameters
01818   // us is the number of random effects parameters
01819   int j;
01820   int i;
01821   ADUNCONST(init_df1b2vector,y)
01822   {
01823     df1b2_gradlist::set_yes_derivatives();
01824     //cout << re_objective_function_value::pobjfun << endl;
01825     //cout << re_objective_function_value::pobjfun->ptr << endl;
01826     (*re_objective_function_value::pobjfun)=0;
01827     df1b2variable pen=0.0;
01828     df1b2variable zz=0.0;
01829     initial_df1b2params::reset(y,pen);
01830     pfmin->user_function();
01831     //pfmin->user_function(y,zz);
01832     (*re_objective_function_value::pobjfun)+=pen;
01833     (*re_objective_function_value::pobjfun)+=zz;
01834     set_dependent_variable(*re_objective_function_value::pobjfun);
01835     df1b2_gradlist::set_no_derivatives();
01836     df1b2variable::passnumber=1;
01837   }
01838 
01839   int mind=y(1).minder;
01840   int jmin=max(mind,xs+1);
01841   int jmax=min(y(1).maxder,xs+us);
01842   if (!initial_df1b2params::separable_flag)
01843   {
01844     df1b2_gradcalc1();
01845     for (i=1;i<=us;i++)
01846       for (j=jmin;j<=jmax;j++)
01847         Hess(i,j-xs)=y(i+xs).u_bar[j-mind];
01848     for (j=jmin;j<=jmax;j++)
01849       grad(j-xs)= re_objective_function_value::pobjfun->u_dot[j-mind];
01850   }
01851 }
01852 
01856 void laplace_approximation_calculator::check_for_need_to_reallocate(int ip)
01857 {
01858 }
01859   //void laplace_approximation_calculator::get_newton_raphson_info
01860   //  (function_minimizer * pfmin)
01861   //{
01862   //  int i,j,ip;
01863   //
01864   //  for (ip=1;ip<=num_der_blocks;ip++)
01865   //  {
01866   //    // do we need to reallocate memory for df1b2variables?
01867   //    check_for_need_to_reallocate(ip);
01868   //    df1b2_gradlist::set_yes_derivatives();
01869   //    //cout << re_objective_function_value::pobjfun << endl;
01870   //    //cout << re_objective_function_value::pobjfun->ptr << endl;
01871   //    (*re_objective_function_value::pobjfun)=0;
01872   //    df1b2variable pen=0.0;
01873   //    df1b2variable zz=0.0;
01874   //    initial_df1b2params::reset(y,pen);
01875   //    if (initial_df1b2params::separable_flag)
01876   //    {
01877   //      Hess.initialize();
01878   //      grad.initialize();
01879   //    }
01880   //    pfmin->user_function();
01881   //    //pfmin->user_function(y,zz);
01882   //    (*re_objective_function_value::pobjfun)+=pen;
01883   //    (*re_objective_function_value::pobjfun)+=zz;
01884   //
01885   //    if (!initial_df1b2params::separable_flag)
01886   //    {
01887   //      set_dependent_variable(*re_objective_function_value::pobjfun);
01888   //      df1b2_gradlist::set_no_derivatives();
01889   //      df1b2variable::passnumber=1;
01890   //      df1b2_gradcalc1();
01891   //      int mind=y(1).minder;
01892   //      int jmin=max(mind,xsize+1);
01893   //      int jmax=min(y(1).maxder,xsize+usize);
01894   //      for (i=1;i<=usize;i++)
01895   //        for (j=jmin;j<=jmax;j++)
01896   //          Hess(i,j-xsize)=y(i+xsize).u_bar[j-mind];
01897   //
01898   //      for (j=jmin;j<=jmax;j++)
01899   //        grad(j-xsize)= re_objective_function_value::pobjfun->u_dot[j-mind];
01900   //    }
01901   //  }
01902   //  {
01903   //    //ofstream ofs("hess.tmp");
01904   //    //ofs << Hess << endl;
01905   //    //ofs << grad << endl;
01906   //  }
01907   //}
01908   //
01909 
01914 double evaluate_function(const dvector& x,function_minimizer * pfmin)
01915 {
01916   int usize=initial_params::nvarcalc();
01917   //double f=0.0;
01918   dvector g(1,usize);
01919   independent_variables u(1,usize);
01920   u=x;
01921   dvariable vf=0.0;
01922   vf=initial_params::reset(dvar_vector(u));
01923   //vf=0.0;
01924   *objective_function_value::pobjfun=0.0;
01925   laplace_approximation_calculator::where_are_we_flag=2;
01926   pfmin->AD_uf_inner();
01927   if ( no_stuff==0 && quadratic_prior::get_num_quadratic_prior()>0)
01928   {
01929     quadratic_prior::get_M_calculations();
01930   }
01931   vf+=*objective_function_value::pobjfun;
01932   laplace_approximation_calculator::where_are_we_flag=0;
01933   initial_df1b2params::cobjfun=value(vf);
01934   gradcalc(usize,g);
01935   double maxg=max(fabs(g));
01936   if (!initial_params::mc_phase)
01937   {
01938     cout << setprecision(16) << " f = " << vf
01939          << " max g = " << maxg << endl;
01940   }
01941   return maxg;
01942 }
01943 
01948 double evaluate_function(double& fval,const dvector& x,
01949   function_minimizer* pfmin)
01950 {
01951   int usize=initial_params::nvarcalc();
01952   //double f=0.0;
01953   dvector g(1,usize);
01954   independent_variables u(1,usize);
01955   u=x;
01956   dvariable vf=0.0;
01957   vf=initial_params::reset(dvar_vector(u));
01958   //vf=0.0;
01959   *objective_function_value::pobjfun=0.0;
01960   laplace_approximation_calculator::where_are_we_flag=2;
01961   pfmin->AD_uf_inner();
01962   if ( no_stuff==0 && quadratic_prior::get_num_quadratic_prior()>0)
01963   {
01964     quadratic_prior::get_M_calculations();
01965   }
01966   vf+=*objective_function_value::pobjfun;
01967   laplace_approximation_calculator::where_are_we_flag=0;
01968   initial_df1b2params::cobjfun=value(vf);
01969   gradcalc(usize,g);
01970   double maxg=max(fabs(g));
01971   fval=value(vf);
01972   if (!initial_params::mc_phase)
01973   {
01974     cout << setprecision(10) << " f = " << vf
01975          << " max g = " << maxg << endl;
01976   }
01977   return maxg;
01978 }
01979 
01984 double evaluate_function(double& fval,const dvector& x,const dvector& g,
01985   function_minimizer * pfmin)
01986 {
01987   int usize=initial_params::nvarcalc();
01988   //double f=0.0;
01989   //dvector g(1,usize);
01990   independent_variables u(1,usize);
01991   u=x;
01992   dvariable vf=0.0;
01993   vf=initial_params::reset(dvar_vector(u));
01994   //vf=0.0;
01995   *objective_function_value::pobjfun=0.0;
01996   laplace_approximation_calculator::where_are_we_flag=2;
01997   pfmin->AD_uf_inner();
01998   if ( no_stuff==0 && quadratic_prior::get_num_quadratic_prior()>0)
01999   {
02000     quadratic_prior::get_M_calculations();
02001   }
02002   vf+=*objective_function_value::pobjfun;
02003   laplace_approximation_calculator::where_are_we_flag=0;
02004   initial_df1b2params::cobjfun=value(vf);
02005   gradcalc(usize,g);
02006   double maxg=max(fabs(g));
02007   fval=value(vf);
02008   if (!initial_params::mc_phase)
02009   {
02010     cout << setprecision(15) << " f = " << vf
02011          << " max g = " << maxg << endl;
02012   }
02013   return maxg;
02014 }
02015 
02020 double evaluate_function_quiet(const dvector& x,function_minimizer * pfmin)
02021 {
02022   int usize=initial_params::nvarcalc();
02023   //double f=0.0;
02024   dvector g(1,usize);
02025   independent_variables u(1,usize);
02026   u=x;
02027   dvariable vf=0.0;
02028   vf=initial_params::reset(dvar_vector(u));
02029   //vf=0.0;
02030   *objective_function_value::pobjfun=0.0;
02031   laplace_approximation_calculator::where_are_we_flag=2;
02032   pfmin->AD_uf_inner();
02033   if ( no_stuff==0 && quadratic_prior::get_num_quadratic_prior()>0)
02034   {
02035     quadratic_prior::get_M_calculations();
02036   }
02037   vf+=*objective_function_value::pobjfun;
02038   laplace_approximation_calculator::where_are_we_flag=0;
02039   initial_df1b2params::cobjfun=value(vf);
02040   gradcalc(usize,g);
02041   double maxg=max(fabs(g));
02042   return maxg;
02043 }
02044 
02049 void evaluate_function_gradient(double& f,const dvector& x,
02050   function_minimizer * pfmin,dvector& xadjoint,dvector& uadjoint)
02051 {
02052   int usize=initial_params::nvarcalc();
02053   dvector g(1,usize);
02054   independent_variables u(1,usize);
02055   u=x;
02056   dvariable vf=0.0;
02057   vf=initial_params::reset(dvar_vector(u));
02058   //vf=0.0;
02059   *objective_function_value::pobjfun=0.0;
02060   pfmin->AD_uf_inner();
02061   if ( no_stuff==0 && quadratic_prior::get_num_quadratic_prior()>0)
02062   {
02063     quadratic_prior::get_M_calculations();
02064   }
02065   vf+=*objective_function_value::pobjfun;
02066   initial_df1b2params::cobjfun=value(vf);
02067   f=value(vf);
02068   gradcalc(usize,g);
02069   int xsize=xadjoint.indexmax();
02070   int us=uadjoint.indexmax();
02071   xadjoint=g(1,xsize);
02072   uadjoint=g(xsize+1,xsize+us).shift(1);
02073 }
02074 
02079 double evaluate_function_no_derivatives(const dvector& x,
02080   function_minimizer* pfmin)
02081 {
02082   double fval;
02083   gradient_structure::set_NO_DERIVATIVES();
02084   int usize=initial_params::nvarcalc();
02085   //double f=0.0;
02086   dvector g(1,usize);
02087   independent_variables u(1,usize);
02088   u=x;
02089   dvariable vf=0.0;
02090   vf=initial_params::reset(dvar_vector(u));
02091   //vf=0.0;
02092   *objective_function_value::pobjfun=0.0;
02093   pfmin->AD_uf_inner();
02094   if ( no_stuff==0 && quadratic_prior::get_num_quadratic_prior()>0)
02095   {
02096     quadratic_prior::get_M_calculations();
02097   }
02098   vf+=*objective_function_value::pobjfun;
02099   initial_df1b2params::cobjfun=value(vf);
02100   fval=value(vf);
02101   gradient_structure::set_YES_DERIVATIVES();
02102   return fval;
02103 }
02104 
02109 void cleanup_laplace_stuff(laplace_approximation_calculator * l)
02110 {
02111   if (l)
02112   {
02113     delete l;
02114     l=0;
02115     df1b2variable::pool->deallocate();
02116 
02117     for (int i=0;i<df1b2variable::adpool_counter;i++)
02118     {
02119       delete df1b2variable::adpool_vector[i];
02120       df1b2variable::adpool_vector[i]=0;
02121       df1b2variable::nvar_vector[i]=0;
02122       df1b2variable::adpool_counter=0;
02123     }
02124   }
02125 }
02126 
02131 gauss_hermite_stuff::gauss_hermite_stuff
02132   (laplace_approximation_calculator * lapprox,
02133   int use_gauss_hermite,int num_separable_calls ,const ivector& itmp) :
02134   x(1,use_gauss_hermite), w(1,use_gauss_hermite), mi(0)
02135 {
02136   // for now we must have the same number of random effects in each
02137   // separable function call"
02138   int i;
02139   for (i=2;i<=num_separable_calls;i++)
02140   {
02141     if (itmp(i)!=itmp(i-1))
02142     {
02143       cerr << " At present for the adaptive gauss_hermite must have the same"
02144            << endl
02145            << " number of random effects in each separable function call"
02146            << endl;
02147       ad_exit(1);
02148     }
02149   }
02150   for (i=1;i<=num_separable_calls;i++)
02151   {
02152     if (itmp(i)>1)
02153     {
02154       lapprox->multi_random_effects=
02155         max(lapprox->multi_random_effects,itmp(i));
02156      /*
02157       cerr << "At present gauss-hermite is only implemented for"
02158         " one random effect per separable function call "
02159        << endl;
02160       ad_exit(1);
02161      */
02162     }
02163   }
02164   if (allocated(gauss_hermite_values))
02165     gauss_hermite_values.deallocate();
02166   if (lapprox->multi_random_effects==0)
02167   {
02168     gauss_hermite_values.allocate(1,num_separable_calls,1,use_gauss_hermite);
02169   }
02170   else
02171   {
02172    ivector indx=pow(use_gauss_hermite,itmp);
02173     gauss_hermite_values.allocate(1,num_separable_calls,1,indx);
02174   }
02175 
02176   normalized_gauss_hermite(x,w);
02177   is=0;
02178 }
02179 
02184 void laplace_approximation_calculator::set_default_hessian_type(void )
02185 {
02186   //if (function_minimizer::hesstype==0)
02187   {
02188     if (initial_df1b2params::separable_flag)
02189     {
02190      // have SEPARABLE_FUNCTION but no object of type quadratic_penalty or
02191      // normal prior -- can't tell if we should use hesstype 2 or 3
02192      // until we run
02193       hesstype=3;
02194     }
02195     else
02196     {
02197       hesstype=1;  // assume block diagonal
02198     }
02199   }
02200  /*
02201   else
02202   {
02203     hesstype=function_minimizer::hesstype;
02204   }
02205  */
02206 }
02207 
02212 double laplace_approximation_calculator::get_fx_fu(function_minimizer * pfmin)
02213 {
02214   initial_params::set_inactive_only_random_effects();
02215   initial_params::set_active_random_effects();
02216 
02217   if (grad_x==0)
02218   {
02219     grad_x=new dvector(1,xsize);
02220   }
02221   grad_x->initialize();
02222 
02223   if (grad_x_u==0)
02224   {
02225     grad_x_u=new dvector(1,xsize+usize);
02226   }
02227   pfmin->inner_opt_flag=0;
02228   independent_variables u(1,xsize+usize);
02229   //gradcalc(0,*grad_x_u);
02230   initial_params::xinit(u);    // get the initial values into the
02231 
02232   dvariable pen=0.0;
02233   dvariable vf=0.0;
02234   pen=initial_params::reset(dvar_vector(u));
02235   *objective_function_value::pobjfun=0.0;
02236   pfmin->AD_uf_inner();
02237   vf+=*objective_function_value::pobjfun;
02238   objective_function_value::fun_without_pen=value(vf);
02239   vf+=pen;
02240   gradcalc(xsize+usize,*grad_x_u);
02241   double f=value(vf);
02242   return f;
02243 }
02244 
02249   void laplace_approximation_calculator::begin_separable_call_stuff(void)
02250   {
02251     separable_call_level++;
02252     //build_up_nested_shape();
02253     //clean(nested_tree_position,separable_call_level);
02254     //nested_separable_calls_counter(separable_call_level)++;
02255     //nested_tree_position(separable_call_level)++;
02256   }
02257 
02262   void laplace_approximation_calculator::end_separable_call_stuff(void)
02263   {
02264     //build_up_nested_shape();
02265     //clean(nested_tree_position,separable_call_level);
02266     separable_call_level--;
02267   }
02268 
02273 void laplace_approximation_calculator::build_up_nested_shape(void)
02274 {
02275   int ll;
02276   ivector& a =nested_separable_calls_counter;
02277   int clean_level=0;
02278   switch(separable_call_level)
02279   {
02280   case 1:
02281     ll=1;
02282     if (nested_separable_calls_counter(ll)>0)
02283     {
02284       if (clean_level==0) clean_level=ll+1;
02285       if (nested_separable_calls_counter(ll+1)>0)
02286       {
02287         nested_shape(a(ll))=a(ll+1);
02288       }
02289       else
02290       {
02291         break;
02292       }
02293     }
02294     else
02295     {
02296       break;
02297     }
02298    case 2:
02299     ll=2;
02300     if (nested_separable_calls_counter(ll)>0)
02301     {
02302       if (clean_level==0) clean_level=ll+1;
02303       if (nested_separable_calls_counter(ll+1)>0)
02304       {
02305         nested_shape(a(1),a(2))=a(ll+1);
02306       }
02307       else
02308       {
02309         break;
02310       }
02311     }
02312     else
02313     {
02314       break;
02315     }
02316    case 3:
02317     ll=3;
02318     if (nested_separable_calls_counter(ll)>0)
02319     {
02320       if (clean_level==0) clean_level=ll+1;
02321       if (nested_separable_calls_counter(ll+1)>0)
02322       {
02323         nested_shape(a(1),a(2),a(3))=a(ll+1);
02324       }
02325       else
02326       {
02327         break;
02328       }
02329     }
02330     else
02331     {
02332       break;
02333     }
02334    case 4:
02335     ll=4;
02336     if (nested_separable_calls_counter(ll)>0)
02337     {
02338       if (clean_level==0) clean_level=ll+1;
02339       if (nested_separable_calls_counter(ll+1)>0)
02340       {
02341         nested_shape(a(1),a(2),a(3),a(4))=a(ll+1);
02342       }
02343       else
02344       {
02345         break;
02346       }
02347     }
02348     else
02349     {
02350       break;
02351     }
02352    default:
02353      cerr << "illegal value in " <<
02354        "laplace_approximation_calculator::build_up_nested_shape"
02355        << endl;
02356   }
02357   if (clean_level>0)
02358   {
02359     int mmax=a.indexmax();
02360     for (int i=clean_level;i<=mmax;i++)
02361     {
02362       a(i)=0;
02363     }
02364   }
02365 }
02366 
02371 ivector * nested_calls_shape::get_ptr1(void)
02372 {
02373   return ptr1;
02374 }
02375 
02380 imatrix * nested_calls_shape::get_ptr2(void)
02381 {
02382   return ptr2;
02383 }
02384 
02389 i3_array * nested_calls_shape::get_ptr3(void)
02390 {
02391   return ptr3;
02392 }
02393 
02398 i4_array * nested_calls_shape::get_ptr4(void)
02399 {
02400   return ptr4;
02401 }
02402 
02407 ostream &  operator << (const ostream& _s,const nested_calls_shape& _m)
02408 {
02409   ADUNCONST(nested_calls_shape,m)
02410   ADUNCONST(ofstream,s)
02411   if (m.get_ptr1())
02412     s<< *(m.get_ptr1()) << endl << endl;
02413 
02414   if (m.get_ptr2())
02415     s<< *(m.get_ptr2()) << endl << endl;
02416 
02417   if (m.get_ptr3())
02418     s<< *(m.get_ptr3()) << endl << endl;
02419 
02420   if (m.get_ptr4())
02421     s<< *(m.get_ptr4()) << endl << endl;
02422 
02423   return s;
02424 }
02425 
02430 void nested_calls_shape::trim(void)
02431 {
02432   int mmax1=0;
02433   if (ptr1)
02434   {
02435     int imin=ptr1->indexmin();
02436     int imax=ptr1->indexmax();
02437     for (int i=imin;i<=imax;i++)
02438     {
02439       if ( (*ptr1)(i)==0) break;
02440       mmax1++;
02441     }
02442   }
02443   if (mmax1==0)
02444   {
02445     delete ptr1;
02446     ptr1=0;
02447   }
02448   else
02449   {
02450     ivector * tmp=new ivector(1,mmax1);
02451     (*tmp)=(*ptr1)(1,mmax1);
02452     delete ptr1;
02453     ptr1=tmp;
02454   }
02455   if (ptr2)
02456   {
02457     if (!ptr1)
02458     {
02459       delete ptr2;
02460       ptr2=0;
02461     }
02462     else
02463     {
02464       imatrix * tmp=new imatrix(1,mmax1);
02465       int imin=tmp->indexmin();
02466       int imax=tmp->indexmax();
02467       for (int i=imin;i<=imax;i++)
02468       {
02469         int jmin=(*ptr2)(i).indexmin();
02470         int jmax=(*ptr2)(i).indexmax();
02471         int mmax2=0;
02472         for (int j=jmin;j<=jmax;j++)
02473         {
02474           if ((*ptr2)(i,j)==0) break;
02475           mmax2++;
02476         }
02477         if (mmax2>0)
02478         {
02479           (*tmp)(i).allocate(1,mmax2);
02480           (*tmp)(i)=(*ptr2)(i)(1,mmax2);
02481         }
02482       }
02483       delete ptr2;
02484       ptr2=tmp;
02485     }
02486   }
02487   if (ptr3)
02488   {
02489     cerr << "warning not dealitn with prt3" << endl;
02490     delete ptr3;
02491     ptr3=0;
02492   }
02493 }
02494 
02499 nested_calls_shape::~nested_calls_shape()
02500 {
02501   if (ptr1)
02502   {
02503     delete ptr1;
02504     ptr1=0;
02505   }
02506   if (ptr2)
02507   {
02508     delete ptr2;
02509     ptr2=0;
02510   }
02511   if (ptr3)
02512   {
02513     delete ptr3;
02514     ptr3=0;
02515   }
02516   if (ptr4)
02517   {
02518     delete ptr4;
02519     ptr4=0;
02520   }
02521 }
02522 
02527 void nested_calls_shape::allocate(int n,int m,int p )
02528 {
02529   if (ptr1)
02530   {
02531     delete ptr1;
02532     ptr1=0;
02533   }
02534   ptr1=new ivector(1,n);
02535 
02536   if (ptr2)
02537   {
02538     delete ptr2;
02539     ptr2=0;
02540   }
02541   ptr2=new imatrix(1,n,1,m);
02542   if (ptr3)
02543   {
02544     delete ptr3;
02545     ptr3=0;
02546   }
02547   ptr3=new i3_array(1,n,1,m,1,p);
02548   /*
02549   if (ptr4)
02550   {
02551     delete ptr4;
02552     ptr4=0;
02553   }
02554   */
02555 }
02556 
02561 void nested_calls_shape::allocate(int n,int m)
02562 {
02563   if (ptr1)
02564   {
02565     delete ptr1;
02566     ptr1=0;
02567   }
02568   ptr1=new ivector(1,n);
02569 
02570   if (ptr2)
02571   {
02572     delete ptr2;
02573     ptr2=0;
02574   }
02575   ptr2=new imatrix(1,n,1,m);
02576   /*
02577   if (ptr3)
02578   {
02579     delete ptr3;
02580     ptr3=0;
02581   }
02582   ptr=new i3_array(1,n,1,m,1,p);
02583   if (ptr4)
02584   {
02585     delete ptr4;
02586     ptr4=0;
02587   }
02588   */
02589 }
02590 
02595 void nested_calls_shape::initialize(void)
02596 {
02597   if (ptr1)
02598   {
02599     ptr1->initialize();
02600   }
02601 
02602   if (ptr2)
02603   {
02604     ptr2->initialize();
02605   }
02606 
02607   if (ptr3)
02608   {
02609     ptr3->initialize();
02610   }
02611 
02612   if (ptr4)
02613   {
02614     ptr4->initialize();
02615   }
02616 }
02617 
02622 void nested_calls_shape::allocate(int n)
02623 {
02624   if (ptr1)
02625   {
02626     delete ptr1;
02627     ptr1=0;
02628   }
02629   ptr1=new ivector(1,n);
02630 
02631   /*
02632   if (ptr2)
02633   {
02634     delete ptr2;
02635     ptr2=0;
02636   }
02637   ptr=new imatrix(1,n,1,m);
02638   if (ptr3)
02639   {
02640     delete ptr3;
02641     ptr3=0;
02642   }
02643   ptr=new i3_array(1,n,1,m,1,p);
02644   if (ptr4)
02645   {
02646     delete ptr4;
02647     ptr4=0;
02648   }
02649   */
02650 }
02651 
02656 void nested_calls_indices::allocate(const nested_calls_shape& _nsc)
02657 {
02658   int mmax=0;
02659   ADUNCONST(nested_calls_shape,nsc)
02660   mmax=nsc.get_ptr1()->indexmax();
02661   if (ptr1)
02662   {
02663     delete ptr1;
02664     ptr1=0;
02665   }
02666   if (nsc.get_ptr1())
02667   {
02668     ptr1=new imatrix(1,mmax,1,*nsc.get_ptr1());
02669   }
02670   if (ptr2)
02671   {
02672     delete ptr2;
02673     ptr2=0;
02674   }
02675   if (nsc.get_ptr2())
02676   {
02677     ptr2=new i3_array(1,mmax,1,*nsc.get_ptr1(),1,*nsc.get_ptr2());
02678   }
02679   if (ptr3)
02680   {
02681     delete ptr3;
02682     ptr3=0;
02683   }
02684   if (nsc.get_ptr3())
02685   {
02686     ptr3=new i4_array(1,mmax,1,*nsc.get_ptr1(),1,*nsc.get_ptr2(),
02687       1,*nsc.get_ptr3());
02688   }
02689 }
02690 
02695 dvector laplace_approximation_calculator::get_uhat_lm_newton2
02696   (const dvector& x,function_minimizer * pfmin)
02697 {
02698   //int on,nopt;
02699   pfmin->inner_opt_flag=1;
02700   double f=0.0;
02701   double fb=1.e+100;
02702   dvector g(1,usize);
02703   dvector ub(1,usize);
02704   independent_variables u(1,usize);
02705   gradcalc(0,g);
02706   fmc1.itn=0;
02707   fmc1.ifn=0;
02708   fmc1.ireturn=0;
02709   initial_params::xinit(u);    // get the initial values into the
02710   initial_params::xinit(ubest);    // get the initial values into the
02711   fmc1.ialph=0;
02712   fmc1.ihang=0;
02713   fmc1.ihflag=0;
02714 
02715   if (init_switch)
02716   {
02717     u.initialize();
02718   }
02719   else
02720   {
02721     u=ubest;
02722   }
02723 
02724   fmc1.dfn=1.e-2;
02725   dvariable pen=0.0;
02726   //cout << "starting  norm(u) = " << norm(u) << endl;
02727   while (fmc1.ireturn>=0)
02728   {
02729    /*
02730     double nu=norm(value(u));
02731     if (nu>400)
02732     {
02733       cout << "U norm(u) = " << nu  << endl;
02734     }
02735     cout << "V norm(u) = " << nu
02736          << " f = " << f  << endl;
02737     */
02738     fmc1.fmin(f,u,g);
02739     //cout << "W norm(u) = " << norm(value(u)) << endl;
02740     if (fmc1.ireturn>0)
02741     {
02742       dvariable vf=0.0;
02743       pen=initial_params::reset(dvar_vector(u));
02744       *objective_function_value::pobjfun=0.0;
02745       pfmin->AD_uf_inner();
02746       if (saddlepointflag)
02747       {
02748         *objective_function_value::pobjfun*=-1.0;
02749       }
02750       if ( no_stuff==0 && quadratic_prior::get_num_quadratic_prior()>0)
02751       {
02752         quadratic_prior::get_M_calculations();
02753       }
02754       vf+=*objective_function_value::pobjfun;
02755 
02756      /*  this is now done in the operator = function
02757       if (quadratic_prior::get_num_quadratic_prior()>0)
02758       {
02759         vf+= quadratic_prior::get_quadratic_priors();
02760       }
02761       */
02762 
02763       objective_function_value::fun_without_pen=value(vf);
02764 
02765       //cout << " pen = " << pen << endl;
02766       if (noboundepen_flag==0)
02767       {
02768         vf+=pen;
02769       }
02770       f=value(vf);
02771       if (f<fb)
02772       {
02773         fb=f;
02774         ub=u;
02775       }
02776       gradcalc(usize,g);
02777       //cout << " f = " << setprecision(17) << f << " " << norm(g)
02778        // << " " << norm(u) << endl;
02779     }
02780     u=ub;
02781   }
02782   cout <<  " inner maxg = " <<  fmc1.gmax;
02783   if (fabs(fmc1.gmax)>1.e+3)
02784     trapper();
02785 
02786   if (fabs(fmc1.gmax)>1.e-4)
02787   {
02788     fmc1.itn=0;
02789     //fmc1.crit=1.e-9;
02790     fmc1.ifn=0;
02791     fmc1.ireturn=0;
02792     fmc1.ihang=0;
02793     fmc1.ihflag=0;
02794     fmc1.ialph=0;
02795     initial_params::xinit(u);    // get the initial values into the
02796     //u.initialize();
02797     while (fmc1.ireturn>=0)
02798     {
02799       fmc1.fmin(f,u,g);
02800       if (fmc1.ireturn>0)
02801       {
02802         dvariable vf=0.0;
02803         pen=initial_params::reset(dvar_vector(u));
02804         *objective_function_value::pobjfun=0.0;
02805         pfmin->AD_uf_inner();
02806         if ( no_stuff==0 && quadratic_prior::get_num_quadratic_prior()>0)
02807         {
02808           quadratic_prior::get_M_calculations();
02809         }
02810         vf+=*objective_function_value::pobjfun;
02811         objective_function_value::fun_without_pen=value(vf);
02812 
02813         vf+=pen;
02814         f=value(vf);
02815         if (f<fb)
02816         {
02817           fb=f;
02818           ub=u;
02819         }
02820         gradcalc(usize,g);
02821         //cout << " f = " << setprecision(15) << f << " " << norm(g) << endl;
02822       }
02823     }
02824     u=ub;
02825     cout <<  "  Inner second time = " << fmc1.gmax;
02826   }
02827   cout << "  Inner f = " << fb << endl;
02828   fmc1.ireturn=0;
02829   fmc1.fbest=fb;
02830   gradient_structure::set_NO_DERIVATIVES();
02831   *objective_function_value::pobjfun=0.0;
02832   pfmin->AD_uf_inner();
02833   if ( no_stuff==0 && quadratic_prior::get_num_quadratic_prior()>0)
02834   {
02835     quadratic_prior::get_M_calculations();
02836   }
02837   gradient_structure::set_YES_DERIVATIVES();
02838   pfmin->inner_opt_flag=0;
02839   return u;
02840 }