CoolProp/include/CPnumerics.h

#ifndef COOLPROP_NUMERICS_H
#define COOLPROP_NUMERICS_H

#include <vector>
#include <set>
#include <cfloat>
#include <stdlib.h>   // For abs
#include <algorithm>  // For max
#include <numeric>
#include <cmath>
#include <array>
#include "PlatformDetermination.h"
#include "CPstrings.h"
#include "Exceptions.h"

#if defined(HUGE_VAL) && !defined(_HUGE)
#    define _HUGE HUGE_VAL
#else
// GCC Version of huge value macro
#    if defined(HUGE) && !defined(_HUGE)
#        define _HUGE HUGE
#    endif
#endif

template <typename T, size_t N>
std::array<T, N> create_filled_array(T value) {
    std::array<T, N> arr;
    arr.fill(value);
    return arr;
}

inline bool ValidNumber(double x) {
    // Idea from http://www.johndcook.com/IEEE_exceptions_in_cpp.html
    return (x <= DBL_MAX && x >= -DBL_MAX);
};

#ifndef M_PI
#    define M_PI 3.14159265358979323846
#endif

#ifndef COOLPROP_OK
#    define COOLPROP_OK 1
#endif

// Undefine these terrible macros defined in windows header
#undef min
#undef max

/* "THE BEER-WARE LICENSE" (Revision 42): Devin Lane wrote this file. As long as you retain
 * this notice you can do whatever you want with this stuff. If we meet some day, and you
 * think this stuff is worth it, you can buy me a beer in return.
 *
 * From http://shiftedbits.org/2011/01/30/cubic-spline-interpolation/
 *
 * IHB(05/01/2016): Removed overload and renamed the interpolate function (cython cannot disambiguate the functions)
 *
 * Templated on type of X, Y. X and Y must have operator +, -, *, /. Y must have defined
 * a constructor that takes a scalar. */
template <typename X, typename Y>
class Spline
{
   public:
    /** An empty, invalid spline */
    Spline() {}

    /** A spline with x and y values */
    Spline(const std::vector<X>& x, const std::vector<Y>& y) {
        if (x.size() != y.size()) {
            std::cerr << "X and Y must be the same size " << std::endl;
            return;
        }

        if (x.size() < 3) {
            std::cerr << "Must have at least three points for interpolation" << std::endl;
            return;
        }

        typedef typename std::vector<X>::difference_type size_type;

        size_type n = y.size() - 1;

        std::vector<Y> b(n), d(n), a(n), c(n + 1), l(n + 1), u(n + 1), z(n + 1);
        std::vector<X> h(n + 1);

        l[0] = Y(1);
        u[0] = Y(0);
        z[0] = Y(0);
        h[0] = x[1] - x[0];

        for (size_type i = 1; i < n; i++) {
            h[i] = x[i + 1] - x[i];
            l[i] = Y(2 * (x[i + 1] - x[i - 1])) - Y(h[i - 1]) * u[i - 1];
            u[i] = Y(h[i]) / l[i];
            a[i] = (Y(3) / Y(h[i])) * (y[i + 1] - y[i]) - (Y(3) / Y(h[i - 1])) * (y[i] - y[i - 1]);
            z[i] = (a[i] - Y(h[i - 1]) * z[i - 1]) / l[i];
        }

        l[n] = Y(1);
        z[n] = c[n] = Y(0);

        for (size_type j = n - 1; j >= 0; j--) {
            c[j] = z[j] - u[j] * c[j + 1];
            b[j] = (y[j + 1] - y[j]) / Y(h[j]) - (Y(h[j]) * (c[j + 1] + Y(2) * c[j])) / Y(3);
            d[j] = (c[j + 1] - c[j]) / Y(3 * h[j]);
        }

        for (size_type i = 0; i < n; i++) {
            mElements.push_back(Element(x[i], y[i], b[i], c[i], d[i]));
        }
    }
    virtual ~Spline() {}

    Y operator[](const X& x) const {
        return interpolate(x);
    }

    Y interpolate(const X& x) const {
        if (mElements.size() == 0) return Y();

        typename std::vector<element_type>::const_iterator it;
        it = std::lower_bound(mElements.begin(), mElements.end(), element_type(x));
        if (it != mElements.begin()) {
            it--;
        }

        return it->eval(x);
    }

    /* Evaluate at multiple locations, assuming xx is sorted ascending */
    std::vector<Y> interpolate_vec(const std::vector<X>& xx) const {
        if (mElements.size() == 0) return std::vector<Y>(xx.size());

        typename std::vector<X>::const_iterator it;
        typename std::vector<element_type>::const_iterator it2;
        it2 = mElements.begin();
        std::vector<Y> ys;
        for (it = xx.begin(); it != xx.end(); it++) {
            it2 = std::lower_bound(it2, mElements.end(), element_type(*it));
            if (it2 != mElements.begin()) {
                it2--;
            }

            ys.push_back(it2->eval(*it));
        }

        return ys;
    }

   protected:
    class Element
    {
       public:
        Element(X _x) : x(_x) {}
        Element(X _x, Y _a, Y _b, Y _c, Y _d) : x(_x), a(_a), b(_b), c(_c), d(_d) {}

        Y eval(const X& xx) const {
            X xix(xx - x);
            return a + b * xix + c * (xix * xix) + d * (xix * xix * xix);
        }

        bool operator<(const Element& e) const {
            return x < e.x;
        }
        bool operator<(const X& xx) const {
            return x < xx;
        }

        X x;
        Y a, b, c, d;
    };

    typedef Element element_type;
    std::vector<element_type> mElements;
};

/// Return the maximum difference between elements in two vectors where comparing z1[i] and z2[i]
template <typename T>
T maxvectordiff(const std::vector<T>& z1, const std::vector<T>& z2) {
    T maxvecdiff = 0;
    for (std::size_t i = 0; i < z1.size(); ++i) {
        T diff = std::abs(z1[i] - z2[i]);
        if (std::abs(diff) > maxvecdiff) {
            maxvecdiff = diff;
        }
    }
    return maxvecdiff;
}

/// Make a linearly spaced vector of points
template <typename T>
std::vector<T> linspace(T xmin, T xmax, std::size_t n) {
    std::vector<T> x(n, 0.0);

    for (std::size_t i = 0; i < n; ++i) {
        x[i] = (xmax - xmin) / (n - 1) * i + xmin;
    }
    return x;
}
/// Make a base-10 logarithmically spaced vector of points
template <typename T>
std::vector<T> log10space(T xmin, T xmax, std::size_t n) {
    std::vector<T> x(n, 0.0);
    T logxmin = log10(xmin), logxmax = log10(xmax);

    for (std::size_t i = 0; i < n; ++i) {
        x[i] = exp((logxmax - logxmin) / (n - 1) * i + logxmin);
    }
    return x;
}
/// Make a base-e logarithmically spaced vector of points
template <typename T>
std::vector<T> logspace(T xmin, T xmax, std::size_t n) {
    std::vector<T> x(n, 0.0);
    T logxmin = log(xmin), logxmax = log(xmax);

    for (std::size_t i = 0; i < n; ++i) {
        x[i] = exp((logxmax - logxmin) / (n - 1) * i + logxmin);
    }
    return x;
}

/**
 * @brief Use bisection to find the inputs that bisect the value you want, the trick
 * here is that this function is allowed to have "holes" where parts of the the array are
 * also filled with invalid numbers for which ValidNumber(x) is false
 * @param vec The vector to be bisected
 * @param val The value to be found
 * @param i The index to the left of the final point; i and i+1 bound the value
 */
template <typename T>
void bisect_vector(const std::vector<T>& vec, T val, std::size_t& i) {
    T rL, rM, rR;
    std::size_t N = vec.size(), L = 0, R = N - 1, M = (L + R) / 2;
    // Move the right limits in until they are good
    while (!ValidNumber(vec[R])) {
        if (R == 1) {
            throw CoolProp::ValueError("All the values in bisection vector are invalid");
        }
        R--;
    }
    // Move the left limits in until they are good
    while (!ValidNumber(vec[L])) {
        if (L == vec.size() - 1) {
            throw CoolProp::ValueError("All the values in bisection vector are invalid");
        }
        L++;
    }
    rL = vec[L] - val;
    rR = vec[R] - val;
    while (R - L > 1) {
        if (!ValidNumber(vec[M])) {
            std::size_t MR = M, ML = M;
            // Move middle-right to the right until it is ok
            while (!ValidNumber(vec[MR])) {
                if (MR == vec.size() - 1) {
                    throw CoolProp::ValueError("All the values in bisection vector are invalid");
                }
                MR++;
            }
            // Move middle-left to the left until it is ok
            while (!ValidNumber(vec[ML])) {
                if (ML == 1) {
                    throw CoolProp::ValueError("All the values in bisection vector are invalid");
                }
                ML--;
            }
            T rML = vec[ML] - val;
            T rMR = vec[MR] - val;
            // Figure out which chunk is the good part
            if (rR * rML > 0 && rL * rML < 0) {
                // solution is between L and ML
                R = ML;
                rR = vec[ML] - val;
            } else if (rR * rMR < 0 && rL * rMR > 0) {
                // solution is between R and MR
                L = MR;
                rL = vec[MR] - val;
            } else {
                throw CoolProp::ValueError(
                  format("Unable to bisect segmented vector; neither chunk contains the solution val:%g left:(%g,%g) right:(%g,%g)", val, vec[L],
                         vec[ML], vec[MR], vec[R]));
            }
            M = (L + R) / 2;
        } else {
            rM = vec[M] - val;
            if (rR * rM > 0 && rL * rM < 0) {
                // solution is between L and M
                R = M;
                rR = vec[R] - val;
            } else {
                // solution is between R and M
                L = M;
                rL = vec[L] - val;
            }
            M = (L + R) / 2;
        }
    }
    i = L;
}

/**
 * @brief Use bisection to find the inputs that bisect the value you want, the trick
 * here is that this function is allowed to have "holes" where parts of the the array are
 * also filled with invalid numbers for which ValidNumber(x) is false
 * @param matrix The vector to be bisected
 * @param j The index of the matric in the off-grain dimension
 * @param val The value to be found
 * @param i The index to the left of the final point; i and i+1 bound the value
 */
template <typename T>
void bisect_segmented_vector_slice(const std::vector<std::vector<T>>& mat, std::size_t j, T val, std::size_t& i) {
    T rL, rM, rR;
    std::size_t N = mat[j].size(), L = 0, R = N - 1, M = (L + R) / 2;
    // Move the right limits in until they are good
    while (!ValidNumber(mat[R][j])) {
        if (R == 1) {
            throw CoolProp::ValueError("All the values in bisection vector are invalid");
        }
        R--;
    }
    rR = mat[R][j] - val;
    // Move the left limits in until they are good
    while (!ValidNumber(mat[L][j])) {
        if (L == mat.size() - 1) {
            throw CoolProp::ValueError("All the values in bisection vector are invalid");
        }
        L++;
    }
    rL = mat[L][j] - val;
    while (R - L > 1) {
        if (!ValidNumber(mat[M][j])) {
            std::size_t MR = M, ML = M;
            // Move middle-right to the right until it is ok
            while (!ValidNumber(mat[MR][j])) {
                if (MR == mat.size() - 1) {
                    throw CoolProp::ValueError("All the values in bisection vector are invalid");
                }
                MR++;
            }
            // Move middle-left to the left until it is ok
            while (!ValidNumber(mat[ML][j])) {
                if (ML == 1) {
                    throw CoolProp::ValueError("All the values in bisection vector are invalid");
                }
                ML--;
            }
            T rML = mat[ML][j] - val;
            T rMR = mat[MR][j] - val;
            // Figure out which chunk is the good part
            if (rR * rMR > 0 && rL * rML < 0) {
                // solution is between L and ML
                R = ML;
                rR = mat[ML][j] - val;
            } else if (rR * rMR < 0 && rL * rML > 0) {
                // solution is between R and MR
                L = MR;
                rL = mat[MR][j] - val;
            } else {
                throw CoolProp::ValueError(
                  format("Unable to bisect segmented vector slice; neither chunk contains the solution %g lef:(%g,%g) right:(%g,%g)", val, mat[L][j],
                         mat[ML][j], mat[MR][j], mat[R][j]));
            }
            M = (L + R) / 2;
        } else {
            rM = mat[M][j] - val;
            if (rR * rM > 0 && rL * rM < 0) {
                // solution is between L and M
                R = M;
                rR = mat[R][j] - val;
            } else {
                // solution is between R and M
                L = M;
                rL = mat[L][j] - val;
            }
            M = (L + R) / 2;
        }
    }
    i = L;
}

// From http://rosettacode.org/wiki/Power_set#C.2B.2B
inline std::size_t powerset_dereference(std::set<std::size_t>::const_iterator v) {
    return *v;
};

// From http://rosettacode.org/wiki/Power_set#C.2B.2B
inline std::set<std::set<std::size_t>> powerset(std::set<std::size_t> const& set) {
    std::set<std::set<std::size_t>> result;
    std::vector<std::set<std::size_t>::const_iterator> elements;
    do {
        std::set<std::size_t> tmp;
        std::transform(elements.begin(), elements.end(), std::inserter(tmp, tmp.end()), powerset_dereference);
        result.insert(tmp);
        if (!elements.empty() && ++elements.back() == set.end()) {
            elements.pop_back();
        } else {
            std::set<std::size_t>::const_iterator iter;
            if (elements.empty()) {
                iter = set.begin();
            } else {
                iter = elements.back();
                ++iter;
            }
            for (; iter != set.end(); ++iter) {
                elements.push_back(iter);
            }
        }
    } while (!elements.empty());

    return result;
}

/// Some functions related to testing and comparison of values
bool inline check_abs(double A, double B, double D) {
    double max = std::abs(A);
    double min = std::abs(B);
    if (min > max) {
        max = min;
        min = std::abs(A);
    }
    if (max > DBL_EPSILON * 1e3)
        return ((1.0 - min / max * 1e0) < D);
    else
        throw CoolProp::ValueError(
          format("Too small numbers: %f cannot be tested with an accepted error of %f for a machine precision of %f. ", max, D, DBL_EPSILON));
};
bool inline check_abs(double A, double B) {
    return check_abs(A, B, 1e5 * DBL_EPSILON);
};

template <class T>
void normalize_vector(std::vector<T>& x) {
    // Sum up all the elements in the vector
    T sumx = std::accumulate(x.begin(), x.end(), static_cast<T>(0));
    // Normalize the components by dividing each by the sum
    for (std::size_t i = 0; i < x.size(); ++i) {
        x[i] /= sumx;
    }
};

/// A spline is a curve given by the form y = ax^3 + bx^2 + c*x + d
/// As there are 4 constants, 4 constraints are needed to create the spline.  These constraints could be the derivative or value at a point
/// Often, the value and derivative of the value are known at two points.
class SplineClass
{
   protected:
    int Nconstraints;
    std::vector<std::vector<double>> A;
    std::vector<double> B;

   public:
    double a, b, c, d;
    SplineClass() : Nconstraints(0), A(4, std::vector<double>(4, 0)), B(4, 0), a(_HUGE), b(_HUGE), c(_HUGE), d(_HUGE) {}
    bool build(void);
    bool add_value_constraint(double x, double y);
    void add_4value_constraints(double x1, double x2, double x3, double x4, double y1, double y2, double y3, double y4);
    bool add_derivative_constraint(double x, double dydx);
    double evaluate(double x);
};

/// from http://stackoverflow.com/a/5721830/1360263
template <class T>
T factorial(T n) {
    if (n == 0) return 1;
    return n * factorial(n - 1);
}
/// see https://proofwiki.org/wiki/Nth_Derivative_of_Mth_Power
/// and https://proofwiki.org/wiki/Definition:Falling_Factorial
template <class T1, class T2>
T1 nth_derivative_of_x_to_m(T1 x, T2 n, T2 m) {
    if (n > m) {
        return 0;
    } else {
        return factorial(m) / factorial(m - n) * pow(x, m - n);
    }
}

void MatInv_2(double A[2][2], double B[2][2]);

double root_sum_square(const std::vector<double>& x);
double interp1d(const std::vector<double>* x, const std::vector<double>* y, double x0);
double powInt(double x, int y);

template <class T>
T POW2(T x) {
    return x * x;
}
template <class T>
T POW3(T x) {
    return POW2(x) * x;
}
template <class T>
T POW4(T x) {
    return POW2(x) * POW2(x);
}
#define POW5(x) ((x) * (x) * (x) * (x) * (x))
#define POW6(x) ((x) * (x) * (x) * (x) * (x) * (x))
#define POW7(x) ((x) * (x) * (x) * (x) * (x) * (x) * (x))

template <class T>
T LinearInterp(T x0, T x1, T y0, T y1, T x) {
    return (y1 - y0) / (x1 - x0) * (x - x0) + y0;
};
template <class T1, class T2>
T2 LinearInterp(const std::vector<T1>& x, const std::vector<T1>& y, std::size_t i0, std::size_t i1, T2 val) {
    return LinearInterp(x[i0], x[i1], y[i0], y[i1], static_cast<T1>(val));
};

template <class T>
T QuadInterp(T x0, T x1, T x2, T f0, T f1, T f2, T x) {
    /* Quadratic interpolation.
    Based on method from Kreyszig,
    Advanced Engineering Mathematics, 9th Edition
    */
    T L0, L1, L2;
    L0 = ((x - x1) * (x - x2)) / ((x0 - x1) * (x0 - x2));
    L1 = ((x - x0) * (x - x2)) / ((x1 - x0) * (x1 - x2));
    L2 = ((x - x0) * (x - x1)) / ((x2 - x0) * (x2 - x1));
    return L0 * f0 + L1 * f1 + L2 * f2;
};
template <class T1, class T2>
T2 QuadInterp(const std::vector<T1>& x, const std::vector<T1>& y, std::size_t i0, std::size_t i1, std::size_t i2, T2 val) {
    return QuadInterp(x[i0], x[i1], x[i2], y[i0], y[i1], y[i2], static_cast<T1>(val));
};

template <class T>
T CubicInterp(T x0, T x1, T x2, T x3, T f0, T f1, T f2, T f3, T x) {
    /*
    Lagrange cubic interpolation as from
    http://nd.edu/~jjwteach/441/PdfNotes/lecture6.pdf
    */
    T L0, L1, L2, L3;
    L0 = ((x - x1) * (x - x2) * (x - x3)) / ((x0 - x1) * (x0 - x2) * (x0 - x3));
    L1 = ((x - x0) * (x - x2) * (x - x3)) / ((x1 - x0) * (x1 - x2) * (x1 - x3));
    L2 = ((x - x0) * (x - x1) * (x - x3)) / ((x2 - x0) * (x2 - x1) * (x2 - x3));
    L3 = ((x - x0) * (x - x1) * (x - x2)) / ((x3 - x0) * (x3 - x1) * (x3 - x2));
    return L0 * f0 + L1 * f1 + L2 * f2 + L3 * f3;
};
/** /brief Calculate the first derivative of the function using a cubic interpolation form
    */
template <class T>
T CubicInterpFirstDeriv(T x0, T x1, T x2, T x3, T f0, T f1, T f2, T f3, T x) {
    // Based on http://math.stackexchange.com/a/809946/66405
    T L0 = ((x - x1) * (x - x2) * (x - x3)) / ((x0 - x1) * (x0 - x2) * (x0 - x3));
    T dL0_dx = (1 / (x - x1) + 1 / (x - x2) + 1 / (x - x3)) * L0;
    T L1 = ((x - x0) * (x - x2) * (x - x3)) / ((x1 - x0) * (x1 - x2) * (x1 - x3));
    T dL1_dx = (1 / (x - x0) + 1 / (x - x2) + 1 / (x - x3)) * L1;
    T L2 = ((x - x0) * (x - x1) * (x - x3)) / ((x2 - x0) * (x2 - x1) * (x2 - x3));
    T dL2_dx = (1 / (x - x0) + 1 / (x - x1) + 1 / (x - x3)) * L2;
    T L3 = ((x - x0) * (x - x1) * (x - x2)) / ((x3 - x0) * (x3 - x1) * (x3 - x2));
    T dL3_dx = (1 / (x - x0) + 1 / (x - x1) + 1 / (x - x2)) * L3;
    return dL0_dx * f0 + dL1_dx * f1 + dL2_dx * f2 + dL3_dx * f3;
};
template <class T1, class T2>
T2 CubicInterp(const std::vector<T1>& x, const std::vector<T1>& y, std::size_t i0, std::size_t i1, std::size_t i2, std::size_t i3, T2 val) {
    return CubicInterp(x[i0], x[i1], x[i2], x[i3], y[i0], y[i1], y[i2], y[i3], static_cast<T1>(val));
};

template <class T>
T is_in_closed_range(T x1, T x2, T x) {
    return (x >= std::min(x1, x2) && x <= std::max(x1, x2));
};

/** \brief Solve a cubic with coefficients in decreasing order
    *
    * 0 = ax^3 + b*x^2 + c*x + d
    *
    * @param a The x^3 coefficient
    * @param b The x^2 coefficient
    * @param c The x^1 coefficient
    * @param d The x^0 coefficient
    * @param N The number of unique real solutions found
    * @param x0 The first solution found
    * @param x1 The second solution found
    * @param x2 The third solution found
    */
void solve_cubic(double a, double b, double c, double d, int& N, double& x0, double& x1, double& x2);

void solve_quartic(double a, double b, double c, double d, double e, int& N, double& x0, double& x1, double& x2, double& x3);

template <class T>
inline T min3(T x1, T x2, T x3) {
    return std::min(std::min(x1, x2), x3);
};
template <class T>
inline T max3(T x1, T x2, T x3) {
    return std::max(std::max(x1, x2), x3);
};
template <class T>
inline T min4(T x1, T x2, T x3, T x4) {
    return std::min(std::min(std::min(x1, x2), x3), x4);
};
template <class T>
inline T max4(T x1, T x2, T x3, T x4) {
    return std::max(std::max(std::max(x1, x2), x3), x4);
};

inline bool double_equal(double a, double b) {
    return std::abs(a - b) <= 1 * DBL_EPSILON * std::max(std::abs(a), std::abs(b));
};

template <class T>
T max_abs_value(const std::vector<T>& x) {
    T max = 0;
    std::size_t N = x.size();
    for (std::size_t i = 0; i < N; ++i) {
        T axi = std::abs(x[i]);
        if (axi > max) {
            max = axi;
        }
    }
    return max;
}

template <class T>
T min_abs_value(const std::vector<T>& x) {
    T min = 1e40;
    std::size_t N = x.size();
    for (std::size_t i = 0; i < N; ++i) {
        T axi = std::abs(x[i]);
        if (axi < min) {
            min = axi;
        }
    }
    return min;
}

inline int Kronecker_delta(std::size_t i, std::size_t j) {
    if (i == j) {
        return static_cast<int>(1);
    } else {
        return static_cast<int>(0);
    }
};
inline int Kronecker_delta(int i, int j) {
    if (i == j) {
        return 1;
    } else {
        return 0;
    }
};

/// Sort three values in place; see http://codereview.stackexchange.com/a/64763
template <typename T>
void sort3(T& a, T& b, T& c) {
    if (a > b) {
        std::swap(a, b);
    }
    if (a > c) {
        std::swap(a, c);
    }
    if (b > c) {
        std::swap(b, c);
    }
}

/**
* Due to the periodicity of angles, you need to handle the case where the
* angles wrap around - suppose theta_d is 6.28 and you are at an angles of 0.1 rad,
* the difference should be around 0.1, not -6.27
*
* This brilliant method is from http://blog.lexique-du-net.com/index.php?post/Calculate-the-real-difference-between-two-angles-keeping-the-sign
* and the comment of user tk
*
* Originally implemented in PDSim
*/
template <class T>
T angle_difference(T angle1, T angle2) {
    return fmod(angle1 - angle2 + M_PI, 2 * M_PI) - M_PI;
}

/// A simple function for use in wrappers where macros cause problems
inline double get_HUGE() {
    return _HUGE;
}

#if defined(_MSC_VER)
// Microsoft version of math.h doesn't include acosh or asinh, so we just define them here.
// It was included from Visual Studio 2013
#    if _MSC_VER < 1800
double acosh(double x) {
    return log(x + sqrt(x * x - 1.0));
}
double asinh(double value) {
    if (value > 0) {
        return log(value + sqrt(value * value + 1));
    } else {
        return -log(-value + sqrt(value * value + 1));
    }
}
#    endif
#endif

#if defined(__ISPOWERPC__)
// PPC version of math.h doesn't include acosh or asinh, so we just define them here
double acosh(double x) {
    return log(x + sqrt(x * x - 1.0));
}
double asinh(double value) {
    if (value > 0) {
        return log(value + sqrt(value * value + 1));
    } else {
        return -log(-value + sqrt(value * value + 1));
    }
}
#endif

#if defined(__ISPOWERPC__)
#    undef EOS
#endif

#endif