QMinimizer.cc

Go to the documentation of this file.

#include "QMinimizer.hh"
#include <gsl/gsl_errno.h>
#include <gsl/gsl_math.h>
 
Q_BEGIN_NAMESPACE
 
QMathFunction::QRoutinePointer QMinimizer::fCurrentFunction; 
 
QMinimizer::QMinimizer()
{
    const gsl_min_fminimizer_type *T;
    T = gsl_min_fminimizer_brent;
    fMinimizer = gsl_min_fminimizer_alloc (T);
    fVerbose = false;
}
 
QMinimizer::~QMinimizer()
{
//    gsl_min_fminimizer_free (fMinimizer);
}
 
double QMinimizer::GetMinimum(QMathFunction func, double lowbound, double upbound, double epsabs, double epsrel)
{
    double min = Q_DOUBLE_DEFAULT;
    double startx = (upbound - lowbound)/2;
    int status;
    int iter = 0, max_iter = 100;
    gsl_function F;
    F = func.GetGSLFunction();
 
    gsl_min_fminimizer_set (fMinimizer, &F, startx, lowbound, upbound);
 
    if(fVerbose) {
        printf ("%5d [%.7f, %.7f] "
                "%.7f\n",
                iter, lowbound, upbound,
                startx);
    }
 
    do
    {
        iter++;
        status = gsl_min_fminimizer_iterate (fMinimizer);
 
        min = gsl_min_fminimizer_x_minimum (fMinimizer);
        lowbound = gsl_min_fminimizer_x_lower (fMinimizer);
        upbound = gsl_min_fminimizer_x_upper (fMinimizer);
 
        status = gsl_min_test_interval (lowbound,upbound, epsabs, epsrel);
 
        if(fVerbose) {
            if (status == GSL_SUCCESS) printf ("Converged:\n");
            printf ("%5d [%.7f, %.7f] "
                    "%.7f\n",
                    iter, lowbound, upbound,
                    min);
        }
    }
    while (status == GSL_CONTINUE && iter < max_iter);
    
    if(status != GSL_SUCCESS) min = Q_DOUBLE_DEFAULT;
 
    return min;
}
double QMinimizer::MinusFunc(double x, double* params) 
{
    return - (*fCurrentFunction)(x,params);
}
 
double QMinimizer::GetMaximum(QMathFunction func,double minx, double maxx, double epsabs, double epsrel)
{
    fCurrentFunction = func.fRoutine;
    QMathFunction minusFunc = func;
    minusFunc.fRoutine = &(MinusFunc);
   return  GetMinimum(minusFunc,minx,maxx,epsabs,epsrel);
}
 
Q_END_NAMESPACE