polyeqv.h

Go to the documentation of this file.

#pragma once
 
#include <algorithm>
#include <iostream>
#include <cmath>
 
#include "../config.h"
#include "../poly.h"
#include "../util/misc.h"
#include "../util/spline.h"
 
namespace alfi::spline {
    template <typename Number = DefaultNumber, template <typename, typename...> class Container = DefaultContainer>
    class PolyEqvSpline {
    public:
        enum class PolynomialType {
            LAGRANGE,
            IMPROVED_LAGRANGE,
            NEWTON,
 
            Default = LAGRANGE,
        };
 
        enum class OptimizationType {
            ACCURACY,
            SPEED,
 
            Default = ACCURACY,
        };
 
        enum class EvaluationType {
            IGNORE_NANS_AND_PREVIOUS,
            IGNORE_NANS,
            NOT_IGNORE_NANS,
 
            Default = IGNORE_NANS_AND_PREVIOUS,
        };
 
        static Container<Number> compute_coeffs(
                const Container<Number>& X,
                const Container<Number>& Y,
                PolynomialType polynomial_type = PolynomialType::Default,
                OptimizationType optimization_type = OptimizationType::Default
        ) {
            if (X.size() != Y.size()) {
                std::cerr << "Error in function " << __FUNCTION__
                          << ": Vectors X (of size " << X.size()
                          << ") and Y (of size " << Y.size()
                          << ") are not the same size. Returning an empty array..." << std::endl;
                return {};
            }
 
            const SizeT n = X.size();
 
            if (n == 0) {
                return {};
            }
 
            if (n <= 4) {
                return util::spline::simple_spline<Number,Container>(X, Y, n - 1);
            }
 
            Container<Number> coeffs(n * (n - 1)); // n coefficients for n-1 segments
 
            switch (optimization_type) {
                case OptimizationType::ACCURACY: {
#if defined(_OPENMP) && !defined(ALFI_DISABLE_OPENMP)
#pragma omp parallel for
#endif
                    for (SizeT i = 0; i < n - 1; ++i) {
                        auto X_i = X;
                        for (auto& x : X_i) {
                            x -= X[i];
                        }
                        Container<Number> P_i;
                        switch (polynomial_type) {
                            case PolynomialType::LAGRANGE: P_i = poly::lagrange(X_i, Y); break;
                            case PolynomialType::IMPROVED_LAGRANGE: P_i = poly::imp_lagrange(X_i, Y); break;
                            case PolynomialType::NEWTON: P_i = poly::newton(X_i, Y); break;
                        }
                        std::move(P_i.begin(), P_i.end(), coeffs.begin() + (i * n));
                        coeffs[i*n+n-1] = Y[i];
                    }
                    break;
                }
                case OptimizationType::SPEED: {
                    Container<Number> P;
                    switch (polynomial_type) {
                        case PolynomialType::LAGRANGE: P = poly::lagrange(X, Y); break;
                        case PolynomialType::IMPROVED_LAGRANGE: P = poly::imp_lagrange(X, Y); break;
                        case PolynomialType::NEWTON: P = poly::newton(X, Y); break;
                    }
#if defined(_OPENMP) && !defined(ALFI_DISABLE_OPENMP)
#pragma omp parallel for
#endif
                    for (SizeT i = 0; i < n - 1; ++i) {
                        for (SizeT j = 0; j < n; ++j) {
                            coeffs[i*n+j] = P[j];
                        }
                        for (SizeT k = 0; k < n - 1; ++k) {
                            for (SizeT j = 1; j < n - k; ++j) {
                                coeffs[i*n+j] += coeffs[i*n+j-1] * X[i];
                            }
                        }
                        coeffs[i*n+n-1] = Y[i];
                    }
                    break;
                }
            }
            return coeffs;
        }
 
        PolyEqvSpline() = default;
 
        template <typename ContainerXType>
        PolyEqvSpline(ContainerXType&& X, const Container<Number>& Y, OptimizationType optimization_type)
            : PolyEqvSpline(std::forward<ContainerXType>(X), Y, PolynomialType::Default, optimization_type) {}
 
        template <typename ContainerXType>
        PolyEqvSpline(ContainerXType&& X,
                      const Container<Number>& Y,
                      PolynomialType polynomial_type = PolynomialType::Default,
                      OptimizationType optimization_type = OptimizationType::Default,
                      EvaluationType evaluation_type = EvaluationType::Default) {
            construct(std::forward<ContainerXType>(X), Y, polynomial_type, optimization_type, evaluation_type);
        }
 
        PolyEqvSpline(const PolyEqvSpline& other) = default;
        PolyEqvSpline(PolyEqvSpline&& other) noexcept = default;
 
        PolyEqvSpline& operator=(const PolyEqvSpline& other) = default;
        PolyEqvSpline& operator=(PolyEqvSpline&& other) noexcept = default;
 
        template <typename ContainerXType>
        void construct(ContainerXType&& X,
                       const Container<Number>& Y,
                       PolynomialType polynomial_type = PolynomialType::Default,
                       OptimizationType optimization_type = OptimizationType::Default,
                       EvaluationType evaluation_type = EvaluationType::Default) {
            auto coeffs = compute_coeffs(X, Y, polynomial_type, optimization_type);
            if (coeffs.empty() && !X.empty()) {
                std::cerr << "Error in function " << __FUNCTION__
                          << ": failed to construct coefficients. Not changing the object..." << std::endl;
                return;
            }
            _polynomial_type = polynomial_type;
            _optimization_type = optimization_type;
            _evaluation_type = evaluation_type;
            _X = std::forward<ContainerXType>(X);
            _coeffs = std::move(coeffs);
        }
 
        Number eval(const Number& x) const {
            return eval(x, std::distance(_X.begin(), util::misc::first_leq_or_begin(_X.begin(), _X.end(), x)));
        }
        Number eval(const Number& x, SizeT segment) const {
            if (_coeffs.empty()) {
                return NAN;
            } else if (_coeffs.size() == 1) {
                return _coeffs[0];
            }
 
            segment = std::clamp(segment, static_cast<SizeT>(0), static_cast<SizeT>(_X.size() - 2));
 
            const Number x_seg = x - _X[segment];
 
            const SizeT n = _X.size();
 
            Number result = _coeffs[segment*n];
 
            if (std::isnan(result)) {
                switch (_evaluation_type) {
                    case EvaluationType::IGNORE_NANS_AND_PREVIOUS:
                    case EvaluationType::IGNORE_NANS:
                        result = 0;
                        break;
                    case EvaluationType::NOT_IGNORE_NANS:
                        return result;
                }
            }
 
            for (SizeT i = 1; i < n; ++i) {
                result *= x_seg;
                Number current = _coeffs[segment*n+i];
                if (std::isnan(current)) {
                    switch (_evaluation_type) {
                        case EvaluationType::IGNORE_NANS_AND_PREVIOUS:
                            result = 0;
                            [[fallthrough]];
                        case EvaluationType::IGNORE_NANS:
                            current = 0;
                            break;
                        case EvaluationType::NOT_IGNORE_NANS:
                            return current;
                    }
                }
                result += current;
            }
 
            return result;
        }
 
        Container<Number> eval(const Container<Number>& xx, bool sorted = true) const {
            Container<Number> result(xx.size());
            if (sorted) {
                for (SizeT i = 0, i_x = 0; i < xx.size(); ++i) {
                    const Number& x = xx[i];
                    while (i_x + 1 < _X.size() && x >= _X[i_x+1])
                        ++i_x;
                    result[i] = eval(x, i_x);
                }
            } else {
                for (SizeT i = 0; i < xx.size(); ++i) {
                    result[i] = eval(xx[i]);
                }
            }
            return result;
        }
 
        Number operator()(const Number& x) const {
            return eval(x);
        }
        Container<Number> operator()(const Container<Number>& xx) const {
            return eval(xx);
        }
 
        PolynomialType polynomial_type() const {
            return _polynomial_type;
        }
 
        OptimizationType optimization_type() const {
            return _optimization_type;
        }
 
        EvaluationType evaluation_type() const {
            return _evaluation_type;
        }
 
        const Container<Number>& X() const & {
            return _X;
        }
        Container<Number>&& X() && {
            return std::move(_X);
        }
 
        const Container<Number>& coeffs() const & {
            return _coeffs;
        }
        Container<Number>&& coeffs() && {
            return std::move(_coeffs);
        }
 
        std::pair<SizeT, SizeT> segment(SizeT index) const {
            return {_X.size()*index, _X.size()*(index+1)};
        }
 
    private:
        PolynomialType _polynomial_type = PolynomialType::Default;
        OptimizationType _optimization_type = OptimizationType::Default;
        EvaluationType _evaluation_type = EvaluationType::Default;
        Container<Number> _X = {};
        Container<Number> _coeffs = {};
    };
}