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