Main Page | Namespace List | Class Hierarchy | Alphabetical List | Class List | Directories | File List | Namespace Members | Class Members | File Members | Related Pages

hjmin.cpp

Go to the documentation of this file.
00001 // This may look like C code, but it is really -*- C++ -*-
00002 /*
00003  ************************************************************************
00004  *
00005  *                        Numerical Math Package
00006  *              Hooke-Jeeves multidimensional minimization
00007  *
00008  * Synopsis
00009  *      double hjmin(b,h,funct)
00010  *      Vector& b                       An initial guess as to the
00011  *                                      location of the minimum.
00012  *                                      On return, it contains the location
00013  *                                      of the minimum the function has found
00014  *      Vector& h                       Initial values for steps along
00015  *                                      each direction
00016  *                                      On exit contains the actual steps
00017  *                                      right before the termination
00018  *      MultivariateFunctor& f          A function being optimized
00019  *
00020  * The hjmin function returns the value of the minimized function f() at the
00021  * point of the found minimum.
00022  *
00023  * An alternative "double hjmin(b,h0,funct)" where h0 is a "double const"
00024  * uses the same initial steps along each direction. No final steps are
00025  * reported, though.
00026  *
00027  * Algorithm
00028  *      Hooke-Jeeves method of direct search for a function minimum
00029  *      The method is of the 0. order (i.e. requiring no gradient computation)
00030  *      See
00031  *      B.Bondi. Methods of optimization. An Introduction - M.,
00032  *      "Radio i sviaz", 1988 - 127 p. (in Russian)
00033  *
00034  * $Id: hjmin.cpp,v 1.1 2004/05/21 21:02:52 maxx Exp $
00035  *
00036  ************************************************************************
00037  */
00038 
00039 
00040 #include "LAStreams.h"
00041 #include "math_num.h"
00042 #include "std.h"
00043 
00044 namespace linalg 
00045 {
00046     using namespace linalg;
00047     
00048 /*
00049  *------------------------------------------------------------------------
00050  *              Class to operate on points in a space (x,f(x))
00051  */
00052 
00053 class FPoint
00054 {
00055   Vector& x;                            // A point in the function's domain
00056   double fval;                          // Function value at the point
00057   MultivariateFunctor& fproc;           // Procedure to compute that value
00058   const bool free_x_on_destructing;     // The flag telling if this FPoint 
00059                                         // "owns" x, and has to dispose of
00060                                         // its dynamic memory on destruction
00061 
00062 public:
00063   FPoint(Vector& b, MultivariateFunctor& f);
00064   FPoint(const FPoint& fp);
00065   ~FPoint(void);
00066 
00067   FPoint& operator = (const FPoint& fp);
00068 
00069   double f(void) const                  { return fval; }
00070 
00071   double fiddle_around(LAStreamIn h);// Examine the function in the
00072                                         // neighborhood of the current point.
00073                                         // h defines the radius of the region
00074 
00075                                         // Proceed in the direction the function
00076                                         // seems to decline
00077   friend void update_in_direction(FPoint& from, FPoint& to);
00078 
00079                                         // Decide whether the region embracing
00080                                         // the local min is small enough
00081   bool is_step_relatively_small(LAStreamIn hs, const double tau);
00082 };
00083 
00084                                 // Construct FPoint from array b
00085 inline FPoint::FPoint(Vector& b, MultivariateFunctor& f)
00086         : x(b), fval(f(b)), fproc(f), free_x_on_destructing(false)
00087 {
00088 }
00089 
00090                                 // Constructor by example
00091 FPoint::FPoint(const FPoint& fp)
00092         : x(*(new Vector(fp.x))), 
00093           fval(fp.fval), fproc(fp.fproc),
00094           free_x_on_destructing(true)
00095 {
00096 }
00097                                 // Destructor
00098 FPoint::~FPoint(void)
00099 {
00100   if( free_x_on_destructing )
00101     delete &x;
00102 }
00103 
00104                                 // Assignment; fproc is assumed the same and
00105                                 // is not copied
00106 inline FPoint& FPoint::operator = (const FPoint& fp)
00107 {
00108   x = fp.x;
00109   fval = fp.fval;
00110   return *this;
00111 }
00112 
00113 /*
00114  * Examine the function f in the vicinity of the current point x
00115  * by making tentative steps fro/back along each coordinate.
00116  * Should the function decrease, x is updated to pinpoint thus
00117  * found new local min.
00118  * The procedure returns the minimal function value found in
00119  * the region.
00120  *
00121  */
00122 double FPoint::fiddle_around(LAStreamIn hs)
00123 {
00124                         // Perform a step along each coordinate
00125   for(LAStreamOut xs(x); !xs.eof(); xs.get() )
00126   {
00127     const double hi = hs.get();
00128     register const double xi_old = xs.peek();   // Old value of x[i]
00129     register double fnew;
00130 
00131     if( xs.peek() = xi_old + hi, (fnew = fproc(x)) < fval )
00132       fval = fnew;                      // Step caused f to decrease, OK
00133     else if( xs.peek() = xi_old - hi, (fnew = fproc(x)) < fval )
00134       fval = fnew;
00135     else                                        // No function decline has been
00136       xs.peek() = xi_old;       // found along this coord, back up
00137   }
00138   return fval;
00139 }                                                
00140 
00141                                 // Proceed in the direction the function
00142                                 // seems to decline
00143                                 // to_new = (to - from) + to
00144                                 // from = to (before modification)
00145 void update_in_direction(FPoint& from, FPoint& to)
00146 {
00147   for(LAStreamOut tos(to.x), froms(from.x); !tos.eof(); )
00148   {
00149     register const double t = tos.peek();
00150     tos.get()  += (t - froms.peek());
00151     froms.get() = t;
00152   }
00153   from.fval = to.fval;
00154   to.fval = (to.fproc)(to.x);
00155 }
00156 
00157                                 // Check to see if the point of minimum has
00158                                 // been located accurately enough
00159 bool FPoint::is_step_relatively_small(LAStreamIn hs, const double tau)
00160 {
00161   for(LAStreamIn xs(x); !hs.eof(); )
00162     if( !(hs.get() /(1 + fabs(xs.get())) < tau) )
00163      return false;
00164   return true;
00165 }
00166 
00167 /*
00168  *------------------------------------------------------------------------
00169  *                          Root module
00170  */
00171 #if 0
00172 double hjmin(Vector& b, Vector& h, MultivariateFunctor& ff)
00173 {
00174                                 // Function Parameters
00175   const double tau = 10*DBL_EPSILON;    // Termination criterion
00176   const double threshold = 1e-8;        // Threshold for the function
00177                                         // decay to be treated as
00178                                         // significant
00179   const double step_reduce_factor = 10;   
00180 
00181   are_compatible(b,h);
00182 
00183   FPoint pmin(b,ff);                    // Point of min
00184   FPoint pbase(pmin);                   // Base point
00185 
00186   for(;;)                       // Main iteration loop
00187   {                             // pmin is the approximation to min so far
00188     if( pbase.fiddle_around(h) < pmin.f() - threshold )
00189     {                                   // Function value dropped significantly
00190       do                                // from pmin to the point pbase
00191             update_in_direction(pmin,pbase);// Keep going in the same direction
00192       while( pbase.fiddle_around(h) < pmin.f() - threshold ); // while it works
00193       pbase = pmin;                     // Save the best approx found
00194     }
00195     else                                // Function didn't fall significantly
00196       if(                               // upon wandering around pbas
00197               h *= 1/step_reduce_factor,        // Try to reduce the step then
00198               pbase.is_step_relatively_small(h,tau) )
00199         return pmin.f();
00200   }
00201 }
00202 
00203                         // The same as above with the only difference
00204                         // initial steps are given to be the same
00205                         // along every direction. The final steps
00206                         // aren't reported back though
00207 double hjmin(Vector& b, const double h0, MultivariateFunctor& f)
00208 {
00209   Vector h(b.q_lwb(),b.q_upb()); h = h0;
00210   return hjmin(b,h,f);
00211 }
00212 #endif // of 0
00213 }

Generated on Wed Dec 15 21:20:29 2004 for vuVolume by  doxygen 1.3.9.1