newton.h

Go to the documentation of this file.

#ifndef CHEMMISOL_NEWTON_H
#define CHEMMISOL_NEWTON_H
 
#include <functional>
#include <cmath>
#include "linear.h"
#include "gauss.h"
 
namespace chemmisol {
    
    template<typename T>
        T identity(const T& x) {
            return x;
        }
 
    template<typename X, typename M, X (&G)(const X&)=identity>
        class Newton {
            X x0;
            std::function<X(const X&)> f;
            std::function<M(const X&)> df;
 
            X solve_eps(const X& x, const X& f_x, float epsilon,
                    X x_min, X f_x_min, double n_f_x_min) const;
            X solve_iter(const X& x, const X& f_x, std::size_t n,
                    X x_min, X f_x_min, double n_f_x_min) const;
 
            public:
            Newton(
                    const X& x0,
                    std::function<X(const X&)> f,
                    std::function<M(const X&)> df)
                : x0(x0), f(f), df(df) {
                }
 
            X solve_eps(float epsilon) const;
 
            X solve_iter(std::size_t n) const;
        };
 
    template<typename X, typename M, X (&G)(const X&)>
        X Newton<X, M, G>::solve_eps(float epsilon) const {
            auto f_x0 = f(x0);
            return solve_eps(x0, f(x0), epsilon, x0, f_x0, norm(f_x0));
        }
 
    template<typename X, typename M, X (&G)(const X&)>
        X Newton<X, M, G>::solve_iter(std::size_t n) const {
            auto f_x0 = f(x0);
            return solve_iter(x0, f_x0, n, x0, f_x0, norm(f_x0));
        }
 
 
    template<typename X, typename M, X (&G)(const X&)>
        X Newton<X, M, G>::solve_eps(const X& x, const X& f_x, float epsilon,
                X x_min, X f_x_min, double n_f_x_min) const {
            CHEM_LOG(TRACE) << "[NEWTON] Current X: " << x;
            CHEM_LOG(TRACE) << "[NEWTON] Current F(X): " << f_x;
            CHEM_LOG(TRACE) << "[NEWTON] (e=" << epsilon << ") Epsilon: " << norm(f_x);
            if(norm(f_x) < epsilon)
                return x;
            X _x = gauss::solve(df(x), -f_x);
            X x1 = G(_x + x);
            auto n_f_x = norm(f_x);
            if (n_f_x < n_f_x_min)
                return solve_eps(x1, f(x1), epsilon, x1, f_x, n_f_x);
            if (n_f_x == n_f_x_min) {
                // Maybe we are in a cycle
                if(x_min == x1) {
                    // We are in a cycle, break the algorithm
                    CHEM_LOG(TRACE) << "[NEWTON] Cycle detected, returns current optimum (X: " << x_min
                        << ", f(X): " << f_x_min << ".";
                    return x_min;
                }
            }
            return solve_eps(x1, f(x1), epsilon, x_min, f_x_min, n_f_x_min);
        }
 
    template<typename X, typename M, X (&G)(const X&)>
        X Newton<X, M, G>::solve_iter(
                const X& x, const X& f_x, std::size_t n,
                X x_min, X f_x_min, double n_f_x_min) const {
            CHEM_LOG(TRACE) << "[NEWTON] Current X: " << x;
            CHEM_LOG(TRACE) << "[NEWTON] Current F(X): " << f_x;
            CHEM_LOG(TRACE) << "[NEWTON] (n=" << n << ") Epsilon: " << norm(f_x);
            if(n == 0 || norm(f_x) == typename X::value_type(0))
                return x;
            X _x = gauss::solve(df(x), -f_x);
            X x1 = G(_x + x);
            auto n_f_x = norm(f_x);
            if (n_f_x < n_f_x_min)
                return solve_iter(x1, f(x1), n-1, x1, f_x, n_f_x);
            if (n_f_x == n_f_x_min) {
                // Maybe we are in a cycle
                if(x_min == x1) {
                    // We are in a cycle, break the algorithm
                    CHEM_LOG(TRACE) << "[NEWTON] Cycle detected at n=" << n
                        << ", returns current optimum (X: " << x_min
                        << ", f(X): " << f_x_min << ".";
                    return x_min;
                }
            }
            return solve_iter(x1, f(x1), n-1, x_min, f_x_min, n_f_x_min);
        }
 
    template<typename X, typename M>
        using AbsoluteNewton = Newton<X, M, abs>;
}
#endif