ADMB Documentation  11.1.2537
 All Classes Files Functions Variables Typedefs Friends Defines
mod_hess.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: mod_hess.cpp 2247 2014-08-21 18:21:43Z johnoel $
00003  *
00004  * Author: David Fournier
00005  * Copyright (c) 2008-2012 Regents of the University of California
00006  */
00007 #  include <df1b2fun.h>
00008 #  include <admodel.h>
00009 //#include <parallel.h>
00010 
00011 void hess_calcreport(int i,int nvar);
00012 void hess_errorreport(void);
00013 void set_labels_for_hess(int);
00014 
00015 // estimate the matrix of second derivatives
00016 void ad_update_hess_stats_report(int i,int nvar);
00017 
00018 void function_minimizer::hess_routine(void)
00019 {
00020   if (random_effects_flag && lapprox !=0 )
00021   {
00022     hess_routine_random_effects();
00023   }
00024   else
00025   {
00026 #if defined(USE_ADPVM)
00027     if (!ad_comm::pvm_manager)
00028     {
00029       hess_routine_noparallel();
00030     }
00031     else
00032     {
00033       switch (ad_comm::pvm_manager->mode)
00034       {
00035       case 1: // master
00036         hess_routine_master();
00037         break;
00038       case 2: // slave
00039         hess_routine_slave();
00040         break;
00041       default:
00042         cerr << "Error: Illegal value for pvm_manager->mode." << endl;
00043         ad_exit(1);
00044       }
00045       cout << "finished hess routine" << endl;
00046     }
00047 #else
00048     hess_routine_noparallel();
00049 #endif
00050   }
00051 }
00052 void function_minimizer::hess_routine_noparallel(void)
00053 {
00054   int nvar=initial_params::nvarcalc(); // get the number of active parameters
00055   //if (adjm_ptr) set_labels_for_hess(nvar);
00056   independent_variables x(1,nvar);
00057   initial_params::xinit(x);        // get the initial values into the x vector
00058   double delta=1.e-5;
00059   dvector g1(1,nvar);
00060   dvector g2(1,nvar);
00061   dvector gbest(1,nvar);
00062   dvector hess(1,nvar);
00063   dvector hess1(1,nvar);
00064   dvector hess2(1,nvar);
00065   double eps=.1;
00066   gradient_structure::set_YES_DERIVATIVES();
00067   gbest.fill_seqadd(1.e+50,0.);
00068 
00069   adstring tmpstring="admodel.hes";
00070   if (ad_comm::wd_flag)
00071      tmpstring = ad_comm::adprogram_name + ".hes";
00072   uostream ofs((char*)tmpstring);
00073 
00074   ofs << nvar;
00075   {
00076     {
00077       dvariable vf=0.0;
00078       vf=initial_params::reset(dvar_vector(x));
00079       *objective_function_value::pobjfun=0.0;
00080       pre_userfunction();
00081       vf+=*objective_function_value::pobjfun;
00082       gradcalc(nvar, g1, vf);
00083     }
00084     double sdelta1;
00085     double sdelta2;
00086     for (int i=1;i<=nvar;i++)
00087     {
00088       hess_calcreport(i,nvar);
00089 
00090       double xsave=x(i);
00091       sdelta1=x(i)+delta;
00092       sdelta1-=x(i);
00093       x(i)=xsave+sdelta1;
00094       dvariable vf=0.0;
00095       vf=initial_params::reset(dvar_vector(x));
00096       *objective_function_value::pobjfun=0.0;
00097       pre_userfunction();
00098       vf+=*objective_function_value::pobjfun;
00099       gradcalc(nvar, g1, vf);
00100 
00101       sdelta2=x(i)-delta;
00102       sdelta2-=x(i);
00103       x(i)=xsave+sdelta2;
00104       vf=0.0;
00105       vf=initial_params::reset(dvar_vector(x));
00106       *objective_function_value::pobjfun=0.0;
00107       pre_userfunction();
00108       vf+=*objective_function_value::pobjfun;
00109       gradcalc(nvar, g2, vf);
00110       x(i)=xsave;
00111       hess1=(g1-g2)/(sdelta1-sdelta2);
00112 
00113       sdelta1=x(i)+eps*delta;
00114       sdelta1-=x(i);
00115       x(i)=xsave+sdelta1;
00116       vf=0.0;
00117       vf=initial_params::reset(dvar_vector(x));
00118       *objective_function_value::pobjfun=0.0;
00119       pre_userfunction();
00120       vf+=*objective_function_value::pobjfun;
00121       gradcalc(nvar, g1, vf);
00122 
00123       x(i)=xsave-eps*delta;
00124       sdelta2=x(i)-eps*delta;
00125       sdelta2-=x(i);
00126       x(i)=xsave+sdelta2;
00127       vf=0.0;
00128       vf=initial_params::reset(dvar_vector(x));
00129       *objective_function_value::pobjfun=0.0;
00130       pre_userfunction();
00131       vf+=*objective_function_value::pobjfun;
00132       gradcalc(nvar, g2, vf);
00133       x(i)=xsave;
00134 
00135       vf=initial_params::reset(dvar_vector(x));
00136       double eps2=eps*eps;
00137       hess2=(g1-g2)/(sdelta1-sdelta2);
00138       hess=(eps2*hess1-hess2) /(eps2-1.);
00139 
00140       ofs << hess;
00141       //if (adjm_ptr) ad_update_hess_stats_report(nvar,i);
00142     }
00143   }
00144   ofs << gradient_structure::Hybrid_bounded_flag;
00145   dvector tscale(1,nvar);   // need to get scale from somewhere
00146   /*int check=*/initial_params::stddev_scale(tscale,x);
00147   ofs << tscale;
00148 }
00149 
00150 void function_minimizer::hess_routine_and_constraint(int iprof,
00151   const dvector& g, dvector& fg)
00152 {
00153   int nvar=initial_params::nvarcalc(); // get the number of active parameters
00154   independent_variables x(1,nvar);
00155   initial_params::xinit(x);        // get the initial values into the x vector
00156   double delta=1.e-6;
00157   dvector g1(1,nvar);
00158   dvector g2(1,nvar);
00159   dvector gbest(1,nvar);
00160   dvector hess(1,nvar);
00161   dvector hess1(1,nvar);
00162   dvector hess2(1,nvar);
00163   //double eps=.1;
00164   gradient_structure::set_YES_DERIVATIVES();
00165   gbest.fill_seqadd(1.e+50,0.);
00166   uostream ofs("admodel.hes");
00167   //ofstream ofs5("tmphess");
00168   double lambda=fg*g/norm2(g);
00169   cout << fg-lambda*g << endl;
00170   cout << norm(fg-lambda*g) << " " << fg*g/(norm(g)*norm(fg)) << endl;
00171   ofs << nvar;
00172   {
00173     {
00174       dvariable vf=0.0;
00175       vf=initial_params::reset(dvar_vector(x));
00176       *objective_function_value::pobjfun=0.0;
00177       pre_userfunction();
00178       vf+=*objective_function_value::pobjfun;
00179       vf-=lambda*likeprof_params::likeprofptr[iprof]->variable();
00180       gradcalc(nvar, g1, vf);
00181     }
00182     double sdelta1;
00183     double sdelta2;
00184 
00185     for (int i=1;i<=nvar;i++)
00186     {
00187       hess_calcreport(i,nvar);
00188 
00189       double xsave=x(i);
00190       sdelta1=x(i)+delta;
00191       sdelta1-=x(i);
00192       x(i)=xsave+sdelta1;
00193       dvariable vf=0.0;
00194       vf=initial_params::reset(dvar_vector(x));
00195       *objective_function_value::pobjfun=0.0;
00196       pre_userfunction();
00197       vf+=*objective_function_value::pobjfun;
00198       vf-=lambda*likeprof_params::likeprofptr[iprof]->variable();
00199       gradcalc(nvar, g1, vf);
00200 
00201       sdelta2=x(i)-delta;
00202       sdelta2-=x(i);
00203       x(i)=xsave+sdelta2;
00204       vf=0.0;
00205       vf=initial_params::reset(dvar_vector(x));
00206       *objective_function_value::pobjfun=0.0;
00207       pre_userfunction();
00208       vf+=*objective_function_value::pobjfun;
00209       vf-=lambda*likeprof_params::likeprofptr[iprof]->variable();
00210       gradcalc(nvar, g2, vf);
00211       x(i)=xsave;
00212       hess1=(g1-g2)/(sdelta1-sdelta2);
00213   /*
00214       sdelta1=x(i)+eps*delta;
00215       sdelta1-=x(i);
00216       x(i)=xsave+sdelta1;
00217       vf=0.0;
00218       vf=initial_params::reset(dvar_vector(x));
00219       *objective_function_value::pobjfun=0.0;
00220       pre_userfunction();
00221       vf+=*objective_function_value::pobjfun;
00222       vf-=lambda*likeprof_params::likeprofptr[iprof]->variable();
00223       f=value(vf);
00224       gradcalc(nvar,g1);
00225 
00226       x(i)=xsave-eps*delta;
00227       sdelta2=x(i)-eps*delta;
00228       sdelta2-=x(i);
00229       x(i)=xsave+sdelta2;
00230       vf=0.0;
00231       vf=0.0;
00232       vf=initial_params::reset(dvar_vector(x));
00233       *objective_function_value::pobjfun=0.0;
00234       pre_userfunction();
00235       vf+=*objective_function_value::pobjfun;
00236       vf-=lambda*likeprof_params::likeprofptr[iprof]->variable();
00237       f=value(vf);
00238       gradcalc(nvar,g2);
00239       x(i)=xsave;
00240 
00241       double eps2=eps*eps;
00242       hess2=(g1-g2)/(sdelta1-sdelta2);
00243       hess=(eps2*hess1-hess2) /(eps2-1.);
00244     */
00245       hess=hess1;
00246       ofs << hess;
00247     }
00248   }
00249   gradient_structure::set_NO_DERIVATIVES();
00250 }
00251 /*
00252 void function_minimizer::hess_routine_and_constraint(int iprof)
00253 {
00254   int nvar=initial_params::nvarcalc(); // get the number of active parameters
00255   independent_variables x(1,nvar);
00256   initial_params::xinit(x);        // get the initial values into the x vector
00257   double f;
00258   double delta=1.e-6;
00259   dvector g1(1,nvar);
00260   dvector g2(1,nvar);
00261   dvector gbest(1,nvar);
00262   dvector hess(1,nvar);
00263   dvector hess1(1,nvar);
00264   dvector hess2(1,nvar);
00265   double eps=.1;
00266   gradient_structure::set_YES_DERIVATIVES();
00267   gbest.fill_seqadd(1.e+50,0.);
00268   uostream ofs("admodel.hes");
00269   //ofstream ofs5("tmphess");
00270   ofs << nvar;
00271   {
00272     {
00273       dvariable vf=0.0;
00274       vf=initial_params::reset(dvar_vector(x));
00275       *objective_function_value::pobjfun=0.0;
00276       pre_userfunction();
00277       vf+=*objective_function_value::pobjfun;
00278       f=value(vf);
00279       gradcalc(nvar,g1);
00280     }
00281     for (int i=1;i<=nvar;i++)
00282     {
00283       cout << "Estimating row " << i << " out of " << nvar
00284            << " for hessian" << endl;
00285 
00286       double f=0.0;
00287       double xsave=x(i);
00288       double sdelta1=x(i)+delta;
00289       sdelta1-=x(i);
00290       x(i)=xsave+sdelta1;
00291       dvariable vf=0.0;
00292       vf=initial_params::reset(dvar_vector(x));
00293       *objective_function_value::pobjfun=0.0;
00294       pre_userfunction();
00295       vf+=*objective_function_value::pobjfun;
00296       f=value(vf);
00297       gradcalc(nvar,g1);
00298 
00299       double sdelta2=x(i)-delta;
00300       sdelta2-=x(i);
00301       x(i)=xsave+sdelta2;
00302       vf=0.0;
00303       vf=initial_params::reset(dvar_vector(x));
00304       *objective_function_value::pobjfun=0.0;
00305       pre_userfunction();
00306       vf+=*objective_function_value::pobjfun;
00307       f=value(vf);
00308       gradcalc(nvar,g2);
00309       x(i)=xsave;
00310       hess1=(g1-g2)/(sdelta1-sdelta2);
00311 
00312       sdelta1=x(i)+eps*delta;
00313       sdelta1-=x(i);
00314       x(i)=xsave+sdelta1;
00315       vf=0.0;
00316       vf=initial_params::reset(dvar_vector(x));
00317       *objective_function_value::pobjfun=0.0;
00318       pre_userfunction();
00319       vf+=*objective_function_value::pobjfun;
00320       f=value(vf);
00321       gradcalc(nvar,g1);
00322 
00323       x(i)=xsave-eps*delta;
00324       sdelta2=x(i)-eps*delta;
00325       sdelta2-=x(i);
00326       x(i)=xsave+sdelta2;
00327       vf=0.0;
00328       vf=0.0;
00329       vf=initial_params::reset(dvar_vector(x));
00330       *objective_function_value::pobjfun=0.0;
00331       pre_userfunction();
00332       vf+=*objective_function_value::pobjfun;
00333       f=value(vf);
00334       gradcalc(nvar,g2);
00335       x(i)=xsave;
00336 
00337       double eps2=eps*eps;
00338       hess2=(g1-g2)/(sdelta1-sdelta2);
00339       hess=(eps2*hess1-hess2) /(eps2-1.);
00340       ofs << hess;
00341     }
00342     // do the hessian for the constraint function
00343     uostream cofs("constrnt.hes");
00344     cofs << nvar;
00345     for (i=1;i<=nvar;i++)
00346     {
00347       cout << "Estimating row " << i << " out of " << nvar
00348            << " for hessian" << endl;
00349 
00350       double f=0.0;
00351       double xsave=x(i);
00352       sdelta1=x(i)+delta;
00353       sdelta1-=x(i);
00354       x(i)=xsave+sdelta1;
00355       dvariable vf=0.0;
00356       vf=initial_params::reset(dvar_vector(x));
00357       *objective_function_value::pobjfun=0.0;
00358       pre_userfunction();
00359       vf=likeprof_params::likeprofptr[iprof]->variable();
00360       f=value(vf);
00361       gradcalc(nvar,g1);
00362 
00363       sdelta2=x(i)-delta;
00364       sdelta2-=x(i);
00365       x(i)=xsave+sdelta2;
00366       vf=0.0;
00367       vf=initial_params::reset(dvar_vector(x));
00368       *objective_function_value::pobjfun=0.0;
00369       pre_userfunction();
00370       vf=likeprof_params::likeprofptr[iprof]->variable();
00371       f=value(vf);
00372       gradcalc(nvar,g2);
00373       x(i)=xsave;
00374       hess1=(g1-g2)/(sdelta1-sdelta2);
00375 
00376       sdelta1=x(i)+eps*delta;
00377       sdelta1-=x(i);
00378       x(i)=xsave+sdelta1;
00379       vf=0.0;
00380       vf=initial_params::reset(dvar_vector(x));
00381       *objective_function_value::pobjfun=0.0;
00382       pre_userfunction();
00383       vf=likeprof_params::likeprofptr[iprof]->variable();
00384       f=value(vf);
00385       gradcalc(nvar,g1);
00386 
00387       x(i)=xsave-eps*delta;
00388       sdelta2=x(i)-eps*delta;
00389       sdelta2-=x(i);
00390       x(i)=xsave+sdelta2;
00391       vf=0.0;
00392       vf=initial_params::reset(dvar_vector(x));
00393       *objective_function_value::pobjfun=0.0;
00394       pre_userfunction();
00395       vf=likeprof_params::likeprofptr[iprof]->variable();
00396       f=value(vf);
00397       gradcalc(nvar,g2);
00398       x(i)=xsave;
00399 
00400       double eps2=eps*eps;
00401       hess2=(g1-g2)/(sdelta1-sdelta2);
00402       hess=(eps2*hess1-hess2) /(eps2-1.);
00403       cofs << hess;
00404     }
00405   }
00406   gradient_structure::set_NO_DERIVATIVES();
00407 }
00408 */
00409 
00414 void function_minimizer::depvars_routine(void)
00415 {
00416   reset_gradient_stack();
00417   dvector ggg(1,1);
00418   gradcalc(0,ggg);
00419   gradient_structure::set_YES_DERIVATIVES();
00420   initial_params::restore_start_phase();
00421   int nvar=initial_params::nvarcalc(); // get the number of active parameters
00422   int ndvar=stddev_params::num_stddev_calc();
00423   independent_variables x(1,nvar);
00424   initial_params::xinit(x);        // get the initial values into the x vector
00425   //double f;
00426   //double delta=1.e-7;
00427   adstring tmpstring="admodel.dep";
00428   if (ad_comm::wd_flag)
00429      tmpstring = ad_comm::adprogram_name + ".dep";
00430   ofstream ofs((char*)tmpstring);
00431   if (lapprox)
00432   {
00433     lapprox->no_function_component_flag=1;
00434   }
00435 
00436   dvariable vf;
00437   vf=initial_params::reset(dvar_vector(x));
00438   *objective_function_value::pobjfun=0.0;
00439   pre_userfunction();
00440   vf+=*objective_function_value::pobjfun;
00441 
00442   ofs << nvar << "  "  << ndvar << endl;
00443   int i;
00444   for (i=0;i< stddev_params::num_stddev_params;i++)
00445   {
00446      stddev_params::stddevptr[i]->set_dependent_variables();
00447   }
00448   gradient_structure::jacobcalc(nvar,ofs);
00449   for (i=0;i< stddev_params::num_stddev_params;i++)
00450   {
00451      ofs << stddev_params::stddevptr[i]->label() << "  ";
00452      ofs << stddev_params::stddevptr[i]->size_count() << endl;
00453   }
00454   if (lapprox)
00455   {
00456     lapprox->no_function_component_flag=0;
00457   }
00458   gradient_structure::set_NO_DERIVATIVES();
00459 }
00463 void function_minimizer::hess_inv(void)
00464 {
00465   initial_params::set_inactive_only_random_effects();
00466   int nvar=initial_params::nvarcalc(); // get the number of active parameters
00467   independent_variables x(1,nvar);
00468 
00469   initial_params::xinit(x);        // get the initial values into the x vector
00470   //double f;
00471   dmatrix hess(1,nvar,1,nvar);
00472   uistream ifs("admodel.hes");
00473   int file_nvar = 0;
00474   ifs >> file_nvar;
00475   if (nvar != file_nvar)
00476   {
00477     cerr << "Number of active variables in file mod_hess.rpt is wrong"
00478          << endl;
00479   }
00480 
00481   for (int i = 1;i <= nvar; i++)
00482   {
00483     ifs >> hess(i);
00484     if (!ifs)
00485     {
00486       cerr << "Error reading line " << i  << " of the hessian"
00487            << " in routine hess_inv()" << endl;
00488       exit(1);
00489     }
00490   }
00491   int hybflag = 0;
00492   ifs >> hybflag;
00493   dvector sscale(1,nvar);
00494   ifs >> sscale;
00495   if (!ifs)
00496   {
00497     cerr << "Error reading sscale"
00498          << " in routine hess_inv()" << endl;
00499   }
00500 
00501   double maxerr=0.0;
00502   for (int i = 1;i <= nvar; i++)
00503   {
00504     for (int j=1;j<i;j++)
00505     {
00506       double tmp=(hess(i,j)+hess(j,i))/2.;
00507       double tmp1=fabs(hess(i,j)-hess(j,i));
00508       tmp1/=(1.e-4+fabs(hess(i,j))+fabs(hess(j,i)));
00509       if (tmp1>maxerr) maxerr=tmp1;
00510       hess(i,j)=tmp;
00511       hess(j,i)=tmp;
00512     }
00513   }
00514   /*
00515   if (maxerr>1.e-2)
00516   {
00517     cerr << "warning -- hessian aprroximation is poor" << endl;
00518   }
00519  */
00520 
00521   for (int i = 1;i <= nvar; i++)
00522   {
00523     int zero_switch=0;
00524     for (int j=1;j<=nvar;j++)
00525     {
00526       if (hess(i,j)!=0.0)
00527       {
00528         zero_switch=1;
00529       }
00530     }
00531     if (!zero_switch)
00532     {
00533       cerr << " Hessian is 0 in row " << i << endl;
00534       cerr << " This means that the derivative if probably identically 0 "
00535               " for this parameter" << endl;
00536     }
00537   }
00538 
00539   int ssggnn;
00540   ln_det(hess,ssggnn);
00541   int on1=0;
00542   {
00543     ofstream ofs3((char*)(ad_comm::adprogram_name + adstring(".eva")));
00544     {
00545       dvector se=eigenvalues(hess);
00546       ofs3 << setshowpoint() << setw(14) << setprecision(10)
00547            << "unsorted:\t" << se << endl;
00548      se=sort(se);
00549      ofs3 << setshowpoint() << setw(14) << setprecision(10)
00550      << "sorted:\t" << se << endl;
00551      if (se(se.indexmin())<=0.0)
00552       {
00553         negative_eigenvalue_flag=1;
00554         cout << "Warning -- Hessian does not appear to be"
00555          " positive definite" << endl;
00556       }
00557     }
00558     ivector negflags(0,hess.indexmax());
00559     int num_negflags=0;
00560     {
00561       int on = option_match(ad_comm::argc,ad_comm::argv,"-eigvec");
00562       on1=option_match(ad_comm::argc,ad_comm::argv,"-spmin");
00563       if (on > -1 || on1 >-1 )
00564       {
00565         ofs3 << setshowpoint() << setw(14) << setprecision(10)
00566           << eigenvalues(hess) << endl;
00567         dmatrix ev=trans(eigenvectors(hess));
00568         ofs3 << setshowpoint() << setw(14) << setprecision(10)
00569           << ev << endl;
00570         for (int i=1;i<=ev.indexmax();i++)
00571         {
00572           double lam=ev(i)*hess*ev(i);
00573           ofs3 << setshowpoint() << setw(14) << setprecision(10)
00574             << lam << "  "  << ev(i)*ev(i) << endl;
00575           if (lam<0.0)
00576           {
00577             num_negflags++;
00578             negflags(num_negflags)=i;
00579           }
00580         }
00581         if ( (on1>-1) && (num_negflags>0))   // we will try to get away from
00582         {                                     // saddle point
00583           negative_eigenvalue_flag=0;
00584           spminflag=1;
00585           if(negdirections)
00586           {
00587             delete negdirections;
00588           }
00589           negdirections = new dmatrix(1,num_negflags);
00590           for (int i=1;i<=num_negflags;i++)
00591           {
00592             (*negdirections)(i)=ev(negflags(i));
00593           }
00594         }
00595         int on2 = option_match(ad_comm::argc,ad_comm::argv,"-cross");
00596         if (on2>-1)
00597         {                                     // saddle point
00598           dmatrix cross(1,ev.indexmax(),1,ev.indexmax());
00599           for (int i = 1;i <= ev.indexmax(); i++)
00600           {
00601             for (int j=1;j<=ev.indexmax();j++)
00602             {
00603               cross(i,j)=ev(i)*ev(j);
00604             }
00605           }
00606           ofs3 <<  endl << "  e(i)*e(j) ";
00607           ofs3 << endl << cross << endl;
00608         }
00609       }
00610     }
00611 
00612     if (spminflag==0)
00613     {
00614       if (num_negflags==0)
00615       {
00616         hess=inv(hess);
00617         int on=0;
00618         if ( (on=option_match(ad_comm::argc,ad_comm::argv,"-eigvec"))>-1)
00619         {
00620           int i;
00621           ofs3 << "choleski decomp of correlation" << endl;
00622           dmatrix ch=choleski_decomp(hess);
00623           for (i=1;i<=ch.indexmax();i++)
00624             ofs3 << ch(i)/norm(ch(i)) << endl;
00625           ofs3 << "parameterization of choleski decomnp of correlation" << endl;
00626           for (i=1;i<=ch.indexmax();i++)
00627           {
00628             dvector tmp=ch(i)/norm(ch(i));
00629             ofs3 << tmp(1,i)/tmp(i) << endl;
00630           }
00631         }
00632       }
00633     }
00634   }
00635   if (spminflag==0)
00636   {
00637     if (on1<0)
00638     {
00639       for (int i = 1;i <= nvar; i++)
00640       {
00641         if (hess(i,i) <= 0.0)
00642         {
00643           hess_errorreport();
00644           ad_exit(1);
00645         }
00646       }
00647     }
00648     {
00649       adstring tmpstring="admodel.cov";
00650       if (ad_comm::wd_flag)
00651         tmpstring = ad_comm::adprogram_name + ".cov";
00652       uostream ofs((char*)tmpstring);
00653       ofs << nvar << hess;
00654       ofs << gradient_structure::Hybrid_bounded_flag;
00655       ofs << sscale;
00656     }
00657   }
00658 }
00659 void hess_calcreport(int i,int nvar)
00660 {
00661   if (ad_printf)
00662     (*ad_printf)("Estimating row %d out of %d for hessian\n",i,nvar);
00663   else
00664     cout << "Estimating row " << i << " out of " << nvar << " for hessian"
00665          << endl;
00666 }
00667 void hess_errorreport(void)
00668 {
00669   if (ad_printf)
00670     (*ad_printf)("Hessian does not appear to be positive definite\n");
00671   else
00672     cout << "Hessian does not appear to be positive definite\n" << endl;
00673 }