ADMB Documentation  11.1.1890
 All Classes Files Functions Variables Typedefs Friends Defines
xfmmtr1.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: xfmmtr1.cpp 1890 2014-04-15 23:00:16Z johnoel $
00003  *
00004  * Author: David Fournier
00005  * Copyright (c) 2008-2012 Regents of the University of California
00006  */
00011 // this is to get UNIX systems to use getchar
00012 // #define UNIXKLUDGE
00013 
00014 #ifdef __ZTC__
00015   #include <conio.h>
00016 #endif
00017 
00018 #include <admodel.h>
00019 extern int ctlc_flag;
00020 
00021 #if defined(__TURBOC__) && !defined(__linux__)
00022   #pragma hdrstop
00023   #include <iostream.h>
00024   #include <conio.h>
00025 #endif
00026 
00027 #if defined (__WAT32__) || defined(_MSC_VER)
00028   #include <conio.h>
00029 #endif
00030 
00031 #ifdef __ZTC__
00032   #include <iostream.hpp>
00033   #include <disp.h>
00034   #define endl "\n"
00035 //  #define if (ad_printf) (*ad_printf) disp_if (ad_printf) (*ad_printf)
00036 #endif
00037 
00038 #ifdef __SUN__
00039   #include <iostream.h>
00040   #include <signal.h>
00041   #define getch getchar
00042 #endif
00043 
00044 #if defined(__GNU__) || defined(UNIXKLUDGE)
00045   #if (__GNUC__ >3)
00046      #include <iostream>
00047      using namespace std;
00048   #else
00049     #include <iostream.h>
00050   #endif
00051   #include <signal.h>
00052   #define getch getchar
00053 #endif
00054 
00055 extern "C" void onintr(int k);
00056 
00057 #ifdef __NDPX__
00058   #include <iostream.hxx>
00059 #endif
00060 
00061 #include <math.h>
00062 #include <stdlib.h>
00063 #include <stdio.h>
00064 #include <ctype.h>
00065 
00070 void do_evaluation(double& f,independent_variables& x,dvector& g,int nvar,
00071   function_minimizer * pmp)
00072 {
00073   *objective_function_value::pobjfun=0.0;
00074   dvariable vf=initial_params::reset(dvar_vector(x));
00075   pmp->userfunction();
00076   vf+=*objective_function_value::pobjfun;
00077   gradcalc(nvar,g);
00078   f=value(vf);
00079 }
00080 
00085 double get_second_derivative(double f,independent_variables& x,
00086   dvector& g,dvector & r,int nvar,function_minimizer * pmp)
00087 {
00088   const double stepsize=1.e-5;
00089   dvector g1(1,nvar);
00090   dvector g2(1,nvar);
00091   x+=stepsize*r;
00092   do_evaluation(f,x,g1,nvar,pmp);
00093   x-=2.*stepsize*r;
00094   do_evaluation(f,x,g2,nvar,pmp);
00095   double scder=r*(g1-g2)/(2.0*stepsize);
00096   cout << " f = " << f << endl;
00097   cout << "  second derivative =  " ;
00098   cout  << " r*(g1-g)/stepsize  = " << scder << endl;
00099   return scder;
00100 }
00101 
00102 
00103 dvector update1(int nvar, int iter, int m, const dvector& g,
00104   const dmatrix& xalpha, dmatrix& y, const dvector& x, const dvector& xold,
00105   const dvector& gold, const dvector& xrho);
00106 
00107 double dafsqrt( double x );
00108 
00113 void fmmt1::fmin2(const double& _f, const independent_variables &_x,
00114   const dvector& _g, function_minimizer *pmp)
00115 {
00116   //int itn=0; int bigbreak=0; int smallbreak=0; int midbreak=0;
00117   int itn=0; int smallbreak=0; int midbreak=0;
00118   int nvar=initial_params::nvarcalc(); // get the number of active
00119   //double a,f, curf, rnorm, stepsize,b,epsilon;
00120   double a,f, curf, stepsize,b,epsilon;
00121   independent_variables x(1,nvar);
00122   initial_params::xinit(x);    // get the initial values into the
00123   dvariable vf=0.0;            // x vector
00124   epsilon=1.e-2;
00125 
00126   dvector curx(1,nvar), g1(1,nvar), xtry(1,nvar), g(1,nvar), r(1,nvar);
00127   do_evaluation(f,x,g,nvar,pmp); // get initial vales for f and g
00128   curf=f; curx=x;
00129   cout << " f = " << f << endl;
00130 
00131 
00132   do
00133   {
00134     r=update1(n,itn,xm,g,xstep,xy,x,xold,gold,xrho); // get search
00135 
00136     cout << "  norm(g) =  " << norm(g) ;
00137     cout << "  r*g/norm(g) =  " << r*g/norm(g) << endl;
00138     do
00139     {
00140       x=curx;
00141 
00142       a=get_second_derivative(f,x,g,r,nvar,pmp);
00143       b=r*g;
00144 
00145       stepsize=-b/a;
00146       do
00147       {
00148         xtry=curx+stepsize*r; x=xtry;
00149 
00150         do_evaluation(f,x,g,nvar,pmp);
00151         cout << " f = " << f << endl;
00152         cout << "  r*g/norm(g) =  " << r*g/norm(g) << endl;
00153 
00154         if (f<curf+1.e-10)
00155         {
00156           curx=x; curf=f;
00157           smallbreak=1;
00158           if (fabs(g*r)/norm(g)<epsilon)
00159           {
00160             midbreak=1;
00161           }
00162         }
00163         else
00164         {
00165           cout << setprecision(10) << f-curf << endl;
00166           stepsize=0.001*stepsize; xtry=curx+stepsize*r;
00167         }
00168       }
00169       while(!smallbreak);
00170       smallbreak=0;
00171     }
00172     while (!midbreak);
00173     midbreak=0;
00174     itn++;
00175   }
00176   while (1);
00177 
00178       exit(1);
00179 }
00180 
00185 dvector update1(int nvar, int iter, int m, const dvector& g, const dmatrix& _s,
00186   dmatrix& y, const dvector& x, const dvector& _xold, const dvector& _gold,
00187   const dvector& _xrho)
00188   {
00189     dvector& xold= (dvector&) _xold;
00190     dmatrix& s= (dmatrix&) _s;
00191     dvector& gold= (dvector&) _gold;
00192     dvector& xrho=(dvector&)_xrho;
00193     dvector beta(1,nvar);
00194     dvector alpha(0,m);
00195     dvector r(1,nvar);
00196     dvector t(1,nvar);
00197     int m1=m+1;
00198     if (iter<1)
00199     {
00200       xold=x;
00201       gold=g;
00202       r=g;
00203     }
00204     else
00205     {
00206       int k=iter-1;
00207       int k1=k%(m1);
00208       y(k1)=g-gold;
00209       s(k1)=x-xold;
00210       xrho(k1)=1./(y(k1)*s(k1));
00211       xold=x;
00212       gold=g;
00213 
00214       int i;
00215       int lb=k-m+1;
00216       if (lb <0) lb=0;
00217       t=g;
00218       for (i=k;i>=lb;i--)
00219       {
00220         int i1=i%(m1);
00221         //int i2=(i+1)%(m1);
00222         //if (i==k)
00223         {
00224           alpha(i-lb)=xrho(i1)*(s(i1)*t);
00225         }
00226         //else
00227         //{
00228         //  alpha(i-lb)=xrho(i1)*( (s(i1)-(s(i1)*s(i2))*s(i2)) * t);
00229         //}
00230         t-=alpha(i-lb)*y(i1);
00231       }
00232       r=t;
00233       for (i=lb;i<=k;i++)
00234       {
00235         int i1=i%(m1);
00236         //r+= (alpha(i1)-xrho(i1)*(y(i1)*r)) * s(i1);
00237         r+= (alpha(i-lb)-xrho(i1)*(y(i1)*r)) * s(i1);
00238       }
00239     }
00240     //cout << r*g/(norm(r)*norm(g)) << endl;
00241     r/=norm(r);
00242     return -1.0*r;
00243   }