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