ADMB Documentation  11.2.2823
 All Classes Files Functions Variables Typedefs Friends Defines
df1b2lap.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: df1b2lap.cpp 2689 2014-11-18 22:04:38Z jsibert $
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 }
01255 
01260 dvector laplace_approximation_calculator::operator()
01261   (const dvector& _x, const double& _f, function_minimizer * pfmin)
01262 {
01263   dvector g;
01264 
01265 #if defined(USE_ADPVM)
01266   if (pfmin->test_trust_flag)
01267   {
01268     return test_trust_region_method(_x,_f,pfmin);
01269   }
01270   if (ad_comm::pvm_manager)
01271   {
01272     switch (ad_comm::pvm_manager->mode)
01273     {
01274     case 1:
01275       return default_calculations_parallel_master(_x,_f,pfmin);
01276       break;
01277     case 2:
01278     {
01279       dvector g(1,1);
01280       default_calculations_parallel_slave(_x,_f,pfmin);
01281       return g;
01282       break;
01283     }
01284     default:
01285       cerr << "illegal value for mode " << endl;
01286       ad_exit(1);
01287     }
01288   }
01289   else
01290 #endif  //# if defined(USE_ADPVM)
01291 
01292   {
01293     //hesstype=1;
01294     //cout << hesstype << endl;
01295     switch (hesstype)
01296     {
01297       case 1:
01298       {
01299         int check_der_flag=0;
01300         int on=-1;
01301         int nopt=0;
01302         if ((on=option_match(ad_comm::argc,ad_comm::argv,"-chkder",nopt))>-1)
01303         {
01304           check_der_flag=1;
01305         }
01306         if (check_der_flag==1)
01307         {
01308           g = default_calculations_check_derivatives(_x,pfmin,_f);
01309         }
01310         else
01311         {
01312           g = default_calculations(_x,_f,pfmin);
01313         }
01314         break;
01315       }
01316       case 2:
01317       {
01318         // Hessian for the random effects is block diagonal
01319         g = block_diagonal_calculations(_x,_f,pfmin);
01320         break;
01321       }
01322       case 3:
01323       case 4:  // not banded but full hessian
01324       {
01325         // Hessian for the random effects is block diagonal
01326         if (laplace_approximation_calculator::variance_components_vector)
01327         {
01328           g = banded_calculations_lme(_x,_f,pfmin);
01329         }
01330         else
01331         {
01332           g = banded_calculations(_x,_f,pfmin);
01333         }
01334         break;
01335       }
01336       default:
01337       {
01338         cerr << "illegal value for hesstype " << endl;
01339         ad_exit(1);
01340       }
01341     }
01342   }
01343 
01344   return g;
01345 }
01346 
01347 void random_effects_userfunction(double f,const dvector& x, const dvector& g);
01348 
01353 void get_second_ders(int xs,int us,const init_df1b2vector _y,dmatrix& Hess,
01354   dmatrix& Dux, df1b2_gradlist * f1b2gradlist,function_minimizer * pfmin,
01355   laplace_approximation_calculator * lpc)
01356 {
01357   // Note that xs is the number of active non random effects
01358   // parameters
01359   // us is the number of random effects parameters
01360   int j;
01361   int i;
01362   ADUNCONST(init_df1b2vector,y)
01363   int num_der_blocks=lpc->num_der_blocks;
01364   int xsize=lpc->xsize;
01365   int usize=lpc->usize;
01366 
01367   for (int ip=1;ip<=num_der_blocks;ip++)
01368   {
01369     df1b2variable::minder=lpc->minder(ip);
01370     df1b2variable::maxder=lpc->maxder(ip);
01371     lpc->set_u_dot(ip);
01372     df1b2_gradlist::set_yes_derivatives();
01373     (*re_objective_function_value::pobjfun)=0;
01374     df1b2variable pen=0.0;
01375     df1b2variable zz=0.0;
01376     initial_params::straight_through_flag=0;
01377     //initial_params::straight_through_flag=1;
01378     initial_df1b2params::reset(y,pen);
01379     initial_params::straight_through_flag=0;
01380     if (initial_df1b2params::separable_flag)
01381     {
01382       initial_df1b2params::separable_calculation_type=2;
01383       Hess.initialize();
01384       Dux.initialize();
01385     }
01386 
01387     //cout << "2D" << endl;
01388     pfmin->user_function();
01389 
01390     //pfmin->user_function(y,zz);
01391     (*re_objective_function_value::pobjfun)+=pen;
01392     (*re_objective_function_value::pobjfun)+=zz;
01393 
01394     if (!initial_df1b2params::separable_flag)
01395     {
01396       set_dependent_variable(*re_objective_function_value::pobjfun);
01397       df1b2_gradlist::set_no_derivatives();
01398       df1b2variable::passnumber=1;
01399       df1b2_gradcalc1();
01400 
01401       int mind=y(1).minder;
01402       int jmin=max(mind,xsize+1);
01403       int jmax=min(y(1).maxder,xsize+usize);
01404       for (i=1;i<=usize;i++)
01405         for (j=jmin;j<=jmax;j++)
01406           Hess(i,j-xsize)=y(i+xsize).u_bar[j-mind];
01407 
01408       jmin=max(mind,1);
01409       jmax=min(y(1).maxder,xsize);
01410       for (i=1;i<=usize;i++)
01411         for (j=jmin;j<=jmax;j++)
01412           Dux(i,j)=y(i+xsize).u_bar[j-1];
01413     }
01414     if (ip<num_der_blocks)
01415       f1b2gradlist->reset();
01416   }
01417 
01418   if  (ad_comm::print_hess_and_exit_flag)
01419   {
01420     cout << "Hess" << endl;
01421     cout << Hess << endl;
01422     ad_exit(0);
01423   }
01424   //cout << "Dux" << endl;
01425   //cout << setprecision(16) << Dux << endl;
01426 }
01427 
01432 double calculate_laplace_approximation(const dvector& x,const dvector& u0,
01433   const dmatrix& _Hess,const dvector& _xadjoint,const dvector& _uadjoint,
01434   const dmatrix& _Hessadjoint,function_minimizer * pmin)
01435 {
01436   ADUNCONST(dvector,xadjoint)
01437   ADUNCONST(dvector,uadjoint)
01438   ADUNCONST(dmatrix,Hessadjoint)
01439   ADUNCONST(dmatrix,Hess)
01440   const int xs=x.size();
01441   const int us=u0.size();
01442   gradient_structure::set_YES_DERIVATIVES();
01443   int nvar;
01444   //dcompressed_triplet & lst2 = *(pmin->lapprox->sparse_triplet);
01445   //hs_symbolic & ssymb=*(pmin->lapprox->sparse_symbolic);
01446   //dcompressed_triplet & xxxt = *(pmin->lapprox->sparse_triplet);
01447   dcompressed_triplet & lst = *(pmin->lapprox->sparse_triplet2);
01448   hs_symbolic & ssymb=*(pmin->lapprox->sparse_symbolic2);
01449   {
01450   /*
01451     cout << norm2(make_dmatrix(lst)-make_dmatrix(xxxt)) << endl;
01452     ofstream ofs1("m1");
01453     ofs1 << setfixed() << setprecision(3) << setw(10)
01454          << make_dmatrix(lst) << endl;
01455     ofstream ofs2("m2");
01456     ofs2 << setfixed() << setprecision(3) << setw(10)
01457          << make_dmatrix(xxxt) << endl;
01458    */
01459   }
01460 
01461   if (pmin->lapprox->sparse_hessian_flag==0)
01462   {
01463     nvar=x.size()+u0.size()+u0.size()*u0.size();
01464   }
01465   else
01466   {
01467     int sz= lst.indexmax()-lst.indexmin()+1;
01468     nvar=x.size()+u0.size()+sz;
01469   }
01470   independent_variables y(1,nvar);
01471 
01472   // need to set random effects active together with whatever
01473   // init parameters should be active in this phase
01474   initial_params::set_inactive_only_random_effects();
01475   initial_params::set_active_random_effects();
01476   /*int onvar=*/initial_params::nvarcalc();
01477   initial_params::xinit(y);    // get the initial values into the
01478   y(1,xs)=x;
01479 
01480   int i,j;
01481   dvar_vector d(1,xs+us);
01482 
01483   dmatrix Hess_save;
01484   // contribution for quadratic prior
01485   if (quadratic_prior::get_num_quadratic_prior()>0)
01486   {
01487     if (allocated(Hess_save)) Hess_save.deallocate();
01488     int mmin=Hess.indexmin();
01489     int mmax=Hess.indexmax();
01490     Hess_save.allocate(mmin,mmax,mmin,mmax);
01491     Hess_save=Hess;
01492     int & vxs = (int&)(xs);
01493     quadratic_prior::get_cHessian_contribution(Hess,vxs);
01494   }
01495  // Here need hooks for sparse matrix structures
01496   int ii=xs+us+1;
01497   if (pmin->lapprox->sparse_hessian_flag==0)
01498   {
01499     for (i=1;i<=us;i++)
01500       for (j=1;j<=us;j++)
01501       y(ii++)=Hess(i,j);
01502   }
01503   else
01504   {
01505     int smin=lst.indexmin();
01506     int smax=lst.indexmax();
01507     for (i=smin;i<=smax;i++)
01508       y(ii++)=lst(i);
01509   }
01510 
01511   if (quadratic_prior::get_num_quadratic_prior()>0)
01512   {
01513     Hess=Hess_save;
01514   }
01515 
01516   dvar_vector vy=dvar_vector(y);
01517   initial_params::stddev_vscale(d,vy);
01518   dvar_matrix vHess;
01519   if (pmin->lapprox->sparse_hessian_flag==0)
01520   {
01521     vHess.allocate(1,us,1,us);
01522   }
01523   initial_params::reset(vy);    // get the initial values into the
01524   ii=xs+us+1;
01525   if (initial_df1b2params::have_bounded_random_effects)
01526   {
01527     for (i=1;i<=us;i++)
01528     {
01529       if (d(i+xs)>0.0)
01530       {
01531         for (j=1;j<=us;j++)
01532         {
01533           if (d(j+xs)>0.0)
01534             vHess(i,j)=vy(ii++)/(d(i+xs)*d(j+xs));
01535           else
01536             vHess(i,j)=vy(ii++)/d(i+xs);
01537         }
01538       }
01539       else
01540       {
01541         for (j=1;j<=us;j++)
01542         {
01543           if (d(j+xs)>0.0)
01544             vHess(i,j)=vy(ii++)/d(j+xs);
01545           else
01546             vHess(i,j)=vy(ii++);
01547         }
01548       }
01549     }
01550   }
01551   else
01552   {
01553     if (laplace_approximation_calculator::sparse_hessian_flag==0)
01554     {
01555       for (i=1;i<=us;i++)
01556       {
01557         for (j=1;j<=us;j++)
01558         {
01559           vHess(i,j)=vy(ii++);
01560         }
01561       }
01562     }
01563     else
01564     {
01565       int mmin=lst.indexmin();
01566       int mmax=lst.indexmax();
01567       dvar_compressed_triplet * vsparse_triplet =
01568         pmin->lapprox->vsparse_triplet;
01569 
01570       if (vsparse_triplet==0)
01571       {
01572         pmin->lapprox->vsparse_triplet=
01573           new dvar_compressed_triplet(mmin,mmax,us,us);
01574         vsparse_triplet = pmin->lapprox->vsparse_triplet;
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       else
01582       {
01583         if (!allocated(*vsparse_triplet))
01584         {
01585           (*vsparse_triplet).allocate(mmin,mmax,us,us);
01586           for (i=mmin;i<=mmax;i++)
01587           {
01588             (*vsparse_triplet)(1,i)=lst(1,i);
01589             (*vsparse_triplet)(2,i)=lst(2,i);
01590           }
01591         }
01592       }
01593       dcompressed_triplet * vsparse_triplet_adjoint =
01594         pmin->lapprox->vsparse_triplet_adjoint;
01595 
01596       if (vsparse_triplet_adjoint==0)
01597       {
01598         pmin->lapprox->vsparse_triplet_adjoint=
01599           new dcompressed_triplet(mmin,mmax,us,us);
01600         vsparse_triplet_adjoint = pmin->lapprox->vsparse_triplet_adjoint;
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       else
01608       {
01609         if (!allocated(*vsparse_triplet_adjoint))
01610         {
01611           (*vsparse_triplet_adjoint).allocate(mmin,mmax,us,us);
01612           for (i=mmin;i<=mmax;i++)
01613           {
01614             (*vsparse_triplet_adjoint)(1,i)=lst(1,i);
01615             (*vsparse_triplet_adjoint)(2,i)=lst(2,i);
01616           }
01617         }
01618       }
01619       vsparse_triplet->get_x()=vy(ii,ii+mmax-mmin).shift(1);
01620     }
01621   }
01622    dvariable vf=0.0;
01623 
01624    *objective_function_value::pobjfun=0.0;
01625    pmin->AD_uf_outer();
01626    if ( no_stuff==0 && quadratic_prior::get_num_quadratic_prior()>0)
01627    {
01628      quadratic_prior::get_M_calculations();
01629    }
01630    vf+=*objective_function_value::pobjfun;
01631    //cout << setprecision(15) << vf << endl;
01632   // *********************************************
01633   // *********************************************
01634   // *********************************************
01635   // dvector tmpg(1,nvar);
01636   // tmpg.initialize();
01637   // gradcalc(nvar,tmpg);
01638   // *********************************************
01639   // *********************************************
01640   // *********************************************
01641 
01642    int sgn=0;
01643    dvariable ld = 0;
01644    if (ad_comm::no_ln_det_choleski_flag)
01645    {
01646      if(laplace_approximation_calculator::saddlepointflag==0)
01647      {
01648        ld=0.5*ln_det(vHess,sgn);
01649      }
01650      else
01651      {
01652        ld=0.5*ln_det(-vHess,sgn);
01653      }
01654    }
01655    else
01656    {
01657      if(laplace_approximation_calculator::saddlepointflag==0)
01658      {
01659        int ierr=0;
01660        if (pmin->lapprox->sparse_hessian_flag==0)
01661        {
01662          ld=0.5*ln_det_choleski_error(vHess,ierr);
01663          if (ierr==1)
01664          {
01665            ofstream ofs("hessian.diag");
01666            ofs << vHess << endl;
01667            ofs <<  x << endl;
01668            ofs <<  u0 << endl;
01669            ofs << "Matrix not positive definite in Ln_det_choleski"
01670                << endl;
01671            cerr << "Matrix not positive definite in Ln_det_choleski\n"
01672                 << "see file hessian.diag for details"
01673                << endl;
01674            ad_exit(1);
01675          }
01676        }
01677        else
01678        {
01679          //double ld1=0.5*ln_det(*(pmin->lapprox->sparse_triplet),
01680          //  *(pmin->lapprox->sparse_symbolic));
01681 
01682          if (write_sparse_flag)
01683          {
01684            //ofstream ofs("sparse");
01685            //ofs << *(pmin->lapprox->vsparse_triplet) << endl;
01686          }
01687          ld=0.5*ln_det(*(pmin->lapprox->vsparse_triplet),
01688            ssymb,*(pmin->lapprox->sparse_triplet2));
01689            //*(pmin->lapprox->sparse_symbolic2),pmin->lapprox);
01690          //cout << ld-ld1 << endl;
01691        }
01692      }
01693      else
01694      {
01695        if (pmin->lapprox->sparse_hessian_flag==0)
01696        {
01697          ld=0.5*ln_det_choleski(-vHess);
01698        }
01699        else
01700        {
01701          cerr << "need to fix this" << endl;
01702          ad_exit(1);
01703          //ld+=ln_det(-*(pmin->lapprox->vsparse_triplet));
01704        }
01705      }
01706    }
01707 
01708    int ps1=0;
01709    if (ps1)
01710    {
01711      dmatrix cHess=value(vHess);
01712      cout << " ln_det = " << ld << " ";
01713      dvector eig=eigenvalues(cHess);
01714      dmatrix r(1,2,1,eig.indexmax());
01715      dvector num(1,eig.indexmax());
01716      num.fill_seqadd(1,1);
01717      r(1)=eig;
01718      r(2)=num;
01719      dmatrix tt=trans(r);
01720      dmatrix ss=trans(sort(tt,1));
01721      cout << ss(1,1)  << " " << ss(1,eig.indexmax()) << " " ;
01722      int nnn=(int)ss(2,1);
01723      dmatrix d=eigenvectors(cHess);
01724      //cout << " d " << d(nnn) << " " << d(nnn)*cHess*d(nnn) << endl;
01725      dmatrix t=trans(d);
01726      r(1)=t(nnn);
01727      r(2)=num;
01728      dmatrix tt2=trans(r);
01729      dmatrix ss2=trans(sort(tt2,1));
01730 
01731      cout << " t " << setprecision(3) << ss2(1)(1,5) << " --- "
01732           << t(nnn)*cHess*t(nnn) << endl;
01733      cout << "   " << setprecision(3) << ss2(2)(1,5) << endl;
01734      //cout << " t " << t(1) << " " << t(1)*cHess*t(2) << endl;
01735    }
01736 
01737    int nx=0;
01738    if (nx==0)
01739    {
01740      vf+=ld;
01741    }
01742    double f=value(vf);
01743    f-=us*0.5*log(2.0*PI);
01744    dvector g(1,nvar);
01745    gradcalc(nvar,g);
01746 
01747   ii=1;
01748   for (i=1;i<=xs;i++)
01749     xadjoint(i)=g(ii++);
01750   for (i=1;i<=us;i++)
01751     uadjoint(i)=g(ii++);
01752   if (pmin->lapprox->sparse_hessian_flag==0)
01753   {
01754     if (allocated(Hessadjoint))
01755     {
01756       for (i=1;i<=us;i++)
01757         for (j=1;j<=us;j++)
01758           Hessadjoint(i,j)=g(ii++);
01759     }
01760   }
01761   else
01762   {
01763     dcompressed_triplet * vsparse_triplet_adjoint =
01764       pmin->lapprox->vsparse_triplet_adjoint;
01765 
01766     int smin=lst.indexmin();
01767     int smax=lst.indexmax();
01768     for (i=smin;i<=smax;i++)
01769     {
01770       (*vsparse_triplet_adjoint)(i)=g(ii);
01771       ii++;
01772     }
01773   }
01774 
01775    /*
01776    if (quadratic_prior::get_num_quadratic_prior()>0)
01777    {
01778      // *******************************************************
01779      // *******************************************************
01780      // for quadratic prior option
01781      // temporarily get der wrt x of x ---> ln_det(F_uu + inv(A(x)))
01782      int nvar=x.size()+u0.size();
01783      independent_variables y(1,nvar);
01784      initial_params::xinit(y);    // get the initial values into the
01785      y(1,xs)=x;
01786      dvar_vector vy=dvar_vector(y);
01787      initial_params::reset(vy);    // get the initial values into the
01788 
01789      pmin->AD_uf_outer();
01790 
01791      dvar_matrix tmpH=quadratic_prior::get_vHessian_contribution();
01792 
01793      tmpH+=Hess;
01794      dvariable ld;
01795      ld=0.5*ln_det(tmpH,sgn);
01796      dvector g(1,nvar);
01797      gradcalc(nvar,g);
01798      int ii=1;
01799      double checksum=0.0;
01800      for (i=1;i<=xs;i++)
01801        xadjoint(i)+=g(ii++);
01802      for (i=1;i<=us;i++)
01803        checksum+=g(ii++);
01804 
01805      if (fabs(checksum)> 1.e-10)
01806      {
01807        cerr << "checksum too big " << checksum << endl;
01808        ad_exit(1);
01809      }
01810 
01811    }
01812   */
01813 
01814 
01815    // *******************************************************
01816    // *******************************************************
01817 
01818 
01819   return f;
01820 }
01821 
01826 void get_newton_raphson_info(int xs,int us,const init_df1b2vector _y,
01827   dmatrix& Hess, dvector& grad, df1b2_gradlist * f1b2gradlist,
01828   function_minimizer * pfmin)
01829 {
01830   // Note that xs is the number of active non random effects
01831   // parameters
01832   // us is the number of random effects parameters
01833   int j;
01834   int i;
01835   ADUNCONST(init_df1b2vector,y)
01836   {
01837     df1b2_gradlist::set_yes_derivatives();
01838     //cout << re_objective_function_value::pobjfun << endl;
01839     //cout << re_objective_function_value::pobjfun->ptr << endl;
01840     (*re_objective_function_value::pobjfun)=0;
01841     df1b2variable pen=0.0;
01842     df1b2variable zz=0.0;
01843     initial_df1b2params::reset(y,pen);
01844     pfmin->user_function();
01845     //pfmin->user_function(y,zz);
01846     (*re_objective_function_value::pobjfun)+=pen;
01847     (*re_objective_function_value::pobjfun)+=zz;
01848     set_dependent_variable(*re_objective_function_value::pobjfun);
01849     df1b2_gradlist::set_no_derivatives();
01850     df1b2variable::passnumber=1;
01851   }
01852 
01853   int mind=y(1).minder;
01854   int jmin=max(mind,xs+1);
01855   int jmax=min(y(1).maxder,xs+us);
01856   if (!initial_df1b2params::separable_flag)
01857   {
01858     df1b2_gradcalc1();
01859     for (i=1;i<=us;i++)
01860       for (j=jmin;j<=jmax;j++)
01861         Hess(i,j-xs)=y(i+xs).u_bar[j-mind];
01862     for (j=jmin;j<=jmax;j++)
01863       grad(j-xs)= re_objective_function_value::pobjfun->u_dot[j-mind];
01864   }
01865 }
01866 
01870 void laplace_approximation_calculator::check_for_need_to_reallocate(int ip)
01871 {
01872 }
01873   //void laplace_approximation_calculator::get_newton_raphson_info
01874   //  (function_minimizer * pfmin)
01875   //{
01876   //  int i,j,ip;
01877   //
01878   //  for (ip=1;ip<=num_der_blocks;ip++)
01879   //  {
01880   //    // do we need to reallocate memory for df1b2variables?
01881   //    check_for_need_to_reallocate(ip);
01882   //    df1b2_gradlist::set_yes_derivatives();
01883   //    //cout << re_objective_function_value::pobjfun << endl;
01884   //    //cout << re_objective_function_value::pobjfun->ptr << endl;
01885   //    (*re_objective_function_value::pobjfun)=0;
01886   //    df1b2variable pen=0.0;
01887   //    df1b2variable zz=0.0;
01888   //    initial_df1b2params::reset(y,pen);
01889   //    if (initial_df1b2params::separable_flag)
01890   //    {
01891   //      Hess.initialize();
01892   //      grad.initialize();
01893   //    }
01894   //    pfmin->user_function();
01895   //    //pfmin->user_function(y,zz);
01896   //    (*re_objective_function_value::pobjfun)+=pen;
01897   //    (*re_objective_function_value::pobjfun)+=zz;
01898   //
01899   //    if (!initial_df1b2params::separable_flag)
01900   //    {
01901   //      set_dependent_variable(*re_objective_function_value::pobjfun);
01902   //      df1b2_gradlist::set_no_derivatives();
01903   //      df1b2variable::passnumber=1;
01904   //      df1b2_gradcalc1();
01905   //      int mind=y(1).minder;
01906   //      int jmin=max(mind,xsize+1);
01907   //      int jmax=min(y(1).maxder,xsize+usize);
01908   //      for (i=1;i<=usize;i++)
01909   //        for (j=jmin;j<=jmax;j++)
01910   //          Hess(i,j-xsize)=y(i+xsize).u_bar[j-mind];
01911   //
01912   //      for (j=jmin;j<=jmax;j++)
01913   //        grad(j-xsize)= re_objective_function_value::pobjfun->u_dot[j-mind];
01914   //    }
01915   //  }
01916   //  {
01917   //    //ofstream ofs("hess.tmp");
01918   //    //ofs << Hess << endl;
01919   //    //ofs << grad << endl;
01920   //  }
01921   //}
01922   //
01923 
01928 double evaluate_function(const dvector& x,function_minimizer * pfmin)
01929 {
01930   int usize=initial_params::nvarcalc();
01931   //double f=0.0;
01932   dvector g(1,usize);
01933   independent_variables u(1,usize);
01934   u=x;
01935   dvariable vf=0.0;
01936   vf=initial_params::reset(dvar_vector(u));
01937   //vf=0.0;
01938   *objective_function_value::pobjfun=0.0;
01939   laplace_approximation_calculator::where_are_we_flag=2;
01940   pfmin->AD_uf_inner();
01941   if ( no_stuff==0 && quadratic_prior::get_num_quadratic_prior()>0)
01942   {
01943     quadratic_prior::get_M_calculations();
01944   }
01945   vf+=*objective_function_value::pobjfun;
01946   laplace_approximation_calculator::where_are_we_flag=0;
01947   initial_df1b2params::cobjfun=value(vf);
01948   gradcalc(usize,g);
01949   double maxg=max(fabs(g));
01950   if (!initial_params::mc_phase)
01951   {
01952     cout << setprecision(16) << " f = " << vf
01953          << " max g = " << maxg << endl;
01954   }
01955   return maxg;
01956 }
01957 
01962 double evaluate_function(double& fval,const dvector& x,
01963   function_minimizer* pfmin)
01964 {
01965   int usize=initial_params::nvarcalc();
01966   //double f=0.0;
01967   dvector g(1,usize);
01968   independent_variables u(1,usize);
01969   u=x;
01970   dvariable vf=0.0;
01971   vf=initial_params::reset(dvar_vector(u));
01972   //vf=0.0;
01973   *objective_function_value::pobjfun=0.0;
01974   laplace_approximation_calculator::where_are_we_flag=2;
01975   pfmin->AD_uf_inner();
01976   if ( no_stuff==0 && quadratic_prior::get_num_quadratic_prior()>0)
01977   {
01978     quadratic_prior::get_M_calculations();
01979   }
01980   vf+=*objective_function_value::pobjfun;
01981   laplace_approximation_calculator::where_are_we_flag=0;
01982   initial_df1b2params::cobjfun=value(vf);
01983   gradcalc(usize,g);
01984   double maxg=max(fabs(g));
01985   fval=value(vf);
01986   if (!initial_params::mc_phase)
01987   {
01988     cout << setprecision(10) << " f = " << vf
01989          << " max g = " << maxg << endl;
01990   }
01991   return maxg;
01992 }
01993 
01998 double evaluate_function(double& fval,const dvector& x,const dvector& g,
01999   function_minimizer * pfmin)
02000 {
02001   int usize=initial_params::nvarcalc();
02002   //double f=0.0;
02003   //dvector g(1,usize);
02004   independent_variables u(1,usize);
02005   u=x;
02006   dvariable vf=0.0;
02007   vf=initial_params::reset(dvar_vector(u));
02008   //vf=0.0;
02009   *objective_function_value::pobjfun=0.0;
02010   laplace_approximation_calculator::where_are_we_flag=2;
02011   pfmin->AD_uf_inner();
02012   if ( no_stuff==0 && quadratic_prior::get_num_quadratic_prior()>0)
02013   {
02014     quadratic_prior::get_M_calculations();
02015   }
02016   vf+=*objective_function_value::pobjfun;
02017   laplace_approximation_calculator::where_are_we_flag=0;
02018   initial_df1b2params::cobjfun=value(vf);
02019   gradcalc(usize,g);
02020   double maxg=max(fabs(g));
02021   fval=value(vf);
02022   if (!initial_params::mc_phase)
02023   {
02024     cout << setprecision(15) << " f = " << vf
02025          << " max g = " << maxg << endl;
02026   }
02027   return maxg;
02028 }
02029 
02034 double evaluate_function_quiet(const dvector& x,function_minimizer * pfmin)
02035 {
02036   int usize=initial_params::nvarcalc();
02037   //double f=0.0;
02038   dvector g(1,usize);
02039   independent_variables u(1,usize);
02040   u=x;
02041   dvariable vf=0.0;
02042   vf=initial_params::reset(dvar_vector(u));
02043   //vf=0.0;
02044   *objective_function_value::pobjfun=0.0;
02045   laplace_approximation_calculator::where_are_we_flag=2;
02046   pfmin->AD_uf_inner();
02047   if ( no_stuff==0 && quadratic_prior::get_num_quadratic_prior()>0)
02048   {
02049     quadratic_prior::get_M_calculations();
02050   }
02051   vf+=*objective_function_value::pobjfun;
02052   laplace_approximation_calculator::where_are_we_flag=0;
02053   initial_df1b2params::cobjfun=value(vf);
02054   gradcalc(usize,g);
02055   double maxg=max(fabs(g));
02056   return maxg;
02057 }
02058 
02063 void evaluate_function_gradient(double& f,const dvector& x,
02064   function_minimizer * pfmin,dvector& xadjoint,dvector& uadjoint)
02065 {
02066   int usize=initial_params::nvarcalc();
02067   dvector g(1,usize);
02068   independent_variables u(1,usize);
02069   u=x;
02070   dvariable vf=0.0;
02071   vf=initial_params::reset(dvar_vector(u));
02072   //vf=0.0;
02073   *objective_function_value::pobjfun=0.0;
02074   pfmin->AD_uf_inner();
02075   if ( no_stuff==0 && quadratic_prior::get_num_quadratic_prior()>0)
02076   {
02077     quadratic_prior::get_M_calculations();
02078   }
02079   vf+=*objective_function_value::pobjfun;
02080   initial_df1b2params::cobjfun=value(vf);
02081   f=value(vf);
02082   gradcalc(usize,g);
02083   int xsize=xadjoint.indexmax();
02084   int us=uadjoint.indexmax();
02085   xadjoint=g(1,xsize);
02086   uadjoint=g(xsize+1,xsize+us).shift(1);
02087 }
02088 
02093 double evaluate_function_no_derivatives(const dvector& x,
02094   function_minimizer* pfmin)
02095 {
02096   double fval;
02097   gradient_structure::set_NO_DERIVATIVES();
02098   int usize=initial_params::nvarcalc();
02099   //double f=0.0;
02100   dvector g(1,usize);
02101   independent_variables u(1,usize);
02102   u=x;
02103   dvariable vf=0.0;
02104   vf=initial_params::reset(dvar_vector(u));
02105   //vf=0.0;
02106   *objective_function_value::pobjfun=0.0;
02107   pfmin->AD_uf_inner();
02108   if ( no_stuff==0 && quadratic_prior::get_num_quadratic_prior()>0)
02109   {
02110     quadratic_prior::get_M_calculations();
02111   }
02112   vf+=*objective_function_value::pobjfun;
02113   initial_df1b2params::cobjfun=value(vf);
02114   fval=value(vf);
02115   gradient_structure::set_YES_DERIVATIVES();
02116   return fval;
02117 }
02118 
02123 void cleanup_laplace_stuff(laplace_approximation_calculator * l)
02124 {
02125   if (l)
02126   {
02127     delete l;
02128     l=0;
02129     df1b2variable::pool->deallocate();
02130 
02131     for (int i=0;i<df1b2variable::adpool_counter;i++)
02132     {
02133       delete df1b2variable::adpool_vector[i];
02134       df1b2variable::adpool_vector[i]=0;
02135       df1b2variable::nvar_vector[i]=0;
02136       df1b2variable::adpool_counter=0;
02137     }
02138   }
02139 }
02140 
02145 gauss_hermite_stuff::gauss_hermite_stuff
02146   (laplace_approximation_calculator * lapprox,
02147   int use_gauss_hermite,int num_separable_calls ,const ivector& itmp) :
02148   x(1,use_gauss_hermite), w(1,use_gauss_hermite), mi(0)
02149 {
02150   // for now we must have the same number of random effects in each
02151   // separable function call"
02152   int i;
02153   for (i=2;i<=num_separable_calls;i++)
02154   {
02155     if (itmp(i)!=itmp(i-1))
02156     {
02157       cerr << " At present for the adaptive gauss_hermite must have the same"
02158            << endl
02159            << " number of random effects in each separable function call"
02160            << endl;
02161       ad_exit(1);
02162     }
02163   }
02164   for (i=1;i<=num_separable_calls;i++)
02165   {
02166     if (itmp(i)>1)
02167     {
02168       lapprox->multi_random_effects=
02169         max(lapprox->multi_random_effects,itmp(i));
02170      /*
02171       cerr << "At present gauss-hermite is only implemented for"
02172         " one random effect per separable function call "
02173        << endl;
02174       ad_exit(1);
02175      */
02176     }
02177   }
02178   if (allocated(gauss_hermite_values))
02179     gauss_hermite_values.deallocate();
02180   if (lapprox->multi_random_effects==0)
02181   {
02182     gauss_hermite_values.allocate(1,num_separable_calls,1,use_gauss_hermite);
02183   }
02184   else
02185   {
02186    ivector indx=pow(use_gauss_hermite,itmp);
02187     gauss_hermite_values.allocate(1,num_separable_calls,1,indx);
02188   }
02189 
02190   normalized_gauss_hermite(x,w);
02191   is=0;
02192 }
02193 
02198 void laplace_approximation_calculator::set_default_hessian_type(void )
02199 {
02200   //if (function_minimizer::hesstype==0)
02201   {
02202     if (initial_df1b2params::separable_flag)
02203     {
02204      // have SEPARABLE_FUNCTION but no object of type quadratic_penalty or
02205      // normal prior -- can't tell if we should use hesstype 2 or 3
02206      // until we run
02207       hesstype=3;
02208     }
02209     else
02210     {
02211       hesstype=1;  // assume block diagonal
02212     }
02213   }
02214  /*
02215   else
02216   {
02217     hesstype=function_minimizer::hesstype;
02218   }
02219  */
02220 }
02221 
02226 double laplace_approximation_calculator::get_fx_fu(function_minimizer * pfmin)
02227 {
02228   initial_params::set_inactive_only_random_effects();
02229   initial_params::set_active_random_effects();
02230 
02231   if (grad_x==0)
02232   {
02233     grad_x=new dvector(1,xsize);
02234   }
02235   grad_x->initialize();
02236 
02237   if (grad_x_u==0)
02238   {
02239     grad_x_u=new dvector(1,xsize+usize);
02240   }
02241   pfmin->inner_opt_flag=0;
02242   independent_variables u(1,xsize+usize);
02243   //gradcalc(0,*grad_x_u);
02244   initial_params::xinit(u);    // get the initial values into the
02245 
02246   dvariable pen=0.0;
02247   dvariable vf=0.0;
02248   pen=initial_params::reset(dvar_vector(u));
02249   *objective_function_value::pobjfun=0.0;
02250   pfmin->AD_uf_inner();
02251   vf+=*objective_function_value::pobjfun;
02252   objective_function_value::fun_without_pen=value(vf);
02253   vf+=pen;
02254   gradcalc(xsize+usize,*grad_x_u);
02255   double f=value(vf);
02256   return f;
02257 }
02258 
02263   void laplace_approximation_calculator::begin_separable_call_stuff(void)
02264   {
02265     separable_call_level++;
02266     //build_up_nested_shape();
02267     //clean(nested_tree_position,separable_call_level);
02268     //nested_separable_calls_counter(separable_call_level)++;
02269     //nested_tree_position(separable_call_level)++;
02270   }
02271 
02276   void laplace_approximation_calculator::end_separable_call_stuff(void)
02277   {
02278     //build_up_nested_shape();
02279     //clean(nested_tree_position,separable_call_level);
02280     separable_call_level--;
02281   }
02282 
02287 void laplace_approximation_calculator::build_up_nested_shape(void)
02288 {
02289   int ll;
02290   ivector& a =nested_separable_calls_counter;
02291   int clean_level=0;
02292   switch(separable_call_level)
02293   {
02294   case 1:
02295     ll=1;
02296     if (nested_separable_calls_counter(ll)>0)
02297     {
02298       if (clean_level==0) clean_level=ll+1;
02299       if (nested_separable_calls_counter(ll+1)>0)
02300       {
02301         nested_shape(a(ll))=a(ll+1);
02302       }
02303       else
02304       {
02305         break;
02306       }
02307     }
02308     else
02309     {
02310       break;
02311     }
02312    case 2:
02313     ll=2;
02314     if (nested_separable_calls_counter(ll)>0)
02315     {
02316       if (clean_level==0) clean_level=ll+1;
02317       if (nested_separable_calls_counter(ll+1)>0)
02318       {
02319         nested_shape(a(1),a(2))=a(ll+1);
02320       }
02321       else
02322       {
02323         break;
02324       }
02325     }
02326     else
02327     {
02328       break;
02329     }
02330    case 3:
02331     ll=3;
02332     if (nested_separable_calls_counter(ll)>0)
02333     {
02334       if (clean_level==0) clean_level=ll+1;
02335       if (nested_separable_calls_counter(ll+1)>0)
02336       {
02337         nested_shape(a(1),a(2),a(3))=a(ll+1);
02338       }
02339       else
02340       {
02341         break;
02342       }
02343     }
02344     else
02345     {
02346       break;
02347     }
02348    case 4:
02349     ll=4;
02350     if (nested_separable_calls_counter(ll)>0)
02351     {
02352       if (clean_level==0) clean_level=ll+1;
02353       if (nested_separable_calls_counter(ll+1)>0)
02354       {
02355         nested_shape(a(1),a(2),a(3),a(4))=a(ll+1);
02356       }
02357       else
02358       {
02359         break;
02360       }
02361     }
02362     else
02363     {
02364       break;
02365     }
02366    default:
02367      cerr << "illegal value in " <<
02368        "laplace_approximation_calculator::build_up_nested_shape"
02369        << endl;
02370   }
02371   if (clean_level>0)
02372   {
02373     int mmax=a.indexmax();
02374     for (int i=clean_level;i<=mmax;i++)
02375     {
02376       a(i)=0;
02377     }
02378   }
02379 }
02380 
02385 ivector * nested_calls_shape::get_ptr1(void)
02386 {
02387   return ptr1;
02388 }
02389 
02394 imatrix * nested_calls_shape::get_ptr2(void)
02395 {
02396   return ptr2;
02397 }
02398 
02403 i3_array * nested_calls_shape::get_ptr3(void)
02404 {
02405   return ptr3;
02406 }
02407 
02412 i4_array * nested_calls_shape::get_ptr4(void)
02413 {
02414   return ptr4;
02415 }
02416 
02421 ostream &  operator << (const ostream& _s,const nested_calls_shape& _m)
02422 {
02423   ADUNCONST(nested_calls_shape,m)
02424   ADUNCONST(ofstream,s)
02425   if (m.get_ptr1())
02426     s<< *(m.get_ptr1()) << endl << endl;
02427 
02428   if (m.get_ptr2())
02429     s<< *(m.get_ptr2()) << endl << endl;
02430 
02431   if (m.get_ptr3())
02432     s<< *(m.get_ptr3()) << endl << endl;
02433 
02434   if (m.get_ptr4())
02435     s<< *(m.get_ptr4()) << endl << endl;
02436 
02437   return s;
02438 }
02439 
02444 void nested_calls_shape::trim(void)
02445 {
02446   int mmax1=0;
02447   if (ptr1)
02448   {
02449     int imin=ptr1->indexmin();
02450     int imax=ptr1->indexmax();
02451     for (int i=imin;i<=imax;i++)
02452     {
02453       if ( (*ptr1)(i)==0) break;
02454       mmax1++;
02455     }
02456   }
02457   if (mmax1==0)
02458   {
02459     delete ptr1;
02460     ptr1=0;
02461   }
02462   else
02463   {
02464     ivector * tmp=new ivector(1,mmax1);
02465     (*tmp)=(*ptr1)(1,mmax1);
02466     delete ptr1;
02467     ptr1=tmp;
02468   }
02469   if (ptr2)
02470   {
02471     if (!ptr1)
02472     {
02473       delete ptr2;
02474       ptr2=0;
02475     }
02476     else
02477     {
02478       imatrix * tmp=new imatrix(1,mmax1);
02479       int imin=tmp->indexmin();
02480       int imax=tmp->indexmax();
02481       for (int i=imin;i<=imax;i++)
02482       {
02483         int jmin=(*ptr2)(i).indexmin();
02484         int jmax=(*ptr2)(i).indexmax();
02485         int mmax2=0;
02486         for (int j=jmin;j<=jmax;j++)
02487         {
02488           if ((*ptr2)(i,j)==0) break;
02489           mmax2++;
02490         }
02491         if (mmax2>0)
02492         {
02493           (*tmp)(i).allocate(1,mmax2);
02494           (*tmp)(i)=(*ptr2)(i)(1,mmax2);
02495         }
02496       }
02497       delete ptr2;
02498       ptr2=tmp;
02499     }
02500   }
02501   if (ptr3)
02502   {
02503     cerr << "warning not dealitn with prt3" << endl;
02504     delete ptr3;
02505     ptr3=0;
02506   }
02507 }
02508 
02513 nested_calls_shape::~nested_calls_shape()
02514 {
02515   if (ptr1)
02516   {
02517     delete ptr1;
02518     ptr1=0;
02519   }
02520   if (ptr2)
02521   {
02522     delete ptr2;
02523     ptr2=0;
02524   }
02525   if (ptr3)
02526   {
02527     delete ptr3;
02528     ptr3=0;
02529   }
02530   if (ptr4)
02531   {
02532     delete ptr4;
02533     ptr4=0;
02534   }
02535 }
02536 
02541 void nested_calls_shape::allocate(int n,int m,int p )
02542 {
02543   if (ptr1)
02544   {
02545     delete ptr1;
02546     ptr1=0;
02547   }
02548   ptr1=new ivector(1,n);
02549 
02550   if (ptr2)
02551   {
02552     delete ptr2;
02553     ptr2=0;
02554   }
02555   ptr2=new imatrix(1,n,1,m);
02556   if (ptr3)
02557   {
02558     delete ptr3;
02559     ptr3=0;
02560   }
02561   ptr3=new i3_array(1,n,1,m,1,p);
02562   /*
02563   if (ptr4)
02564   {
02565     delete ptr4;
02566     ptr4=0;
02567   }
02568   */
02569 }
02570 
02575 void nested_calls_shape::allocate(int n,int m)
02576 {
02577   if (ptr1)
02578   {
02579     delete ptr1;
02580     ptr1=0;
02581   }
02582   ptr1=new ivector(1,n);
02583 
02584   if (ptr2)
02585   {
02586     delete ptr2;
02587     ptr2=0;
02588   }
02589   ptr2=new imatrix(1,n,1,m);
02590   /*
02591   if (ptr3)
02592   {
02593     delete ptr3;
02594     ptr3=0;
02595   }
02596   ptr=new i3_array(1,n,1,m,1,p);
02597   if (ptr4)
02598   {
02599     delete ptr4;
02600     ptr4=0;
02601   }
02602   */
02603 }
02604 
02609 void nested_calls_shape::initialize(void)
02610 {
02611   if (ptr1)
02612   {
02613     ptr1->initialize();
02614   }
02615 
02616   if (ptr2)
02617   {
02618     ptr2->initialize();
02619   }
02620 
02621   if (ptr3)
02622   {
02623     ptr3->initialize();
02624   }
02625 
02626   if (ptr4)
02627   {
02628     ptr4->initialize();
02629   }
02630 }
02631 
02636 void nested_calls_shape::allocate(int n)
02637 {
02638   if (ptr1)
02639   {
02640     delete ptr1;
02641     ptr1=0;
02642   }
02643   ptr1=new ivector(1,n);
02644 
02645   /*
02646   if (ptr2)
02647   {
02648     delete ptr2;
02649     ptr2=0;
02650   }
02651   ptr=new imatrix(1,n,1,m);
02652   if (ptr3)
02653   {
02654     delete ptr3;
02655     ptr3=0;
02656   }
02657   ptr=new i3_array(1,n,1,m,1,p);
02658   if (ptr4)
02659   {
02660     delete ptr4;
02661     ptr4=0;
02662   }
02663   */
02664 }
02665 
02670 void nested_calls_indices::allocate(const nested_calls_shape& _nsc)
02671 {
02672   int mmax=0;
02673   ADUNCONST(nested_calls_shape,nsc)
02674   mmax=nsc.get_ptr1()->indexmax();
02675   if (ptr1)
02676   {
02677     delete ptr1;
02678     ptr1=0;
02679   }
02680   if (nsc.get_ptr1())
02681   {
02682     ptr1=new imatrix(1,mmax,1,*nsc.get_ptr1());
02683   }
02684   if (ptr2)
02685   {
02686     delete ptr2;
02687     ptr2=0;
02688   }
02689   if (nsc.get_ptr2())
02690   {
02691     ptr2=new i3_array(1,mmax,1,*nsc.get_ptr1(),1,*nsc.get_ptr2());
02692   }
02693   if (ptr3)
02694   {
02695     delete ptr3;
02696     ptr3=0;
02697   }
02698   if (nsc.get_ptr3())
02699   {
02700     ptr3=new i4_array(1,mmax,1,*nsc.get_ptr1(),1,*nsc.get_ptr2(),
02701       1,*nsc.get_ptr3());
02702   }
02703 }
02704 
02709 dvector laplace_approximation_calculator::get_uhat_lm_newton2
02710   (const dvector& x,function_minimizer * pfmin)
02711 {
02712   //int on,nopt;
02713   pfmin->inner_opt_flag=1;
02714   double f=0.0;
02715   double fb=1.e+100;
02716   dvector g(1,usize);
02717   dvector ub(1,usize);
02718   independent_variables u(1,usize);
02719   gradcalc(0,g);
02720   fmc1.itn=0;
02721   fmc1.ifn=0;
02722   fmc1.ireturn=0;
02723   initial_params::xinit(u);    // get the initial values into the
02724   initial_params::xinit(ubest);    // get the initial values into the
02725   fmc1.ialph=0;
02726   fmc1.ihang=0;
02727   fmc1.ihflag=0;
02728 
02729   if (init_switch)
02730   {
02731     u.initialize();
02732   }
02733   else
02734   {
02735     u=ubest;
02736   }
02737 
02738   fmc1.dfn=1.e-2;
02739   dvariable pen=0.0;
02740   //cout << "starting  norm(u) = " << norm(u) << endl;
02741   while (fmc1.ireturn>=0)
02742   {
02743    /*
02744     double nu=norm(value(u));
02745     if (nu>400)
02746     {
02747       cout << "U norm(u) = " << nu  << endl;
02748     }
02749     cout << "V norm(u) = " << nu
02750          << " f = " << f  << endl;
02751     */
02752     fmc1.fmin(f,u,g);
02753     //cout << "W norm(u) = " << norm(value(u)) << endl;
02754     if (fmc1.ireturn>0)
02755     {
02756       dvariable vf=0.0;
02757       pen=initial_params::reset(dvar_vector(u));
02758       *objective_function_value::pobjfun=0.0;
02759       pfmin->AD_uf_inner();
02760       if (saddlepointflag)
02761       {
02762         *objective_function_value::pobjfun*=-1.0;
02763       }
02764       if ( no_stuff==0 && quadratic_prior::get_num_quadratic_prior()>0)
02765       {
02766         quadratic_prior::get_M_calculations();
02767       }
02768       vf+=*objective_function_value::pobjfun;
02769 
02770      /*  this is now done in the operator = function
02771       if (quadratic_prior::get_num_quadratic_prior()>0)
02772       {
02773         vf+= quadratic_prior::get_quadratic_priors();
02774       }
02775       */
02776 
02777       objective_function_value::fun_without_pen=value(vf);
02778 
02779       //cout << " pen = " << pen << endl;
02780       if (noboundepen_flag==0)
02781       {
02782         vf+=pen;
02783       }
02784       f=value(vf);
02785       if (f<fb)
02786       {
02787         fb=f;
02788         ub=u;
02789       }
02790       gradcalc(usize,g);
02791       //cout << " f = " << setprecision(17) << f << " " << norm(g)
02792        // << " " << norm(u) << endl;
02793     }
02794     u=ub;
02795   }
02796   cout <<  " inner maxg = " <<  fmc1.gmax;
02797   if (fabs(fmc1.gmax)>1.e+3)
02798     trapper();
02799 
02800   if (fabs(fmc1.gmax)>1.e-4)
02801   {
02802     fmc1.itn=0;
02803     //fmc1.crit=1.e-9;
02804     fmc1.ifn=0;
02805     fmc1.ireturn=0;
02806     fmc1.ihang=0;
02807     fmc1.ihflag=0;
02808     fmc1.ialph=0;
02809     initial_params::xinit(u);    // get the initial values into the
02810     //u.initialize();
02811     while (fmc1.ireturn>=0)
02812     {
02813       fmc1.fmin(f,u,g);
02814       if (fmc1.ireturn>0)
02815       {
02816         dvariable vf=0.0;
02817         pen=initial_params::reset(dvar_vector(u));
02818         *objective_function_value::pobjfun=0.0;
02819         pfmin->AD_uf_inner();
02820         if ( no_stuff==0 && quadratic_prior::get_num_quadratic_prior()>0)
02821         {
02822           quadratic_prior::get_M_calculations();
02823         }
02824         vf+=*objective_function_value::pobjfun;
02825         objective_function_value::fun_without_pen=value(vf);
02826 
02827         vf+=pen;
02828         f=value(vf);
02829         if (f<fb)
02830         {
02831           fb=f;
02832           ub=u;
02833         }
02834         gradcalc(usize,g);
02835         //cout << " f = " << setprecision(15) << f << " " << norm(g) << endl;
02836       }
02837     }
02838     u=ub;
02839     cout <<  "  Inner second time = " << fmc1.gmax;
02840   }
02841   cout << "  Inner f = " << fb << endl;
02842   fmc1.ireturn=0;
02843   fmc1.fbest=fb;
02844   gradient_structure::set_NO_DERIVATIVES();
02845   *objective_function_value::pobjfun=0.0;
02846   pfmin->AD_uf_inner();
02847   if ( no_stuff==0 && quadratic_prior::get_num_quadratic_prior()>0)
02848   {
02849     quadratic_prior::get_M_calculations();
02850   }
02851   gradient_structure::set_YES_DERIVATIVES();
02852   pfmin->inner_opt_flag=0;
02853   return u;
02854 }