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