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