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