cern_minimize.h

/*
  -------------------------------------------------------------------
  
  Copyright (C) 2006, 2007, 2008, Andrew W. Steiner
  
  This file is part of O2scl.
  
  O2scl is free software; you can redistribute it and/or modify
  it under the terms of the GNU General Public License as published by
  the Free Software Foundation; either version 3 of the License, or
  (at your option) any later version.
  
  O2scl is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  GNU General Public License for more details.
  
  You should have received a copy of the GNU General Public License
  along with O2scl. If not, see <http://www.gnu.org/licenses/>.

  -------------------------------------------------------------------
*/
#ifndef O2SCL_CERN_MINIMIZE_H
#define O2SCL_CERN_MINIMIZE_H

#include <o2scl/minimize.h>

#ifndef DOXYGENP
namespace o2scl {
#endif

  /** 
      \brief One-dimensional minimization (CERNLIB)

      The golden section search is applied in the interval  \f$ (a,b) \f$ 
      using a fixed number \f$ n \f$ of function evaluations where 
      \f[
      n=\left[ 2.08 \ln(|a-b|/\mathrm{tolx})+\frac{1}{2}\right]+1 
      \f]
      
      The accuracy depends on the function. A choice of \f$
      \mathrm{tolx}>10^{-8} \f$ usually results in a relative error of
      $x$ which is smaller than or of the order of \f$ \mathrm{tolx}
      \f$ .
      
      This routine strictly searches the interval \f$ (a,b) \f$ . If
      the function is nowhere flat in this interval, then min_bkt()
      will return either \f$ a \f$ or \f$ b \f$ and \ref min_type
      is set to 1.
      
      \note The number of function evaluations can be quite large if
      multi_min::tolx is sufficiently small. If multi_min::tolx is
      exactly zero, then \f$ 10^{-8} \f$ will be used instead.

      Based on the CERNLIB routines RMINFC and DMINFC, which was 
      based on \ref Fletcher87, and \ref Krabs83 and is documented at
      http://wwwasdoc.web.cern.ch/wwwasdoc/shortwrupsdir/d503/top.html
  */
#ifdef DOXYGENP
  template<class param_t, class func_t=funct<param_t> > class cern_minimize : 
  public minimize_bkt
#else
  template<class param_t, class func_t=funct<param_t> > class cern_minimize : 
  public minimize_bkt<param_t,func_t> 
#endif
    {
      
#ifndef DOXGYEN_INTERNAL

      protected:

      /// C analog of Fortran's "Nearest integer" function 
      inline int nint(double x) {
        if (x<0.0) return ((int)(x-0.5));
        else return ((int)(x+0.5));
      }
    
#endif      

    public:
      
      cern_minimize() {
        delta_set=false;
        min_type=0;
      }
    
      /** 
          \brief Calculate the minimum \c min of \c func between
          \c a and \c b
          
          The initial value of \c x is ignored.
          
          If there is no minimum in the given interval, then on exit
          \c x will be equal to either \c a or \c b and \ref min_type
          will be set to 1 instead of zero. The error handler is not
          called, as this need not be interpreted as an error.
      */
      virtual int min_bkt(double &x, double a, double b, double &y,
                          param_t &pa, func_t &func) {
      
        // The square root of 5
        static const double w5=2.23606797749979;
        static const double hv=(3.0-w5)/2.0, hw=(w5-1.0)/2.0;
        double c, d, v=0.0, fv, w=0.0, fw, h;
      
        if (!delta_set) delta=10.0*this->tolx;
        
        int n=-1;
        if (this->tolx==0.0) {
          n=nint(2.0*log(fabs((a-b)/1.0e-8)));
        } else {
          if (a!=b) n=nint(2.0*log(fabs((a-b)/this->tolx)));
        }
        this->last_ntrial=n;
        c=a;
        d=b;
        if (a>b) {
          c=b;
          d=a;
        }
        bool llt=true, lge=true;
      
        // This is never an infinite loop so long as 'n' is finite. We
        // should maybe check somehow that 'n' is not too large
        while (true) {
          h=d-c;
          if (this->verbose>0) {
            x=(c+d)/2.0;
            func(x,y,pa);
            this->print_iter(x,y,this->last_ntrial-n,0.0,this->tolx);
          }
          if (n<0) {
            x=(c+d)/2.0;
            func(x,y,pa);
            if ((fabs(x-a)>delta && fabs(x-b)>delta)) min_type=0;
            min_type=1;
            return 0;
          }
          if (llt) {
            v=c+hv*h;
            func(v,fv,pa);
          }
          if (lge) {
            w=c+hw*h;
            func(w,fw,pa);
          }
          if (fv<fw) {
            llt=true;
            lge=false;
            d=w;
            w=v;
            fw=fv;
          } else {
            llt=false;
            lge=true;
            c=v;
            v=w;
            fv=fw;
          }
          n--;
        }

        O2SCL_ERR("Failed sanity check in cern_minimize::min_bkt().",
                  gsl_esanity);
        return gsl_esanity;
      }
    
      /** 
          \brief Set the value of  \f$ \delta \f$ 
          
          If this is not called before min_bkt() is used, then the
          suggested value \f$ \delta=10 \mathrm{tolx} \f$ is used.
      */
      int set_delta(double d) {
        if (d>0.0) {
          delta=d;
          delta_set=true;
        }
        return 0;
      }
    
      /// Return string denoting type ("cern_minimize")
      virtual const char *type() { return "cern_minimize"; }

      /// Type of minimum found
      int min_type;

    protected:

#ifndef DOXYGEN_INTERNAL

      /// The value of delta as specified by the user
      double delta;

      /// True if the value of delta has been set
      bool delta_set;

#endif

    };

#ifndef DOXYGENP
}
#endif

#endif

---

Documentation generated with Doxygen and provided under the GNU Free Documentation License. See License Information for details.

Project hosting provided by , O2scl Sourceforge Project Page