dep_hess.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  */
00012
00013 // estimate the matrix of second derivatives
00014
00019 dmatrix function_minimizer::dep_hess_routine(const dvariable& dep)
00020 {
00021   int nvar=initial_params::nvarcalc(); // get the number of active parameters
00022   independent_variables x(1,nvar);
00023   initial_params::xinit(x);        // get the initial values into the x vector
00024   dvector scale(1,nvar);   // need to get scale from somewhere
00025   dvector tscale(1,nvar);   // need to get scale from somewhere
00026   initial_params::stddev_scale(scale,x);
00027   //check=initial_params::stddev_curvscale(curv,x);
00028   double delta=1.e-7;
00029   dvector g1(1,nvar);
00030   dvector depg(1,nvar);
00031   dvector g2(1,nvar);
00032   dvector gbest(1,nvar);
00033   dmatrix hess(1,nvar,1,nvar);
00034   dvector hess1(1,nvar);
00035   dvector hess2(1,nvar);
00036   double eps=.1;
00039   {
00041     dvariable vf=0.0;
00042     vf=initial_params::reset(dvar_vector(x));
00043     userfunction();
00044     vf=dep;
00046   }
00047   double sdelta1;
00048   double sdelta2;
00049   int i;
00050   for (i=1;i<=nvar;i++)
00051   {
00052     cout << "Estimating row " << i << " out of " << nvar
00053          << " for dependent variable hessian" << endl;
00054
00055     double xsave=x(i);
00056     sdelta1=x(i)+delta;
00057     sdelta1-=x(i);
00058     x(i)=xsave+sdelta1;
00059     dvariable vf=0.0;
00060     initial_params::stddev_scale(tscale,x);
00062     vf=initial_params::reset(dvar_vector(x));
00063     userfunction();
00064     vf=dep;
00066     g1=elem_div(g1,tscale);
00067
00068     sdelta2=x(i)-delta;
00069     sdelta2-=x(i);
00070     x(i)=xsave+sdelta2;
00071     initial_params::stddev_scale(tscale,x);
00073     vf=0.0;
00074     vf=initial_params::reset(dvar_vector(x));
00075     userfunction();
00076     vf=dep;
00078     g2=elem_div(g2,tscale);
00079     x(i)=xsave;
00080     hess1=(g1-g2)/(sdelta1-sdelta2)/scale(i);
00081
00082     sdelta1=x(i)+eps*delta;
00083     sdelta1-=x(i);
00084     x(i)=xsave+sdelta1;
00085     initial_params::stddev_scale(tscale,x);
00087     vf=0.0;
00088     vf=initial_params::reset(dvar_vector(x));
00089     userfunction();
00090     vf=dep;
00092     g1=elem_div(g1,tscale);
00093
00094     x(i)=xsave-eps*delta;
00095     sdelta2=x(i)-eps*delta;
00096     sdelta2-=x(i);
00097     x(i)=xsave+sdelta2;
00098     initial_params::stddev_scale(tscale,x);
00100     vf=0.0;
00101     vf=initial_params::reset(dvar_vector(x));
00102     userfunction();
00103     vf=dep;
00105     g2=elem_div(g2,tscale);
00106     x(i)=xsave;
00107
00108     hess2=(g1-g2)/(sdelta1-sdelta2)/scale(i);
00109     double eps2=eps*eps;
00110     hess(i)=(eps2*hess1-hess2) /(eps2-1.);
00111   }
00112   /*
00113   for (i=1;i<=nvar;i++)
00114   {
00115     hess(i,i)/=(scale(i)*scale(i));
00116     for (int j=1;j<i;j++)
00117     {
00118       double tmp=(hess(i,j)+hess(j,i))/(2.*scale(i)*scale(j));
00119       hess(i,j)=tmp;
00120       hess(j,i)=tmp;
00121     }
00122     hess(i,i)-=depg(i)*curv(i)/(scale(i)*scale(i)*scale(i));
00123   }
00124  */
00125   return hess;
00126 }
```