ADMB Documentation  11.1.2274
 All Classes Files Functions Variables Typedefs Friends Defines
df1b2lp6.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: df1b2lp6.cpp 2258 2014-08-27 22:37:49Z 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   do
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   while(1);
00652 
00653   return xadjoint;
00654 }
00655   //int check_pool_flag1=0;
00656 
00661 void laplace_approximation_calculator::
00662   do_separable_stuff_newton_raphson_banded(df1b2variable& ff)
00663 {
00664   //***********************************************************
00665   //***********************************************************
00666   set_dependent_variable(ff);
00667   df1b2_gradlist::set_no_derivatives();
00668   df1b2variable::passnumber=1;
00669   //if (check_pool_flag1)
00670    // check_pool_depths();
00671   df1b2_gradcalc1();
00672   //if (check_pool_flag1)
00673   //  check_pool_depths();
00674 
00675   init_df1b2vector & locy= *funnel_init_var::py;
00676   imatrix& list=*funnel_init_var::plist;
00677   int num_local_re=0;
00678   int num_fixed_effects=0;
00679 
00680   int i;
00681   ivector lre_index(1,funnel_init_var::num_active_parameters);
00682   ivector lfe_index(1,funnel_init_var::num_active_parameters);
00683 
00684   for (i=1;i<=funnel_init_var::num_active_parameters;i++)
00685   {
00686     if (list(i,1)>xsize)
00687     {
00688       lre_index(++num_local_re)=i;
00689     }
00690     else if (list(i,1)>0)
00691     {
00692       lfe_index(++num_fixed_effects)=i;
00693     }
00694   }
00695 
00696   if (num_local_re > 0)
00697   {
00698     int j;
00699     switch(hesstype)
00700     {
00701     case 3:
00702       for (i=1;i<=num_local_re;i++)
00703       {
00704         int lrei=lre_index(i);
00705         for (j=1;j<=num_local_re;j++)
00706         {
00707           int lrej=lre_index(j);
00708           int i1=list(lrei,1)-xsize;
00709           int i2=list(lrei,2);
00710           int j1=list(lrej,1)-xsize;
00711           int j2=list(lrej,2);
00712           if (i1>=j1) (*bHess)(i1,j1)+=locy(i2).u_bar[j2-1];
00713         }
00714       }
00715       break;
00716     case 4:
00717       if (sparse_hessian_flag==0)
00718       {
00719         for (i=1;i<=num_local_re;i++)
00720         {
00721           int lrei=lre_index(i);
00722           for (j=1;j<=num_local_re;j++)
00723           {
00724             int lrej=lre_index(j);
00725             int i1=list(lrei,1)-xsize;
00726             int i2=list(lrei,2);
00727             int j1=list(lrej,1)-xsize;
00728             int j2=list(lrej,2);
00729             Hess(i1,j1)+=locy(i2).u_bar[j2-1];
00730           }
00731         }
00732       }
00733       else
00734       {
00735         for (i=1;i<=num_local_re;i++)
00736         {
00737           int lrei=lre_index(i);
00738           for (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 
00746             if (i1<=j1)
00747             {
00748               sparse_count++;
00749               (*sparse_triplet2)((*sparse_iterator)(sparse_count))
00750                 +=locy(i2).u_bar[j2-1];
00751             }
00752           }
00753         }
00754       }
00755       break;
00756     default:
00757       cerr << "illegal value for hesstype" << endl;
00758       ad_exit(1);
00759     }
00760 
00761     for (i=1;i<=num_local_re;i++)
00762     {
00763       int lrei=lre_index(i);
00764       int i1=list(lrei,1);
00765       int i2=list(lrei,2);
00766       //grad(i1-xsize)= re_objective_function_value::pobjfun->u_dot[i2-1];
00767       grad(i1-xsize)+=ff.u_dot[i2-1];
00768     }
00769   }
00770 
00771   f1b2gradlist->reset();
00772   f1b2gradlist->list.initialize();
00773   f1b2gradlist->list2.initialize();
00774   f1b2gradlist->list3.initialize();
00775   f1b2gradlist->nlist.initialize();
00776   f1b2gradlist->nlist2.initialize();
00777   f1b2gradlist->nlist3.initialize();
00778   funnel_init_var::num_vars=0;
00779   funnel_init_var::num_active_parameters=0;
00780   funnel_init_var::num_inactive_vars=0;
00781 }
00782 //int tmp_testcount=0;
00783 df1b2variable * tmp_pen=00;
00784 
00789 dvector laplace_approximation_calculator::
00790   get_newton_raphson_info_banded (function_minimizer * pfmin)
00791 {
00792   int nv=initial_df1b2params::set_index();
00793   if (allocated(used_flags))
00794   {
00795     if (used_flags.indexmax() != nv)
00796     {
00797       used_flags.safe_deallocate();
00798     }
00799   }
00800   if (!allocated(used_flags))
00801   {
00802     used_flags.safe_allocate(1,nv);
00803   }
00804 
00805   for (int ip=1;ip<=num_der_blocks;ip++)
00806   {
00807     if (ip>1)   // change to combine sparse matrix stuff with num der blocks
00808     {           // df  3-4-09
00809       sparse_count=0;
00810     }
00811     used_flags.initialize();
00812     // do we need to reallocate memory for df1b2variables?
00813     check_for_need_to_reallocate(ip);
00814     df1b2_gradlist::set_no_derivatives();
00815     //cout << re_objective_function_value::pobjfun << endl;
00816     //cout << re_objective_function_value::pobjfun->ptr << endl;
00817     (*re_objective_function_value::pobjfun)=0;
00818     df1b2variable pen=0.0;
00819     tmp_pen=&pen;
00820     df1b2variable zz=0.0;
00821 
00822     initial_df1b2params::reset(y,pen);
00823 
00824     // call function to do block diagonal newton-raphson
00825     // the step vector from the newton-raphson is in the vector step
00826     df1b2_gradlist::set_yes_derivatives();
00827 
00828     funnel_init_var::lapprox=this;
00829     //cout << funnel_init_var::lapprox << endl;
00830     block_diagonal_flag=1;
00831    /*
00832     tmp_testcount++;
00833     if (tmp_testcount>=9)
00834     {
00835       pen.deallocate();
00836     }
00837     */
00838 
00839     if (ip==1)
00840     {
00841       if (sparse_triplet2)
00842         sparse_triplet2->initialize();
00843     }
00844 
00845     pfmin->user_function();
00846     /*
00847     if (tmp_testcount>=9)
00848     {
00849       pen.deallocate();
00850     }
00851      */
00852     funnel_init_var::lapprox=0;
00853     block_diagonal_flag=0;
00854     pen.deallocate();
00855   }
00856 
00857   return step;
00858 }
00859 
00864 void laplace_approximation_calculator::
00865   do_separable_stuff_laplace_approximation_banded(df1b2variable& ff)
00866 {
00867   set_dependent_variable(ff);
00868   //df1b2_gradlist::set_no_derivatives();
00869   df1b2variable::passnumber=1;
00870   df1b2_gradcalc1();
00871 
00872   init_df1b2vector & locy= *funnel_init_var::py;
00873   imatrix& list=*funnel_init_var::plist;
00874 
00875   int i; int j; int us=0; int xs=0;
00876   ivector lre_index(1,funnel_init_var::num_active_parameters);
00877   ivector lfe_index(1,funnel_init_var::num_active_parameters);
00878 
00879   for (i=1;i<=funnel_init_var::num_active_parameters;i++)
00880   {
00881     if (list(i,1)>xsize)
00882     {
00883       lre_index(++us)=i;
00884     }
00885     else if (list(i,1)>0)
00886     {
00887       lfe_index(++xs)=i;
00888     }
00889   }
00890 
00891   for (j=1;j<=xs;j++)
00892   {
00893     int j1=list(lfe_index(j),1);
00894     int j2=list(lfe_index(j),2);
00895     xadjoint(j1)+=ff.u_dot[j2-1];
00896   }
00897 
00898   if (us>0)
00899   {
00900     if (hesstype==3)
00901     {
00902       for (i=1;i<=us;i++)
00903       {
00904         for (j=1;j<=us;j++)
00905         {
00906           int i1=list(lre_index(i),1)-xsize;
00907           int i2=list(lre_index(i),2);
00908           int j1=list(lre_index(j),1)-xsize;
00909           int j2=list(lre_index(j),2);
00910           if (i1>=j1) (*bHess)(i1,j1)+=locy(i2).u_bar[j2-1];
00911         }
00912       }
00913     }
00914     else if (hesstype==4)
00915     {
00916       if (sparse_hessian_flag==0)
00917       {
00918         for (i=1;i<=us;i++)
00919         {
00920           for (j=1;j<=us;j++)
00921           {
00922             int i1=list(lre_index(i),1)-xsize;
00923             int i2=list(lre_index(i),2);
00924             int j1=list(lre_index(j),1)-xsize;
00925             int j2=list(lre_index(j),2);
00926             Hess(i1,j1)+=locy(i2).u_bar[j2-1];
00927           }
00928         }
00929       }
00930       else
00931       {
00932         for (i=1;i<=us;i++)
00933         {
00934           for (j=1;j<=us;j++)
00935           {
00936             int i1=list(lre_index(i),1)-xsize;
00937             int i2=list(lre_index(i),2);
00938             int j1=list(lre_index(j),1)-xsize;
00939             int j2=list(lre_index(j),2);
00940 
00941             if (i1<=j1)
00942             {
00943               sparse_count++;
00944               (*sparse_triplet2)((*sparse_iterator)(sparse_count))
00945                 +=locy(i2).u_bar[j2-1];
00946             }
00947           }
00948         }
00949       }
00950     }
00951 
00952     for (i=1;i<=us;i++)
00953     {
00954       int i1=list(lre_index(i),1)-xsize;
00955       int i2=list(lre_index(i),2);
00956       uadjoint(i1)+=ff.u_dot[i2-1];
00957     }
00958 
00959     for (i=1;i<=us;i++)
00960     {
00961       for (j=1;j<=xs;j++)
00962       {
00963         int i1=list(lre_index(i),1)-xsize;
00964         int i2=list(lre_index(i),2);
00965         int j1=list(lfe_index(j),1);
00966         int j2=list(lfe_index(j),2);
00967         Dux(i1,j1)+=locy(i2).u_bar[j2-1];
00968       }
00969     }
00970   }
00971   f1b2gradlist->reset();
00972   f1b2gradlist->list.initialize();
00973   f1b2gradlist->list2.initialize();
00974   f1b2gradlist->list3.initialize();
00975   f1b2gradlist->nlist.initialize();
00976   f1b2gradlist->nlist2.initialize();
00977   f1b2gradlist->nlist3.initialize();
00978   funnel_init_var::num_vars=0;
00979   funnel_init_var::num_active_parameters=0;
00980   funnel_init_var::num_inactive_vars=0;
00981 }
00982 
00987 double calculate_laplace_approximation(const dvector& x,const dvector& u0,
00988   const banded_symmetric_dmatrix& bHess,const dvector& _xadjoint,
00989   const dvector& _uadjoint,
00990   const banded_symmetric_dmatrix& _bHessadjoint,function_minimizer * pmin)
00991 {
00992   ADUNCONST(dvector,xadjoint)
00993   ADUNCONST(dvector,uadjoint)
00994   ADUNCONST(banded_symmetric_dmatrix,bHessadjoint)
00995   int bw=bHess.bandwidth();
00996   const int xs=x.size();
00997   const int us=u0.size();
00998   gradient_structure::set_YES_DERIVATIVES();
00999   int nvar=x.size()+u0.size()+((bw+1)*(2*u0.size()-bw))/2;
01000   independent_variables y(1,nvar);
01001 
01002   // need to set random effects active together with whatever
01003   // init parameters should be active in this phase
01004   initial_params::set_inactive_only_random_effects();
01005   initial_params::set_active_random_effects();
01006   /*int onvar=*/initial_params::nvarcalc();
01007   initial_params::xinit(y);    // get the initial values into the
01008   y(1,xs)=x;
01009 
01010   int i,j;
01011 
01012  // Here need hooks for sparse matrix structures
01013   int ii=xs+us+1;
01014   for (i=1;i<=us;i++)
01015   {
01016     int jmin=admax(1,i-bw+1);
01017     for (j=jmin;j<=i;j++)
01018     y(ii++)=bHess(i,j);
01019   }
01020 
01021   dvar_vector vy=dvar_vector(y);
01022   banded_symmetric_dvar_matrix vHess(1,us,bw);
01023   initial_params::reset(vy);    // get the initial values into the
01024   ii=xs+us+1;
01025   for (i=1;i<=us;i++)
01026   {
01027     int jmin=admax(1,i-bw+1);
01028     for (j=jmin;j<=i;j++)
01029       vHess(i,j)=vy(ii++);
01030   }
01031 
01032    dvariable vf=0.0;
01033 
01034    *objective_function_value::pobjfun=0.0;
01035    pmin->AD_uf_outer();
01036    vf+=*objective_function_value::pobjfun;
01037 
01038    int sgn=0;
01039    dvariable ld;
01040 
01041    int eigswitch=0;
01042    if (eigswitch)
01043    {
01044      ofstream ofs("ee");
01045      dvector ev(bHess.indexmin(),bHess.indexmax());
01046      dmatrix evecs=eigenvectors(dmatrix(bHess),ev);
01047      ofs << setprecision(3) << setw(12) << setscientific() << dmatrix(bHess)
01048          << endl << endl;
01049      ofs << ev << endl << endl << evecs << endl;
01050    }
01051    ld=0.5*ln_det_choleski(vHess,sgn);
01052    if (sgn==1)
01053    {
01054      cout << "Choleski failed" << endl;
01055      pmin->lapprox->bHess_pd_flag=1;
01056    }
01057 
01058    vf+=ld;
01059    const double ltp=0.5*log(2.0*PI);
01060    vf-=us*ltp;
01061    double f=value(vf);
01062    dvector g(1,nvar);
01063    gradcalc(nvar,g);
01064 
01065   ii=1;
01066   for (i=1;i<=xs;i++)
01067     xadjoint(i)=g(ii++);
01068   for (i=1;i<=us;i++)
01069     uadjoint(i)=g(ii++);
01070   for (i=1;i<=us;i++)
01071   {
01072     int jmin=admax(1,i-bw+1);
01073     for (j=jmin;j<=i;j++)
01074       bHessadjoint(i,j)=g(ii++);
01075   }
01076   return f;
01077 }
01078 
01083 dvector laplace_approximation_calculator::
01084   banded_calculations_trust_region_approach(const dvector& _uhat,
01085   function_minimizer * pfmin)
01086 {
01087   dvector& uhat=(dvector&) _uhat;
01088   dvector uhat_old(uhat.indexmin(),uhat.indexmax());
01089   dvector uhat_new(uhat.indexmin(),uhat.indexmax());
01090   dvector uhat_best(uhat.indexmin(),uhat.indexmax());
01091 
01092   double wght=0.0;
01093   double delta=5.e-5;
01094   //do
01095   dvector values(1,300);
01096   double oldfbest=pmin->lapprox->fmc1.fbest;
01097   double newfbest = 0.0;
01098   int have_value=0;
01099   //for (int jj=1;jj<=300;jj++)
01100   int jj=1;
01101   double lastval=oldfbest;
01102   do
01103   {
01104     jj++;
01105     wght+=delta;
01106     //double wght=0.0;
01107     double newval=0.0;
01108     //cout << "Enter weight size " << endl;
01109     //cin >> wght;
01110     if (wght<0.0)
01111       break;
01112     int mmin=bHess->indexmin();
01113     int mmax=bHess->indexmax();
01114     banded_symmetric_dmatrix tmp(mmin,mmax,bHess->bandwidth());
01115     tmp=*bHess;
01116     uhat_old=uhat;
01117     int ierr=0;
01118     for (int i=mmin;i<=mmax;i++)
01119     {
01120       tmp(i,i)+=wght;
01121     }
01122     banded_lower_triangular_dmatrix bltd=choleski_decomp(tmp,ierr);
01123     if (!ierr)
01124     {
01125       dvector v=solve(bltd,grad);
01126       step=-solve_trans(bltd,v);
01127 
01128       uhat_old=uhat;
01129       uhat+=step;
01130       //cout << "norm(uhat_old) = " << norm(uhat_old)
01131        //    << "   norm(uhat) = " << norm(uhat)  << endl;
01132 
01133       /*double maxg=*/evaluate_function(newval,uhat,pfmin);
01134       if (have_value && newval>newfbest)
01135       {
01136         break;
01137       }
01138       if (jj>1)
01139       {
01140         if (newval<lastval)  // we are doing better so increasse step size
01141         {
01142           delta*=2;
01143         }
01144         if (newval>lastval && !have_value)  // we have gone to far go back
01145         {
01146           wght-=delta;
01147           delta/=16;
01148         }
01149       }
01150       lastval=newval;
01151 
01152       if (newval<newfbest)
01153       {
01154         newfbest=newval;
01155         uhat_best=uhat;
01156         have_value=jj;
01157       }
01158       uhat_new=uhat;
01159       uhat=uhat_old;
01160     }
01161     else
01162     {
01163       delta*=2;
01164     }
01165   }
01166   while(jj<10);
01167   if (!have_value)
01168   {
01169     cerr << "can't improve function value in trust region calculations"
01170          << endl;
01171     //ad_exit(1);
01172   }
01173   return uhat_best;
01174   initial_params::set_active_only_random_effects();
01175   if (!inner_lmnflag)
01176   {
01177     if (!ADqd_flag)
01178     {
01179       uhat=get_uhat_quasi_newton(uhat_new,pfmin);
01180       //double maxg=fabs(fmc1.gmax);
01181       //double f_from_1=fmc1.fbest;
01182     }
01183     else
01184     {
01185       uhat=get_uhat_quasi_newton_qd(uhat_new,pfmin);
01186     }
01187   }
01188   else
01189   {
01190     uhat=get_uhat_lm_newton(uhat_new,pfmin);
01191   }
01192   return uhat;
01193 }