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