ADMB Documentation  11.2.2817
 All Classes Files Functions Variables Typedefs Friends Defines
dflocmin.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: dflocmin.cpp 1935 2014-04-26 02:02:58Z johnoel $
00003  *
00004  * Author: David Fournier
00005  * Copyright (c) 2008-2012 Regents of the University of California
00006  */
00011 #  include <admodel.h>
00012 #  include <df1b2fun.h>
00013 #  include <adrndeff.h>
00014 /*
00015         int fcount =0;
00016 static int no_stuff=0;
00017 
00018 double evaluate_function(const dvector& x,function_minimizer * pfmin);
00019 void get_newton_raphson_info(int xs,int us,const init_df1b2vector _y,dmatrix& Hess,
00020   dvector& grad, df1b2_gradlist * f1b2gradlist,function_minimizer * pfmin);
00021 
00022 //dvariable AD_uf_inner(const dvector& x,const dvar_vector& u);
00023 void get_second_ders(int xs,int us,const init_df1b2vector y,dmatrix& Hess,
00024   dmatrix& Dux, df1b2_gradlist * f1b2gradlist,function_minimizer * pfmin,
00025   laplace_approximation_calculator* lap);
00026 
00027   double re_objective_function_value::fun_without_pen=0;
00028 
00029 int laplace_approximation_calculator::saddlepointflag=0;
00030 int laplace_approximation_calculator::print_importance_sampling_weights_flag=0;
00031 
00032 int laplace_approximation_calculator::where_are_we_flag=0;
00033 dvar_vector *
00034   laplace_approximation_calculator::variance_components_vector=0;
00035 */
00036 
00041 dvector laplace_approximation_calculator::local_minimization
00042 (dvector& s,dmatrix& H,dvector& grad,double lambda)
00043 {
00044   dvector vbest(1,usize);
00045   vbest.initialize();
00046 
00047   int better_flag=0;
00048   int counter=0;
00049   s.initialize();
00050   s(1)=1.0;
00051   double fbest=evaluate_function_no_derivatives(uhat,pmin);
00052   do
00053   {
00054     dvector v=local_minimization_routine(s,H,grad,lambda);
00055     dvector xx=uhat+v;
00056     double f2=evaluate_function_no_derivatives(xx,pmin);
00057     cout << endl << fbest-f2 << endl;
00058     if(f2<fbest)
00059     {
00060       better_flag=1;
00061       fbest=f2;
00062       lambda*=5.0;
00063       vbest=v;
00064       s=v;
00065     }
00066     else
00067     {
00068       if (better_flag==1)
00069       {
00070         // we have a better value so go with it
00071         break;
00072       }
00073       else
00074       {
00075         // try a smaller trust region
00076         lambda/=5;
00077         s=vbest;
00078       }
00079     }
00080   }
00081   while (counter<20);
00082 
00083   if (better_flag==0)
00084   {
00085     cerr << "Error cannot find better value to try and get a"
00086       " positive definite hessian" << endl;
00087   }
00088 
00089   return vbest;
00090 }
00091 
00096 dvector laplace_approximation_calculator::local_minimization_routine
00097 (dvector& s,dmatrix& H,dvector& grad,double lambda)
00098 {
00099   double f=0.0;
00100   double fb=1.e+100;
00101   dvector g(1,usize);
00102   dvector ub(1,usize);
00103   fmc1.itn=0;
00104   fmc1.ifn=0;
00105   fmc1.ireturn=0;
00106   fmc1.ialph=0;
00107   fmc1.ihang=0;
00108   fmc1.ihflag=0;
00109   fmc1.crit=1.e-12;
00110   long fmsave = fmc1.maxfn;
00111   fmc1.maxfn=1000;;
00112 
00113   fmc1.dfn=1.e-2;
00114   while (fmc1.ireturn>=0)
00115   {
00116     fmc1.fmin(f,s,g);
00117     if (fmc1.ireturn>0)
00118     {
00119       double ns=norm(s);
00120       double ns2=square(ns);
00121       dvector v=s/ns;
00122 
00123       dvector z=H*v;
00124       double vHv=v*z;
00125 
00126       double gradv=grad*v;
00127       f=lambda*gradv+0.5*lambda*lambda*vHv+ square(ns2-1.0);
00128       //f=0.5*lambda*lambda*s*H*s;
00129       if (f<fb)
00130       {
00131         fb=f;
00132         ub=s;
00133       }
00134       g=lambda*grad/ns -lambda * gradv*s/ns2
00135            + lambda * lambda * z/ns
00136            - lambda * lambda * vHv*s/ns2 + 4.0*(ns2-1.0)*s;
00137     }
00138   }
00139   s=ub;
00140   cout <<  " inner maxg = " <<  fmc1.gmax;
00141 
00142   fmc1.maxfn=fmsave;
00143   fmc1.ireturn=0;
00144   fmc1.fbest=fb;
00145   return ub;
00146 }
00147 
00148  //
00149  //
00150  //dvector laplace_approximation_calculator::local_minimization_routine
00151  //(dvector& s,dmatrix& H,dvector& grad,double lambda)
00152  //{
00153  //  double f=0.0;
00154  //  double fb=1.e+100;
00155  //  dvector g(1,usize);
00156  //  dvector ub(1,usize);
00157  //  fmc1.itn=0;
00158  //  fmc1.ifn=0;
00159  //  fmc1.ireturn=0;
00160  //  fmc1.ialph=0;
00161  //  fmc1.ihang=0;
00162  //  fmc1.ihflag=0;
00163  //  fmc1.crit=1.e-12;
00164  //  double beta=.1;
00165  //
00166  //  s.initialize();
00167  //
00168  //  fmc1.dfn=1.e-2;
00169  //  while (fmc1.ireturn>=0)
00170  //  {
00171  //    fmc1.fmin(f,s,g);
00172  //    if (fmc1.ireturn>0)
00173  //    {
00174  //      double den=lambda-norm2(s);
00175  //
00176  //      if (den<=0)
00177  //      {
00178  //        f=1.e+60;
00179  //      }
00180  //      else
00181  //      {
00182  //        f=grad*s+0.5*(s*(H*s))+0.5*beta/den;
00183  //        if (f<fb)
00184  //        {
00185  //          fb=f;
00186  //          ub=s;
00187  //        }
00188  //        g=grad + H*s + lambda*s +beta*s/square(den);
00189  //      }
00190  //    }
00191  //  }
00192  //  s=ub;
00193  //  cout <<  " inner maxg = " <<  fmc1.gmax;
00194  //
00195  //  fmc1.ireturn=0;
00196  //  fmc1.fbest=fb;
00197  //  return ub;
00198  //}