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