ADMB Documentation: para3.cpp (* $Id: para3.cpp 542 2012-07-10 21:04:06Z johnoel $) Source File

Go to the documentation of this file.
00001 /*
00002  * $Id: para3.cpp 542 2012-07-10 21:04:06Z johnoel $
00003  *
00004  * Author: David Fournier
00005  * Copyright (c) 2008-2012 Regents of the University of California 
00006  */
00007 #include <admodel.h>
00008 //#include <iomanip.h>
00009 #ifdef __GNUDOS__
00010   #include <gccmanip.h>
00011 #endif
00012 
00013 void hess_calcreport(int i,int nvar);
00014 void hess_errorreport(void);
00015 void set_labels_for_hess(int);
00016 
00017 class admb_javapointers;
00018 extern admb_javapointers * adjm_ptr;
00019 void useless(const double& sdelta2);
00020 // estimate the matrix of second derivatives
00021 void ad_update_hess_stats_report(int i,int nvar);
00022 
00023 #if defined(USE_ADPVM)
00024 void function_minimizer::hess_routine_master()
00025 {
00026 
00027   int nvar=initial_params::nvarcalc(); // get the number of active parameters
00028   //if (adjm_ptr) set_labels_for_hess(nvar);
00029   independent_variables x(1,nvar);
00030   initial_params::xinit(x);        // get the initial values into the x vector
00031   double f;
00032   double delta=1.e-6;
00033   dvector g1(1,nvar);
00034   dvector g2(1,nvar);
00035   dvector gbest(1,nvar);
00036   dvector hess(1,nvar);
00037   dvector hess1(1,nvar);
00038   dvector hess2(1,nvar);
00039   double eps=.1;
00040   gradient_structure::set_YES_DERIVATIVES();
00041   gbest.fill_seqadd(1.e+50,0.);
00042 
00043   adstring tmpstring="admodel.hes";
00044   if (ad_comm::wd_flag)
00045      tmpstring = ad_comm::adprogram_name + ".hes";
00046   uostream ofs((char*)tmpstring);
00047     
00048   ofs << nvar;
00049   {
00050     pvm_master_function_evaluation(f,x,g1,nvar);
00051     double sdelta1;
00052     double sdelta2;
00053     for (int i=1;i<=nvar;i++)
00054     {
00055 #if defined (__SPDLL__)
00056       hess_calcreport(i,nvar);
00057 #else
00058       cout << "Estimating row " << i << " out of " << nvar
00059            << " for hessian" << endl;
00060 #endif
00061 
00062       double f=0.0;
00063       double xsave=x(i);
00064       sdelta1=x(i)+delta;
00065       useless(sdelta1);
00066       sdelta1-=x(i);
00067       x(i)=xsave+sdelta1;
00068 
00069       pvm_master_function_evaluation(f,x,g1,nvar);
00070 
00071       sdelta2=x(i)-delta;
00072       useless(sdelta2);
00073       sdelta2-=x(i);
00074       x(i)=xsave+sdelta2;
00075 
00076       pvm_master_function_evaluation(f,x,g2,nvar);
00077 
00078       x(i)=xsave;
00079       hess1=(g1-g2)/(sdelta1-sdelta2);
00080  
00081       sdelta1=x(i)+eps*delta;
00082       useless(sdelta1);
00083       sdelta1-=x(i);
00084       x(i)=xsave+sdelta1;
00085 
00086       pvm_master_function_evaluation(f,x,g1,nvar);
00087 
00088       x(i)=xsave-eps*delta;
00089       sdelta2=x(i)-eps*delta;
00090       useless(sdelta2);
00091       sdelta2-=x(i);
00092       x(i)=xsave+sdelta2;
00093 
00094       pvm_master_function_evaluation(f,x,g2,nvar);
00095 
00096       x(i)=xsave;
00097 
00098       dvariable vf=initial_params::reset(dvar_vector(x));
00099       double eps2=eps*eps;
00100       hess2=(g1-g2)/(sdelta1-sdelta2);
00101       hess=(eps2*hess1-hess2) /(eps2-1.);
00102    
00103       ofs << hess;
00104       //if (adjm_ptr) ad_update_hess_stats_report(nvar,i);
00105     }
00106   }
00107   gradient_structure::set_NO_DERIVATIVES();
00108 }
00109 
00110 void function_minimizer::hess_routine_slave()
00111 {
00112   int nvar=initial_params::nvarcalc(); // get the number of active parameters
00113   gradient_structure::set_YES_DERIVATIVES();
00114 
00115   pvm_slave_function_evaluation();
00116   for (int i=1;i<=nvar;i++)
00117   {
00118     for (int it=1;it<=4;it++) pvm_slave_function_evaluation();
00119   }
00120   gradient_structure::set_NO_DERIVATIVES();
00121 }
00122 #endif // #if defined(USE_ADPVM)
00123