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