ADMB Documentation  11.1.2501
 All Classes Files Functions Variables Typedefs Friends Defines
mod_rhes.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: mod_rhes.cpp 1935 2014-04-26 02:02:58Z johnoel $
00003  *
00004  * Author: David Fournier
00005  * Copyright (c) 2008-2012 Regents of the University of California
00006  */
00011 #include <sstream>
00012 using std::istringstream;
00013 
00014 #include <admodel.h>
00015 #include <df1b2fun.h>
00016 #include <adrndeff.h>
00017 
00018 void get_inverse_sparse_hessian(dcompressed_triplet & st, hs_symbolic& S,
00019   uostream& ofs1,ofstream& ofs,int usize,int xsize,dvector& u);
00020 
00025 banded_lower_triangular_dmatrix quiet_choleski_decomp(
00026   const banded_symmetric_dmatrix& _M,const int& _ierr)
00027 {
00028   int & ierr = (int &) _ierr;
00029   ADUNCONST(banded_symmetric_dmatrix,M)
00030   int minsave=M.indexmin();
00031   M.shift(1);
00032   int n=M.indexmax();
00033 
00034   int bw=M.bandwidth();
00035   banded_lower_triangular_dmatrix L(1,n,bw);
00036 #ifndef SAFE_INITIALIZE
00037     L.initialize();
00038 #endif
00039 
00040   int i,j,k;
00041   double tmp;
00042     if (M(1,1)<=0)
00043     {
00044       ierr=1;
00045       return L;
00046     }
00047   L(1,1)=sqrt(M(1,1));
00048   for (i=2;i<=bw;i++)
00049   {
00050     L(i,1)=M(i,1)/L(1,1);
00051   }
00052 
00053   for (i=2;i<=n;i++)
00054   {
00055     for (j=i-bw+1;j<=i-1;j++)
00056     {
00057       if (j>1)
00058       {
00059         tmp=M(i,j);
00060         for (k=i-bw+1;k<=j-1;k++)
00061         {
00062           if (k>0 && k>j-bw)
00063             tmp-=L(i,k)*L(j,k);
00064         }
00065         L(i,j)=tmp/L(j,j);
00066       }
00067     }
00068     tmp=M(i,i);
00069     for (k=i-bw+1;k<=i-1;k++)
00070     {
00071       if (k>0)
00072         tmp-=L(i,k)*L(i,k);
00073     }
00074     if (tmp<=0)
00075     {
00076       ierr=1;
00077       return L;
00078     }
00079     L(i,i)=sqrt(tmp);
00080   }
00081   M.shift(minsave);
00082   L.shift(minsave);
00083 
00084   return L;
00085 }
00086 
00091 void function_minimizer::hess_routine_random_effects(void)
00092 {
00093 #if defined(USE_ADPVM)
00094   if (ad_comm::pvm_manager)
00095   {
00096     switch (ad_comm::pvm_manager->mode)
00097     {
00098     case 1: //master
00099       hess_routine_noparallel_random_effects();
00100       break;
00101     case 2: //slave
00102       hess_routine_slave_random_effects();
00103       break;
00104     default:
00105       cerr << "Illegal value for mode" << endl;
00106       ad_exit(1);
00107     }
00108   }
00109   else
00110 #endif
00111   {
00112       hess_routine_noparallel_random_effects();
00113   }
00114 }
00115 dvector get_solution_vector(int npts);
00116 
00121 void function_minimizer::hess_routine_noparallel_random_effects(void)
00122 {
00123   int nvar=initial_params::nvarcalc(); // get the number of active parameters
00124   //if (adjm_ptr) set_labels_for_hess(nvar);
00125   independent_variables x(1,nvar);
00126   initial_params::xinit(x);        // get the initial values into the x vector
00127   double f;
00128   double delta=1.e-4;
00129   dvector g1(1,nvar);
00130   dvector g0(1,nvar);
00131   dvector g2(1,nvar);
00132   dvector gbest(1,nvar);
00133   dvector hess(1,nvar);
00134   dvector hess1(1,nvar);
00135   dvector hess2(1,nvar);
00136   //double eps=.1;
00137   gradient_structure::set_YES_DERIVATIVES();
00138   gbest.fill_seqadd(1.e+50,0.);
00139 
00140   dvector ddd(1,nvar);
00141   gradcalc(0,ddd);
00142   adstring tmpstring;
00143 
00144   {
00145     first_hessian_flag=1;
00146     {
00147       g1=(*lapprox)(x,f,this);
00148       g0=g1;
00149     }
00150     // modify so thaqt we have l_uu and dux for delta method
00151     // DF feb 15 05
00152     //if (lapprox->hesstype==2 || lapprox->hesstype==3)
00153     if (lapprox->hesstype==2 )
00154     {
00155       if (lapprox->block_diagonal_hessian)
00156       {
00157         //if (ad_comm::wd_flag)
00158         tmpstring = ad_comm::adprogram_name + ".rhes";
00159         ofstream ofs((char*)(tmpstring));
00160             ofs << "   value      std.dev" << endl;
00161         int mmin=lapprox->block_diagonal_hessian->indexmin();
00162         int mmax=lapprox->block_diagonal_hessian->indexmax();
00163         int i,j;
00164         int ii=1;
00165         dvector & u= lapprox->uhat;
00166         for (i=mmin;i<=mmax;i++)
00167         {
00168           if (allocated((*(lapprox->block_diagonal_hessian))(i)))
00169           {
00170             dmatrix m= inv((*(lapprox->block_diagonal_hessian))(i));
00171             dvector d=sqrt(diagonal(m));
00172             int jmin=d.indexmin();
00173             int jmax=d.indexmax();
00174             for (j=jmin;j<=jmax;j++)
00175             {
00176               //if (ii<=u.indexmax())
00177               {
00178                 ofs << setprecision(5) << setscientific()
00179                     << setw(14) << u(ii++) << " " << d(j) << endl;;
00180               }
00181             }
00182           }
00183         }
00184       }
00185       else if (lapprox->bHess)
00186       {
00187         //if (ad_comm::wd_flag)
00188         tmpstring = ad_comm::adprogram_name + ".rhes";
00189         ofstream ofs((char*)(tmpstring));
00190             ofs << "   value      std.dev" << endl;
00191         int mmin=lapprox->bHess->indexmin();
00192         int mmax=lapprox->bHess->indexmax();
00193         //int i,j;
00194         int i;
00195         //int ii=1;
00196         dvector & u= lapprox->uhat;
00197         dvector e(mmin,mmax);
00198         //choleski_decomp(*lapprox->bHess);
00199         int ierr;
00200 
00201         banded_lower_triangular_dmatrix tmp=choleski_decomp(*lapprox->bHess,
00202           ierr);
00203         e.initialize();
00204         for (i=mmin;i<=mmax;i++)
00205         {
00206           e(i)=1.0;
00207           dvector v=solve(tmp,e);
00208           e(i)=0;
00209 
00210           double d=sqrt(v*v);
00211             ofs << setprecision(5) << setscientific()
00212                 << setw(14) << u(i) << " " << d << endl;;
00213         }
00214       }
00215     }
00216     else
00217     {
00218       int i;
00219       //if (ad_comm::wd_flag)
00220       dmatrix m;
00221       tmpstring = ad_comm::adprogram_name + ".rhes";
00222       ofstream ofs((char*)(tmpstring));
00223           ofs << "   value      std.dev" << endl;
00224       //int ii=1;
00225       tmpstring = ad_comm::adprogram_name + ".luu";
00226       uostream ofs1((char*)(tmpstring));
00227       dvector & u= lapprox->uhat;
00228       if (lapprox->hesstype !=3)
00229       {
00230         if (allocated(lapprox->Hess))
00231         {
00232           m= inv(lapprox->Hess);
00233           int mmin=m.indexmin();
00234           int mmax=m.indexmax();
00235           for (i=mmin;i<=mmax;i++)
00236           {
00237             ofs << setprecision(5) << setscientific()
00238                 << setw(14) << u(i) << " " << sqrt(m(i,i)) << endl;;
00239           }
00240           // save l_uu and l_xu for covariance calculations
00241           ofs1 << lapprox->usize << lapprox->xsize;
00242           ofs1 << m;
00243         }
00244         else if (lapprox->sparse_triplet2)
00245         {
00246           dcompressed_triplet & st= *(lapprox->sparse_triplet2);
00247           hs_symbolic& S= *(lapprox->sparse_symbolic2);
00248           get_inverse_sparse_hessian(st,S,ofs1,ofs,lapprox->usize,
00249             lapprox->xsize,u);
00250           // save l_uu and l_xu for covariance calculations
00251         }
00252       }
00253       else
00254       {
00255         if (lapprox->bHess)
00256         {
00257           int ierr=0;
00258           int mmin=lapprox->bHess->indexmin();
00259           int mmax=lapprox->bHess->indexmax();
00260           const banded_lower_triangular_dmatrix& C=
00261             quiet_choleski_decomp(*lapprox->bHess,ierr);
00262           ivector e(mmin,mmax);
00263           e.initialize();
00264           if (ierr==0)
00265           {
00266             ofs1 << lapprox->usize << lapprox->xsize;
00267             for (int i=mmin;i<=mmax;i++)
00268             {
00269               if (i>1) e(i-1)=0;
00270               e(i)=1;
00271               dvector w=solve_trans(C,solve(C,e));
00272               ofs << setprecision(5) << setscientific()
00273                   << setw(14) << u(i) << " " << sqrt(w(i)) << endl;;
00274               ofs1 << w;
00275             }
00276           }
00277           else
00278           {
00279           }
00280         }
00281       }
00282       if (!ofs)
00283       {
00284         cerr << "Error writing to file " << tmpstring << endl;
00285         ad_exit(1);
00286       }
00287       // save l_uu and l_xu for covariance calculations
00288       ofs1 << lapprox->Dux;
00289       if (!ofs1)
00290       {
00291         cerr << "Error writing to file " << tmpstring << endl;
00292         ad_exit(1);
00293       }
00294       ofs1.close();
00295     }
00296 
00297     {
00298       tmpstring = ad_comm::adprogram_name + ".luu";
00299       uistream uis1((char*)(tmpstring));
00300       int i = 0, j = 0;
00301       uis1 >> i >> j;
00302       cout << i << " " << j << endl;
00303     }
00304 
00305     int npts=2;
00306     int on,nopt = 0;
00307     if ( (on=option_match(ad_comm::argc,ad_comm::argv,"-hpts",nopt))>-1)
00308     {
00309       if (nopt !=1)
00310       {
00311         cerr << "Usage -hpts option needs non-negative integer  -- ignored.\n";
00312       }
00313       else
00314       {
00315         npts=atoi(ad_comm::argv[on+1]);
00316       }
00317     }
00318 
00319     double _delta=0.0;
00320 
00321     if ( (on=option_match(ad_comm::argc,ad_comm::argv,"-hsize",nopt))>-1)
00322     {
00323       if (!nopt)
00324       {
00325         cerr << "Usage -hsize option needs number  -- ignored" << endl;
00326       }
00327       else
00328       {
00329         istringstream ist(ad_comm::argv[on+1]);
00330         ist >> _delta;
00331 
00332         if (_delta<=0)
00333         {
00334           cerr << "Usage -hsize option needs positive number  -- ignored.\n";
00335           _delta=0.0;
00336         }
00337       }
00338       if (_delta>0)
00339       {
00340         delta=_delta;
00341       }
00342     }
00343 
00344     // get a number which is exactly representable
00345     double sdelta=1.0+delta;
00346     sdelta-=1.0;
00347     {
00348       //
00349       uostream uos("hessian.bin");
00350       uos << npts;
00351       for (int i=1;i<=nvar;i++)
00352       {
00353         cout << "Estimating row " << i << " out of " << nvar
00354              << " for hessian" << endl;
00355 
00356         for (int j=-npts;j<=npts;j++)
00357         {
00358           if (j !=0)
00359           {
00360             double f=0.0;
00361             double xsave=x(i);
00362             x(i)=xsave+j*sdelta;
00363             g1=(*lapprox)(x,f,this);
00364             x(i)=xsave;
00365             uos << i << j << sdelta << g1;
00366           }
00367           else
00368           {
00369             uos << i << j << sdelta << g0;
00370           }
00371         }
00372       }
00373     }
00374     // check for accuracy
00375     {
00376       uistream uis("hessian.bin");
00377       uis >> npts;
00378       dvector v=get_solution_vector(npts);
00379       v.shift(-npts);
00380       dmatrix tmp(-npts,npts,1,nvar);
00381       dmatrix hess(1,nvar,1,nvar);
00382       ivector iind(-npts,npts);
00383       ivector jind(-npts,npts);
00384       double sd = 0;
00385       int i;
00386       for (i=1;i<=nvar;i++)
00387       {
00388         for (int j=-npts;j<=npts;j++)
00389         {
00390           uis >> iind(j) >> jind(j) >> sd >> tmp(j);
00391         }
00392         hess(i)=(v*tmp).shift(1);
00393         hess(i)/=sd;
00394       }
00395       {
00396         adstring tmpstring="admodel.hes";
00397         if (ad_comm::wd_flag)
00398         {
00399           tmpstring = ad_comm::adprogram_name + ".hes";
00400         }
00401         uostream ofs((char*)tmpstring);
00402         ofs << nvar;
00403         dmatrix shess(1,nvar,1,nvar);
00404         double maxerr=0.0;
00405         for (i=1;i<=nvar;i++)
00406         {
00407           for (int j=1;j<=nvar;j++)
00408           {
00409             shess(i,j)=(hess(i,j)-hess(j,i))/
00410              (1.e-3+sfabs(hess(i,j))+fabs(hess(j,i)));
00411             if (shess(i,j)>maxerr) maxerr=shess(i,j);
00412           }
00413         }
00414         ofstream ofs1("hesscheck");
00415         ofs1 << "maxerr = " << maxerr << endl << endl;
00416         ofs1 << setw(10) << hess << endl << endl;
00417         ofs1 << setw(10) << shess << endl;
00418         ofs << hess;
00419         ofs << gradient_structure::Hybrid_bounded_flag;
00420         initial_params::set_inactive_only_random_effects();
00421         dvector tscale(1,nvar);   // need to get scale from somewhere
00422         /*int check=*/initial_params::stddev_scale(tscale,x);
00423         ofs << tscale;
00424       }
00425     }
00426    /*
00427     first_hessian_flag=0;
00428     double sdelta1;
00429     double sdelta2;
00430     lapprox->fmc1.fringe=1.e-9;
00431     for (int i=1;i<=nvar;i++)
00432     {
00433       hess_calcreport(i,nvar);
00434 
00435       double f=0.0;
00436       double xsave=x(i);
00437       sdelta1=x(i)+delta;
00438       sdelta1-=x(i);
00439       x(i)=xsave+sdelta1;
00440       g1=(*lapprox)(x,f,this);
00441 
00442       sdelta2=x(i)-delta;
00443       sdelta2-=x(i);
00444       x(i)=xsave+sdelta2;
00445       g2=(*lapprox)(x,f,this);
00446       x(i)=xsave;
00447       hess1=(g1-g2)/(sdelta1-sdelta2);
00448 
00449       sdelta1=x(i)+eps*delta;
00450       sdelta1-=x(i);
00451       x(i)=xsave+sdelta1;
00452       g1=(*lapprox)(x,f,this);
00453 
00454       x(i)=xsave-eps*delta;
00455       sdelta2=x(i)-eps*delta;
00456       sdelta2-=x(i);
00457       x(i)=xsave+sdelta2;
00458       g2=(*lapprox)(x,f,this);
00459       x(i)=xsave;
00460 
00461       initial_params::set_inactive_only_random_effects();
00462       initial_params::reset(dvar_vector(x));
00463       double eps2=eps*eps;
00464       hess2=(g1-g2)/(sdelta1-sdelta2);
00465       hess=(eps2*hess1-hess2) /(eps2-1.);
00466 
00467       ofs << hess;
00468     }
00469    */
00470   }
00471 }
00472 
00473 #if defined(USE_ADPVM)
00474 
00479 void function_minimizer::hess_routine_slave_random_effects(void)
00480 {
00481   int nvar=initial_params::nvarcalc(); // get the number of active parameters
00482   //if (adjm_ptr) set_labels_for_hess(nvar);
00483   independent_variables x(1,nvar);
00484   initial_params::xinit(x);        // get the initial values into the x vector
00485   double f;
00486   double delta=1.e-6;
00487   double eps=.1;
00488   gradient_structure::set_YES_DERIVATIVES();
00489 
00490   dvector ddd(1,nvar);
00491   gradcalc(0,ddd);
00492   {
00493     {
00494       (*lapprox)(x,f,this);
00495     }
00496     double sdelta1;
00497     double sdelta2;
00498     lapprox->fmc1.fringe=1.e-9;
00499     for (int i=1;i<=nvar;i++)
00500     {
00501       double f=0.0;
00502       double xsave=x(i);
00503       sdelta1=x(i)+delta;
00504       sdelta1-=x(i);
00505       x(i)=xsave+sdelta1;
00506       (*lapprox)(x,f,this);
00507 
00508       sdelta2=x(i)-delta;
00509       sdelta2-=x(i);
00510       x(i)=xsave+sdelta2;
00511       (*lapprox)(x,f,this);
00512       x(i)=xsave;
00513 
00514       sdelta1=x(i)+eps*delta;
00515       sdelta1-=x(i);
00516       x(i)=xsave+sdelta1;
00517       (*lapprox)(x,f,this);
00518 
00519       x(i)=xsave-eps*delta;
00520       sdelta2=x(i)-eps*delta;
00521       sdelta2-=x(i);
00522       x(i)=xsave+sdelta2;
00523       (*lapprox)(x,f,this);
00524       x(i)=xsave;
00525 
00526       initial_params::set_inactive_only_random_effects();
00527       initial_params::reset(dvar_vector(x));
00528     }
00529   }
00530 }
00531 #endif // #if defined(USE_ADPVM)
00532 
00537 dvector get_solution_vector(int n)
00538 {
00539   int i;
00540   int n1=2*n+1;
00541   dmatrix tmp(1,n1,1,n1);
00542   dvector v(1,n1);
00543   v.initialize();
00544   v(2)=1.0;
00545   dvector c(1,n1);
00546   c.fill_seqadd(-n,1);
00547   tmp.initialize();
00548   tmp(1)=1;
00549   tmp(2)=c;
00550   for (i=3;i<=n1;i++)
00551   {
00552     tmp(i)=elem_prod(tmp(i-1),c);
00553   }
00554   dmatrix tmp1=inv(tmp);
00555   return tmp1*v;
00556 }