ADMB Documentation  11.1x.2730
 All Classes Files Functions Variables Typedefs Friends Defines
nnewmod2.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: nnewmod2.cpp 2534 2014-10-31 09:32:16Z johnoel $
00003  *
00004  * Author: David Fournier
00005  * Copyright (c) 2008-2012 Regents of the University of California
00006  */
00007 #include <admodel.h>
00008 
00009 double function_minimizer::projected_hess_determinant(const dvector& g,
00010   const int underflow_flag, const dvector& xscale,
00011   const double& _ln_det_proj_jac)
00012 {
00013  double& ln_det_proj_jac=(double&) _ln_det_proj_jac;
00014  int ibreak=-1;
00015  int sgn=0;
00016  double lndet=0.0;
00017  //char ch;
00018  if (!underflow_flag)
00019  {
00020   uistream ifs("admodel.hes");
00021   if (!ifs)
00022   {
00023     cerr << "Error opening file admodel.hes" << endl;
00024   }
00025   int nvar = 0;
00026   ifs >> nvar;
00027   //dmatrix S(1,nvar,1,nvar);
00028   if (nvar > 0)
00029   {
00030     if (nvar != initial_params::nvarcalc())
00031     {
00032       cout << "the number of independent variables is wrong in admodel.hes"
00033          << endl;
00034     }
00035     dmatrix p(1,nvar,1,nvar);
00036     dmatrix p1(1,nvar-1,1,nvar);
00037     dmatrix h(1,nvar,1,nvar);
00038     dvector ss(1,nvar);
00039     ifs >> h;
00040     if (!ifs)
00041     {
00042       cerr << "Error reading the hessian from file admodel.hes" << endl;
00043     }
00044     dvector n=g/norm(g);
00045     // project the standard basis vectors onto the tangent space
00046     int i;
00047     for (i=1;i<=nvar;i++)
00048     {
00049       p(i)=-n(i)*n;
00050       p(i,i)+=1;
00051     }
00052 
00053     for (i=1;i<=nvar;i++)
00054     {
00055       ss(i)=norm(p(i));
00056     }
00057     double minsize = min(ss);
00058 
00059     for (i=1;i<=nvar;i++)
00060     {
00061       if (ss(i) == minsize)
00062       {
00063         ibreak = i;
00064         break;
00065       }
00066       p1(i)=p(i);
00067     }
00068 
00069     int ii;
00070     for (ii=i+1;ii<=nvar;ii++)
00071     {
00072       p1(ii-1)=p(ii);
00073     }
00074 
00075     dmatrix tmpS(1,nvar-1,1,nvar-1);
00076 
00077     //for (ii=1;ii<=nvar-1;ii++)
00078     //{
00079       //for (i=1;i<=nvar;i++)
00080       //{
00081         //p1(ii,i)*=xscale(i);
00082       //}
00083     //}
00084 
00085     for (i=1;i<=nvar-1;i++)
00086     {
00087       tmpS(i,i)=p1(i)*p1(i);
00088       for (int j=1;j<i;j++)
00089       {
00090         tmpS(i,j)=p1(i)*p1(j);
00091         tmpS(j,i)=tmpS(i,j);
00092       }
00093     }
00094     ln_det_proj_jac=ln_det(tmpS,sgn);
00095 
00096     // reset the p1 basis
00097     for (i=1;i<=nvar;i++)
00098     {
00099       if (i==ibreak) break;
00100       p1(i)=p(i);
00101     }
00102 
00103     for (ii=i+1;ii<=nvar;ii++)
00104     {
00105       p1(ii-1)=p(ii);
00106     }
00107 
00108     for (i=1;i<=nvar;i++)
00109     {
00110       for (int j=1;j<i;j++)
00111       {
00112         double tmp=(h(i,j)+h(j,i))/2.;
00113         h(i,j)=tmp;
00114         h(j,i)=tmp;
00115       }
00116     }
00117 
00118     // move to "model space"
00119     for (i=1;i<=nvar;i++)
00120     {
00121       for (int j=1;j<=nvar;j++)
00122       {
00123         h(i,j)/=(xscale(i)*xscale(j));
00124       }
00125     }
00126 
00127     for (i=1;i<nvar;i++)
00128     {
00129       dvector tmp = h*p1(i);
00130       tmpS(i,i)=tmp*p1(i);
00131       for (int j=1;j<i;j++)
00132       {
00133         tmpS(i,j)=tmp*p1(j);
00134         tmpS(j,i)=tmpS(i,j);
00135       }
00136     }
00137     lndet=ln_det(tmpS,sgn);
00138   }
00139   if (sgn <= 0)
00140   {
00141     cerr << "Error restricted Hessian is not positive definite" << endl;
00142   }
00143  }
00144  else
00145  {
00146    lndet=-50.0;
00147  }
00148  return lndet;
00149 }