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