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