ADMB Documentation  11.2.2853
 All Classes Files Functions Variables Typedefs Friends Defines
df1b2lp6.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: df1b2lp6.cpp 2796 2014-12-13 20:09:55Z 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)
00424         {
00425           cerr << "Block diagonal Hessian is unallocated" << endl;
00426           ad_exit(1);
00427         }
00428         else
00429         {
00430           f+=0.5*ln_det_choleski(*bHess,sgn);
00431         }
00432       }
00433       else
00434       {
00435         //hs_symbolic & ssymb=*(pmin->lapprox->sparse_symbolic2);
00436         //dvariable tmp=0.5*ln_det(*(pmin->lapprox->vsparse_triplet),
00437         //  ssymb,*(pmin->lapprox->sparse_triplet2));
00438         f+=0.5*ln_det(*(pmin->lapprox->sparse_triplet2),ssymb);
00439       }
00440     }
00441     else
00442     {
00443       xadjoint.initialize();
00444       uadjoint.initialize();
00445       Dux.initialize();
00446 
00447       if (hesstype==3)
00448         bHess->initialize();
00449       else if (hesstype==4)
00450         Hess.initialize();
00451 
00452       block_diagonal_flag=2;
00453       used_flags.initialize();
00454       funnel_init_var::lapprox=this;
00455       sparse_count = 0;
00456 
00457       initial_params::straight_through_flag=1;
00458 
00459       if (sparse_triplet2)
00460         sparse_triplet2->initialize();
00461 
00462       pfmin->user_function();
00463       initial_params::straight_through_flag=0;
00464 
00465       int ierr=0;
00466 
00467       laplace_approximation_calculator::where_are_we_flag=3;
00468       if (!ierr)
00469       {
00470         if (num_importance_samples==0)
00471         {
00472           if (hesstype==3)
00473           {
00474             f=calculate_laplace_approximation(x,uhat,*bHess,xadjoint,uadjoint,
00475               *bHessadjoint,pfmin);
00476           }
00477           else if (hesstype==4)
00478           {
00479             //cout << "Hess" << endl << Hess << endl;
00480             f=calculate_laplace_approximation(x,uhat,Hess,xadjoint,uadjoint,
00481               Hessadjoint,pfmin);
00482           }
00483           else
00484           {
00485             cerr << "Error in hesstype" << endl;
00486             ad_exit(1);
00487           }
00488         }
00489         else
00490         {
00491           if (hesstype==3)
00492           {
00493             //cerr << "Error not implemented yet" << endl;
00494             //ad_exit(1);
00495             f=calculate_importance_sample(x,uhat,*bHess,xadjoint,uadjoint,
00496               *bHessadjoint,pfmin);
00497           }
00498           else if (hesstype==4)
00499           {
00500             if (pmin->lapprox->sparse_hessian_flag==0)
00501               f=calculate_importance_sample(x,uhat,Hess,xadjoint,uadjoint,
00502                 Hessadjoint,pfmin);
00503             else
00504               f=calculate_importance_sample_shess(x,uhat,Hess,xadjoint,uadjoint,
00505                 Hessadjoint,pfmin);
00506           }
00507           else
00508           {
00509             cerr << "Error in hesstype" << endl;
00510             ad_exit(1);
00511           }
00512         }
00513       }
00514       else
00515       {
00516         f=1.e+30;
00517       }
00518 
00519       // set flag for thrid erivatvies and call function again because
00520       // stack is wiped out
00521 
00522       if (hesstype==3)
00523       {
00524         bHess->initialize();
00525       }
00526       else if (hesstype==4)
00527       {
00528         if (sparse_hessian_flag==0)
00529         {
00530           Hess.initialize();
00531         }
00532         else
00533         {
00534           sparse_triplet2->initialize();
00535         }
00536       }
00537       else
00538       {
00539         cerr << "Illegal value for hesstype" << endl;
00540         ad_exit(1);
00541       }
00542       initial_params::straight_through_flag=1;
00543       block_diagonal_flag=3;
00544       local_dtemp.initialize();
00545 
00546       // *****  Note for quadratic prior code: I don't think that this
00547       // part gets added to the Hessian here.
00548       sparse_count=0;
00549       sparse_count_adjoint=0;
00550       pfmin->user_function();
00551 
00552       // *** Hessian calculated just above did not have quadratic prior
00553       // in it so can save this part for quadratci prioer adjoint calculations
00554       if (quadratic_prior::get_num_quadratic_prior()>0)
00555       {
00556         if (pHess_non_quadprior_part)
00557         {
00558           if (pHess_non_quadprior_part->indexmax() != Hess.indexmax())
00559           {
00560             delete pHess_non_quadprior_part;
00561             pHess_non_quadprior_part=0;
00562           }
00563         }
00564         if (!pHess_non_quadprior_part)
00565         {
00566           pHess_non_quadprior_part=new dmatrix(1,usize,1,usize);
00567           if (!pHess_non_quadprior_part)
00568           {
00569             cerr << "Error allocating memory for Hesssian part" << endl;
00570             ad_exit(1);
00571           }
00572         }
00573         (*pHess_non_quadprior_part)=Hess;
00574       }
00575 
00576       block_diagonal_flag=0;
00577       initial_params::straight_through_flag=1;
00578 
00579       //dmatrix tHess=dmatrix(*bHess);
00580       initial_params::straight_through_flag=0;
00581       funnel_init_var::lapprox=0;
00582       //cout << initial_df1b2params::cobjfun << endl;
00583       //f=initial_df1b2params::cobjfun;
00584       block_diagonal_flag=0;
00585 #ifndef OPT_LIB
00586       assert(nvar <= INT_MAX);
00587 #endif
00588       dvector scale1(1,(int)nvar);   // need to get scale from somewhere
00589       initial_params::set_inactive_only_random_effects();
00590       /*int check=*/initial_params::stddev_scale(scale1,x);
00591       //for (i=1;i<=xadjoint.indexmax();i++)
00592       //  xadjoint(i)*=scale1(i);
00593       laplace_approximation_calculator::where_are_we_flag=0;
00594 
00595       if (df1b2quadratic_prior::get_num_quadratic_prior()>0)
00596       {
00597        // !!!! need to fix this!!!!!!!!!!!!!!!!!!!!!!!
00598         laplace_approximation_calculator::where_are_we_flag=3;
00599         quadratic_prior::in_qp_calculations=1;
00600         funnel_init_var::lapprox=this;
00601         df1b2_gradlist::set_no_derivatives();
00602         df1b2quadratic_prior::get_Lxu_contribution(Dux);
00603         quadratic_prior::in_qp_calculations=0;
00604         funnel_init_var::lapprox=0;
00605         laplace_approximation_calculator::where_are_we_flag=0;
00606       }
00607       if (initial_df1b2params::separable_flag)
00608       {
00609         dvector scale(1,(int)nvar);   // need to get scale from somewhere
00610         /*int check=*/initial_params::stddev_scale(scale,x);
00611         dvector sscale=scale(1,Dux(1).indexmax());
00612         for (i=1;i<=usize;i++)
00613         {
00614           Dux(i)=elem_prod(Dux(i),sscale);
00615         }
00616         local_dtemp=elem_prod(local_dtemp,sscale);
00617       }
00618       //cout << trans(Dux)(1) << endl;
00619       //cout << trans(Dux)(3) << endl;
00620       if (quadratic_prior::get_num_quadratic_prior()>0)
00621       {
00622         dvector tmp=evaluate_function_with_quadprior(x,usize,pfmin);
00623         local_dtemp+=tmp;
00624       }
00625 
00626       for (i=1;i<=xsize;i++)
00627       {
00628         xadjoint(i)+=local_dtemp(i);
00629       }
00630       if (df1b2quadratic_prior::get_num_quadratic_prior()>0)
00631       {
00632        // !!!! need to fix this!!!!!!!!!!!!!!!!!!!!!!!
00633         quadratic_prior::get_cHessian_contribution_from_vHessian(Hess,xsize);
00634       }
00635 
00636       if (hesstype==3)
00637       {
00638         //xadjoint -= uadjoint*solve(*bHess,Dux);
00639         if (bHess_pd_flag==0)
00640         {
00641           xadjoint -= solve(*bHess,uadjoint)*Dux;
00642         }
00643       }
00644       else if (hesstype==4)
00645       {
00646         if (sparse_hessian_flag)
00647         {
00648           if (!sparse_triplet2 || !sparse_symbolic2)
00649           {
00650             throw std::bad_alloc();
00651           }
00652           else
00653           {
00654        //dvector tmp=solve(*sparse_triplet,Hess,uadjoint,*sparse_symbolic)*Dux;
00655             dvector tmp=solve(*sparse_triplet2,uadjoint,*sparse_symbolic2)*Dux;
00656             xadjoint -= tmp;
00657           }
00658         }
00659         else
00660         {
00661           xadjoint -= solve(Hess,uadjoint)*Dux;
00662         }
00663       }
00664     }
00665     if (bHess_pd_flag==0) break;
00666   }
00667 
00668   return xadjoint;
00669 }
00670   //int check_pool_flag1=0;
00671 
00676 void laplace_approximation_calculator::
00677   do_separable_stuff_newton_raphson_banded(df1b2variable& ff)
00678 {
00679   //***********************************************************
00680   //***********************************************************
00681   set_dependent_variable(ff);
00682   df1b2_gradlist::set_no_derivatives();
00683   df1b2variable::passnumber=1;
00684   //if (check_pool_flag1)
00685    // check_pool_depths();
00686   df1b2_gradcalc1();
00687   //if (check_pool_flag1)
00688   //  check_pool_depths();
00689 
00690   init_df1b2vector & locy= *funnel_init_var::py;
00691   imatrix& list=*funnel_init_var::plist;
00692   int num_local_re=0;
00693   int num_fixed_effects=0;
00694 
00695 #ifndef OPT_LIB
00696   assert(funnel_init_var::num_active_parameters <= INT_MAX);
00697 #endif
00698   ivector lre_index(1, (int)funnel_init_var::num_active_parameters);
00699   ivector lfe_index(1, (int)funnel_init_var::num_active_parameters);
00700 
00701   for (int i=1;i<=(int)funnel_init_var::num_active_parameters;i++)
00702   {
00703     if (list(i,1)>xsize)
00704     {
00705       lre_index(++num_local_re)=i;
00706     }
00707     else if (list(i,1)>0)
00708     {
00709       lfe_index(++num_fixed_effects)=i;
00710     }
00711   }
00712 
00713   if (num_local_re > 0)
00714   {
00715     switch(hesstype)
00716     {
00717     case 3:
00718       for (int i=1;i<=num_local_re;i++)
00719       {
00720         int lrei=lre_index(i);
00721         for (int j=1;j<=num_local_re;j++)
00722         {
00723           int lrej=lre_index(j);
00724           int i1=list(lrei,1)-xsize;
00725           int i2=list(lrei,2);
00726           int j1=list(lrej,1)-xsize;
00727           int j2=list(lrej,2);
00728           if (i1>=j1) (*bHess)(i1,j1)+=locy(i2).u_bar[j2-1];
00729         }
00730       }
00731       break;
00732     case 4:
00733       if (sparse_hessian_flag==0)
00734       {
00735         for (int i=1;i<=num_local_re;i++)
00736         {
00737           int lrei=lre_index(i);
00738           for (int j=1;j<=num_local_re;j++)
00739           {
00740             int lrej=lre_index(j);
00741             int i1=list(lrei,1)-xsize;
00742             int i2=list(lrei,2);
00743             int j1=list(lrej,1)-xsize;
00744             int j2=list(lrej,2);
00745             Hess(i1,j1)+=locy(i2).u_bar[j2-1];
00746           }
00747         }
00748       }
00749       else
00750       {
00751         for (int i=1;i<=num_local_re;i++)
00752         {
00753           int lrei=lre_index(i);
00754           for (int j=1;j<=num_local_re;j++)
00755           {
00756             int lrej=lre_index(j);
00757             int i1=list(lrei,1)-xsize;
00758             int i2=list(lrei,2);
00759             int j1=list(lrej,1)-xsize;
00760             int j2=list(lrej,2);
00761 
00762             if (i1<=j1)
00763             {
00764               sparse_count++;
00765               (*sparse_triplet2)((*sparse_iterator)(sparse_count))
00766                 +=locy(i2).u_bar[j2-1];
00767             }
00768           }
00769         }
00770       }
00771       break;
00772     default:
00773       cerr << "illegal value for hesstype" << endl;
00774       ad_exit(1);
00775     }
00776 
00777     for (int i=1;i<=num_local_re;i++)
00778     {
00779       int lrei=lre_index(i);
00780       int i1=list(lrei,1);
00781       int i2=list(lrei,2);
00782       //grad(i1-xsize)= re_objective_function_value::pobjfun->u_dot[i2-1];
00783       grad(i1-xsize)+=ff.u_dot[i2-1];
00784     }
00785   }
00786 
00787   f1b2gradlist->reset();
00788   f1b2gradlist->list.initialize();
00789   f1b2gradlist->list2.initialize();
00790   f1b2gradlist->list3.initialize();
00791   f1b2gradlist->nlist.initialize();
00792   f1b2gradlist->nlist2.initialize();
00793   f1b2gradlist->nlist3.initialize();
00794   funnel_init_var::num_vars=0;
00795   funnel_init_var::num_active_parameters=0;
00796   funnel_init_var::num_inactive_vars=0;
00797 }
00798 //int tmp_testcount=0;
00799 df1b2variable * tmp_pen=00;
00800 
00805 dvector laplace_approximation_calculator::
00806   get_newton_raphson_info_banded (function_minimizer * pfmin)
00807 {
00808   int nv=initial_df1b2params::set_index();
00809   if (allocated(used_flags))
00810   {
00811     if (used_flags.indexmax() != nv)
00812     {
00813       used_flags.safe_deallocate();
00814     }
00815   }
00816   if (!allocated(used_flags))
00817   {
00818     used_flags.safe_allocate(1,nv);
00819   }
00820 
00821   for (int ip=1;ip<=num_der_blocks;ip++)
00822   {
00823     if (ip>1)   // change to combine sparse matrix stuff with num der blocks
00824     {           // df  3-4-09
00825       sparse_count=0;
00826     }
00827     used_flags.initialize();
00828     // do we need to reallocate memory for df1b2variables?
00829     check_for_need_to_reallocate(ip);
00830     df1b2_gradlist::set_no_derivatives();
00831     //cout << re_objective_function_value::pobjfun << endl;
00832     //cout << re_objective_function_value::pobjfun->ptr << endl;
00833     (*re_objective_function_value::pobjfun)=0;
00834     df1b2variable pen=0.0;
00835     tmp_pen=&pen;
00836     df1b2variable zz=0.0;
00837 
00838     initial_df1b2params::reset(y,pen);
00839 
00840     // call function to do block diagonal newton-raphson
00841     // the step vector from the newton-raphson is in the vector step
00842     df1b2_gradlist::set_yes_derivatives();
00843 
00844     funnel_init_var::lapprox=this;
00845     //cout << funnel_init_var::lapprox << endl;
00846     block_diagonal_flag=1;
00847    /*
00848     tmp_testcount++;
00849     if (tmp_testcount>=9)
00850     {
00851       pen.deallocate();
00852     }
00853     */
00854 
00855     if (ip==1)
00856     {
00857       if (sparse_triplet2)
00858         sparse_triplet2->initialize();
00859     }
00860 
00861     pfmin->user_function();
00862     /*
00863     if (tmp_testcount>=9)
00864     {
00865       pen.deallocate();
00866     }
00867      */
00868     funnel_init_var::lapprox=0;
00869     block_diagonal_flag=0;
00870     pen.deallocate();
00871   }
00872 
00873   return step;
00874 }
00875 
00880 void laplace_approximation_calculator::
00881   do_separable_stuff_laplace_approximation_banded(df1b2variable& ff)
00882 {
00883   set_dependent_variable(ff);
00884   //df1b2_gradlist::set_no_derivatives();
00885   df1b2variable::passnumber=1;
00886   df1b2_gradcalc1();
00887 
00888   init_df1b2vector & locy= *funnel_init_var::py;
00889   imatrix& list=*funnel_init_var::plist;
00890 
00891   int us=0; int xs=0;
00892 #ifndef OPT_LIB
00893   assert(funnel_init_var::num_active_parameters <= INT_MAX);
00894 #endif
00895   ivector lre_index(1,(int)funnel_init_var::num_active_parameters);
00896   ivector lfe_index(1,(int)funnel_init_var::num_active_parameters);
00897 
00898   for (int i=1;i<=(int)funnel_init_var::num_active_parameters;i++)
00899   {
00900     if (list(i,1)>xsize)
00901     {
00902       lre_index(++us)=i;
00903     }
00904     else if (list(i,1)>0)
00905     {
00906       lfe_index(++xs)=i;
00907     }
00908   }
00909 
00910   for (int j=1;j<=xs;j++)
00911   {
00912     int j1=list(lfe_index(j),1);
00913     int j2=list(lfe_index(j),2);
00914     xadjoint(j1)+=ff.u_dot[j2-1];
00915   }
00916 
00917   if (us>0)
00918   {
00919     if (hesstype==3)
00920     {
00921       for (int i=1;i<=us;i++)
00922       {
00923         for (int j=1;j<=us;j++)
00924         {
00925           int i1=list(lre_index(i),1)-xsize;
00926           int i2=list(lre_index(i),2);
00927           int j1=list(lre_index(j),1)-xsize;
00928           int j2=list(lre_index(j),2);
00929           if (i1>=j1) (*bHess)(i1,j1)+=locy(i2).u_bar[j2-1];
00930         }
00931       }
00932     }
00933     else if (hesstype==4)
00934     {
00935       if (sparse_hessian_flag==0)
00936       {
00937         for (int i=1;i<=us;i++)
00938         {
00939           for (int j=1;j<=us;j++)
00940           {
00941             int i1=list(lre_index(i),1)-xsize;
00942             int i2=list(lre_index(i),2);
00943             int j1=list(lre_index(j),1)-xsize;
00944             int j2=list(lre_index(j),2);
00945             Hess(i1,j1)+=locy(i2).u_bar[j2-1];
00946           }
00947         }
00948       }
00949       else
00950       {
00951         for (int i=1;i<=us;i++)
00952         {
00953           for (int j=1;j<=us;j++)
00954           {
00955             int i1=list(lre_index(i),1)-xsize;
00956             int i2=list(lre_index(i),2);
00957             int j1=list(lre_index(j),1)-xsize;
00958             int j2=list(lre_index(j),2);
00959 
00960             if (i1<=j1)
00961             {
00962               sparse_count++;
00963               (*sparse_triplet2)((*sparse_iterator)(sparse_count))
00964                 +=locy(i2).u_bar[j2-1];
00965             }
00966           }
00967         }
00968       }
00969     }
00970 
00971     for (int i=1;i<=us;i++)
00972     {
00973       int i1=list(lre_index(i),1)-xsize;
00974       int i2=list(lre_index(i),2);
00975       uadjoint(i1)+=ff.u_dot[i2-1];
00976     }
00977 
00978     for (int i=1;i<=us;i++)
00979     {
00980       for (int j=1;j<=xs;j++)
00981       {
00982         int i1=list(lre_index(i),1)-xsize;
00983         int i2=list(lre_index(i),2);
00984         int j1=list(lfe_index(j),1);
00985         int j2=list(lfe_index(j),2);
00986         Dux(i1,j1)+=locy(i2).u_bar[j2-1];
00987       }
00988     }
00989   }
00990   f1b2gradlist->reset();
00991   f1b2gradlist->list.initialize();
00992   f1b2gradlist->list2.initialize();
00993   f1b2gradlist->list3.initialize();
00994   f1b2gradlist->nlist.initialize();
00995   f1b2gradlist->nlist2.initialize();
00996   f1b2gradlist->nlist3.initialize();
00997   funnel_init_var::num_vars=0;
00998   funnel_init_var::num_active_parameters=0;
00999   funnel_init_var::num_inactive_vars=0;
01000 }
01001 
01006 double calculate_laplace_approximation(const dvector& x,const dvector& u0,
01007   const banded_symmetric_dmatrix& bHess,const dvector& _xadjoint,
01008   const dvector& _uadjoint,
01009   const banded_symmetric_dmatrix& _bHessadjoint,function_minimizer * pmin)
01010 {
01011   ADUNCONST(dvector,xadjoint)
01012   ADUNCONST(dvector,uadjoint)
01013   ADUNCONST(banded_symmetric_dmatrix,bHessadjoint)
01014   int bw=bHess.bandwidth();
01015   const int xs=x.size();
01016   const int us=u0.size();
01017   gradient_structure::set_YES_DERIVATIVES();
01018   int nvar=x.size()+u0.size()+((bw+1)*(2*u0.size()-bw))/2;
01019   independent_variables y(1,nvar);
01020 
01021   // need to set random effects active together with whatever
01022   // init parameters should be active in this phase
01023   initial_params::set_inactive_only_random_effects();
01024   initial_params::set_active_random_effects();
01025   /*int onvar=*/initial_params::nvarcalc();
01026   initial_params::xinit(y);    // get the initial values into the
01027   y(1,xs)=x;
01028 
01029  // Here need hooks for sparse matrix structures
01030   int ii=xs+us+1;
01031   for (int i=1;i<=us;i++)
01032   {
01033     int jmin=admax(1,i-bw+1);
01034     for (int j=jmin;j<=i;j++)
01035       y(ii++)=bHess(i,j);
01036   }
01037 
01038   dvar_vector vy=dvar_vector(y);
01039   banded_symmetric_dvar_matrix vHess(1,us,bw);
01040   initial_params::reset(vy);    // get the initial values into the
01041   ii=xs+us+1;
01042   for (int i=1;i<=us;i++)
01043   {
01044     int jmin=admax(1,i-bw+1);
01045     for (int j=jmin;j<=i;j++)
01046       vHess(i,j)=vy(ii++);
01047   }
01048 
01049    dvariable vf=0.0;
01050 
01051    *objective_function_value::pobjfun=0.0;
01052    pmin->AD_uf_outer();
01053    vf+=*objective_function_value::pobjfun;
01054 
01055    int sgn=0;
01056    dvariable ld;
01057 
01058    int eigswitch=0;
01059    if (eigswitch)
01060    {
01061      ofstream ofs("ee");
01062      dvector ev(bHess.indexmin(),bHess.indexmax());
01063      dmatrix evecs=eigenvectors(dmatrix(bHess),ev);
01064      ofs << setprecision(3) << setw(12) << setscientific() << dmatrix(bHess)
01065          << endl << endl;
01066      ofs << ev << endl << endl << evecs << endl;
01067    }
01068    ld=0.5*ln_det_choleski(vHess,sgn);
01069    if (sgn==1)
01070    {
01071      cout << "Choleski failed" << endl;
01072      pmin->lapprox->bHess_pd_flag=1;
01073    }
01074 
01075    vf+=ld;
01076    const double ltp=0.5*log(2.0*PI);
01077    vf-=us*ltp;
01078    double f=value(vf);
01079    dvector g(1,nvar);
01080    gradcalc(nvar,g);
01081 
01082   ii=1;
01083   for (int i=1;i<=xs;i++)
01084     xadjoint(i)=g(ii++);
01085   for (int i=1;i<=us;i++)
01086     uadjoint(i)=g(ii++);
01087   for (int i=1;i<=us;i++)
01088   {
01089     int jmin=admax(1,i-bw+1);
01090     for (int j=jmin;j<=i;j++)
01091       bHessadjoint(i,j)=g(ii++);
01092   }
01093   return f;
01094 }
01095 
01100 dvector laplace_approximation_calculator::
01101   banded_calculations_trust_region_approach(const dvector& _uhat,
01102   function_minimizer * pfmin)
01103 {
01104   dvector& uhat=(dvector&) _uhat;
01105   dvector uhat_old(uhat.indexmin(),uhat.indexmax());
01106   dvector uhat_new(uhat.indexmin(),uhat.indexmax());
01107   dvector uhat_best(uhat.indexmin(),uhat.indexmax());
01108 
01109   double wght=0.0;
01110   double delta=5.e-5;
01111   //do
01112   dvector values(1,300);
01113   double oldfbest=pmin->lapprox->fmc1.fbest;
01114   double newfbest = 0.0;
01115   int have_value=0;
01116   //for (int jj=1;jj<=300;jj++)
01117   int jj=1;
01118   double lastval=oldfbest;
01119   do
01120   {
01121     jj++;
01122     wght+=delta;
01123     //double wght=0.0;
01124     double newval=0.0;
01125     //cout << "Enter weight size " << endl;
01126     //cin >> wght;
01127     if (wght<0.0)
01128       break;
01129     int mmin=bHess->indexmin();
01130     int mmax=bHess->indexmax();
01131     banded_symmetric_dmatrix tmp(mmin,mmax,bHess->bandwidth());
01132     tmp=*bHess;
01133     uhat_old=uhat;
01134     int ierr=0;
01135     for (int i=mmin;i<=mmax;i++)
01136     {
01137       tmp(i,i)+=wght;
01138     }
01139     banded_lower_triangular_dmatrix bltd=choleski_decomp(tmp,ierr);
01140     if (!ierr)
01141     {
01142       dvector v=solve(bltd,grad);
01143       step=-solve_trans(bltd,v);
01144 
01145       uhat_old=uhat;
01146       uhat+=step;
01147       //cout << "norm(uhat_old) = " << norm(uhat_old)
01148        //    << "   norm(uhat) = " << norm(uhat)  << endl;
01149 
01150       /*double maxg=*/evaluate_function(newval,uhat,pfmin);
01151       if (have_value && newval>newfbest)
01152       {
01153         break;
01154       }
01155       if (jj>1)
01156       {
01157         if (newval<lastval)  // we are doing better so increasse step size
01158         {
01159           delta*=2;
01160         }
01161         if (newval>lastval && !have_value)  // we have gone to far go back
01162         {
01163           wght-=delta;
01164           delta/=16;
01165         }
01166       }
01167       lastval=newval;
01168 
01169       if (newval<newfbest)
01170       {
01171         newfbest=newval;
01172         uhat_best=uhat;
01173         have_value=jj;
01174       }
01175       uhat_new=uhat;
01176       uhat=uhat_old;
01177     }
01178     else
01179     {
01180       delta*=2;
01181     }
01182   }
01183   while(jj<10);
01184   if (!have_value)
01185   {
01186     cerr << "can't improve function value in trust region calculations"
01187          << endl;
01188     //ad_exit(1);
01189   }
01190   return uhat_best;
01213 }