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