extrema.h

Go to the documentation of this file.

 
 
#ifndef THEORETICA_EXTREMA_H
#define THEORETICA_EXTREMA_H
 
#include "../core/constants.h"
#include "./roots.h"
 
 
namespace theoretica {
 
 
    template<typename RealFunction>
    inline iter_result<real> maximize_golden(
        RealFunction f, real a, real b,
        real tolerance = OPTIMIZATION_TOL,
        unsigned int max_iter = OPTIMIZATION_GOLDENSECTION_ITER) {
 
        if(a > b) {
            TH_MATH_ERROR("maximize_golden", b, MathError::InvalidArgument);
            return iter_result<real>(ConvergenceStatus::InvalidInput);
        }
 
        real x1 = a;
        real x2 = b;
        real x3 = b - (b - a) / PHI;
        real x4 = a + (b - a) / PHI;
 
        unsigned int iter = 0;
 
        // Keep iterating until the bracketing interval is closer than the tolerance,
        // or until the maximum number of iterations is reached.
        while(abs(x2 - x1) > tolerance && iter <= max_iter) {
 
            if(f(x3) > f(x4)) {
                x2 = x4;
            } else {
                x1 = x3;
            }
 
            x3 = x2 - (x2 - x1) / PHI;
            x4 = x1 + (x2 - x1) / PHI;
 
            iter++;
        }
 
        if(iter > max_iter) {
            TH_MATH_ERROR("maximize_golden", iter, MathError::NoConvergence);
            return iter_result<real>(ConvergenceStatus::MaxIterations, iter, abs(x2 - x1) / 2.0);
        }
 
        return iter_result<real>((x2 + x1) / 2.0, iter, abs(x2 - x1) / 2.0);
    }
 
 
    template<typename RealFunction>
    inline iter_result<real> minimize_golden(
        RealFunction f, real a, real b,
        real tolerance = OPTIMIZATION_TOL,
        unsigned int max_iter = OPTIMIZATION_GOLDENSECTION_ITER) {
 
        if(a > b) {
            TH_MATH_ERROR("minimize_golden", b, MathError::InvalidArgument);
            return iter_result<real>(ConvergenceStatus::InvalidInput);
        }
 
        real x1 = a;
        real x2 = b;
        real x3 = b - (b - a) / PHI;
        real x4 = a + (b - a) / PHI;
 
        unsigned int iter = 0;
 
        while(abs(x2 - x1) > tolerance && iter <= max_iter) {
 
            if(f(x3) < f(x4)) {
                x2 = x4;
            } else {
                x1 = x3;
            }
 
            x3 = x2 - (x2 - x1) / PHI;
            x4 = x1 + (x2 - x1) / PHI;
 
            iter++;
        }
 
        if(iter > max_iter) {
            TH_MATH_ERROR("minimize_golden", iter, MathError::NoConvergence);
            return iter_result<real>(ConvergenceStatus::MaxIterations, iter, abs(x2 - x1) / 2.0);
        }
 
        return iter_result<real>((x2 + x1) / 2.0, iter, abs(x2 - x1) / 2.0);
    }
 
 
    template<typename RealFunction>
    inline iter_result<real> maximize_newton(
        RealFunction f, RealFunction Df, RealFunction D2f,
        real guess = 0.0, real tolerance = OPTIMIZATION_TOL,
        unsigned int max_iter = OPTIMIZATION_NEWTON_ITER) {
 
        iter_result<real> z = root_newton(Df, D2f, guess, tolerance, max_iter);
 
        if(D2f(z) > 0) {
            
            TH_MATH_ERROR("maximize_newton", z, MathError::NoConvergence);
 
            // If the root finding algorithm converged but found a solution
            // with the wrong concavity (minimum instead of maximum),
            // mark the status as diverged.
            return iter_result<real>(
                (z.status == ConvergenceStatus::Success) ? ConvergenceStatus::Diverged : z.status,
                z.iterations,
                z.residual
            );
        }
 
        return iter_result<real>(z, z.iterations, z.residual);
    }
 
 
    template<typename RealFunction>
    inline iter_result<real> minimize_newton(
        RealFunction f, RealFunction Df,
        RealFunction D2f, real guess = 0,
        real tolerance = OPTIMIZATION_TOL,
        unsigned int max_iter = OPTIMIZATION_NEWTON_ITER) {
 
        iter_result<real> z = root_newton(Df, D2f, guess, tolerance, max_iter);
 
        if(D2f(z) < 0) {
 
            TH_MATH_ERROR("minimize_newton", z, MathError::NoConvergence);
 
            // If the root finding algorithm converged but found a solution
            // with the wrong concavity (maximum instead of minimum),
            // mark the status as diverged.
            return iter_result<real>(
                (z.status == ConvergenceStatus::Success) ? ConvergenceStatus::Diverged : z.status,
                z.iterations,
                z.residual
            );
        }
 
        return iter_result<real>(z, z.iterations, z.residual);
    }
 
 
    template<typename RealFunction>
    inline iter_result<real> maximize_bisect(
        RealFunction f, RealFunction Df,
        real a, real b,
        real tolerance = OPTIMIZATION_TOL,
        unsigned int max_iter = OPTIMIZATION_BISECTION_ITER) {
 
        iter_result<real> z = root_bisect(Df, a, b, tolerance, max_iter);
 
        // Approximate the function concavity
        if(deriv_central(Df, z) > 0) {
 
            TH_MATH_ERROR("maximize_bisect", z, MathError::NoConvergence);
 
            // If the root finding algorithm converged but found a solution
            // with the wrong concavity (minimum instead of maximum),
            // mark the status as diverged.
            return iter_result<real>(
                (z.status == ConvergenceStatus::Success) ? ConvergenceStatus::Diverged : z.status,
                z.iterations,
                z.residual
            );
        }
 
        return iter_result<real>(z, z.iterations, z.residual);
    }
 
 
    template<typename RealFunction>
    inline iter_result<real> minimize_bisect(
        RealFunction f, RealFunction Df, real a, real b,
        real tolerance = OPTIMIZATION_TOL,
        unsigned int max_iter = OPTIMIZATION_BISECTION_ITER) {
 
        iter_result<real> z = root_bisect(Df, a, b, tolerance, max_iter);
 
        // Approximate the function concavity
        if(deriv_central(Df, z) < 0) {
 
            TH_MATH_ERROR("minimize_bisect", z, MathError::NoConvergence);
 
            // If the root finding algorithm converged but found a solution
            // with the wrong concavity (maximum instead of minimum),
            // mark the status as diverged.
            return iter_result<real>(
                (z.status == ConvergenceStatus::Success) ? ConvergenceStatus::Diverged : z.status,
                z.iterations,
                z.residual
            );
        }
 
        return iter_result<real>(z, z.iterations, z.residual);
    }
 
}
 
#endif