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