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