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