ADMB Documentation  11.1.1927
 All Classes Files Functions Variables Typedefs Friends Defines
df1b2lp1.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: df1b2lp1.cpp 1919 2014-04-22 22:02:01Z johnoel $
00003  *
00004  * Author: David Fournier
00005  * Copyright (c) 2008-2012 Regents of the University of California
00006  */
00011 #if defined(USE_LAPLACE)
00012 #  include <fvar.hpp>
00013 #  include <admodel.h>
00014 #  include <df1b2fun.h>
00015 #  include <adrndeff.h>
00016 void evaluate_function_gradient(double& f,const dvector& x,
00017   function_minimizer * pfmin,dvector&,dvector&);
00018 double evaluate_function(const dvector& x,function_minimizer * pfmin);
00019 void get_second_ders(int xs,int us,const init_df1b2vector y,dmatrix& Hess,
00020   dmatrix& Dux, df1b2_gradlist * f1b2gradlist,function_minimizer * pfmin,
00021   laplace_approximation_calculator* lap);
00022 double calculate_laplace_approximation(const dvector& x,const dvector& u0,
00023   const dmatrix& Hess,const dvector& _xadjoint,const dvector& _uadjoint,
00024   const dmatrix& _Hessadjoint,function_minimizer * pmin);
00025 
00026 double calculate_importance_sample(const dvector& x,const dvector& u0,
00027   const dmatrix& Hess,const dvector& _xadjoint,const dvector& _uadjoint,
00028   const dmatrix& _Hessadjoint,function_minimizer * pmin);
00029 
00030 double calculate_importance_sample_funnel(const dvector& x,const dvector& u0,
00031   const dmatrix& Hess,const dvector& _xadjoint,const dvector& _uadjoint,
00032   const dmatrix& _Hessadjoint,function_minimizer * pmin);
00033 
00034 dmatrix choleski_decomp_positive(const dmatrix& M,double b);
00035 
00040 dvector laplace_approximation_calculator::default_calculations
00041   (const dvector& _x,const double& _f,function_minimizer * pfmin)
00042 {
00043   // for use when there is no separability
00044   ADUNCONST(dvector,x)
00045   ADUNCONST(double,f)
00046   int i,j;
00047 
00048   initial_params::set_inactive_only_random_effects();
00049   gradient_structure::set_NO_DERIVATIVES();
00050   initial_params::reset(x);    // get current x values into the model
00051   gradient_structure::set_YES_DERIVATIVES();
00052 
00053   initial_params::set_active_only_random_effects();
00054   //int lmn_flag=0;
00055   if (ad_comm::time_flag)
00056   {
00057     if (ad_comm::ptm1)
00058     {
00059       ad_comm::ptm1->get_elapsed_time_and_reset();
00060     }
00061     if (ad_comm::ptm)
00062     {
00063       ad_comm::ptm->get_elapsed_time_and_reset();
00064     }
00065   }
00066   if (ad_comm::time_flag)
00067   {
00068     if (ad_comm::ptm)
00069     {
00070       double time=ad_comm::ptm->get_elapsed_time();
00071       if (ad_comm::global_logfile)
00072       {
00073         (*ad_comm::global_logfile) << " Time pos 0 "
00074           << time << endl;
00075       }
00076     }
00077   }
00078 
00079   double maxg=1.e+200;
00080   //double maxg_save;
00081   dvector uhat_old(1,usize);
00082   if (inner_maxfn>0)
00083   {
00084     if (!inner_lmnflag)
00085     {
00086       if (!ADqd_flag)
00087       {
00088         uhat=get_uhat_quasi_newton(x,pfmin);
00089         maxg=fabs(fmc1.gmax);
00090         //double f_from_1=fmc1.fbest;
00091       }
00092       else
00093         uhat=get_uhat_quasi_newton_qd(x,pfmin);
00094     }
00095     else
00096     {
00097       uhat=get_uhat_lm_newton(x,pfmin);
00098     }
00099 
00100     if (ad_comm::time_flag)
00101     {
00102       if (ad_comm::ptm)
00103       {
00104         double time=ad_comm::ptm->get_elapsed_time_and_reset();
00105         if (ad_comm::global_logfile)
00106         {
00107           (*ad_comm::global_logfile) << " Time in inner optimization "
00108             << time << endl;
00109         }
00110       }
00111     }
00112   }
00113 
00114   for (i=1;i<=xsize;i++)
00115   {
00116     y(i)=x(i);
00117   }
00118   for (i=1;i<=usize;i++)
00119   {
00120     y(i+xsize)=uhat(i);
00121   }
00122 
00123   int ierr=0;
00124   int niters=0;
00125   if (function_minimizer::first_hessian_flag)
00126     niters=num_nr_iters+1;
00127   else
00128     niters=num_nr_iters;
00129 
00130   int nv=0;
00131   if (quadratic_prior::get_num_quadratic_prior()>0)
00132   {
00133     nv=initial_df1b2params::set_index();
00134     if (allocated(used_flags))
00135     {
00136       if (used_flags.indexmax() != nv)
00137       {
00138         used_flags.safe_deallocate();
00139       }
00140     }
00141     if (!allocated(used_flags))
00142     {
00143       used_flags.safe_allocate(1,nv);
00144     }
00145   }
00146 
00147   for(int ii=1;ii<=niters;ii++)
00148   {
00149     if (quadratic_prior::get_num_quadratic_prior()>0)
00150     {
00151       check_pool_size();
00152     }
00153     {
00154       // test newton raphson
00155       Hess.initialize();
00156      cout << "Newton raphson " << ii << endl;
00157       pmin->inner_opt_flag=1;
00158       get_newton_raphson_info(pfmin);
00159       pmin->inner_opt_flag=0;
00160 
00161       if (quadratic_prior::get_num_quadratic_prior()>0)
00162       {
00163         laplace_approximation_calculator::where_are_we_flag=2;
00164         evaluate_function_quiet(uhat,pfmin);
00165         laplace_approximation_calculator::where_are_we_flag=0;
00166         quadratic_prior::get_cHessian_contribution(Hess,xsize);
00167         quadratic_prior::get_cgradient_contribution(grad,xsize);
00168       }
00169 
00170       if (ii==1)
00171       {
00172         /*
00173         double diff = fabs(re_objective_function_value::fun_without_pen
00174                       - objective_function_value::fun_without_pen);
00175         if (diff>1.e-7)
00176         {
00177           cout << "there is a difference in the the user_functions "
00178             << setprecision(15) <<  re_objective_function_value::fun_without_pen
00179             <<  " from autodif  we have "
00180             << setprecision(15) << objective_function_value::fun_without_pen
00181             << " diff = "
00182             << diff  << endl;
00183           //ad_exit(1);
00184         }
00185        */
00186       }
00187 
00188       if (ad_comm::time_flag)
00189       {
00190         if (ad_comm::ptm)
00191         {
00192           double time=ad_comm::ptm->get_elapsed_time_and_reset();
00193           if (ad_comm::global_logfile)
00194           {
00195             (*ad_comm::global_logfile) << " Time in newton-raphson " <<  ii
00196               << "  " << time << endl;
00197           }
00198         }
00199       }
00200 
00201       dvector step;
00202       int print_hess_in_newton_raphson_flag=0;
00203       if (print_hess_in_newton_raphson_flag)
00204       {
00205         cout << norm2(Hess-trans(Hess)) << endl;
00206         if (ad_comm::global_logfile)
00207         {
00208           (*ad_comm::global_logfile) << setprecision(4) << setscientific()
00209             << setw(12) << sort(eigenvalues(Hess)) << endl;
00210           (*ad_comm::global_logfile) << setprecision(4) << setscientific()
00211             << setw(12) << Hess << endl;
00212         }
00213       }
00214 
00215 #if defined(USE_ATLAS)
00216       if (!ad_comm::no_atlas_flag)
00217       {
00218         step=-atlas_solve_spd(Hess,grad,ierr);
00219       }
00220       else
00221       {
00222         dmatrix A=choleski_decomp_positive(Hess,ierr);
00223         if (!ierr)
00224         {
00225           step=-solve(Hess,grad);
00226           //step=-solve(A*trans(A),grad);
00227         }
00228       }
00229       if (ierr)
00230       {
00231         f1b2gradlist->reset();
00232         f1b2gradlist->list.initialize();
00233         f1b2gradlist->list2.initialize();
00234         f1b2gradlist->list3.initialize();
00235         f1b2gradlist->nlist.initialize();
00236         f1b2gradlist->nlist2.initialize();
00237         f1b2gradlist->nlist3.initialize();
00238         break;
00239       }
00240 #else
00241       //step=-solve(Hess,grad);
00242       int ierror=0;
00243       int icount=0;
00244       int trust_update_flag=0;
00245       do
00246       {
00247         icount++;
00248         if (saddlepointflag==1 || saddlepointflag==2)
00249         {
00250           step=choleski_solve_neghess_error(Hess,grad,ierror);
00251         }
00252         else
00253         {
00254           step=-choleski_solve_error(Hess,grad,ierror);
00255         }
00256         if (ierror==1)
00257         {
00258           trust_update_flag=1;
00259           uhat_old=uhat;
00260           dvector vv=sort(eigenvalues(Hess));
00261           // matrix is not positive definite
00262           // do minimization
00263           dvector s(grad.indexmin(),grad.indexmax());
00264           double lambda=0.01;
00265           if (saddlepointflag==2) // max not min
00266           {
00267             const dmatrix  & cmHess=-Hess;
00268             const dvector & cmgrad = -grad;
00269             dmatrix  & mHess = (dmatrix  &) (cmHess);
00270             dvector & mgrad = (dvector &) (cmgrad);
00271             step=local_minimization(s,mHess,mgrad,lambda);
00272           }
00273           else
00274           {
00275             step=local_minimization(s,Hess,grad,lambda);
00276           }
00277           uhat+=step;
00278           for (i=1;i<=usize;i++)
00279           {
00280             y(i+xsize)=uhat(i);
00281           }
00282 
00283           f1b2gradlist->reset();
00284           f1b2gradlist->list.initialize();
00285           f1b2gradlist->list2.initialize();
00286           f1b2gradlist->list3.initialize();
00287           f1b2gradlist->nlist.initialize();
00288           f1b2gradlist->nlist2.initialize();
00289           f1b2gradlist->nlist3.initialize();
00290 
00291           get_newton_raphson_info(pfmin);
00292 
00293           if (quadratic_prior::get_num_quadratic_prior()>0)
00294           {
00295             laplace_approximation_calculator::where_are_we_flag=2;
00296             evaluate_function_quiet(uhat,pfmin);
00297             laplace_approximation_calculator::where_are_we_flag=0;
00298             quadratic_prior::get_cHessian_contribution(Hess,xsize);
00299             quadratic_prior::get_cgradient_contribution(grad,xsize);
00300           }
00301         }
00302         else if (trust_update_flag==1)
00303         {
00304           cout << "Found positive definite Hessian with trust region method"
00305                << endl;
00306         }
00307       }
00308       while (ierror==1 && icount<100);
00309       if (ierror==1)
00310       {
00311         cerr << "Can't make a positive definite Hessian" << endl;
00312         ad_exit(1);
00313       }
00314 #endif
00315 
00316       if (ad_comm::time_flag)
00317       {
00318         if (ad_comm::ptm)
00319         {
00320           double time=ad_comm::ptm->get_elapsed_time_and_reset();
00321           if (ad_comm::global_logfile)
00322           {
00323             (*ad_comm::global_logfile) << " time_in solve " <<  ii << "  "
00324               << time << endl;
00325           }
00326         }
00327       }
00328 
00329       f1b2gradlist->reset();
00330       f1b2gradlist->list.initialize();
00331       f1b2gradlist->list2.initialize();
00332       f1b2gradlist->list3.initialize();
00333       f1b2gradlist->nlist.initialize();
00334       f1b2gradlist->nlist2.initialize();
00335       f1b2gradlist->nlist3.initialize();
00336 
00337       if (trust_update_flag==0)
00338       {
00339         uhat_old=uhat;
00340         uhat+=step;
00341       }
00342       else
00343       {
00344         trust_update_flag=0;
00345       }
00346 
00347       double maxg_old=maxg;
00348       pmin->inner_opt_flag=1;
00349       maxg=fabs(evaluate_function(uhat,pfmin));
00350       if (maxg>maxg_old)
00351       {
00352         uhat=uhat_old;
00353         evaluate_function(uhat,pfmin);
00354         break;
00355       }
00356       if (maxg < 1.e-13)
00357       {
00358         break;
00359       }
00360     }
00361     for (i=1;i<=usize;i++)
00362     {
00363       y(i+xsize)=uhat(i);
00364     }
00365   }
00366 
00367   if (num_nr_iters<=0)
00368   {
00369     evaluate_function(uhat,pfmin);
00370   }
00371 
00372   for (i=1;i<=usize;i++)
00373   {
00374     y(i+xsize)=uhat(i);
00375   }
00376 
00377 
00378   if (ad_comm::time_flag)
00379   {
00380     if (ad_comm::ptm)
00381     {
00382       double time=ad_comm::ptm->get_elapsed_time_and_reset();
00383       if (ad_comm::global_logfile)
00384       {
00385         (*ad_comm::global_logfile) << " Time in reset and evaluate function"
00386           << time << endl;
00387       }
00388     }
00389   }
00390   pmin->inner_opt_flag=0;
00391 
00392 
00393   if (saddlepointflag==2)
00394   {
00395     dvector localx(1,xsize+usize);
00396     for (int i=1;i<=xsize+usize;i++)
00397     {
00398       localx(i)=value(y(i));
00399     }
00400     initial_params::set_inactive_only_random_effects();
00401     initial_params::set_active_random_effects();
00402     /*int nnn=*/initial_params::nvarcalc();
00403     evaluate_function_gradient(f,localx,pfmin,xadjoint,uadjoint);
00404     pmin->inner_opt_flag=1;
00405     get_second_ders(xsize,usize,y,Hess,Dux,f1b2gradlist,pfmin,this);
00406     pmin->inner_opt_flag=0;
00407     xadjoint -= solve(Hess,uadjoint)*Dux;
00408     f1b2gradlist->reset();
00409     f1b2gradlist->list.initialize();
00410     f1b2gradlist->list2.initialize();
00411     f1b2gradlist->list3.initialize();
00412     f1b2gradlist->nlist.initialize();
00413     f1b2gradlist->nlist2.initialize();
00414     f1b2gradlist->nlist3.initialize();
00415   }
00416   else // don't need this for minimax calculations
00417   {
00418     get_second_ders(xsize,usize,y,Hess,Dux,f1b2gradlist,pfmin,this);
00419     //int sgn=0;
00420 
00421     if (ad_comm::time_flag)
00422     {
00423       if (ad_comm::ptm)
00424       {
00425         double time=ad_comm::ptm->get_elapsed_time_and_reset();
00426         if (ad_comm::global_logfile)
00427         {
00428           (*ad_comm::global_logfile) << " Time in dget second ders "
00429             << time << endl;
00430         }
00431       }
00432     }
00433     if (!ierr)
00434     {
00435       if (num_importance_samples==0)
00436       {
00437         //cout << "Hess " << endl << Hess << endl;
00438         f=calculate_laplace_approximation(x,uhat,Hess,xadjoint,uadjoint,
00439           Hessadjoint,pfmin);
00440       }
00441       else
00442       {
00443         if (isfunnel_flag==0)
00444         {
00445           f=calculate_importance_sample(x,uhat,Hess,xadjoint,uadjoint,
00446             Hessadjoint,pfmin);
00447         }
00448         else
00449         {
00450           f=calculate_importance_sample_funnel(x,uhat,Hess,xadjoint,uadjoint,
00451             Hessadjoint,pfmin);
00452         }
00453       }
00454     }
00455     else
00456     {
00457       f=1.e+30;
00458     }
00459 
00460     if (ad_comm::time_flag)
00461     {
00462       if (ad_comm::ptm)
00463       {
00464         double time=ad_comm::ptm->get_elapsed_time_and_reset();
00465         if (ad_comm::global_logfile)
00466         {
00467           (*ad_comm::global_logfile) <<
00468             " Time in calculate laplace approximation " << time << endl;
00469         }
00470       }
00471     }
00472 
00473     for (int ip=num_der_blocks;ip>=1;ip--)
00474     {
00475       df1b2variable::minder=minder(ip);
00476       df1b2variable::maxder=maxder(ip);
00477       int mind=y(1).minder;
00478       int jmin=max(mind,xsize+1);
00479       int jmax=min(y(1).maxder,xsize+usize);
00480       for (i=1;i<=usize;i++)
00481       {
00482         for (j=jmin;j<=jmax;j++)
00483         {
00484           //Hess(i,j-xsize)=y(i+xsize).u_bar[j-mind];
00485           y(i+xsize).get_u_bar_tilde()[j-mind]=Hessadjoint(i,j-xsize);
00486         }
00487       }
00488 
00489       if (initial_df1b2params::separable_flag)
00490       {
00491         for (j=1;j<=xsize+usize;j++)
00492         {
00493           *y(j).get_u_tilde()=0;
00494         }
00495         Hess.initialize();
00496         initial_df1b2params::separable_calculation_type=3;
00497         pfmin->user_function();
00498       }
00499       else
00500       {
00501         if (ip<num_der_blocks)
00502         {
00503           f1b2gradlist->reset();
00504           set_u_dot(ip);
00505           df1b2_gradlist::set_yes_derivatives();
00506           (*re_objective_function_value::pobjfun)=0;
00507           df1b2variable pen=0.0;
00508           df1b2variable zz=0.0;
00509 
00510           initial_df1b2params::reset(y,pen);
00511           pfmin->user_function();
00512 
00513           re_objective_function_value::fun_without_pen=
00514             value(*re_objective_function_value::pobjfun);
00515 
00516           (*re_objective_function_value::pobjfun)+=pen;
00517           (*re_objective_function_value::pobjfun)+=zz;
00518 
00519           set_dependent_variable(*re_objective_function_value::pobjfun);
00520           df1b2_gradlist::set_no_derivatives();
00521           df1b2variable::passnumber=1;
00522           df1b2_gradcalc1();
00523         }
00524 
00525         for (i=1;i<=usize;i++)
00526         {
00527           for (j=jmin;j<=jmax;j++)
00528           {
00529             //Hess(i,j-xsize)=y(i+xsize).u_bar[j-mind];
00530             y(i+xsize).get_u_bar_tilde()[j-mind]=Hessadjoint(i,j-xsize);
00531           }
00532         }
00533 
00534         //int mind=y(1).minder;
00535         df1b2variable::passnumber=2;
00536         df1b2_gradcalc1();
00537 
00538         df1b2variable::passnumber=3;
00539         df1b2_gradcalc1();
00540 
00541         f1b2gradlist->reset();
00542         f1b2gradlist->list.initialize();
00543         f1b2gradlist->list2.initialize();
00544         f1b2gradlist->list3.initialize();
00545         f1b2gradlist->nlist.initialize();
00546         f1b2gradlist->nlist2.initialize();
00547         f1b2gradlist->nlist3.initialize();
00548       }
00549 
00550       if (ad_comm::time_flag)
00551       {
00552         if (ad_comm::ptm)
00553         {
00554           double time=ad_comm::ptm->get_elapsed_time_and_reset();
00555           if (ad_comm::global_logfile)
00556           {
00557             (*ad_comm::global_logfile) << " time for 3rd derivatives "
00558               << time << endl;
00559           }
00560         }
00561       }
00562 
00563       dvector dtmp(1,xsize);
00564       for (i=1;i<=xsize;i++)
00565       {
00566         dtmp(i)=*y(i).get_u_tilde();
00567       }
00568       if (initial_df1b2params::separable_flag)
00569       {
00570         dvector scale(1,nvar);   // need to get scale from somewhere
00571         /*int check=*/initial_params::stddev_scale(scale,x);
00572         dvector sscale=scale(1,Dux(1).indexmax());
00573         for (i=1;i<=usize;i++)
00574         {
00575           Dux(i)=elem_prod(Dux(i),sscale);
00576         }
00577         dtmp=elem_prod(dtmp,sscale);
00578       }
00579 
00580       for (i=1;i<=xsize;i++)
00581       {
00582         xadjoint(i)+=dtmp(i);
00583       }
00584       for (i=1;i<=usize;i++)
00585         uadjoint(i)+=*y(xsize+i).get_u_tilde();
00586     }
00587    // *****************************************************************
00588    // new stuff to deal with quadraticprior
00589    // *****************************************************************
00590 
00591       int xstuff=3;
00592       if (xstuff && df1b2quadratic_prior::get_num_quadratic_prior()>0)
00593       {
00594         initial_params::straight_through_flag=0;
00595         funnel_init_var::lapprox=0;
00596         block_diagonal_flag=0;
00597         dvector scale1(1,nvar);   // need to get scale from somewhere
00598         initial_params::set_inactive_only_random_effects();
00599         /*int check=*/initial_params::stddev_scale(scale1,x);
00600 
00601         laplace_approximation_calculator::where_are_we_flag=3;
00602         quadratic_prior::in_qp_calculations=1;
00603         funnel_init_var::lapprox=this;
00604         df1b2_gradlist::set_no_derivatives();
00605         dvector scale(1,nvar);   // need to get scale from somewhere
00606         /*check=*/initial_params::stddev_scale(scale,x);
00607         dvector sscale=scale(1,Dux(1).indexmax());
00608 
00609         for (i=1;i<=usize;i++)
00610         {
00611           Dux(i)=elem_div(Dux(i),sscale);
00612         }
00613 
00614         if (xstuff>1)
00615         {
00616           df1b2quadratic_prior::get_Lxu_contribution(Dux);
00617         }
00618         quadratic_prior::in_qp_calculations=0;
00619         funnel_init_var::lapprox=0;
00620         laplace_approximation_calculator::where_are_we_flag=0;
00621 
00622         for (i=1;i<=usize;i++)
00623         {
00624           Dux(i)=elem_prod(Dux(i),sscale);
00625         }
00626         //local_dtemp=elem_prod(local_dtemp,sscale);
00627 
00628         if (xstuff>2)
00629         {
00630           dvector tmp=evaluate_function_with_quadprior(x,usize,pfmin);
00631           for (i=1;i<=xsize;i++)
00632           {
00633             xadjoint(i)+=tmp(i);
00634           }
00635         }
00636 
00637         if (xstuff>2)
00638         {
00639           quadratic_prior::get_cHessian_contribution_from_vHessian(Hess,xsize);
00640         }
00641       }
00642 
00643    // *****************************************************************
00644    // new stuff to deal with quadraticprior
00645    // *****************************************************************
00646     if (ad_comm::ptm)
00647     {
00648       /*double time=*/ad_comm::ptm->get_elapsed_time_and_reset();
00649     }
00650 
00651   #if defined(USE_ATLAS)
00652         if (!ad_comm::no_atlas_flag)
00653         {
00654           //xadjoint -= uadjoint*atlas_solve_spd_trans(Hess,Dux);
00655           xadjoint -= atlas_solve_spd_trans(Hess,uadjoint)*Dux;
00656         }
00657         else
00658         {
00659           //xadjoint -= uadjoint*solve(Hess,Dux);
00660           xadjoint -= solve(Hess,uadjoint)*Dux;
00661         }
00662   #else
00663         //xadjoint -= uadjoint*solve(Hess,Dux);
00664         xadjoint -= solve(Hess,uadjoint)*Dux;
00665   #endif
00666 
00667     if (ad_comm::ptm)
00668     {
00669       double time=ad_comm::ptm->get_elapsed_time_and_reset();
00670       if (ad_comm::global_logfile)
00671       {
00672         (*ad_comm::global_logfile) << " Time in second solve "
00673           << time << endl;
00674       }
00675     }
00676     if (ad_comm::ptm1)
00677     {
00678       double time=ad_comm::ptm1->get_elapsed_time_and_reset();
00679       if (ad_comm::global_logfile)
00680       {
00681         (*ad_comm::global_logfile) << " Total time in function evaluation "
00682           << time << endl << endl;
00683       }
00684     }
00685   }
00686   return xadjoint;
00687 }
00688 
00693 void laplace_approximation_calculator::get_newton_raphson_info
00694   (function_minimizer * pfmin)
00695 {
00696   int i,j,ip;
00697 
00698   if (ad_comm::time_flag)
00699   {
00700     if (ad_comm::ptm)
00701     {
00702         (*ad_comm::global_logfile) << " Starting Newton-Raphson "
00703           <<  endl;
00704     }
00705   }
00706 
00707   for (ip=1;ip<=num_der_blocks;ip++)
00708   {
00709     df1b2variable::minder=minder(ip);
00710     df1b2variable::maxder=maxder(ip);
00711 
00712     set_u_dot(ip);
00713 
00714     // do we need to reallocate memory for df1b2variables?
00715     check_for_need_to_reallocate(ip);
00716     df1b2_gradlist::set_yes_derivatives();
00717     //cout << re_objective_function_value::pobjfun << endl;
00718     //cout << re_objective_function_value::pobjfun->ptr << endl;
00719     (*re_objective_function_value::pobjfun)=0;
00720     df1b2variable pen=0.0;
00721     df1b2variable zz=0.0;
00722     initial_df1b2params::reset(y,pen);
00723      //cout << setprecision(15) << y << endl;
00724     if (initial_df1b2params::separable_flag)
00725     {
00726       Hess.initialize();
00727       grad.initialize();
00728     }
00729 
00730     double time1 = 0;
00731     if (ad_comm::time_flag)
00732     {
00733       if (ad_comm::ptm)
00734       {
00735         time1 = ad_comm::ptm->get_elapsed_time();
00736       }
00737     }
00738     pfmin->user_function();
00739     if (ad_comm::time_flag)
00740     {
00741       if (ad_comm::ptm)
00742       {
00743         if (ad_comm::global_logfile)
00744         {
00745           double time=ad_comm::ptm->get_elapsed_time();
00746           (*ad_comm::global_logfile) <<
00747             "       Time in user_function() " <<  ip << "  " << time-time1
00748             << endl;
00749         }
00750       }
00751     }
00752 
00753     re_objective_function_value::fun_without_pen
00754       =value(*re_objective_function_value::pobjfun);
00755 
00756     (*re_objective_function_value::pobjfun)+=pen;
00757     (*re_objective_function_value::pobjfun)+=zz;
00758 
00759     //cout << setprecision(15) << *re_objective_function_value::pobjfun << endl;
00760     //cout << setprecision(15) << pen << endl;
00761     if (!initial_df1b2params::separable_flag)
00762     {
00763       set_dependent_variable(*re_objective_function_value::pobjfun);
00764       df1b2_gradlist::set_no_derivatives();
00765       df1b2variable::passnumber=1;
00766       df1b2_gradcalc1();
00767       int mind=y(1).minder;
00768       int jmin=max(mind,xsize+1);
00769       int jmax=min(y(1).maxder,xsize+usize);
00770       for (i=1;i<=usize;i++)
00771         for (j=jmin;j<=jmax;j++)
00772           Hess(i,j-xsize)=y(i+xsize).u_bar[j-mind];
00773      /*
00774       {
00775         ofstream ofs("hessreport");
00776         ofs << setw(6) << Hess << endl << endl;
00777         ofs << setw(10) << setfixed() << setprecision(3) << choleski_decomp(Hess) << endl << endl;
00778         ofs << setw(10) << setfixed() << setprecision(3) << inv(Hess) << endl << endl;
00779         ofs << setw(10) << setfixed() << setprecision(3) << choleski_decomp(inv(Hess)) << endl << endl;
00780       }
00781       */
00782     // ****************************************************************
00783     // ****************************************************************
00784     // ****************************************************************
00785     // temporary shit
00786      /*
00787       for (i=1;i<=usize;i++)
00788       {
00789         for (j=jmin;j<=jmax;j++)
00790         {
00791           //Hess(i,j-xsize)=y(i+xsize).u_bar[j-mind];
00792           y(i+xsize).get_u_bar_tilde()[j-mind]=Hessadjoint(i,j-xsize);
00793         }
00794       }
00795       cout << Hess << endl;
00796       df1b2variable::passnumber=2;
00797       df1b2_gradcalc1();
00798 
00799       df1b2variable::passnumber=3;
00800       df1b2_gradcalc1();
00801      */
00802     // ****************************************************************
00803     // ****************************************************************
00804     // ****************************************************************
00805     // ****************************************************************
00806       for (j=jmin;j<=jmax;j++)
00807         grad(j-xsize)= re_objective_function_value::pobjfun->u_dot[j-mind];
00808     }
00809     if (ip<num_der_blocks)
00810       f1b2gradlist->reset();
00811   }
00812 
00813 
00814   // just to match master pvm routine
00815   if (ad_comm::time_flag)
00816   {
00817     if (ad_comm::ptm)
00818     {
00819       /*double time=*/ad_comm::ptm->get_elapsed_time();
00820     }
00821   }
00822 }
00823 
00828 void laplace_approximation_calculator::set_u_dot(int ip)
00829 {
00830   int mmin=y.indexmin();
00831   int mmax=y.indexmax();
00832 
00833   for (int i=mmin;i<=mmax;i++)
00834   {
00835     y(i).set_u_dot();
00836   }
00837 }
00838 
00843 void laplace_approximation_calculator::check_pool_size(void)
00844 {
00845   int num_active_parameters=nvar;
00846   f1b2gradlist->reset();
00847 
00848   adpool * tmppool=df1b2variable::pool;
00849   if (tmppool)
00850   {
00851     //cout << tmppool << endl;
00852     // check if current pool is the right size
00853     if (tmppool->nvar != num_active_parameters)
00854     {
00855       // check sizes of other pools
00856       int found_pool_flag=0;
00857       for (int i=0;i<df1b2variable::adpool_counter;i++)
00858       {
00859         if (df1b2variable::adpool_vector[i]->nvar == num_active_parameters)
00860         {
00861           adpool * tmp = df1b2variable::pool;
00862           //cout << "Old Depth = " << df1b2variable::pool->depth_check()
00863            //    << "  "  << df1b2variable::pool  << "  ";
00864           df1b2variable::pool=df1b2variable::adpool_vector[i];
00865           //cout << "Depth = " << df1b2variable::pool->depth_check()
00866            //    << "  "  << df1b2variable::pool  << endl;
00867           df1b2variable::adpool_vector[i]=tmp;
00868           df1b2variable::nvar_vector[i]=df1b2variable::pool->nvar;
00869           found_pool_flag=1;
00870           break;
00871         }
00872       }
00873       if (!found_pool_flag)
00874       {
00875         if (df1b2variable::adpool_counter>=df1b2variable::adpool_vectorsize)
00876         {
00877            cerr << "Need to increase adpool_vectorsize" << endl;
00878            ad_exit(1);
00879         }
00880         df1b2variable::adpool_vector[df1b2variable::adpool_counter]=
00881           df1b2variable::pool;
00882         df1b2variable::nvar_vector[df1b2variable::adpool_counter]=
00883           df1b2variable::pool->nvar;
00884         //df1b2variable::adpool_counter++;
00885         df1b2variable::increment_adpool_counter();
00886         df1b2variable::pool=new adpool();
00887         if (!df1b2variable::pool)
00888         {
00889           cerr << "Memory allocation error" << endl;
00890           ad_exit(1);
00891         }
00892       }
00893     }
00894   }
00895   else
00896   {
00897     df1b2variable::pool=new adpool();
00898     if (!df1b2variable::pool)
00899     {
00900       cerr << "Memory allocation error" << endl;
00901       ad_exit(1);
00902     }
00903   }
00904   df1b2variable::nvar=num_active_parameters;
00905   df1b2variable::set_blocksize();
00906 }
00907 #endif //#if defined(USE_LAPLACE)