ADMB Documentation  11.1x.2704
 All Classes Files Functions Variables Typedefs Friends Defines
df1b2lp6.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: df1b2lp6.cpp 2637 2014-11-13 05:58:08Z johnoel $
00003  *
00004  * Author: David Fournier
00005  * Copyright (c) 2008-2012 Regents of the University of California
00006  */
00011 //#define USE_DD_STUFF
00012 //#define USE_DD
00013 #include <admodel.h>
00014 #include <df1b2fun.h>
00015 #include <adrndeff.h>
00016 #ifndef OPT_LIB
00017   #include <cassert>
00018   #include <climits>
00019 #endif
00020 
00021 double calculate_laplace_approximation(const dvector& x,const dvector& u0,
00022   const banded_symmetric_dmatrix& bHess,const dvector& _xadjoint,
00023   const dvector& _uadjoint,
00024   const banded_symmetric_dmatrix& _bHessadjoint,function_minimizer * pmin);
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_shess(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 int use_dd_nr=0;
00034 int admb_ssflag=0;
00035 dvector solve(const dmatrix & st,const dmatrix & Hess,
00036   const dvector& grad);
00037 
00038 #if defined(USE_DD_STUFF)
00039 #  if defined(_MSC_VER)
00040     extern "C" _export  void dd_newton_raphson(int n,double * v,double * diag,
00041       double * ldiag, double * yy);
00042 #  else
00043     extern "C" void dd_newton_raphson(int n,double * v,double * diag,
00044       double * ldiag, double * yy);
00045 #  endif
00046 #endif
00047 
00052 void laplace_approximation_calculator::
00053   do_newton_raphson_banded(function_minimizer * pfmin,double f_from_1,
00054   int& no_converge_flag)
00055 {
00056   //quadratic_prior * tmpptr=quadratic_prior::ptr[0];
00057   //cout << tmpptr << endl;
00058 
00059 
00060   laplace_approximation_calculator::where_are_we_flag=2;
00061   double maxg=fabs(evaluate_function(uhat,pfmin));
00062 
00063 
00064   laplace_approximation_calculator::where_are_we_flag=0;
00065   dvector uhat_old(1,usize);
00066   for(int ii=1;ii<=num_nr_iters;ii++)
00067   {
00068     // test newton raphson
00069     switch(hesstype)
00070     {
00071     case 3:
00072       bHess->initialize();
00073       break;
00074     case 4:
00075       Hess.initialize();
00076       break;
00077     default:
00078       cerr << "Illegal value for hesstype here" << endl;
00079       ad_exit(1);
00080     }
00081 
00082     grad.initialize();
00083     //int check=initial_params::stddev_scale(scale,uhat);
00084     //check=initial_params::stddev_curvscale(curv,uhat);
00085     //max_separable_g=0.0;
00086     sparse_count = 0;
00087 
00088     step=get_newton_raphson_info_banded(pfmin);
00089     //if (bHess)
00090      // cout << "norm(*bHess) = " << norm(*bHess) << endl;
00091     //cout << "norm(Hess) = " << norm(Hess) << endl;
00092     //cout << grad << endl;
00093     //check_pool_depths();
00094     if (!initial_params::mc_phase)
00095       cout << "Newton raphson " << ii << "  ";
00096     if (quadratic_prior::get_num_quadratic_prior()>0)
00097     {
00098       quadratic_prior::get_cHessian_contribution(Hess,xsize);
00099       quadratic_prior::get_cgradient_contribution(grad,xsize);
00100     }
00101 
00102     int ierr=0;
00103     if (hesstype==3)
00104     {
00105       if (use_dd_nr==0)
00106       {
00107         banded_lower_triangular_dmatrix bltd=choleski_decomp(*bHess,ierr);
00108         if (ierr && no_converge_flag ==0)
00109         {
00110           no_converge_flag=1;
00111           //break;
00112         }
00113         if (ierr)
00114         {
00115           double oldval;
00116           evaluate_function(oldval,uhat,pfmin);
00117           uhat=banded_calculations_trust_region_approach(uhat,pfmin);
00118         }
00119         else
00120         {
00121           if (dd_nr_flag==0)
00122           {
00123             dvector v=solve(bltd,grad);
00124             step=-solve_trans(bltd,v);
00125             //uhat_old=uhat;
00126             uhat+=step;
00127           }
00128           else
00129           {
00130 #if defined(USE_DD_STUFF)
00131             int n=grad.indexmax();
00132             maxg=fabs(evaluate_function(uhat,pfmin));
00133             uhat=dd_newton_raphson2(grad,*bHess,uhat);
00134 #else
00135             cerr << "high precision Newton Raphson not implemented" << endl;
00136             ad_exit(1);
00137 #endif
00138           }
00139           maxg=fabs(evaluate_function(uhat,pfmin));
00140           if (f_from_1< pfmin->lapprox->fmc1.fbest)
00141           {
00142             uhat=banded_calculations_trust_region_approach(uhat,pfmin);
00143             maxg=fabs(evaluate_function(uhat,pfmin));
00144           }
00145         }
00146       }
00147       else
00148       {
00149         cout << "error not used" << endl;
00150         ad_exit(1);
00151        /*
00152         banded_symmetric_ddmatrix bHessdd=banded_symmetric_ddmatrix(*bHess);
00153         ddvector gradd=ddvector(grad);
00154         //banded_lower_triangular_ddmatrix bltdd=choleski_decomp(bHessdd,ierr);
00155         if (ierr && no_converge_flag ==0)
00156         {
00157           no_converge_flag=1;
00158           break;
00159         }
00160         if (ierr)
00161         {
00162           double oldval;
00163           evaluate_function(oldval,uhat,pfmin);
00164           uhat=banded_calculations_trust_region_approach(uhat,pfmin);
00165           maxg=fabs(evaluate_function(uhat,pfmin));
00166         }
00167         else
00168         {
00169           ddvector v=solve(bHessdd,gradd);
00170           step=-make_dvector(v);
00171           //uhat_old=uhat;
00172           uhat=make_dvector(ddvector(uhat)+step);
00173           maxg=fabs(evaluate_function(uhat,pfmin));
00174           if (f_from_1< pfmin->lapprox->fmc1.fbest)
00175           {
00176             uhat=banded_calculations_trust_region_approach(uhat,pfmin);
00177             maxg=fabs(evaluate_function(uhat,pfmin));
00178           }
00179         }
00180         */
00181       }
00182 
00183       if (maxg < 1.e-13)
00184       {
00185         break;
00186       }
00187     }
00188     else if (hesstype==4)
00189     {
00190       dvector step;
00191 
00192 #     if defined(USE_ATLAS)
00193         if (!ad_comm::no_atlas_flag)
00194         {
00195           step=-atlas_solve_spd(Hess,grad,ierr);
00196         }
00197         else
00198         {
00199           dmatrix A=choleski_decomp_positive(Hess,ierr);
00200           if (!ierr)
00201           {
00202             step=-solve(Hess,grad);
00203             //step=-solve(A*trans(A),grad);
00204           }
00205         }
00206         if (!ierr) break;
00207 #     else
00208         if (sparse_hessian_flag)
00209         {
00210           //step=-solve(*sparse_triplet,Hess,grad,*sparse_symbolic);
00211           dvector temp=solve(*sparse_triplet2,grad,*sparse_symbolic2,ierr);
00212           if (ierr)
00213           {
00214             step=-temp;
00215           }
00216           else
00217           {
00218             cerr << "matrix not pos definite in sparse choleski"  << endl;
00219             pfmin->bad_step_flag=1;
00220 
00221             int on;
00222             int nopt;
00223             if ((on=option_match(ad_comm::argc,ad_comm::argv,"-ieigvec",nopt))
00224               >-1)
00225             {
00226               dmatrix M=make_dmatrix(*sparse_triplet2);
00227 
00228               ofstream ofs3("inner-eigvectors");
00229               ofs3 << "eigenvalues and eigenvectors " << endl;
00230               dvector v=eigenvalues(M);
00231               dmatrix ev=trans(eigenvectors(M));
00232               ofs3 << "eigenvectors" << endl;
00233               int i;
00234               for (i=1;i<=ev.indexmax();i++)
00235                {
00236                   ofs3 << setw(4) << i  << " "
00237                    << setshowpoint() << setw(14) << setprecision(10) << v(i)
00238                    << " "
00239                    << setshowpoint() << setw(14) << setprecision(10)
00240                    << ev(i) << endl;
00241                }
00242             }
00243           }
00244           //cout << norm2(step-tmpstep) << endl;
00245           //dvector step1=-solve(Hess,grad);
00246           //cout << norm2(step-step1) << endl;
00247         }
00248         else
00249         {
00250           step=-solve(Hess,grad);
00251         }
00252 #     endif
00253       if (pmin->bad_step_flag)
00254         break;
00255       uhat_old=uhat;
00256       uhat+=step;
00257 
00258       double maxg_old=maxg;
00259       maxg=fabs(evaluate_function(uhat,pfmin));
00260       if (maxg>maxg_old)
00261       {
00262         uhat=uhat_old;
00263         evaluate_function(uhat,pfmin);
00264         break;
00265       }
00266       if (maxg < 1.e-13)
00267       {
00268         break;
00269       }
00270     }
00271 
00272     if (sparse_hessian_flag==0)
00273     {
00274       for (int i=1;i<=usize;i++)
00275       {
00276         y(i+xsize)=uhat(i);
00277       }
00278     }
00279     else
00280     {
00281       for (int i=1;i<=usize;i++)
00282       {
00283         value(y(i+xsize))=uhat(i);
00284       }
00285     }
00286   }
00287 }
00288 
00293 double laplace_approximation_calculator::
00294   inner_optimization_banded(/*dvector& uhat,*/ dvector& x,
00295   function_minimizer * pfmin,int& no_converge_flag)
00296 {
00297   if (no_converge_flag)
00298   {
00299     no_converge_flag=0;
00300   }
00301 
00302   if (!inner_lmnflag)
00303   {
00304     if (!ADqd_flag)
00305     {
00306       uhat=get_uhat_quasi_newton(x,pfmin);
00307       double maxg=fabs(fmc1.gmax);
00308       if (maxg>1.0)
00309       {
00310         uhat=get_uhat_quasi_newton(x,pfmin);
00311       }
00312     }
00313     else
00314     {
00315       uhat=get_uhat_quasi_newton_qd(x,pfmin);
00316     }
00317   }
00318   else
00319   {
00320     uhat=get_uhat_lm_newton(x,pfmin);
00321     //uhat=get_uhat_lm_newton2(x,pfmin);
00322     //maxg=objective_function_value::gmax;
00323   }
00324   return fmc1.fbest;
00325 }
00326 
00331 dvector laplace_approximation_calculator::banded_calculations
00332   (const dvector& _x,const double& _f,function_minimizer * pfmin)
00333 {
00334   // for use when there is no separability
00335   ADUNCONST(dvector,x)
00336   ADUNCONST(double,f)
00337   //int i,j;
00338   int i;
00339 
00340   initial_params::set_inactive_only_random_effects();
00341   gradient_structure::set_NO_DERIVATIVES();
00342   initial_params::reset(x);    // get current x values into the model
00343   gradient_structure::set_YES_DERIVATIVES();
00344 
00345   initial_params::set_active_only_random_effects();
00346   if (init_switch==0)
00347   {
00348     gradient_structure::set_NO_DERIVATIVES();
00349     initial_params::xinit(ubest);
00350     gradient_structure::set_YES_DERIVATIVES();
00351   }
00352   //double maxg;
00353   //double maxg_save;
00354   //dvector uhat(1,usize);
00355   double f_from_1=0.0;
00356 
00357   int no_converge_flag=0;
00358 
00359   // this is the main loop to do inner optimization
00360   for (;;)
00361   {
00362     int icount=0;
00363     do
00364     {
00365       icount++;
00366       // do the inner optimziation
00367       if (inner_maxfn>0)
00368       {
00369         f_from_1=inner_optimization_banded(/*uhat,*/ x,pfmin,
00370           no_converge_flag);
00371       }
00372 
00373       if (sparse_hessian_flag==0)
00374       {
00375         for (i=1;i<=xsize;i++) { y(i)=x(i); }
00376         for (i=1;i<=usize;i++) { y(i+xsize)=uhat(i); }
00377       }
00378       else
00379       {
00380         for (i=1;i<=xsize;i++) { value(y(i))=x(i); }
00381         for (i=1;i<=usize;i++) { value(y(i+xsize))=uhat(i); }
00382       }
00383 
00384       laplace_approximation_calculator::where_are_we_flag=2;
00385       if (admb_ssflag==0)
00386       {
00387         do_newton_raphson_banded(pfmin,f_from_1,no_converge_flag);
00388       }
00389       else
00390       {
00391         do_newton_raphson_state_space(pfmin,f_from_1,no_converge_flag);
00392       }
00393       laplace_approximation_calculator::where_are_we_flag=0;
00394 
00395       if (num_nr_iters<=0) { evaluate_function(uhat,pfmin); }
00396 
00397       if (sparse_hessian_flag==0)
00398       {
00399         for (i=1;i<=usize;i++) { y(i+xsize)=uhat(i); }
00400       }
00401       else
00402       {
00403         for (i=1;i<=usize;i++) { value(y(i+xsize))=uhat(i); }
00404       }
00405       if (icount>2) pfmin->bad_step_flag=1;
00406       if (pfmin->bad_step_flag)
00407         return xadjoint;
00408     }
00409     while(no_converge_flag);
00410 
00411     /* If we are in mcmc phase we just need to calcualte the
00412        ln_det(Hess) and return
00413     */
00414     hs_symbolic & ssymb=*(pmin->lapprox->sparse_symbolic2);
00415     if (initial_params::mc_phase)
00416     {
00417       do_newton_raphson_banded(pfmin,f_from_1,no_converge_flag);
00418       int sgn=0;
00419       double& f = (double&) _f;
00420       f=initial_df1b2params::cobjfun;
00421      if (pmin->lapprox->sparse_hessian_flag==0)
00422      {
00423         if (bHess==0)
00424         {
00425           cerr << "Block diagonal Hessian is unallocated" << endl;
00426           ad_exit(1);
00427         }
00428         f+=0.5*ln_det_choleski(*bHess,sgn);
00429       }
00430       else
00431       {
00432         //hs_symbolic & ssymb=*(pmin->lapprox->sparse_symbolic2);
00433         //dvariable tmp=0.5*ln_det(*(pmin->lapprox->vsparse_triplet),
00434         //  ssymb,*(pmin->lapprox->sparse_triplet2));
00435         f+=0.5*ln_det(*(pmin->lapprox->sparse_triplet2),ssymb);
00436       }
00437     }
00438     else
00439     {
00440       xadjoint.initialize();
00441       uadjoint.initialize();
00442       Dux.initialize();
00443 
00444       if (hesstype==3)
00445         bHess->initialize();
00446       else if (hesstype==4)
00447         Hess.initialize();
00448 
00449       block_diagonal_flag=2;
00450       used_flags.initialize();
00451       funnel_init_var::lapprox=this;
00452       sparse_count = 0;
00453 
00454       initial_params::straight_through_flag=1;
00455 
00456       if (sparse_triplet2)
00457         sparse_triplet2->initialize();
00458 
00459       pfmin->user_function();
00460       initial_params::straight_through_flag=0;
00461 
00462       int ierr=0;
00463 
00464       laplace_approximation_calculator::where_are_we_flag=3;
00465       if (!ierr)
00466       {
00467         if (num_importance_samples==0)
00468         {
00469           if (hesstype==3)
00470           {
00471             f=calculate_laplace_approximation(x,uhat,*bHess,xadjoint,uadjoint,
00472               *bHessadjoint,pfmin);
00473           }
00474           else if (hesstype==4)
00475           {
00476             //cout << "Hess" << endl << Hess << endl;
00477             f=calculate_laplace_approximation(x,uhat,Hess,xadjoint,uadjoint,
00478               Hessadjoint,pfmin);
00479           }
00480           else
00481           {
00482             cerr << "Error in hesstype" << endl;
00483             ad_exit(1);
00484           }
00485         }
00486         else
00487         {
00488           if (hesstype==3)
00489           {
00490             //cerr << "Error not implemented yet" << endl;
00491             //ad_exit(1);
00492             f=calculate_importance_sample(x,uhat,*bHess,xadjoint,uadjoint,
00493               *bHessadjoint,pfmin);
00494           }
00495           else if (hesstype==4)
00496           {
00497             if (pmin->lapprox->sparse_hessian_flag==0)
00498               f=calculate_importance_sample(x,uhat,Hess,xadjoint,uadjoint,
00499                 Hessadjoint,pfmin);
00500             else
00501               f=calculate_importance_sample_shess(x,uhat,Hess,xadjoint,uadjoint,
00502                 Hessadjoint,pfmin);
00503           }
00504           else
00505           {
00506             cerr << "Error in hesstype" << endl;
00507             ad_exit(1);
00508           }
00509         }
00510       }
00511       else
00512       {
00513         f=1.e+30;
00514       }
00515 
00516       // set flag for thrid erivatvies and call function again because
00517       // stack is wiped out
00518 
00519       if (hesstype==3)
00520       {
00521         bHess->initialize();
00522       }
00523       else if (hesstype==4)
00524       {
00525         if (sparse_hessian_flag==0)
00526         {
00527           Hess.initialize();
00528         }
00529         else
00530         {
00531           sparse_triplet2->initialize();
00532         }
00533       }
00534       else
00535       {
00536         cerr << "Illegal value for hesstype" << endl;
00537         ad_exit(1);
00538       }
00539       initial_params::straight_through_flag=1;
00540       block_diagonal_flag=3;
00541       local_dtemp.initialize();
00542 
00543       // *****  Note for quadratic prior code: I don't think that this
00544       // part gets added to the Hessian here.
00545       sparse_count=0;
00546       sparse_count_adjoint=0;
00547       pfmin->user_function();
00548 
00549       // *** Hessian calculated just above did not have quadratic prior
00550       // in it so can save this part for quadratci prioer adjoint calculations
00551       if (quadratic_prior::get_num_quadratic_prior()>0)
00552       {
00553         if (pHess_non_quadprior_part)
00554         {
00555           if (pHess_non_quadprior_part->indexmax() != Hess.indexmax())
00556           {
00557             delete pHess_non_quadprior_part;
00558             pHess_non_quadprior_part=0;
00559           }
00560         }
00561         if (!pHess_non_quadprior_part)
00562         {
00563           pHess_non_quadprior_part=new dmatrix(1,usize,1,usize);
00564           if (!pHess_non_quadprior_part)
00565           {
00566             cerr << "Error allocating memory for Hesssian part" << endl;
00567             ad_exit(1);
00568           }
00569         }
00570         (*pHess_non_quadprior_part)=Hess;
00571       }
00572 
00573       block_diagonal_flag=0;
00574       initial_params::straight_through_flag=1;
00575 
00576       //dmatrix tHess=dmatrix(*bHess);
00577       initial_params::straight_through_flag=0;
00578       funnel_init_var::lapprox=0;
00579       //cout << initial_df1b2params::cobjfun << endl;
00580       //f=initial_df1b2params::cobjfun;
00581       block_diagonal_flag=0;
00582 #ifndef OPT_LIB
00583       assert(nvar <= INT_MAX);
00584 #endif
00585       dvector scale1(1,(int)nvar);   // need to get scale from somewhere
00586       initial_params::set_inactive_only_random_effects();
00587       /*int check=*/initial_params::stddev_scale(scale1,x);
00588       //for (i=1;i<=xadjoint.indexmax();i++)
00589       //  xadjoint(i)*=scale1(i);
00590       laplace_approximation_calculator::where_are_we_flag=0;
00591 
00592       if (df1b2quadratic_prior::get_num_quadratic_prior()>0)
00593       {
00594        // !!!! need to fix this!!!!!!!!!!!!!!!!!!!!!!!
00595         laplace_approximation_calculator::where_are_we_flag=3;
00596         quadratic_prior::in_qp_calculations=1;
00597         funnel_init_var::lapprox=this;
00598         df1b2_gradlist::set_no_derivatives();
00599         df1b2quadratic_prior::get_Lxu_contribution(Dux);
00600         quadratic_prior::in_qp_calculations=0;
00601         funnel_init_var::lapprox=0;
00602         laplace_approximation_calculator::where_are_we_flag=0;
00603       }
00604       if (initial_df1b2params::separable_flag)
00605       {
00606         dvector scale(1,(int)nvar);   // need to get scale from somewhere
00607         /*int check=*/initial_params::stddev_scale(scale,x);
00608         dvector sscale=scale(1,Dux(1).indexmax());
00609         for (i=1;i<=usize;i++)
00610         {
00611           Dux(i)=elem_prod(Dux(i),sscale);
00612         }
00613         local_dtemp=elem_prod(local_dtemp,sscale);
00614       }
00615       //cout << trans(Dux)(1) << endl;
00616       //cout << trans(Dux)(3) << endl;
00617       if (quadratic_prior::get_num_quadratic_prior()>0)
00618       {
00619         dvector tmp=evaluate_function_with_quadprior(x,usize,pfmin);
00620         local_dtemp+=tmp;
00621       }
00622 
00623       for (i=1;i<=xsize;i++)
00624       {
00625         xadjoint(i)+=local_dtemp(i);
00626       }
00627       if (df1b2quadratic_prior::get_num_quadratic_prior()>0)
00628       {
00629        // !!!! need to fix this!!!!!!!!!!!!!!!!!!!!!!!
00630         quadratic_prior::get_cHessian_contribution_from_vHessian(Hess,xsize);
00631       }
00632 
00633       if (hesstype==3)
00634       {
00635         //xadjoint -= uadjoint*solve(*bHess,Dux);
00636         if (bHess_pd_flag==0)
00637         {
00638           xadjoint -= solve(*bHess,uadjoint)*Dux;
00639         }
00640       }
00641       else if (hesstype==4)
00642       {
00643         if (sparse_hessian_flag)
00644         {
00645 //dvector tmp=solve(*sparse_triplet,Hess,uadjoint,*sparse_symbolic)*Dux;
00646           dvector tmp=solve(*sparse_triplet2,uadjoint,*sparse_symbolic2)*Dux;
00647           xadjoint -= tmp;
00648         }
00649         else
00650         {
00651           xadjoint -= solve(Hess,uadjoint)*Dux;
00652         }
00653       }
00654     }
00655     if (bHess_pd_flag==0) break;
00656   }
00657 
00658   return xadjoint;
00659 }
00660   //int check_pool_flag1=0;
00661 
00666 void laplace_approximation_calculator::
00667   do_separable_stuff_newton_raphson_banded(df1b2variable& ff)
00668 {
00669   //***********************************************************
00670   //***********************************************************
00671   set_dependent_variable(ff);
00672   df1b2_gradlist::set_no_derivatives();
00673   df1b2variable::passnumber=1;
00674   //if (check_pool_flag1)
00675    // check_pool_depths();
00676   df1b2_gradcalc1();
00677   //if (check_pool_flag1)
00678   //  check_pool_depths();
00679 
00680   init_df1b2vector & locy= *funnel_init_var::py;
00681   imatrix& list=*funnel_init_var::plist;
00682   int num_local_re=0;
00683   int num_fixed_effects=0;
00684 
00685 #ifndef OPT_LIB
00686   assert(funnel_init_var::num_active_parameters <= INT_MAX);
00687 #endif
00688   ivector lre_index(1, (int)funnel_init_var::num_active_parameters);
00689   ivector lfe_index(1, (int)funnel_init_var::num_active_parameters);
00690 
00691   for (int i=1;i<=(int)funnel_init_var::num_active_parameters;i++)
00692   {
00693     if (list(i,1)>xsize)
00694     {
00695       lre_index(++num_local_re)=i;
00696     }
00697     else if (list(i,1)>0)
00698     {
00699       lfe_index(++num_fixed_effects)=i;
00700     }
00701   }
00702 
00703   if (num_local_re > 0)
00704   {
00705     switch(hesstype)
00706     {
00707     case 3:
00708       for (int i=1;i<=num_local_re;i++)
00709       {
00710         int lrei=lre_index(i);
00711         for (int j=1;j<=num_local_re;j++)
00712         {
00713           int lrej=lre_index(j);
00714           int i1=list(lrei,1)-xsize;
00715           int i2=list(lrei,2);
00716           int j1=list(lrej,1)-xsize;
00717           int j2=list(lrej,2);
00718           if (i1>=j1) (*bHess)(i1,j1)+=locy(i2).u_bar[j2-1];
00719         }
00720       }
00721       break;
00722     case 4:
00723       if (sparse_hessian_flag==0)
00724       {
00725         for (int i=1;i<=num_local_re;i++)
00726         {
00727           int lrei=lre_index(i);
00728           for (int j=1;j<=num_local_re;j++)
00729           {
00730             int lrej=lre_index(j);
00731             int i1=list(lrei,1)-xsize;
00732             int i2=list(lrei,2);
00733             int j1=list(lrej,1)-xsize;
00734             int j2=list(lrej,2);
00735             Hess(i1,j1)+=locy(i2).u_bar[j2-1];
00736           }
00737         }
00738       }
00739       else
00740       {
00741         for (int i=1;i<=num_local_re;i++)
00742         {
00743           int lrei=lre_index(i);
00744           for (int j=1;j<=num_local_re;j++)
00745           {
00746             int lrej=lre_index(j);
00747             int i1=list(lrei,1)-xsize;
00748             int i2=list(lrei,2);
00749             int j1=list(lrej,1)-xsize;
00750             int j2=list(lrej,2);
00751 
00752             if (i1<=j1)
00753             {
00754               sparse_count++;
00755               (*sparse_triplet2)((*sparse_iterator)(sparse_count))
00756                 +=locy(i2).u_bar[j2-1];
00757             }
00758           }
00759         }
00760       }
00761       break;
00762     default:
00763       cerr << "illegal value for hesstype" << endl;
00764       ad_exit(1);
00765     }
00766 
00767     for (int i=1;i<=num_local_re;i++)
00768     {
00769       int lrei=lre_index(i);
00770       int i1=list(lrei,1);
00771       int i2=list(lrei,2);
00772       //grad(i1-xsize)= re_objective_function_value::pobjfun->u_dot[i2-1];
00773       grad(i1-xsize)+=ff.u_dot[i2-1];
00774     }
00775   }
00776 
00777   f1b2gradlist->reset();
00778   f1b2gradlist->list.initialize();
00779   f1b2gradlist->list2.initialize();
00780   f1b2gradlist->list3.initialize();
00781   f1b2gradlist->nlist.initialize();
00782   f1b2gradlist->nlist2.initialize();
00783   f1b2gradlist->nlist3.initialize();
00784   funnel_init_var::num_vars=0;
00785   funnel_init_var::num_active_parameters=0;
00786   funnel_init_var::num_inactive_vars=0;
00787 }
00788 //int tmp_testcount=0;
00789 df1b2variable * tmp_pen=00;
00790 
00795 dvector laplace_approximation_calculator::
00796   get_newton_raphson_info_banded (function_minimizer * pfmin)
00797 {
00798   int nv=initial_df1b2params::set_index();
00799   if (allocated(used_flags))
00800   {
00801     if (used_flags.indexmax() != nv)
00802     {
00803       used_flags.safe_deallocate();
00804     }
00805   }
00806   if (!allocated(used_flags))
00807   {
00808     used_flags.safe_allocate(1,nv);
00809   }
00810 
00811   for (int ip=1;ip<=num_der_blocks;ip++)
00812   {
00813     if (ip>1)   // change to combine sparse matrix stuff with num der blocks
00814     {           // df  3-4-09
00815       sparse_count=0;
00816     }
00817     used_flags.initialize();
00818     // do we need to reallocate memory for df1b2variables?
00819     check_for_need_to_reallocate(ip);
00820     df1b2_gradlist::set_no_derivatives();
00821     //cout << re_objective_function_value::pobjfun << endl;
00822     //cout << re_objective_function_value::pobjfun->ptr << endl;
00823     (*re_objective_function_value::pobjfun)=0;
00824     df1b2variable pen=0.0;
00825     tmp_pen=&pen;
00826     df1b2variable zz=0.0;
00827 
00828     initial_df1b2params::reset(y,pen);
00829 
00830     // call function to do block diagonal newton-raphson
00831     // the step vector from the newton-raphson is in the vector step
00832     df1b2_gradlist::set_yes_derivatives();
00833 
00834     funnel_init_var::lapprox=this;
00835     //cout << funnel_init_var::lapprox << endl;
00836     block_diagonal_flag=1;
00837    /*
00838     tmp_testcount++;
00839     if (tmp_testcount>=9)
00840     {
00841       pen.deallocate();
00842     }
00843     */
00844 
00845     if (ip==1)
00846     {
00847       if (sparse_triplet2)
00848         sparse_triplet2->initialize();
00849     }
00850 
00851     pfmin->user_function();
00852     /*
00853     if (tmp_testcount>=9)
00854     {
00855       pen.deallocate();
00856     }
00857      */
00858     funnel_init_var::lapprox=0;
00859     block_diagonal_flag=0;
00860     pen.deallocate();
00861   }
00862 
00863   return step;
00864 }
00865 
00870 void laplace_approximation_calculator::
00871   do_separable_stuff_laplace_approximation_banded(df1b2variable& ff)
00872 {
00873   set_dependent_variable(ff);
00874   //df1b2_gradlist::set_no_derivatives();
00875   df1b2variable::passnumber=1;
00876   df1b2_gradcalc1();
00877 
00878   init_df1b2vector & locy= *funnel_init_var::py;
00879   imatrix& list=*funnel_init_var::plist;
00880 
00881   int us=0; int xs=0;
00882 #ifndef OPT_LIB
00883   assert(funnel_init_var::num_active_parameters <= INT_MAX);
00884 #endif
00885   ivector lre_index(1,(int)funnel_init_var::num_active_parameters);
00886   ivector lfe_index(1,(int)funnel_init_var::num_active_parameters);
00887 
00888   for (int i=1;i<=(int)funnel_init_var::num_active_parameters;i++)
00889   {
00890     if (list(i,1)>xsize)
00891     {
00892       lre_index(++us)=i;
00893     }
00894     else if (list(i,1)>0)
00895     {
00896       lfe_index(++xs)=i;
00897     }
00898   }
00899 
00900   for (int j=1;j<=xs;j++)
00901   {
00902     int j1=list(lfe_index(j),1);
00903     int j2=list(lfe_index(j),2);
00904     xadjoint(j1)+=ff.u_dot[j2-1];
00905   }
00906 
00907   if (us>0)
00908   {
00909     if (hesstype==3)
00910     {
00911       for (int i=1;i<=us;i++)
00912       {
00913         for (int j=1;j<=us;j++)
00914         {
00915           int i1=list(lre_index(i),1)-xsize;
00916           int i2=list(lre_index(i),2);
00917           int j1=list(lre_index(j),1)-xsize;
00918           int j2=list(lre_index(j),2);
00919           if (i1>=j1) (*bHess)(i1,j1)+=locy(i2).u_bar[j2-1];
00920         }
00921       }
00922     }
00923     else if (hesstype==4)
00924     {
00925       if (sparse_hessian_flag==0)
00926       {
00927         for (int i=1;i<=us;i++)
00928         {
00929           for (int j=1;j<=us;j++)
00930           {
00931             int i1=list(lre_index(i),1)-xsize;
00932             int i2=list(lre_index(i),2);
00933             int j1=list(lre_index(j),1)-xsize;
00934             int j2=list(lre_index(j),2);
00935             Hess(i1,j1)+=locy(i2).u_bar[j2-1];
00936           }
00937         }
00938       }
00939       else
00940       {
00941         for (int i=1;i<=us;i++)
00942         {
00943           for (int j=1;j<=us;j++)
00944           {
00945             int i1=list(lre_index(i),1)-xsize;
00946             int i2=list(lre_index(i),2);
00947             int j1=list(lre_index(j),1)-xsize;
00948             int j2=list(lre_index(j),2);
00949 
00950             if (i1<=j1)
00951             {
00952               sparse_count++;
00953               (*sparse_triplet2)((*sparse_iterator)(sparse_count))
00954                 +=locy(i2).u_bar[j2-1];
00955             }
00956           }
00957         }
00958       }
00959     }
00960 
00961     for (int i=1;i<=us;i++)
00962     {
00963       int i1=list(lre_index(i),1)-xsize;
00964       int i2=list(lre_index(i),2);
00965       uadjoint(i1)+=ff.u_dot[i2-1];
00966     }
00967 
00968     for (int i=1;i<=us;i++)
00969     {
00970       for (int j=1;j<=xs;j++)
00971       {
00972         int i1=list(lre_index(i),1)-xsize;
00973         int i2=list(lre_index(i),2);
00974         int j1=list(lfe_index(j),1);
00975         int j2=list(lfe_index(j),2);
00976         Dux(i1,j1)+=locy(i2).u_bar[j2-1];
00977       }
00978     }
00979   }
00980   f1b2gradlist->reset();
00981   f1b2gradlist->list.initialize();
00982   f1b2gradlist->list2.initialize();
00983   f1b2gradlist->list3.initialize();
00984   f1b2gradlist->nlist.initialize();
00985   f1b2gradlist->nlist2.initialize();
00986   f1b2gradlist->nlist3.initialize();
00987   funnel_init_var::num_vars=0;
00988   funnel_init_var::num_active_parameters=0;
00989   funnel_init_var::num_inactive_vars=0;
00990 }
00991 
00996 double calculate_laplace_approximation(const dvector& x,const dvector& u0,
00997   const banded_symmetric_dmatrix& bHess,const dvector& _xadjoint,
00998   const dvector& _uadjoint,
00999   const banded_symmetric_dmatrix& _bHessadjoint,function_minimizer * pmin)
01000 {
01001   ADUNCONST(dvector,xadjoint)
01002   ADUNCONST(dvector,uadjoint)
01003   ADUNCONST(banded_symmetric_dmatrix,bHessadjoint)
01004   int bw=bHess.bandwidth();
01005   const int xs=x.size();
01006   const int us=u0.size();
01007   gradient_structure::set_YES_DERIVATIVES();
01008   int nvar=x.size()+u0.size()+((bw+1)*(2*u0.size()-bw))/2;
01009   independent_variables y(1,nvar);
01010 
01011   // need to set random effects active together with whatever
01012   // init parameters should be active in this phase
01013   initial_params::set_inactive_only_random_effects();
01014   initial_params::set_active_random_effects();
01015   /*int onvar=*/initial_params::nvarcalc();
01016   initial_params::xinit(y);    // get the initial values into the
01017   y(1,xs)=x;
01018 
01019  // Here need hooks for sparse matrix structures
01020   int ii=xs+us+1;
01021   for (int i=1;i<=us;i++)
01022   {
01023     int jmin=admax(1,i-bw+1);
01024     for (int j=jmin;j<=i;j++)
01025       y(ii++)=bHess(i,j);
01026   }
01027 
01028   dvar_vector vy=dvar_vector(y);
01029   banded_symmetric_dvar_matrix vHess(1,us,bw);
01030   initial_params::reset(vy);    // get the initial values into the
01031   ii=xs+us+1;
01032   for (int i=1;i<=us;i++)
01033   {
01034     int jmin=admax(1,i-bw+1);
01035     for (int j=jmin;j<=i;j++)
01036       vHess(i,j)=vy(ii++);
01037   }
01038 
01039    dvariable vf=0.0;
01040 
01041    *objective_function_value::pobjfun=0.0;
01042    pmin->AD_uf_outer();
01043    vf+=*objective_function_value::pobjfun;
01044 
01045    int sgn=0;
01046    dvariable ld;
01047 
01048    int eigswitch=0;
01049    if (eigswitch)
01050    {
01051      ofstream ofs("ee");
01052      dvector ev(bHess.indexmin(),bHess.indexmax());
01053      dmatrix evecs=eigenvectors(dmatrix(bHess),ev);
01054      ofs << setprecision(3) << setw(12) << setscientific() << dmatrix(bHess)
01055          << endl << endl;
01056      ofs << ev << endl << endl << evecs << endl;
01057    }
01058    ld=0.5*ln_det_choleski(vHess,sgn);
01059    if (sgn==1)
01060    {
01061      cout << "Choleski failed" << endl;
01062      pmin->lapprox->bHess_pd_flag=1;
01063    }
01064 
01065    vf+=ld;
01066    const double ltp=0.5*log(2.0*PI);
01067    vf-=us*ltp;
01068    double f=value(vf);
01069    dvector g(1,nvar);
01070    gradcalc(nvar,g);
01071 
01072   ii=1;
01073   for (int i=1;i<=xs;i++)
01074     xadjoint(i)=g(ii++);
01075   for (int i=1;i<=us;i++)
01076     uadjoint(i)=g(ii++);
01077   for (int i=1;i<=us;i++)
01078   {
01079     int jmin=admax(1,i-bw+1);
01080     for (int j=jmin;j<=i;j++)
01081       bHessadjoint(i,j)=g(ii++);
01082   }
01083   return f;
01084 }
01085 
01090 dvector laplace_approximation_calculator::
01091   banded_calculations_trust_region_approach(const dvector& _uhat,
01092   function_minimizer * pfmin)
01093 {
01094   dvector& uhat=(dvector&) _uhat;
01095   dvector uhat_old(uhat.indexmin(),uhat.indexmax());
01096   dvector uhat_new(uhat.indexmin(),uhat.indexmax());
01097   dvector uhat_best(uhat.indexmin(),uhat.indexmax());
01098 
01099   double wght=0.0;
01100   double delta=5.e-5;
01101   //do
01102   dvector values(1,300);
01103   double oldfbest=pmin->lapprox->fmc1.fbest;
01104   double newfbest = 0.0;
01105   int have_value=0;
01106   //for (int jj=1;jj<=300;jj++)
01107   int jj=1;
01108   double lastval=oldfbest;
01109   do
01110   {
01111     jj++;
01112     wght+=delta;
01113     //double wght=0.0;
01114     double newval=0.0;
01115     //cout << "Enter weight size " << endl;
01116     //cin >> wght;
01117     if (wght<0.0)
01118       break;
01119     int mmin=bHess->indexmin();
01120     int mmax=bHess->indexmax();
01121     banded_symmetric_dmatrix tmp(mmin,mmax,bHess->bandwidth());
01122     tmp=*bHess;
01123     uhat_old=uhat;
01124     int ierr=0;
01125     for (int i=mmin;i<=mmax;i++)
01126     {
01127       tmp(i,i)+=wght;
01128     }
01129     banded_lower_triangular_dmatrix bltd=choleski_decomp(tmp,ierr);
01130     if (!ierr)
01131     {
01132       dvector v=solve(bltd,grad);
01133       step=-solve_trans(bltd,v);
01134 
01135       uhat_old=uhat;
01136       uhat+=step;
01137       //cout << "norm(uhat_old) = " << norm(uhat_old)
01138        //    << "   norm(uhat) = " << norm(uhat)  << endl;
01139 
01140       /*double maxg=*/evaluate_function(newval,uhat,pfmin);
01141       if (have_value && newval>newfbest)
01142       {
01143         break;
01144       }
01145       if (jj>1)
01146       {
01147         if (newval<lastval)  // we are doing better so increasse step size
01148         {
01149           delta*=2;
01150         }
01151         if (newval>lastval && !have_value)  // we have gone to far go back
01152         {
01153           wght-=delta;
01154           delta/=16;
01155         }
01156       }
01157       lastval=newval;
01158 
01159       if (newval<newfbest)
01160       {
01161         newfbest=newval;
01162         uhat_best=uhat;
01163         have_value=jj;
01164       }
01165       uhat_new=uhat;
01166       uhat=uhat_old;
01167     }
01168     else
01169     {
01170       delta*=2;
01171     }
01172   }
01173   while(jj<10);
01174   if (!have_value)
01175   {
01176     cerr << "can't improve function value in trust region calculations"
01177          << endl;
01178     //ad_exit(1);
01179   }
01180   return uhat_best;
01181   initial_params::set_active_only_random_effects();
01182   if (!inner_lmnflag)
01183   {
01184     if (!ADqd_flag)
01185     {
01186       uhat=get_uhat_quasi_newton(uhat_new,pfmin);
01187       //double maxg=fabs(fmc1.gmax);
01188       //double f_from_1=fmc1.fbest;
01189     }
01190     else
01191     {
01192       uhat=get_uhat_quasi_newton_qd(uhat_new,pfmin);
01193     }
01194   }
01195   else
01196   {
01197     uhat=get_uhat_lm_newton(uhat_new,pfmin);
01198   }
01199   return uhat;
01200 }