ADMB Documentation  11.1.2192
 All Classes Files Functions Variables Typedefs Friends Defines
df1b2lp1.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: df1b2lp1.cpp 1935 2014-04-26 02:02:58Z 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       else
00342       {
00343         trust_update_flag=0;
00344       }
00345 
00346       double maxg_old=maxg;
00347       pmin->inner_opt_flag=1;
00348       maxg=fabs(evaluate_function(uhat,pfmin));
00349       if (maxg>maxg_old)
00350       {
00351         uhat=uhat_old;
00352         evaluate_function(uhat,pfmin);
00353         break;
00354       }
00355       if (maxg < 1.e-13)
00356       {
00357         break;
00358       }
00359     }
00360     for (i=1;i<=usize;i++)
00361     {
00362       y(i+xsize)=uhat(i);
00363     }
00364   }
00365 
00366   if (num_nr_iters<=0)
00367   {
00368     evaluate_function(uhat,pfmin);
00369   }
00370 
00371   for (i=1;i<=usize;i++)
00372   {
00373     y(i+xsize)=uhat(i);
00374   }
00375 
00376 
00377   if (ad_comm::time_flag)
00378   {
00379     if (ad_comm::ptm)
00380     {
00381       double time=ad_comm::ptm->get_elapsed_time_and_reset();
00382       if (ad_comm::global_logfile)
00383       {
00384         (*ad_comm::global_logfile) << " Time in reset and evaluate function"
00385           << time << endl;
00386       }
00387     }
00388   }
00389   pmin->inner_opt_flag=0;
00390 
00391 
00392   if (saddlepointflag==2)
00393   {
00394     dvector localx(1,xsize+usize);
00395     for (int i=1;i<=xsize+usize;i++)
00396     {
00397       localx(i)=value(y(i));
00398     }
00399     initial_params::set_inactive_only_random_effects();
00400     initial_params::set_active_random_effects();
00401     /*int nnn=*/initial_params::nvarcalc();
00402     evaluate_function_gradient(f,localx,pfmin,xadjoint,uadjoint);
00403     pmin->inner_opt_flag=1;
00404     get_second_ders(xsize,usize,y,Hess,Dux,f1b2gradlist,pfmin,this);
00405     pmin->inner_opt_flag=0;
00406     xadjoint -= solve(Hess,uadjoint)*Dux;
00407     f1b2gradlist->reset();
00408     f1b2gradlist->list.initialize();
00409     f1b2gradlist->list2.initialize();
00410     f1b2gradlist->list3.initialize();
00411     f1b2gradlist->nlist.initialize();
00412     f1b2gradlist->nlist2.initialize();
00413     f1b2gradlist->nlist3.initialize();
00414   }
00415   else // don't need this for minimax calculations
00416   {
00417     get_second_ders(xsize,usize,y,Hess,Dux,f1b2gradlist,pfmin,this);
00418     //int sgn=0;
00419 
00420     if (ad_comm::time_flag)
00421     {
00422       if (ad_comm::ptm)
00423       {
00424         double time=ad_comm::ptm->get_elapsed_time_and_reset();
00425         if (ad_comm::global_logfile)
00426         {
00427           (*ad_comm::global_logfile) << " Time in dget second ders "
00428             << time << endl;
00429         }
00430       }
00431     }
00432     if (!ierr)
00433     {
00434       if (num_importance_samples==0)
00435       {
00436         //cout << "Hess " << endl << Hess << endl;
00437         f=calculate_laplace_approximation(x,uhat,Hess,xadjoint,uadjoint,
00438           Hessadjoint,pfmin);
00439       }
00440       else
00441       {
00442         if (isfunnel_flag==0)
00443         {
00444           f=calculate_importance_sample(x,uhat,Hess,xadjoint,uadjoint,
00445             Hessadjoint,pfmin);
00446         }
00447         else
00448         {
00449           f=calculate_importance_sample_funnel(x,uhat,Hess,xadjoint,uadjoint,
00450             Hessadjoint,pfmin);
00451         }
00452       }
00453     }
00454     else
00455     {
00456       f=1.e+30;
00457     }
00458 
00459     if (ad_comm::time_flag)
00460     {
00461       if (ad_comm::ptm)
00462       {
00463         double time=ad_comm::ptm->get_elapsed_time_and_reset();
00464         if (ad_comm::global_logfile)
00465         {
00466           (*ad_comm::global_logfile) <<
00467             " Time in calculate laplace approximation " << time << endl;
00468         }
00469       }
00470     }
00471 
00472     for (int ip=num_der_blocks;ip>=1;ip--)
00473     {
00474       df1b2variable::minder=minder(ip);
00475       df1b2variable::maxder=maxder(ip);
00476       int mind=y(1).minder;
00477       int jmin=max(mind,xsize+1);
00478       int jmax=min(y(1).maxder,xsize+usize);
00479       for (i=1;i<=usize;i++)
00480       {
00481         for (j=jmin;j<=jmax;j++)
00482         {
00483           //Hess(i,j-xsize)=y(i+xsize).u_bar[j-mind];
00484           y(i+xsize).get_u_bar_tilde()[j-mind]=Hessadjoint(i,j-xsize);
00485         }
00486       }
00487 
00488       if (initial_df1b2params::separable_flag)
00489       {
00490         for (j=1;j<=xsize+usize;j++)
00491         {
00492           *y(j).get_u_tilde()=0;
00493         }
00494         Hess.initialize();
00495         initial_df1b2params::separable_calculation_type=3;
00496         pfmin->user_function();
00497       }
00498       else
00499       {
00500         if (ip<num_der_blocks)
00501         {
00502           f1b2gradlist->reset();
00503           set_u_dot(ip);
00504           df1b2_gradlist::set_yes_derivatives();
00505           (*re_objective_function_value::pobjfun)=0;
00506           df1b2variable pen=0.0;
00507           df1b2variable zz=0.0;
00508 
00509           initial_df1b2params::reset(y,pen);
00510           pfmin->user_function();
00511 
00512           re_objective_function_value::fun_without_pen=
00513             value(*re_objective_function_value::pobjfun);
00514 
00515           (*re_objective_function_value::pobjfun)+=pen;
00516           (*re_objective_function_value::pobjfun)+=zz;
00517 
00518           set_dependent_variable(*re_objective_function_value::pobjfun);
00519           df1b2_gradlist::set_no_derivatives();
00520           df1b2variable::passnumber=1;
00521           df1b2_gradcalc1();
00522         }
00523 
00524         for (i=1;i<=usize;i++)
00525         {
00526           for (j=jmin;j<=jmax;j++)
00527           {
00528             //Hess(i,j-xsize)=y(i+xsize).u_bar[j-mind];
00529             y(i+xsize).get_u_bar_tilde()[j-mind]=Hessadjoint(i,j-xsize);
00530           }
00531         }
00532 
00533         //int mind=y(1).minder;
00534         df1b2variable::passnumber=2;
00535         df1b2_gradcalc1();
00536 
00537         df1b2variable::passnumber=3;
00538         df1b2_gradcalc1();
00539 
00540         f1b2gradlist->reset();
00541         f1b2gradlist->list.initialize();
00542         f1b2gradlist->list2.initialize();
00543         f1b2gradlist->list3.initialize();
00544         f1b2gradlist->nlist.initialize();
00545         f1b2gradlist->nlist2.initialize();
00546         f1b2gradlist->nlist3.initialize();
00547       }
00548 
00549       if (ad_comm::time_flag)
00550       {
00551         if (ad_comm::ptm)
00552         {
00553           double time=ad_comm::ptm->get_elapsed_time_and_reset();
00554           if (ad_comm::global_logfile)
00555           {
00556             (*ad_comm::global_logfile) << " time for 3rd derivatives "
00557               << time << endl;
00558           }
00559         }
00560       }
00561 
00562       dvector dtmp(1,xsize);
00563       for (i=1;i<=xsize;i++)
00564       {
00565         dtmp(i)=*y(i).get_u_tilde();
00566       }
00567       if (initial_df1b2params::separable_flag)
00568       {
00569         dvector scale(1,nvar);   // need to get scale from somewhere
00570         /*int check=*/initial_params::stddev_scale(scale,x);
00571         dvector sscale=scale(1,Dux(1).indexmax());
00572         for (i=1;i<=usize;i++)
00573         {
00574           Dux(i)=elem_prod(Dux(i),sscale);
00575         }
00576         dtmp=elem_prod(dtmp,sscale);
00577       }
00578 
00579       for (i=1;i<=xsize;i++)
00580       {
00581         xadjoint(i)+=dtmp(i);
00582       }
00583       for (i=1;i<=usize;i++)
00584         uadjoint(i)+=*y(xsize+i).get_u_tilde();
00585     }
00586    // *****************************************************************
00587    // new stuff to deal with quadraticprior
00588    // *****************************************************************
00589 
00590       int xstuff=3;
00591       if (xstuff && df1b2quadratic_prior::get_num_quadratic_prior()>0)
00592       {
00593         initial_params::straight_through_flag=0;
00594         funnel_init_var::lapprox=0;
00595         block_diagonal_flag=0;
00596         dvector scale1(1,nvar);   // need to get scale from somewhere
00597         initial_params::set_inactive_only_random_effects();
00598         /*int check=*/initial_params::stddev_scale(scale1,x);
00599 
00600         laplace_approximation_calculator::where_are_we_flag=3;
00601         quadratic_prior::in_qp_calculations=1;
00602         funnel_init_var::lapprox=this;
00603         df1b2_gradlist::set_no_derivatives();
00604         dvector scale(1,nvar);   // need to get scale from somewhere
00605         /*check=*/initial_params::stddev_scale(scale,x);
00606         dvector sscale=scale(1,Dux(1).indexmax());
00607 
00608         for (i=1;i<=usize;i++)
00609         {
00610           Dux(i)=elem_div(Dux(i),sscale);
00611         }
00612 
00613         if (xstuff>1)
00614         {
00615           df1b2quadratic_prior::get_Lxu_contribution(Dux);
00616         }
00617         quadratic_prior::in_qp_calculations=0;
00618         funnel_init_var::lapprox=0;
00619         laplace_approximation_calculator::where_are_we_flag=0;
00620 
00621         for (i=1;i<=usize;i++)
00622         {
00623           Dux(i)=elem_prod(Dux(i),sscale);
00624         }
00625         //local_dtemp=elem_prod(local_dtemp,sscale);
00626 
00627         if (xstuff>2)
00628         {
00629           dvector tmp=evaluate_function_with_quadprior(x,usize,pfmin);
00630           for (i=1;i<=xsize;i++)
00631           {
00632             xadjoint(i)+=tmp(i);
00633           }
00634         }
00635 
00636         if (xstuff>2)
00637         {
00638           quadratic_prior::get_cHessian_contribution_from_vHessian(Hess,xsize);
00639         }
00640       }
00641 
00642    // *****************************************************************
00643    // new stuff to deal with quadraticprior
00644    // *****************************************************************
00645     if (ad_comm::ptm)
00646     {
00647       /*double time=*/ad_comm::ptm->get_elapsed_time_and_reset();
00648     }
00649 
00650   #if defined(USE_ATLAS)
00651         if (!ad_comm::no_atlas_flag)
00652         {
00653           //xadjoint -= uadjoint*atlas_solve_spd_trans(Hess,Dux);
00654           xadjoint -= atlas_solve_spd_trans(Hess,uadjoint)*Dux;
00655         }
00656         else
00657         {
00658           //xadjoint -= uadjoint*solve(Hess,Dux);
00659           xadjoint -= solve(Hess,uadjoint)*Dux;
00660         }
00661   #else
00662         //xadjoint -= uadjoint*solve(Hess,Dux);
00663         xadjoint -= solve(Hess,uadjoint)*Dux;
00664   #endif
00665 
00666     if (ad_comm::ptm)
00667     {
00668       double time=ad_comm::ptm->get_elapsed_time_and_reset();
00669       if (ad_comm::global_logfile)
00670       {
00671         (*ad_comm::global_logfile) << " Time in second solve "
00672           << time << endl;
00673       }
00674     }
00675     if (ad_comm::ptm1)
00676     {
00677       double time=ad_comm::ptm1->get_elapsed_time_and_reset();
00678       if (ad_comm::global_logfile)
00679       {
00680         (*ad_comm::global_logfile) << " Total time in function evaluation "
00681           << time << endl << endl;
00682       }
00683     }
00684   }
00685   return xadjoint;
00686 }
00687 
00692 void laplace_approximation_calculator::get_newton_raphson_info
00693   (function_minimizer * pfmin)
00694 {
00695   int i,j,ip;
00696 
00697   if (ad_comm::time_flag)
00698   {
00699     if (ad_comm::ptm)
00700     {
00701         (*ad_comm::global_logfile) << " Starting Newton-Raphson "
00702           <<  endl;
00703     }
00704   }
00705 
00706   for (ip=1;ip<=num_der_blocks;ip++)
00707   {
00708     df1b2variable::minder=minder(ip);
00709     df1b2variable::maxder=maxder(ip);
00710 
00711     set_u_dot(ip);
00712 
00713     // do we need to reallocate memory for df1b2variables?
00714     check_for_need_to_reallocate(ip);
00715     df1b2_gradlist::set_yes_derivatives();
00716     //cout << re_objective_function_value::pobjfun << endl;
00717     //cout << re_objective_function_value::pobjfun->ptr << endl;
00718     (*re_objective_function_value::pobjfun)=0;
00719     df1b2variable pen=0.0;
00720     df1b2variable zz=0.0;
00721     initial_df1b2params::reset(y,pen);
00722      //cout << setprecision(15) << y << endl;
00723     if (initial_df1b2params::separable_flag)
00724     {
00725       Hess.initialize();
00726       grad.initialize();
00727     }
00728 
00729     double time1 = 0;
00730     if (ad_comm::time_flag)
00731     {
00732       if (ad_comm::ptm)
00733       {
00734         time1 = ad_comm::ptm->get_elapsed_time();
00735       }
00736     }
00737     pfmin->user_function();
00738     if (ad_comm::time_flag)
00739     {
00740       if (ad_comm::ptm)
00741       {
00742         if (ad_comm::global_logfile)
00743         {
00744           double time=ad_comm::ptm->get_elapsed_time();
00745           (*ad_comm::global_logfile) <<
00746             "       Time in user_function() " <<  ip << "  " << time-time1
00747             << endl;
00748         }
00749       }
00750     }
00751 
00752     re_objective_function_value::fun_without_pen
00753       =value(*re_objective_function_value::pobjfun);
00754 
00755     (*re_objective_function_value::pobjfun)+=pen;
00756     (*re_objective_function_value::pobjfun)+=zz;
00757 
00758     //cout << setprecision(15) << *re_objective_function_value::pobjfun << endl;
00759     //cout << setprecision(15) << pen << endl;
00760     if (!initial_df1b2params::separable_flag)
00761     {
00762       set_dependent_variable(*re_objective_function_value::pobjfun);
00763       df1b2_gradlist::set_no_derivatives();
00764       df1b2variable::passnumber=1;
00765       df1b2_gradcalc1();
00766       int mind=y(1).minder;
00767       int jmin=max(mind,xsize+1);
00768       int jmax=min(y(1).maxder,xsize+usize);
00769       for (i=1;i<=usize;i++)
00770         for (j=jmin;j<=jmax;j++)
00771           Hess(i,j-xsize)=y(i+xsize).u_bar[j-mind];
00772      /*
00773       {
00774         ofstream ofs("hessreport");
00775         ofs << setw(6) << Hess << endl << endl;
00776         ofs << setw(10) << setfixed() << setprecision(3) << choleski_decomp(Hess) << endl << endl;
00777         ofs << setw(10) << setfixed() << setprecision(3) << inv(Hess) << endl << endl;
00778         ofs << setw(10) << setfixed() << setprecision(3) << choleski_decomp(inv(Hess)) << endl << endl;
00779       }
00780       */
00781     // ****************************************************************
00782     // ****************************************************************
00783     // ****************************************************************
00784     // temporary shit
00785      /*
00786       for (i=1;i<=usize;i++)
00787       {
00788         for (j=jmin;j<=jmax;j++)
00789         {
00790           //Hess(i,j-xsize)=y(i+xsize).u_bar[j-mind];
00791           y(i+xsize).get_u_bar_tilde()[j-mind]=Hessadjoint(i,j-xsize);
00792         }
00793       }
00794       cout << Hess << endl;
00795       df1b2variable::passnumber=2;
00796       df1b2_gradcalc1();
00797 
00798       df1b2variable::passnumber=3;
00799       df1b2_gradcalc1();
00800      */
00801     // ****************************************************************
00802     // ****************************************************************
00803     // ****************************************************************
00804     // ****************************************************************
00805       for (j=jmin;j<=jmax;j++)
00806         grad(j-xsize)= re_objective_function_value::pobjfun->u_dot[j-mind];
00807     }
00808     if (ip<num_der_blocks)
00809       f1b2gradlist->reset();
00810   }
00811 
00812 
00813   // just to match master pvm routine
00814   if (ad_comm::time_flag)
00815   {
00816     if (ad_comm::ptm)
00817     {
00818       /*double time=*/ad_comm::ptm->get_elapsed_time();
00819     }
00820   }
00821 }
00822 
00827 void laplace_approximation_calculator::set_u_dot(int ip)
00828 {
00829   int mmin=y.indexmin();
00830   int mmax=y.indexmax();
00831 
00832   for (int i=mmin;i<=mmax;i++)
00833   {
00834     y(i).set_u_dot();
00835   }
00836 }
00837 
00842 void laplace_approximation_calculator::check_pool_size(void)
00843 {
00844   int num_active_parameters=nvar;
00845   f1b2gradlist->reset();
00846 
00847   adpool * tmppool=df1b2variable::pool;
00848   if (tmppool)
00849   {
00850     //cout << tmppool << endl;
00851     // check if current pool is the right size
00852     if (tmppool->nvar != num_active_parameters)
00853     {
00854       // check sizes of other pools
00855       int found_pool_flag=0;
00856       for (int i=0;i<df1b2variable::adpool_counter;i++)
00857       {
00858         if (df1b2variable::adpool_vector[i]->nvar == num_active_parameters)
00859         {
00860           adpool * tmp = df1b2variable::pool;
00861           //cout << "Old Depth = " << df1b2variable::pool->depth_check()
00862            //    << "  "  << df1b2variable::pool  << "  ";
00863           df1b2variable::pool=df1b2variable::adpool_vector[i];
00864           //cout << "Depth = " << df1b2variable::pool->depth_check()
00865            //    << "  "  << df1b2variable::pool  << endl;
00866           df1b2variable::adpool_vector[i]=tmp;
00867           df1b2variable::nvar_vector[i]=df1b2variable::pool->nvar;
00868           found_pool_flag=1;
00869           break;
00870         }
00871       }
00872       if (!found_pool_flag)
00873       {
00874         if (df1b2variable::adpool_counter>=df1b2variable::adpool_vectorsize)
00875         {
00876            cerr << "Need to increase adpool_vectorsize" << endl;
00877            ad_exit(1);
00878         }
00879         df1b2variable::adpool_vector[df1b2variable::adpool_counter]=
00880           df1b2variable::pool;
00881         df1b2variable::nvar_vector[df1b2variable::adpool_counter]=
00882           df1b2variable::pool->nvar;
00883         //df1b2variable::adpool_counter++;
00884         df1b2variable::increment_adpool_counter();
00885         df1b2variable::pool=new adpool();
00886         if (!df1b2variable::pool)
00887         {
00888           cerr << "Memory allocation error" << endl;
00889           ad_exit(1);
00890         }
00891       }
00892     }
00893   }
00894   else
00895   {
00896     df1b2variable::pool=new adpool();
00897     if (!df1b2variable::pool)
00898     {
00899       cerr << "Memory allocation error" << endl;
00900       ad_exit(1);
00901     }
00902   }
00903   df1b2variable::nvar=num_active_parameters;
00904   df1b2variable::set_blocksize();
00905 }