Spline.hpp

Go to the documentation of this file.

// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*
  This file is part of the Open Porous Media project (OPM).

  OPM is free software: you can redistribute it and/or modify
  it under the terms of the GNU General Public License as published by
  the Free Software Foundation, either version 2 of the License, or
  (at your option) any later version.

  OPM is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  GNU General Public License for more details.

  You should have received a copy of the GNU General Public License
  along with OPM.  If not, see <http://www.gnu.org/licenses/>.

  Consult the COPYING file in the top-level source directory of this
  module for the precise wording of the license and the list of
  copyright holders.
*/
#ifndef OPM_SPLINE_HPP
#define OPM_SPLINE_HPP

#include <opm/common/Exceptions.hpp>

#include <opm/material/common/TridiagonalMatrix.hpp>
#include <opm/material/common/PolynomialUtils.hpp>

#include <cstddef>
#include <iosfwd>
#include <vector>

namespace Opm
{
template<class Scalar>
class Spline
{
    typedef TridiagonalMatrix<Scalar> Matrix;
    typedef std::vector<Scalar> Vector;

public:
    enum SplineType {
        Full,
        Natural,
        Periodic,
        Monotonic
    };

    Spline()
    { }

    Spline(Scalar x0, Scalar x1,
           Scalar y0, Scalar y1,
           Scalar m0, Scalar m1)
    { set(x0, x1, y0, y1, m0, m1); }

    template <class ScalarArrayX, class ScalarArrayY>
    Spline(std::size_t nSamples,
           const ScalarArrayX& x,
           const ScalarArrayY& y,
           SplineType splineType = Natural,
           bool sortInputs = true)
    { this->setXYArrays(nSamples, x, y, splineType, sortInputs); }

    template <class PointArray>
    Spline(std::size_t nSamples,
           const PointArray& points,
           SplineType splineType = Natural,
           bool sortInputs = true)
    { this->setArrayOfPoints(nSamples, points, splineType, sortInputs); }

    template <class ScalarContainer>
    Spline(const ScalarContainer& x,
           const ScalarContainer& y,
           SplineType splineType = Natural,
           bool sortInputs = true)
    { this->setXYContainers(x, y, splineType, sortInputs); }

    template <class PointContainer>
    explicit Spline(const PointContainer& points,
                    SplineType splineType = Natural,
                    bool sortInputs = true)
    { this->setContainerOfPoints(points, splineType, sortInputs); }

    template <class ScalarArray>
    Spline(std::size_t nSamples,
           const ScalarArray& x,
           const ScalarArray& y,
           Scalar m0,
           Scalar m1,
           bool sortInputs = true)
    { this->setXYArrays(nSamples, x, y, m0, m1, sortInputs); }

    template <class PointArray>
    Spline(std::size_t nSamples,
           const PointArray& points,
           Scalar m0,
           Scalar m1,
           bool sortInputs = true)
    { this->setArrayOfPoints(nSamples, points, m0, m1, sortInputs); }

    template <class ScalarContainerX, class ScalarContainerY>
    Spline(const ScalarContainerX& x,
           const ScalarContainerY& y,
           Scalar m0,
           Scalar m1,
           bool sortInputs = true)
    { this->setXYContainers(x, y, m0, m1, sortInputs); }

    template <class PointContainer>
    Spline(const PointContainer& points,
           Scalar m0,
           Scalar m1,
           bool sortInputs = true)
    { this->setContainerOfPoints(points, m0, m1, sortInputs); }

    void set(Scalar x0, Scalar x1,
             Scalar y0, Scalar y1,
             Scalar m0, Scalar m1)
    {
        slopeVec_.resize(2);
        xPos_.resize(2);
        yPos_.resize(2);

        if (x0 > x1) {
            xPos_[0] = x1;
            xPos_[1] = x0;
            yPos_[0] = y1;
            yPos_[1] = y0;
        }
        else {
            xPos_[0] = x0;
            xPos_[1] = x1;
            yPos_[0] = y0;
            yPos_[1] = y1;
        }

        slopeVec_[0] = m0;
        slopeVec_[1] = m1;

        Matrix M(numSamples());
        Vector d(numSamples());
        Vector moments(numSamples());
        this->makeFullSystem_(M, d, m0, m1);

        // solve for the moments
        M.solve(moments, d);

        this->setSlopesFromMoments_(slopeVec_, moments);
    }

    // Full splines                      //

    template <class ScalarArrayX, class ScalarArrayY>
    void setXYArrays(std::size_t nSamples,
                     const ScalarArrayX& x,
                     const ScalarArrayY& y,
                     Scalar m0, Scalar m1,
                     bool sortInputs = true)
    {
        assert(nSamples > 1);

        setNumSamples_(nSamples);
        for (std::size_t i = 0; i < nSamples; ++i) {
            xPos_[i] = x[i];
            yPos_[i] = y[i];
        }

        if (sortInputs)
            sortInput_();
        else if (xPos_[0] > xPos_[numSamples() - 1])
            reverseSamplingPoints_();

        makeFullSpline_(m0, m1);
    }

    template <class ScalarContainerX, class ScalarContainerY>
    void setXYContainers(const ScalarContainerX& x,
                         const ScalarContainerY& y,
                         Scalar m0, Scalar m1,
                         bool sortInputs = true)
    {
        assert(x.size() == y.size());
        assert(x.size() > 1);

        setNumSamples_(x.size());

        std::ranges::copy(x, xPos_.begin());
        std::ranges::copy(y, yPos_.begin());

        if (sortInputs)
            sortInput_();
        else if (xPos_[0] > xPos_[numSamples() - 1])
            reverseSamplingPoints_();

        makeFullSpline_(m0, m1);
    }

    template <class PointArray>
    void setArrayOfPoints(std::size_t nSamples,
                          const PointArray& points,
                          Scalar m0,
                          Scalar m1,
                          bool sortInputs = true)
    {
        // a spline with no or just one sampling points? what an
        // incredible bad idea!
        assert(nSamples > 1);

        setNumSamples_(nSamples);
        for (std::size_t i = 0; i < nSamples; ++i) {
            xPos_[i] = points[i][0];
            yPos_[i] = points[i][1];
        }

        if (sortInputs)
            sortInput_();
        else if (xPos_[0] > xPos_[numSamples() - 1])
            reverseSamplingPoints_();

        makeFullSpline_(m0, m1);
    }

    template <class XYContainer>
    void setContainerOfPoints(const XYContainer& points,
                              Scalar m0,
                              Scalar m1,
                              bool sortInputs = true)
    {
        // a spline with no or just one sampling points? what an
        // incredible bad idea!
        assert(points.size() > 1);

        setNumSamples_(points.size());
        typename XYContainer::const_iterator it = points.begin();
        typename XYContainer::const_iterator endIt = points.end();
        for (std::size_t i = 0; it != endIt; ++i, ++it) {
            xPos_[i] = (*it)[0];
            yPos_[i] = (*it)[1];
        }

        if (sortInputs)
            sortInput_();
        else if (xPos_[0] > xPos_[numSamples() - 1])
            reverseSamplingPoints_();

        // make a full spline
        makeFullSpline_(m0, m1);
    }

    template <class XYContainer>
    void setContainerOfTuples(const XYContainer& points,
                              Scalar m0,
                              Scalar m1,
                              bool sortInputs = true)
    {
        // resize internal arrays
        setNumSamples_(points.size());
        typename XYContainer::const_iterator it = points.begin();
        typename XYContainer::const_iterator endIt = points.end();
        for (unsigned i = 0; it != endIt; ++i, ++it) {
            xPos_[i] = std::get<0>(*it);
            yPos_[i] = std::get<1>(*it);
        }

        if (sortInputs)
            sortInput_();
        else if (xPos_[0] > xPos_[numSamples() - 1])
            reverseSamplingPoints_();

        // make a full spline
        makeFullSpline_(m0, m1);
    }

    // Natural/Periodic splines          //

    template <class ScalarArrayX, class ScalarArrayY>
    void setXYArrays(std::size_t nSamples,
                     const ScalarArrayX& x,
                     const ScalarArrayY& y,
                     SplineType splineType = Natural,
                     bool sortInputs = true)
    {
        assert(nSamples > 1);

        setNumSamples_(nSamples);
        for (std::size_t i = 0; i < nSamples; ++i) {
            xPos_[i] = x[i];
            yPos_[i] = y[i];
        }

        if (sortInputs)
            sortInput_();
        else if (xPos_[0] > xPos_[numSamples() - 1])
            reverseSamplingPoints_();

        if (splineType == Periodic)
            makePeriodicSpline_();
        else if (splineType == Natural)
            makeNaturalSpline_();
        else if (splineType == Monotonic)
            this->makeMonotonicSpline_(slopeVec_);
        else
            throw std::runtime_error("Spline type "+std::to_string(int(splineType))+" not supported at this place");
    }

    template <class ScalarContainerX, class ScalarContainerY>
    void setXYContainers(const ScalarContainerX& x,
                         const ScalarContainerY& y,
                         SplineType splineType = Natural,
                         bool sortInputs = true)
    {
        assert(x.size() == y.size());
        assert(x.size() > 1);

        setNumSamples_(x.size());
        std::ranges::copy(x, xPos_.begin());
        std::ranges::copy(y, yPos_.begin());

        if (sortInputs)
            sortInput_();
        else if (xPos_[0] > xPos_[numSamples() - 1])
            reverseSamplingPoints_();

        if (splineType == Periodic)
            makePeriodicSpline_();
        else if (splineType == Natural)
            makeNaturalSpline_();
        else if (splineType == Monotonic)
            this->makeMonotonicSpline_(slopeVec_);
        else
            throw std::runtime_error("Spline type "+std::to_string(int(splineType))+" not supported at this place");
    }

    template <class PointArray>
    void setArrayOfPoints(std::size_t nSamples,
                          const PointArray& points,
                          SplineType splineType = Natural,
                          bool sortInputs = true)
    {
        // a spline with no or just one sampling points? what an
        // incredible bad idea!
        assert(nSamples > 1);

        setNumSamples_(nSamples);
        for (std::size_t i = 0; i < nSamples; ++i) {
            xPos_[i] = points[i][0];
            yPos_[i] = points[i][1];
        }

        if (sortInputs)
            sortInput_();
        else if (xPos_[0] > xPos_[numSamples() - 1])
            reverseSamplingPoints_();

        if (splineType == Periodic)
            makePeriodicSpline_();
        else if (splineType == Natural)
            makeNaturalSpline_();
        else if (splineType == Monotonic)
            this->makeMonotonicSpline_(slopeVec_);
        else
            throw std::runtime_error("Spline type "+std::to_string(int(splineType))+" not supported at this place");
    }

    template <class XYContainer>
    void setContainerOfPoints(const XYContainer& points,
                              SplineType splineType = Natural,
                              bool sortInputs = true)
    {
        // a spline with no or just one sampling points? what an
        // incredible bad idea!
        assert(points.size() > 1);

        setNumSamples_(points.size());
        typename XYContainer::const_iterator it = points.begin();
        typename XYContainer::const_iterator endIt = points.end();
        for (std::size_t i = 0; it != endIt; ++ i, ++it) {
            xPos_[i] = (*it)[0];
            yPos_[i] = (*it)[1];
        }

        if (sortInputs)
            sortInput_();
        else if (xPos_[0] > xPos_[numSamples() - 1])
            reverseSamplingPoints_();

        if (splineType == Periodic)
            makePeriodicSpline_();
        else if (splineType == Natural)
            makeNaturalSpline_();
        else if (splineType == Monotonic)
            this->makeMonotonicSpline_(slopeVec_);
        else
            throw std::runtime_error("Spline type "+std::to_string(int(splineType))+" not supported at this place");
    }

    template <class XYContainer>
    void setContainerOfTuples(const XYContainer& points,
                              SplineType splineType = Natural,
                              bool sortInputs = true)
    {
        // resize internal arrays
        setNumSamples_(points.size());
        typename XYContainer::const_iterator it = points.begin();
        typename XYContainer::const_iterator endIt = points.end();
        for (unsigned i = 0; it != endIt; ++i, ++it) {
            xPos_[i] = std::get<0>(*it);
            yPos_[i] = std::get<1>(*it);
        }

        if (sortInputs)
            sortInput_();
        else if (xPos_[0] > xPos_[numSamples() - 1])
            reverseSamplingPoints_();

        if (splineType == Periodic)
            makePeriodicSpline_();
        else if (splineType == Natural)
            makeNaturalSpline_();
        else if (splineType == Monotonic)
            this->makeMonotonicSpline_(slopeVec_);
        else
            throw std::runtime_error("Spline type "+std::to_string(int(splineType))+" not supported at this place");
    }

    template <class Evaluation>
    bool applies(const Evaluation& x) const
    { return x_(0) <= x && x <= x_(numSamples() - 1); }

    std::size_t numSamples() const
    { return xPos_.size(); }

    Scalar xAt(std::size_t sampleIdx) const
    { return x_(sampleIdx); }

    Scalar valueAt(std::size_t sampleIdx) const
    { return y_(sampleIdx); }

    void printCSV(Scalar xi0, Scalar xi1, std::size_t k, std::ostream& os) const;

    template <class Evaluation>
    Evaluation eval(const Evaluation& x, bool extrapolate = false) const
    {
        if (!extrapolate && !applies(x))
            throw NumericalProblem("Tried to evaluate a spline outside of its range");

        // handle extrapolation
        if (extrapolate) {
            if (x < xAt(0)) {
                Scalar m = evalDerivative_(xAt(0), /*segmentIdx=*/0);
                Scalar y0 = y_(0);
                return y0 + m*(x - xAt(0));
            }
            else if (x > xAt(static_cast<std::size_t>(static_cast<long int>(numSamples()) - 1))) {
                Scalar m = evalDerivative_(xAt(static_cast<std::size_t>(numSamples() - 1)),
                                           /*segmentIdx=*/static_cast<std::size_t>(numSamples()-2));
                Scalar y0 = y_(static_cast<std::size_t>(numSamples() - 1));
                return y0 + m*(x - xAt(static_cast<std::size_t>(numSamples() - 1)));
            }
        }

        return eval_(x, segmentIdx_(scalarValue(x)));
    }

    template <class Evaluation>
    Evaluation evalDerivative(const Evaluation& x, bool extrapolate = false) const
    {
        if (!extrapolate && !applies(x))
            throw NumericalProblem("Tried to evaluate the derivative of a spline outside of its range");

        // handle extrapolation
        if (extrapolate) {
            if (x < xAt(0))
                return evalDerivative_(xAt(0), /*segmentIdx=*/0);
            else if (x > xAt(numSamples() - 1))
                return evalDerivative_(xAt(numSamples() - 1), /*segmentIdx=*/numSamples() - 2);
        }

        return evalDerivative_(x, segmentIdx_(scalarValue(x)));
    }

    template <class Evaluation>
    Evaluation evalSecondDerivative(const Evaluation& x, bool extrapolate = false) const
    {
        if (!extrapolate && !applies(x))
            throw NumericalProblem("Tried to evaluate the second derivative of a spline outside of its range");
        else if (extrapolate)
            return 0.0;

        return evalDerivative2_(x, segmentIdx_(scalarValue(x)));
    }

    template <class Evaluation>
    Evaluation evalThirdDerivative(const Evaluation& x, bool extrapolate = false) const
    {
        if (!extrapolate && !applies(x))
            throw NumericalProblem("Tried to evaluate the third derivative of a spline outside of its range");
        else if (extrapolate)
            return 0.0;

        return evalDerivative3_(x, segmentIdx_(scalarValue(x)));
    }

    template <class Evaluation>
    Evaluation intersect(const Evaluation& a,
                         const Evaluation& b,
                         const Evaluation& c,
                         const Evaluation& d) const
    {
        return intersectInterval(xAt(0), xAt(numSamples() - 1), a, b, c, d);
    }

    template <class Evaluation>
    Evaluation intersectInterval(Scalar x0, Scalar x1,
                                 const Evaluation& a,
                                 const Evaluation& b,
                                 const Evaluation& c,
                                 const Evaluation& d) const
    {
        assert(applies(x0) && applies(x1));

        Evaluation tmpSol[3], sol = 0;
        std::size_t nSol = 0;
        std::size_t iFirst = segmentIdx_(x0);
        std::size_t iLast = segmentIdx_(x1);
        for (std::size_t i = iFirst; i <= iLast; ++i)
        {
            std::size_t nCur = intersectSegment_(tmpSol, i, a, b, c, d, x0, x1);
            if (nCur == 1)
                sol = tmpSol[0];

            nSol += nCur;
            if (nSol > 1) {
                throw std::runtime_error("Spline has more than one intersection"); //<<a<<"x^3 + "<<b<"x^2 + "<<c<"x + "<<d);
            }
        }

        if (nSol != 1)
            throw std::runtime_error("Spline has no intersection"); //<<a<"x^3 + " <<b<"x^2 + "<<c<"x + "<<d<<"!");

        return sol;
    }

    int monotonic(Scalar x0, Scalar x1,
                  [[maybe_unused]] bool extrapolate = false) const
    {
        assert(std::abs(x0 - x1) > 1e-30);

        // make sure that x0 is smaller than x1
        if (x0 > x1)
            std::swap(x0, x1);

        assert(x0 < x1);

        int r = 3;
        if (x0 < xAt(0)) {
            assert(extrapolate);
            Scalar m = evalDerivative_(xAt(0), /*segmentIdx=*/0);
            if (std::abs(m) < 1e-20)
                r = (m < 0)?-1:1;
            x0 = xAt(0);
        };

        std::size_t i = segmentIdx_(x0);
        if (x_(i + 1) >= x1) {
            // interval is fully contained within a single spline
            // segment
            monotonic_(i, x0, x1, r);
            return r;
        }

        // the first segment overlaps with the specified interval
        // partially
        monotonic_(i, x0, x_(i+1), r);
        ++ i;

        // make sure that the segments which are completly in the
        // interval [x0, x1] all exhibit the same monotonicity.
        std::size_t iEnd = segmentIdx_(x1);
        for (; i < iEnd - 1; ++i) {
            monotonic_(i, x_(i), x_(i + 1), r);
            if (!r)
                return 0;
        }

        // if the user asked for a part of the spline which is
        // extrapolated, we need to check the slope at the spline's
        // endpoint
        if (x1 > xAt(numSamples() - 1)) {
            assert(extrapolate);

            Scalar m = evalDerivative_(xAt(numSamples() - 1), /*segmentIdx=*/numSamples() - 2);
            if (m < 0)
                return (r < 0 || r==3) ? -1 : 0;
            else if (m > 0)
                return r > 0  ? 1 : 0;

            return r;
        }

        // check for the last segment
        monotonic_(iEnd, x_(iEnd), x1, r);

        return r;
    }

    int monotonic() const
    { return monotonic(xAt(0), xAt(numSamples() - 1)); }

protected:
    struct ComparatorX_
    {
        explicit ComparatorX_(const std::vector<Scalar>& x)
            : x_(x)
        {}

        bool operator ()(unsigned idxA, unsigned idxB) const
        { return x_.at(idxA) < x_.at(idxB); }

        const std::vector<Scalar>& x_;
    };

    void sortInput_()
    {
        std::size_t n = numSamples();

        // create a vector containing 0...n-1
        std::vector<unsigned> idxVector(n);
        for (unsigned i = 0; i < n; ++i)
            idxVector[i] = i;

        // sort the indices according to the x values of the sample
        // points
        ComparatorX_ cmp(xPos_);
        std::ranges::sort(idxVector, cmp);

        // reorder the sample points
        std::vector<Scalar> tmpX(n), tmpY(n);
        for (std::size_t i = 0; i < idxVector.size(); ++ i) {
            tmpX[i] = xPos_[idxVector[i]];
            tmpY[i] = yPos_[idxVector[i]];
        }
        xPos_ = tmpX;
        yPos_ = tmpY;
    }

    void reverseSamplingPoints_()
    {
        // reverse the arrays
        std::size_t n = numSamples();
        for (unsigned i = 0; i <= (n - 1)/2; ++i) {
            std::swap(xPos_[i], xPos_[n - i - 1]);
            std::swap(yPos_[i], yPos_[n - i - 1]);
        }
    }

    void setNumSamples_(std::size_t nSamples)
    {
        xPos_.resize(nSamples);
        yPos_.resize(nSamples);
        slopeVec_.resize(nSamples);
    }

    void makeFullSpline_(Scalar m0, Scalar m1)
    {
        Matrix M(numSamples());
        std::vector<Scalar> d(numSamples());
        std::vector<Scalar> moments(numSamples());

        // create linear system of equations
        this->makeFullSystem_(M, d, m0, m1);

        // solve for the moments (-> second derivatives)
        M.solve(moments, d);

        // convert the moments to slopes at the sample points
        this->setSlopesFromMoments_(slopeVec_, moments);
    }

    void makeNaturalSpline_()
    {
        Matrix M(numSamples(), numSamples());
        Vector d(numSamples());
        Vector moments(numSamples());

        // create linear system of equations
        this->makeNaturalSystem_(M, d);

        // solve for the moments (-> second derivatives)
        M.solve(moments, d);

        // convert the moments to slopes at the sample points
        this->setSlopesFromMoments_(slopeVec_, moments);
    }

    void makePeriodicSpline_()
    {
        Matrix M(numSamples() - 1);
        Vector d(numSamples() - 1);
        Vector moments(numSamples() - 1);

        // create linear system of equations. This is a bit hacky,
        // because it assumes that std::vector internally stores its
        // data as a big C-style array, but it saves us from yet
        // another copy operation
        this->makePeriodicSystem_(M, d);

        // solve for the moments (-> second derivatives)
        M.solve(moments, d);

        moments.resize(numSamples());
        for (int i = static_cast<int>(numSamples()) - 2; i >= 0; --i) {
            unsigned ui = static_cast<unsigned>(i);
            moments[ui+1] = moments[ui];
        }
        moments[0] = moments[numSamples() - 1];

        // convert the moments to slopes at the sample points
        this->setSlopesFromMoments_(slopeVec_, moments);
    }

    template <class DestVector, class SourceVector>
    void assignSamplingPoints_(DestVector& destX,
                               DestVector& destY,
                               const SourceVector& srcX,
                               const SourceVector& srcY,
                               unsigned nSamples)
    {
        assert(nSamples >= 2);

        // copy sample points, make sure that the first x value is
        // smaller than the last one
        for (unsigned i = 0; i < nSamples; ++i) {
            unsigned idx = i;
            if (srcX[0] > srcX[nSamples - 1])
                idx = nSamples - i - 1;
            destX[i] = srcX[idx];
            destY[i] = srcY[idx];
        }
    }

    template <class DestVector, class ListIterator>
    void assignFromArrayList_(DestVector& destX,
                              DestVector& destY,
                              const ListIterator& srcBegin,
                              const ListIterator& srcEnd,
                              unsigned nSamples)
    {
        assert(nSamples >= 2);

        // find out wether the x values are in reverse order
        ListIterator it = srcBegin;
        ++it;
        bool reverse = false;
        if ((*srcBegin)[0] > (*it)[0])
            reverse = true;
        --it;

        // loop over all sampling points
        for (unsigned i = 0; it != srcEnd; ++i, ++it) {
            unsigned idx = i;
            if (reverse)
                idx = nSamples - i - 1;
            destX[idx] = (*it)[0];
            destY[idx] = (*it)[1];
        }
    }

    template <class DestVector, class ListIterator>
    void assignFromTupleList_(DestVector& destX,
                              DestVector& destY,
                              ListIterator srcBegin,
                              ListIterator srcEnd,
                              unsigned nSamples)
    {
        assert(nSamples >= 2);

        // copy sample points, make sure that the first x value is
        // smaller than the last one

        // find out wether the x values are in reverse order
        ListIterator it = srcBegin;
        ++it;
        bool reverse = false;
        if (std::get<0>(*srcBegin) > std::get<0>(*it))
            reverse = true;
        --it;

        // loop over all sampling points
        for (unsigned i = 0; it != srcEnd; ++i, ++it) {
            unsigned idx = i;
            if (reverse)
                idx = nSamples - i - 1;
            destX[idx] = std::get<0>(*it);
            destY[idx] = std::get<1>(*it);
        }
    }

    template <class Vector, class Matrix>
    void makeFullSystem_(Matrix& M, Vector& d, Scalar m0, Scalar m1)
    {
        makeNaturalSystem_(M, d);

        std::size_t n = numSamples() - 1;
        // first row
        M[0][1] = 1;
        d[0] = 6/h_(1) * ( (y_(1) - y_(0))/h_(1) - m0);

        // last row
        M[n][n - 1] = 1;

        // right hand side
        d[n] =
            6/h_(n)
            *
            (m1 - (y_(n) - y_(n - 1))/h_(n));
    }

    template <class Vector, class Matrix>
    void makeNaturalSystem_(Matrix& M, Vector& d)
    {
        M = 0.0;

        // See: J. Stoer: "Numerische Mathematik 1", 9th edition,
        // Springer, 2005, p. 111
        std::size_t n = numSamples() - 1;

        // second to next to last rows
        for (std::size_t i = 1; i < n; ++i) {
            Scalar lambda_i = h_(i + 1) / (h_(i) + h_(i + 1));
            Scalar mu_i = 1 - lambda_i;
            Scalar d_i =
                6 / (h_(i) + h_(i + 1))
                *
                ( (y_(i + 1) - y_(i))/h_(i + 1) - (y_(i) - y_(i - 1))/h_(i));

            M[i][i-1] = mu_i;
            M[i][i] = 2;
            M[i][i + 1] = lambda_i;
            d[i] = d_i;
        };

        // See Stroer, equation (2.5.2.7)
        Scalar lambda_0 = 0;
        Scalar d_0 = 0;

        Scalar mu_n = 0;
        Scalar d_n = 0;

        // first row
        M[0][0] = 2;
        M[0][1] = lambda_0;
        d[0] = d_0;

        // last row
        M[n][n-1] = mu_n;
        M[n][n] = 2;
        d[n] = d_n;
    }

    template <class Matrix, class Vector>
    void makePeriodicSystem_(Matrix& M, Vector& d)
    {
        M = 0.0;

        // See: J. Stoer: "Numerische Mathematik 1", 9th edition,
        // Springer, 2005, p. 111
        std::size_t n = numSamples() - 1;

        assert(M.rows() == n);

        // second to next to last rows
        for (std::size_t i = 2; i < n; ++i) {
            Scalar lambda_i = h_(i + 1) / (h_(i) + h_(i + 1));
            Scalar mu_i = 1 - lambda_i;
            Scalar d_i =
                6 / (h_(i) + h_(i + 1))
                *
                ( (y_(i + 1) - y_(i))/h_(i + 1) - (y_(i) - y_(i - 1))/h_(i));

            M[i-1][i-2] = mu_i;
            M[i-1][i-1] = 2;
            M[i-1][i] = lambda_i;
            d[i-1] = d_i;
        };

        Scalar lambda_n = h_(1) / (h_(n) + h_(1));
        Scalar lambda_1 = h_(2) / (h_(1) + h_(2));;
        Scalar mu_1 = 1 - lambda_1;
        Scalar mu_n = 1 - lambda_n;

        Scalar d_1 =
            6 / (h_(1) + h_(2))
            *
            ( (y_(2) - y_(1))/h_(2) - (y_(1) - y_(0))/h_(1));
        Scalar d_n =
            6 / (h_(n) + h_(1))
            *
            ( (y_(1) - y_(n))/h_(1) - (y_(n) - y_(n-1))/h_(n));

        // first row
        M[0][0] = 2;
        M[0][1] = lambda_1;
        M[0][n-1] = mu_1;
        d[0] = d_1;

        // last row
        M[n-1][0] = lambda_n;
        M[n-1][n-2] = mu_n;
        M[n-1][n-1] = 2;
        d[n-1] = d_n;
    }

    template <class Vector>
    void makeMonotonicSpline_(Vector& slopes)
    {
        auto n = numSamples();

        // calculate the slopes of the secant lines
        std::vector<Scalar> delta(n);
        for (std::size_t k = 0; k < n - 1; ++k)
            delta[k] = (y_(k + 1) - y_(k))/(x_(k + 1) - x_(k));

        // calculate the "raw" slopes at the sample points
        for (std::size_t k = 1; k < n - 1; ++k)
            slopes[k] = (delta[k - 1] + delta[k])/2;
        slopes[0] = delta[0];
        slopes[n - 1] = delta[n - 2];

        // post-process the "raw" slopes at the sample points
        for (std::size_t k = 0; k < n - 1; ++k) {
            if (std::abs(delta[k]) < 1e-50) {
                // make the spline flat if the inputs are equal
                slopes[k] = 0;
                slopes[k + 1] = 0;
                ++ k;
                continue;
            }
            else {
                Scalar alpha = slopes[k] / delta[k];
                Scalar beta = slopes[k + 1] / delta[k];

                if (alpha < 0 || (k > 0 && slopes[k] / delta[k - 1] < 0)) {
                    slopes[k] = 0;
                }
                // limit (alpha, beta) to a circle of radius 3
                else if (alpha*alpha + beta*beta > 3*3) {
                    Scalar tau = 3.0/std::sqrt(alpha*alpha + beta*beta);
                    slopes[k] = tau*alpha*delta[k];
                    slopes[k + 1] = tau*beta*delta[k];
                }
            }
        }
    }

    template <class MomentsVector, class SlopeVector>
    void setSlopesFromMoments_(SlopeVector& slopes, const MomentsVector& moments)
    {
        std::size_t n = numSamples();

        // evaluate slope at the rightmost point.
        // See: J. Stoer: "Numerische Mathematik 1", 9th edition,
        // Springer, 2005, p. 109
        Scalar mRight;

        {
            Scalar h = this->h_(n - 1);
            Scalar x = h;
            //Scalar x_1 = 0;

            Scalar A =
                (y_(n - 1) - y_(n - 2))/h
                -
                h/6*(moments[n-1] - moments[n - 2]);

            mRight =
                //- moments[n - 2] * x_1*x_1 / (2 * h)
                //+
                moments[n - 1] * x*x / (2 * h)
                +
                A;
        }

        // evaluate the slope for the first n-1 sample points
        for (std::size_t i = 0; i < n - 1; ++ i) {
            // See: J. Stoer: "Numerische Mathematik 1", 9th edition,
            // Springer, 2005, p. 109
            Scalar h_i = this->h_(i + 1);
            //Scalar x_i = 0;
            Scalar x_i1 = h_i;

            Scalar A_i =
                (y_(i+1) - y_(i))/h_i
                -
                h_i/6*(moments[i+1] - moments[i]);

            slopes[i] =
                - moments[i] * x_i1*x_i1 / (2 * h_i)
                +
                //moments[i + 1] * x_i*x_i / (2 * h_i)
                //+
                A_i;

        }
        slopes[n - 1] = mRight;
    }

    // evaluate the spline at a given the position and given the
    // segment index
    template <class Evaluation>
    Evaluation eval_(const Evaluation& x, std::size_t i) const
    {
        // See http://en.wikipedia.org/wiki/Cubic_Hermite_spline
        Scalar delta = h_(i + 1);
        Evaluation t = (x - x_(i))/delta;

        return
            h00_(t) * y_(i)
            + h10_(t) * slope_(i)*delta
            + h01_(t) * y_(i + 1)
            + h11_(t) * slope_(i + 1)*delta;
    }

    // evaluate the derivative of a spline given the actual position
    // and the segment index
    template <class Evaluation>
    Evaluation evalDerivative_(const Evaluation& x, std::size_t i) const
    {
        // See http://en.wikipedia.org/wiki/Cubic_Hermite_spline
        Scalar delta = h_(i + 1);
        Evaluation t = (x - x_(i))/delta;
        Evaluation alpha = 1 / delta;

        return
            alpha *
            (h00_prime_(t) * y_(i)
             + h10_prime_(t) * slope_(i)*delta
             + h01_prime_(t) * y_(i + 1)
             + h11_prime_(t) * slope_(i + 1)*delta);
    }

    // evaluate the second derivative of a spline given the actual
    // position and the segment index
    template <class Evaluation>
    Evaluation evalDerivative2_(const Evaluation& x, std::size_t i) const
    {
        // See http://en.wikipedia.org/wiki/Cubic_Hermite_spline
        Scalar delta = h_(i + 1);
        Evaluation t = (x - x_(i))/delta;
        Evaluation alpha = 1 / delta;

        return
            alpha*alpha
            *(h00_prime2_(t) * y_(i)
              + h10_prime2_(t) * slope_(i)*delta
              + h01_prime2_(t) * y_(i + 1)
              + h11_prime2_(t) * slope_(i + 1)*delta);
    }

    // evaluate the third derivative of a spline given the actual
    // position and the segment index
    template <class Evaluation>
    Evaluation evalDerivative3_(const Evaluation& x, std::size_t i) const
    {
        // See http://en.wikipedia.org/wiki/Cubic_Hermite_spline
        Scalar delta = h_(i + 1);
        Evaluation t = (x - x_(i))/delta;
        Evaluation alpha = 1 / delta;

        return
            alpha*alpha*alpha
            *(h00_prime3_(t)*y_(i)
              + h10_prime3_(t)*slope_(i)*delta
              + h01_prime3_(t)*y_(i + 1)
              + h11_prime3_(t)*slope_(i + 1)*delta);
    }

    // hermite basis functions
    template <class Evaluation>
    Evaluation h00_(const Evaluation& t) const
    { return (2*t - 3)*t*t + 1; }

    template <class Evaluation>
    Evaluation h10_(const Evaluation& t) const
    { return ((t - 2)*t + 1)*t; }

    template <class Evaluation>
    Evaluation h01_(const Evaluation& t) const
    { return (-2*t + 3)*t*t; }

    template <class Evaluation>
    Evaluation h11_(const Evaluation& t) const
    { return (t - 1)*t*t; }

    // first derivative of the hermite basis functions
    template <class Evaluation>
    Evaluation h00_prime_(const Evaluation& t) const
    { return (3*2*t - 2*3)*t; }

    template <class Evaluation>
    Evaluation h10_prime_(const Evaluation& t) const
    { return (3*t - 2*2)*t + 1; }

    template <class Evaluation>
    Evaluation h01_prime_(const Evaluation& t) const
    { return (-3*2*t + 2*3)*t; }

    template <class Evaluation>
    Evaluation h11_prime_(const Evaluation& t) const
    { return (3*t - 2)*t; }

    // second derivative of the hermite basis functions
    template <class Evaluation>
    Evaluation h00_prime2_(const Evaluation& t) const
    { return 2*3*2*t - 2*3; }

    template <class Evaluation>
    Evaluation h10_prime2_(const Evaluation& t) const
    { return 2*3*t - 2*2; }

    template <class Evaluation>
    Evaluation h01_prime2_(const Evaluation& t) const
    { return -2*3*2*t + 2*3; }

    template <class Evaluation>
    Evaluation h11_prime2_(const Evaluation& t) const
    { return 2*3*t - 2; }

    // third derivative of the hermite basis functions
    template <class Evaluation>
    Scalar h00_prime3_(const Evaluation&) const
    { return 2*3*2; }

    template <class Evaluation>
    Scalar h10_prime3_(const Evaluation&) const
    { return 2*3; }

    template <class Evaluation>
    Scalar h01_prime3_(const Evaluation&) const
    { return -2*3*2; }

    template <class Evaluation>
    Scalar h11_prime3_(const Evaluation&) const
    { return 2*3; }

    // returns the monotonicality of an interval of a spline segment
    //
    // The return value have the following meaning:
    //
    // 3: spline is constant within interval [x0, x1]
    // 1: spline is monotonously increasing in the specified interval
    // 0: spline is not monotonic (or constant) in the specified interval
    // -1: spline is monotonously decreasing in the specified interval
    int monotonic_(std::size_t i, Scalar x0, Scalar x1, int& r) const
    {
        // coefficients of derivative in monomial basis
        Scalar a = 3*a_(i);
        Scalar b = 2*b_(i);
        Scalar c = c_(i);

        if (std::abs(a) < 1e-20 && std::abs(b) < 1e-20 && std::abs(c) < 1e-20)
            return 3; // constant in interval, r stays unchanged!

        Scalar disc = b*b - 4*a*c;
        if (disc < 0) {
            // discriminant of derivative is smaller than 0, i.e. the
            // segment's derivative does not exhibit any extrema.
            if (x0*(x0*a + b) + c > 0) {
                r = (r==3 || r == 1)?1:0;
                return 1;
            }
            else {
                r = (r==3 || r == -1)?-1:0;
                return -1;
            }
        }
        disc = std::sqrt(disc);
        Scalar xE1 = (-b + disc)/(2*a);
        Scalar xE2 = (-b - disc)/(2*a);

        if (std::abs(disc) < 1e-30) {
            // saddle point -> no extrema
            if (std::abs(xE1 - x0) < 1e-30)
                // make sure that we're not picking the saddle point
                // to determine whether we're monotonically increasing
                // or decreasing
                x0 = x1;
            if (x0*(x0*a + b) + c > 0) {
                r = (r==3 || r == 1)?1:0;
                return 1;
            }
            else {
                r = (r==3 || r == -1)?-1:0;
                return -1;
            }
        };
        if ((x0 < xE1 && xE1 < x1) ||
            (x0 < xE2 && xE2 < x1))
        {
            // there is an extremum in the range (x0, x1)
            r = 0;
            return 0;
        }
        // no extremum in range (x0, x1)
        x0 = (x0 + x1)/2; // pick point in the middle of the interval
                          // to avoid extrema on the boundaries
        if (x0*(x0*a + b) + c > 0) {
            r = (r==3 || r == 1)?1:0;
            return 1;
        }
        else {
            r = (r==3 || r == -1)?-1:0;
            return -1;
        }
    }

    template <class Evaluation>
    std::size_t intersectSegment_(Evaluation* sol,
                             std::size_t segIdx,
                             const Evaluation& a,
                             const Evaluation& b,
                             const Evaluation& c,
                             const Evaluation& d,
                             Scalar x0 = -1e30, Scalar x1 = 1e30) const
    {
        unsigned n =
            invertCubicPolynomial(sol,
                                       a_(segIdx) - a,
                                       b_(segIdx) - b,
                                       c_(segIdx) - c,
                                       d_(segIdx) - d);
        x0 = std::max(x_(segIdx), x0);
        x1 = std::min(x_(segIdx+1), x1);

        // filter the intersections outside of the specified interval
        std::size_t k = 0;
        for (unsigned j = 0; j < n; ++j) {
            if (x0 <= sol[j] && sol[j] <= x1) {
                sol[k] = sol[j];
                ++k;
            }
        }
        return k;
    }

    // find the segment index for a given x coordinate
    std::size_t segmentIdx_(Scalar x) const
    {
        // bisection
        std::size_t iLow = 0;
        std::size_t iHigh = numSamples() - 1;

        while (iLow + 1 < iHigh) {
            std::size_t i = (iLow + iHigh) / 2;
            if (x_(i) > x)
                iHigh = i;
            else
                iLow = i;
        };
        return iLow;
    }

    Scalar h_(std::size_t i) const
    {
        assert(x_(i) > x_(i-1)); // the sampling points must be given
                                 // in ascending order
        return x_(i) - x_(i - 1);
    }

    Scalar x_(std::size_t i) const
    { return xPos_[i]; }

    Scalar y_(std::size_t i) const
    { return yPos_[i]; }

    Scalar slope_(std::size_t i) const
    { return slopeVec_[i]; }

    // returns the coefficient in front of the x^3 term. In Stoer this
    // is delta.
    Scalar a_(std::size_t i) const
    { return evalDerivative3_(/*x=*/Scalar(0.0), i)/6.0; }

    // returns the coefficient in front of the x^2 term In Stoer this
    // is gamma.
    Scalar b_(std::size_t i) const
    { return evalDerivative2_(/*x=*/Scalar(0.0), i)/2.0; }

    // returns the coefficient in front of the x^1 term. In Stoer this
    // is beta.
    Scalar c_(std::size_t i) const
    { return evalDerivative_(/*x=*/Scalar(0.0), i); }

    // returns the coefficient in front of the x^0 term. In Stoer this
    // is alpha.
    Scalar d_(std::size_t i) const
    { return eval_(/*x=*/Scalar(0.0), i); }

    Vector xPos_;
    Vector yPos_;
    Vector slopeVec_;
};
}

#endif