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