GeneralFunctionMinimizer.h

//===========================================================================
// GoTools - SINTEF Geometry Tools version 1.1
//
// GoTools module: CORE
//
// Copyright (C) 2000-2007 SINTEF ICT, Applied Mathematics, Norway.
//
// This program 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 version 2 of the License. 
//
// This program 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 this program; if not, write to the Free Software            
// Foundation, Inc.,                                                      
// 59 Temple Place - Suite 330,                                           
// Boston, MA  02111-1307, USA.                                           
//
// Contact information: E-mail: tor.dokken@sintef.no                      
// SINTEF ICT, Department of Applied Mathematics,                         
// P.O. Box 124 Blindern,                                                 
// 0314 Oslo, Norway.                                                     
//
// Other licenses are also available for this software, notably licenses
// for:
// - Building commercial software.                                        
// - Building software whose source code you wish to keep private.        
//===========================================================================
#ifndef _GENERAL_FUNCTIONMINIMIZER_H
#define _GENERAL_FUNCTIONMINIMIZER_H

#include "Point.h"
#include <vector>
#include <limits>

namespace Go {


template<class Functor> class FunctionMinimizer;


//===========================================================================
template<class Functor>
void minimise_conjugated_gradient(FunctionMinimizer<Functor>& dfmin);
//===========================================================================

template<class Functor>
class FunctionMinimizer
{
 public:
    FunctionMinimizer(int num_param, const Functor& fun, 
                      const double* const seed, 
                      double tol = std::sqrt(std::numeric_limits<double>::epsilon()));

    ~FunctionMinimizer() 
    {}
    
    double minimize(const Point& dir, bool& hit_domain_edge); 


    double linminBrent(const Point& dir,
                       const double* bracket,
                       const double* fval_brak);


    inline double fval() const; 

    inline double fval(const Point& param) const;

    inline void grad(Point& result) const;

    inline void grad(const Point& param, Point& result) const;

    inline void moveUV(const Point& dir, double multiplier);

    bool atMin(int param_ix) const {return at_min_[param_ix];}
    bool atMax(int param_ix) const {return at_max_[param_ix];}

    double getPar(int param_ix) const {return par_[param_ix];}

    const double* getPar() const {return par_.begin();}

    int numPars() const {return par_.size();}
                              
 private:
    const Functor fun_;

    Point par_;
    std::vector<bool> at_min_;
    std::vector<bool> at_max_;

    mutable bool cached_value_updated_;
    mutable bool cached_grad_updated_;

    mutable double cached_value_;
    mutable Point cached_grad_;

    const double minimization_tol_;
    std::vector<double> param_tol_;

    static const double root_machine_precision_;
    static const double rel_tol_;
    static const double perturbation_;
    static const double golden_ratio_;
    static const double default_partition_;
    static const double tiny_;
    static const int max_iter_;


    // update at_umin_, at_umax_, at_vmin_ and at_vmax_.  Should be called whenever the
    // current parameters ('cur_uv_') are changed.
    inline void checkBorder();

    // clear cached values.  Should be called whenever the current parameters ('cur_uv_') 
    // are changed.
    inline void resetCache() const;

    // determine how far we are allowed to step from the current parameter point in a certain
    // direction ('dir') before going out of the valid parameter domain.
    double determineMaxSteplength(const Point& dir) const;

    // Bracket the function minimum from the current parameter point along a certain 
    // direction ('dir'). The maximum allowed steplength is given in 'max_step'.  
    // Upon completion, two things might happen:
    // Possibility 1 : the function was able to bracket a minimum.  The position of the 
    //                 bracketing points (from current parameters, along 'dir') is 
    //                 returned in 'bracket' (3 values), and the corresponding distance
    //                 function values are returned in fval_brak.  The function returns 'true'.
    // Possibility 2 : the function hit the 'max_step' before it was able to bracket any
    //                 minimum.  The function returns 'false'.
    bool bracketInterval(const Point& dir,  
                         const double max_step,
                         //const double working_tolerance,
                         double* bracket,    // points to 3-array
                         double* fval_brak); // points to 3-array

    // given a direction 'dir', this function generates a direction 'result' that is a scaled
    // version of 'dir', with a possible sign change to make sure that the direction is 
    // "downhill" (negative scalar product with function gradient).  Returns 'true' on success,
    // and 'false' if the direction was practically perpendicular to the function gradient.
    bool orientDirection(const Point& dir, Point& result);

    // Returns -1 if the scalar product of p1 and p2 is negative, +1 if the scalar product
    // is positive and 0 if it close to 0.  In this case, 'close to 0' means that it is 
    // impossible to judge the derivatives sign when taking into account the numerical
    // noise associated with the magnitude of the points' components.
    int scalarProductSign(const Point& p1, const Point& p2);

    // Estimate the minimum of a parabola where two points and one tangent is known.  The 
    // two points are the current parameter point and the point obtained by marching from this
    // point a distance of 'dir' multiplied by 'p2'.  The function value in this second point
    // must also be given by 'fp2'.  The tangent known is the one for the current parameter 
    // point.  The abscissa for the minimum point of the parabola is returned.
    double parabolicEstimate(double p2, double fp2, const Point& dir);

    // Estimate the minimum of a parabola where three points are known.  If the minimum is found
    // to lie within the distance 'max_from_b' from b, the point is returned as 'u', and the
    // function returns 'true'.  Otherwise, it returns 'false', and 'u' is not calculated.
    static bool parabolicFit(double a, double fa, 
                             double b, double fb, 
                             double c, double fc, 
                             double max_from_b, double&u);

    // determine the smallest detectable change around the value 'c'
    static double numerical_tolerance(double c);

    // Determine the numerical tolerance around a given point in a given direction
    static double numerical_tolerance(const Point& x, const Point& dir);

    // Determine the numerical tolerance around a given point in a given direction, also
    // considering the change in function value (fx is the function at x, dfx is the derivative)
    static double numerical_tolerance(const Point& x, const Point& dir,
                                      const double fx, const double dfx);


    // A helper function for linminBrent.  A new point has been evaluated (u, fu), and we must
    // update the other points used by the algorithm.
    static void adjustBrackets(double& a, double& fa,
                               double& b, double& fb,
                               double& c, double& fc,
                               double& b2, double& fb2,
                               double& b3, double& fb3,
                               double u, double fu);

    static void shift3(double& a, double& b, double& c, double d) 
        { a = b; b = c; c = d;}
    static void shift2(double& a, double& b, double c) 
        { a = b; b = c;}
};

}; // end namespace Go

#include "GeneralFunctionMinimizer_implementation.h"

#endif // _GENERAL_FUNCTIONMINIMIZER_H

---