ADMB Documentation  11.6rc.3405
 All Classes Files Functions Variables Typedefs Friends Defines
df1b2lp1.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id$
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 #ifdef DIAG
00205       int print_hess_in_newton_raphson_flag=0;
00206       if (print_hess_in_newton_raphson_flag)
00207       {
00208         cout << norm2(Hess-trans(Hess)) << endl;
00209         if (ad_comm::global_logfile)
00210         {
00211           (*ad_comm::global_logfile) << setprecision(4) << setscientific()
00212             << setw(12) << sort(eigenvalues(Hess)) << endl;
00213           (*ad_comm::global_logfile) << setprecision(4) << setscientific()
00214             << setw(12) << Hess << endl;
00215         }
00216       }
00217 #endif
00218 
00219 #if defined(USE_ATLAS)
00220       if (!ad_comm::no_atlas_flag)
00221       {
00222         step=-atlas_solve_spd(Hess,grad,ierr);
00223       }
00224       else
00225       {
00226         dmatrix A=choleski_decomp_positive(Hess,ierr);
00227         if (!ierr)
00228         {
00229           step=-solve(Hess,grad);
00230           //step=-solve(A*trans(A),grad);
00231         }
00232       }
00233       if (ierr)
00234       {
00235         f1b2gradlist->reset();
00236         f1b2gradlist->list.initialize();
00237         f1b2gradlist->list2.initialize();
00238         f1b2gradlist->list3.initialize();
00239         f1b2gradlist->nlist.initialize();
00240         f1b2gradlist->nlist2.initialize();
00241         f1b2gradlist->nlist3.initialize();
00242         break;
00243       }
00244 #else
00245       //step=-solve(Hess,grad);
00246       int ierror=0;
00247       int icount=0;
00248       int trust_update_flag=0;
00249       do
00250       {
00251         icount++;
00252         if (saddlepointflag==1 || saddlepointflag==2)
00253         {
00254           step=choleski_solve_neghess_error(Hess,grad,ierror);
00255         }
00256         else
00257         {
00258           step=-choleski_solve_error(Hess,grad,ierror);
00259         }
00260         if (ierror==1)
00261         {
00262           trust_update_flag=1;
00263           uhat_old=uhat;
00264           dvector vv=sort(eigenvalues(Hess));
00265           // matrix is not positive definite
00266           // do minimization
00267           dvector s(grad.indexmin(),grad.indexmax());
00268           double lambda=0.01;
00269           if (saddlepointflag==2) // max not min
00270           {
00271             const dmatrix  & cmHess=-Hess;
00272             const dvector & cmgrad = -grad;
00273             dmatrix  & mHess = (dmatrix  &) (cmHess);
00274             dvector & mgrad = (dvector &) (cmgrad);
00275             step=local_minimization(s,mHess,mgrad,lambda);
00276           }
00277           else
00278           {
00279             step=local_minimization(s,Hess,grad,lambda);
00280           }
00281           if(!ad_comm::print_hess_and_exit_flag){
00282             uhat+=step;
00283           }
00284           for (int i=1;i<=usize;i++)
00285           {
00286             y(i+xsize)=uhat(i);
00287           }
00288 
00289           f1b2gradlist->reset();
00290           f1b2gradlist->list.initialize();
00291           f1b2gradlist->list2.initialize();
00292           f1b2gradlist->list3.initialize();
00293           f1b2gradlist->nlist.initialize();
00294           f1b2gradlist->nlist2.initialize();
00295           f1b2gradlist->nlist3.initialize();
00296 
00297           get_newton_raphson_info(pfmin);
00298 
00299           if (quadratic_prior::get_num_quadratic_prior()>0)
00300           {
00301             laplace_approximation_calculator::where_are_we_flag=2;
00302             evaluate_function_quiet(uhat,pfmin);
00303             laplace_approximation_calculator::where_are_we_flag=0;
00304             quadratic_prior::get_cHessian_contribution(Hess,xsize);
00305             quadratic_prior::get_cgradient_contribution(grad,xsize);
00306           }
00307         }
00308         else if (trust_update_flag==1)
00309         {
00310           cout << "Found positive definite Hessian with trust region method"
00311                << endl;
00312         }
00313       }
00314       while (ierror==1 && icount<100);
00315       if (ierror==1)
00316       {
00317         cerr << "Can't make a positive definite Hessian" << endl;
00318         ad_exit(1);
00319       }
00320 #endif
00321 
00322       if (ad_comm::time_flag)
00323       {
00324         if (ad_comm::ptm)
00325         {
00326           double time=ad_comm::ptm->get_elapsed_time_and_reset();
00327           if (ad_comm::global_logfile)
00328           {
00329             (*ad_comm::global_logfile) << " time_in solve " <<  ii << "  "
00330               << time << endl;
00331           }
00332         }
00333       }
00334 
00335       f1b2gradlist->reset();
00336       f1b2gradlist->list.initialize();
00337       f1b2gradlist->list2.initialize();
00338       f1b2gradlist->list3.initialize();
00339       f1b2gradlist->nlist.initialize();
00340       f1b2gradlist->nlist2.initialize();
00341       f1b2gradlist->nlist3.initialize();
00342 
00343       if (trust_update_flag==0)
00344       {
00345         uhat_old=uhat;
00346         if(!ad_comm::print_hess_and_exit_flag){
00347           uhat+=step;
00348         }
00349       }
00350 
00351       double maxg_old=maxg;
00352       pmin->inner_opt_flag=1;
00353       maxg=fabs(evaluate_function(uhat,pfmin));
00354       if (maxg>maxg_old)
00355       {
00356         if(!ad_comm::print_hess_and_exit_flag){
00357           uhat=uhat_old;
00358         }
00359         evaluate_function(uhat,pfmin);
00360         break;
00361       }
00362       if (maxg < 1.e-13)
00363       {
00364         break;
00365       }
00366     }
00367     for (int i=1;i<=usize;i++)
00368     {
00369       y(i+xsize)=uhat(i);
00370     }
00371   }
00372 
00373   if (num_nr_iters<=0)
00374   {
00375     evaluate_function(uhat,pfmin);
00376   }
00377 
00378   for (int i=1;i<=usize;i++)
00379   {
00380     y(i+xsize)=uhat(i);
00381   }
00382 
00383 
00384   if (ad_comm::time_flag)
00385   {
00386     if (ad_comm::ptm)
00387     {
00388       double time=ad_comm::ptm->get_elapsed_time_and_reset();
00389       if (ad_comm::global_logfile)
00390       {
00391         (*ad_comm::global_logfile) << " Time in reset and evaluate function"
00392           << time << endl;
00393       }
00394     }
00395   }
00396   pmin->inner_opt_flag=0;
00397 
00398 
00399   if (saddlepointflag==2)
00400   {
00401     dvector localx(1,xsize+usize);
00402     for (int i=1;i<=xsize+usize;i++)
00403     {
00404       localx(i)=value(y(i));
00405     }
00406     initial_params::set_inactive_only_random_effects();
00407     initial_params::set_active_random_effects();
00408     /*int nnn=*/initial_params::nvarcalc();
00409     evaluate_function_gradient(f,localx,pfmin,xadjoint,uadjoint);
00410     pmin->inner_opt_flag=1;
00411     get_second_ders(xsize,usize,y,Hess,Dux,f1b2gradlist,pfmin,this);
00412     pmin->inner_opt_flag=0;
00413     xadjoint -= solve(Hess,uadjoint)*Dux;
00414     f1b2gradlist->reset();
00415     f1b2gradlist->list.initialize();
00416     f1b2gradlist->list2.initialize();
00417     f1b2gradlist->list3.initialize();
00418     f1b2gradlist->nlist.initialize();
00419     f1b2gradlist->nlist2.initialize();
00420     f1b2gradlist->nlist3.initialize();
00421   }
00422   else // don't need this for minimax calculations
00423   {
00424     get_second_ders(xsize,usize,y,Hess,Dux,f1b2gradlist,pfmin,this);
00425     //int sgn=0;
00426 
00427     if (ad_comm::time_flag)
00428     {
00429       if (ad_comm::ptm)
00430       {
00431         double time=ad_comm::ptm->get_elapsed_time_and_reset();
00432         if (ad_comm::global_logfile)
00433         {
00434           (*ad_comm::global_logfile) << " Time in dget second ders "
00435             << time << endl;
00436         }
00437       }
00438     }
00439     if (!ierr)
00440     {
00441       if (num_importance_samples==0)
00442       {
00443         //cout << "Hess " << endl << Hess << endl;
00444         f=calculate_laplace_approximation(x,uhat,Hess,xadjoint,uadjoint,
00445           Hessadjoint,pfmin);
00446       }
00447       else
00448       {
00449         if (isfunnel_flag==0)
00450         {
00451           f=calculate_importance_sample(x,uhat,Hess,xadjoint,uadjoint,
00452             Hessadjoint,pfmin);
00453         }
00454         else
00455         {
00456           f=calculate_importance_sample_funnel(x,uhat,Hess,xadjoint,uadjoint,
00457             Hessadjoint,pfmin);
00458         }
00459       }
00460     }
00461     else
00462     {
00463       f=1.e+30;
00464     }
00465 
00466     if (ad_comm::time_flag)
00467     {
00468       if (ad_comm::ptm)
00469       {
00470         double time=ad_comm::ptm->get_elapsed_time_and_reset();
00471         if (ad_comm::global_logfile)
00472         {
00473           (*ad_comm::global_logfile) <<
00474             " Time in calculate laplace approximation " << time << endl;
00475         }
00476       }
00477     }
00478 
00479     for (int ip=num_der_blocks;ip>=1;ip--)
00480     {
00481       df1b2variable::minder=minder(ip);
00482       df1b2variable::maxder=maxder(ip);
00483       int mind=y(1).minder;
00484       int jmin=max(mind,xsize+1);
00485       int jmax=min(y(1).maxder,xsize+usize);
00486       for (int i=1;i<=usize;i++)
00487       {
00488         for (int j=jmin;j<=jmax;j++)
00489         {
00490           //Hess(i,j-xsize)=y(i+xsize).u_bar[j-mind];
00491           y(i+xsize).get_u_bar_tilde()[j-mind]=Hessadjoint(i,j-xsize);
00492         }
00493       }
00494 
00495       if (initial_df1b2params::separable_flag)
00496       {
00497         for (int j=1;j<=xsize+usize;j++)
00498         {
00499           *y(j).get_u_tilde()=0;
00500         }
00501         Hess.initialize();
00502         initial_df1b2params::separable_calculation_type=3;
00503         pfmin->user_function();
00504       }
00505       else
00506       {
00507         if (ip<num_der_blocks)
00508         {
00509           f1b2gradlist->reset();
00510           set_u_dot(ip);
00511           df1b2_gradlist::set_yes_derivatives();
00512           (*re_objective_function_value::pobjfun)=0;
00513           df1b2variable pen=0.0;
00514           df1b2variable zz=0.0;
00515 
00516           initial_df1b2params::reset(y,pen);
00517           pfmin->user_function();
00518 
00519           re_objective_function_value::fun_without_pen=
00520             value(*re_objective_function_value::pobjfun);
00521 
00522           (*re_objective_function_value::pobjfun)+=pen;
00523           (*re_objective_function_value::pobjfun)+=zz;
00524 
00525           set_dependent_variable(*re_objective_function_value::pobjfun);
00526           df1b2_gradlist::set_no_derivatives();
00527           df1b2variable::passnumber=1;
00528           df1b2_gradcalc1();
00529         }
00530 
00531         for (int i=1;i<=usize;i++)
00532         {
00533           for (int j=jmin;j<=jmax;j++)
00534           {
00535             //Hess(i,j-xsize)=y(i+xsize).u_bar[j-mind];
00536             y(i+xsize).get_u_bar_tilde()[j-mind]=Hessadjoint(i,j-xsize);
00537           }
00538         }
00539 
00540         //int mind=y(1).minder;
00541         df1b2variable::passnumber=2;
00542         df1b2_gradcalc1();
00543 
00544         df1b2variable::passnumber=3;
00545         df1b2_gradcalc1();
00546 
00547         f1b2gradlist->reset();
00548         f1b2gradlist->list.initialize();
00549         f1b2gradlist->list2.initialize();
00550         f1b2gradlist->list3.initialize();
00551         f1b2gradlist->nlist.initialize();
00552         f1b2gradlist->nlist2.initialize();
00553         f1b2gradlist->nlist3.initialize();
00554       }
00555 
00556       if (ad_comm::time_flag)
00557       {
00558         if (ad_comm::ptm)
00559         {
00560           double time=ad_comm::ptm->get_elapsed_time_and_reset();
00561           if (ad_comm::global_logfile)
00562           {
00563             (*ad_comm::global_logfile) << " time for 3rd derivatives "
00564               << time << endl;
00565           }
00566         }
00567       }
00568 
00569       dvector dtmp(1,xsize);
00570       for (int i=1;i<=xsize;i++)
00571       {
00572         dtmp(i)=*y(i).get_u_tilde();
00573       }
00574       if (initial_df1b2params::separable_flag)
00575       {
00576 #ifndef OPT_LIB
00577         assert(nvar <= INT_MAX);
00578 #endif
00579         dvector scale(1,(int)nvar);   // need to get scale from somewhere
00580         /*int check=*/initial_params::stddev_scale(scale,x);
00581         dvector sscale=scale(1,Dux(1).indexmax());
00582         for (int i=1;i<=usize;i++)
00583         {
00584           Dux(i)=elem_prod(Dux(i),sscale);
00585         }
00586         dtmp=elem_prod(dtmp,sscale);
00587       }
00588 
00589       for (int i=1;i<=xsize;i++)
00590       {
00591         xadjoint(i)+=dtmp(i);
00592       }
00593       for (int i=1;i<=usize;i++)
00594         uadjoint(i)+=*y(xsize+i).get_u_tilde();
00595     }
00596    // *****************************************************************
00597    // new stuff to deal with quadraticprior
00598    // *****************************************************************
00599 
00600       int xstuff=3;
00601       if (xstuff && df1b2quadratic_prior::get_num_quadratic_prior()>0)
00602       {
00603         initial_params::straight_through_flag=0;
00604         funnel_init_var::lapprox=0;
00605         block_diagonal_flag=0;
00606 #ifndef OPT_LIB
00607         assert(nvar <= INT_MAX);
00608 #endif
00609         dvector scale1(1,(int)nvar);   // need to get scale from somewhere
00610         initial_params::set_inactive_only_random_effects();
00611         /*int check=*/initial_params::stddev_scale(scale1,x);
00612 
00613         laplace_approximation_calculator::where_are_we_flag=3;
00614         quadratic_prior::in_qp_calculations=1;
00615         funnel_init_var::lapprox=this;
00616         df1b2_gradlist::set_no_derivatives();
00617         dvector scale(1,(int)nvar);   // need to get scale from somewhere
00618         /*check=*/initial_params::stddev_scale(scale,x);
00619         dvector sscale=scale(1,Dux(1).indexmax());
00620 
00621         for (int i=1;i<=usize;i++)
00622         {
00623           Dux(i)=elem_div(Dux(i),sscale);
00624         }
00625 
00626         if (xstuff>1)
00627         {
00628           df1b2quadratic_prior::get_Lxu_contribution(Dux);
00629         }
00630         quadratic_prior::in_qp_calculations=0;
00631         funnel_init_var::lapprox=0;
00632         laplace_approximation_calculator::where_are_we_flag=0;
00633 
00634         for (int i=1;i<=usize;i++)
00635         {
00636           Dux(i)=elem_prod(Dux(i),sscale);
00637         }
00638         //local_dtemp=elem_prod(local_dtemp,sscale);
00639 
00640         if (xstuff>2)
00641         {
00642           dvector tmp=evaluate_function_with_quadprior(x,usize,pfmin);
00643           for (int i=1;i<=xsize;i++)
00644           {
00645             xadjoint(i)+=tmp(i);
00646           }
00647         }
00648 
00649         if (xstuff>2)
00650         {
00651           quadratic_prior::get_cHessian_contribution_from_vHessian(Hess,xsize);
00652         }
00653       }
00654 
00655    // *****************************************************************
00656    // new stuff to deal with quadraticprior
00657    // *****************************************************************
00658     if (ad_comm::ptm)
00659     {
00660       /*double time=*/ad_comm::ptm->get_elapsed_time_and_reset();
00661     }
00662 
00663   #if defined(USE_ATLAS)
00664         if (!ad_comm::no_atlas_flag)
00665         {
00666           //xadjoint -= uadjoint*atlas_solve_spd_trans(Hess,Dux);
00667           xadjoint -= atlas_solve_spd_trans(Hess,uadjoint)*Dux;
00668         }
00669         else
00670         {
00671           //xadjoint -= uadjoint*solve(Hess,Dux);
00672           xadjoint -= solve(Hess,uadjoint)*Dux;
00673         }
00674   #else
00675         //xadjoint -= uadjoint*solve(Hess,Dux);
00676         xadjoint -= solve(Hess,uadjoint)*Dux;
00677   #endif
00678 
00679     if (ad_comm::ptm)
00680     {
00681       double time=ad_comm::ptm->get_elapsed_time_and_reset();
00682       if (ad_comm::global_logfile)
00683       {
00684         (*ad_comm::global_logfile) << " Time in second solve "
00685           << time << endl;
00686       }
00687     }
00688     if (ad_comm::ptm1)
00689     {
00690       double time=ad_comm::ptm1->get_elapsed_time_and_reset();
00691       if (ad_comm::global_logfile)
00692       {
00693         (*ad_comm::global_logfile) << " Total time in function evaluation "
00694           << time << endl << endl;
00695       }
00696     }
00697   }
00698   return xadjoint;
00699 }
00700 
00705 void laplace_approximation_calculator::get_newton_raphson_info
00706   (function_minimizer * pfmin)
00707 {
00708   if (ad_comm::time_flag)
00709   {
00710     if (ad_comm::ptm)
00711     {
00712         (*ad_comm::global_logfile) << " Starting Newton-Raphson "
00713           <<  endl;
00714     }
00715   }
00716 
00717   for (int ip=1;ip<=num_der_blocks;ip++)
00718   {
00719     df1b2variable::minder=minder(ip);
00720     df1b2variable::maxder=maxder(ip);
00721 
00722     set_u_dot(ip);
00723 
00724     // do we need to reallocate memory for df1b2variables?
00725     check_for_need_to_reallocate(ip);
00726     df1b2_gradlist::set_yes_derivatives();
00727     //cout << re_objective_function_value::pobjfun << endl;
00728     //cout << re_objective_function_value::pobjfun->ptr << endl;
00729     (*re_objective_function_value::pobjfun)=0;
00730     df1b2variable pen=0.0;
00731     df1b2variable zz=0.0;
00732     initial_df1b2params::reset(y,pen);
00733      //cout << setprecision(15) << y << endl;
00734     if (initial_df1b2params::separable_flag)
00735     {
00736       Hess.initialize();
00737       grad.initialize();
00738     }
00739 
00740     double time1 = 0;
00741     if (ad_comm::time_flag)
00742     {
00743       if (ad_comm::ptm)
00744       {
00745         time1 = ad_comm::ptm->get_elapsed_time();
00746       }
00747     }
00748     pfmin->user_function();
00749     if (ad_comm::time_flag)
00750     {
00751       if (ad_comm::ptm)
00752       {
00753         if (ad_comm::global_logfile)
00754         {
00755           double time=ad_comm::ptm->get_elapsed_time();
00756           (*ad_comm::global_logfile) <<
00757             "       Time in user_function() " <<  ip << "  " << time-time1
00758             << endl;
00759         }
00760       }
00761     }
00762 
00763     re_objective_function_value::fun_without_pen
00764       =value(*re_objective_function_value::pobjfun);
00765 
00766     (*re_objective_function_value::pobjfun)+=pen;
00767     (*re_objective_function_value::pobjfun)+=zz;
00768 
00769     //cout << setprecision(15) << *re_objective_function_value::pobjfun << endl;
00770     //cout << setprecision(15) << pen << endl;
00771     if (!initial_df1b2params::separable_flag)
00772     {
00773       set_dependent_variable(*re_objective_function_value::pobjfun);
00774       df1b2_gradlist::set_no_derivatives();
00775       df1b2variable::passnumber=1;
00776       df1b2_gradcalc1();
00777       int mind=y(1).minder;
00778       int jmin=max(mind,xsize+1);
00779       int jmax=min(y(1).maxder,xsize+usize);
00780       for (int i=1;i<=usize;i++)
00781         for (int j=jmin;j<=jmax;j++)
00782           Hess(i,j-xsize)=y(i+xsize).u_bar[j-mind];
00783      /*
00784       {
00785         ofstream ofs("hessreport");
00786         ofs << setw(6) << Hess << endl << endl;
00787         ofs << setw(10) << setfixed() << setprecision(3) << choleski_decomp(Hess) << endl << endl;
00788         ofs << setw(10) << setfixed() << setprecision(3) << inv(Hess) << endl << endl;
00789         ofs << setw(10) << setfixed() << setprecision(3) << choleski_decomp(inv(Hess)) << endl << endl;
00790       }
00791       */
00792     // ****************************************************************
00793     // ****************************************************************
00794     // ****************************************************************
00795     // temporary shit
00796      /*
00797       for (i=1;i<=usize;i++)
00798       {
00799         for (j=jmin;j<=jmax;j++)
00800         {
00801           //Hess(i,j-xsize)=y(i+xsize).u_bar[j-mind];
00802           y(i+xsize).get_u_bar_tilde()[j-mind]=Hessadjoint(i,j-xsize);
00803         }
00804       }
00805       cout << Hess << endl;
00806       df1b2variable::passnumber=2;
00807       df1b2_gradcalc1();
00808 
00809       df1b2variable::passnumber=3;
00810       df1b2_gradcalc1();
00811      */
00812     // ****************************************************************
00813     // ****************************************************************
00814     // ****************************************************************
00815     // ****************************************************************
00816       for (int j=jmin;j<=jmax;j++)
00817         grad(j-xsize)= re_objective_function_value::pobjfun->u_dot[j-mind];
00818     }
00819     if (ip<num_der_blocks)
00820       f1b2gradlist->reset();
00821   }
00822 
00823 
00824   // just to match master pvm routine
00825   if (ad_comm::time_flag)
00826   {
00827     if (ad_comm::ptm)
00828     {
00829       /*double time=*/ad_comm::ptm->get_elapsed_time();
00830     }
00831   }
00832 }
00833 
00838 void laplace_approximation_calculator::set_u_dot(int ip)
00839 {
00840   int mmin=y.indexmin();
00841   int mmax=y.indexmax();
00842 
00843   for (int i=mmin;i<=mmax;i++)
00844   {
00845     y(i).set_u_dot();
00846   }
00847 }
00848 
00853 void laplace_approximation_calculator::check_pool_size(void)
00854 {
00855   unsigned int num_active_parameters = nvar;
00856   f1b2gradlist->reset();
00857 
00858   adpool * tmppool=df1b2variable::pool;
00859   if (tmppool)
00860   {
00861     //cout << tmppool << endl;
00862     // check if current pool is the right size
00863     if (tmppool->nvar != num_active_parameters)
00864     {
00865       // check sizes of other pools
00866       int found_pool_flag=0;
00867       for (int i=0;i<df1b2variable::adpool_counter;i++)
00868       {
00869         if (df1b2variable::adpool_vector[i]->nvar == num_active_parameters)
00870         {
00871           adpool * tmp = df1b2variable::pool;
00872           //cout << "Old Depth = " << df1b2variable::pool->depth_check()
00873            //    << "  "  << df1b2variable::pool  << "  ";
00874           df1b2variable::pool=df1b2variable::adpool_vector[i];
00875           //cout << "Depth = " << df1b2variable::pool->depth_check()
00876            //    << "  "  << df1b2variable::pool  << endl;
00877           df1b2variable::adpool_vector[i]=tmp;
00878           df1b2variable::nvar_vector[i]=df1b2variable::pool->nvar;
00879           found_pool_flag=1;
00880           break;
00881         }
00882       }
00883       if (!found_pool_flag)
00884       {
00885         if (df1b2variable::adpool_counter>=df1b2variable::adpool_vectorsize)
00886         {
00887            cerr << "Need to increase adpool_vectorsize" << endl;
00888            ad_exit(1);
00889         }
00890         df1b2variable::adpool_vector[df1b2variable::adpool_counter]=
00891           df1b2variable::pool;
00892         df1b2variable::nvar_vector[df1b2variable::adpool_counter]=
00893           df1b2variable::pool->nvar;
00894         //df1b2variable::adpool_counter++;
00895         df1b2variable::increment_adpool_counter();
00896         df1b2variable::pool=new adpool();
00897         if (!df1b2variable::pool)
00898         {
00899           cerr << "Memory allocation error" << endl;
00900           ad_exit(1);
00901         }
00902       }
00903     }
00904   }
00905   else
00906   {
00907     df1b2variable::pool=new adpool();
00908     if (!df1b2variable::pool)
00909     {
00910       cerr << "Memory allocation error" << endl;
00911       ad_exit(1);
00912     }
00913   }
00914   df1b2variable::nvar = num_active_parameters;
00915   df1b2variable::set_blocksize();
00916 }