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