Tabulated1DFunction.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:
/*
  Copyright (C) 2015 by Andreas Lauser

  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/>.
*/
#ifndef OPM_TABULATED_1D_FUNCTION_HPP
#define OPM_TABULATED_1D_FUNCTION_HPP

#include <opm/common/ErrorMacros.hpp>
#include <opm/common/Exceptions.hpp>
#include <opm/material/common/Unused.hpp>

#include <algorithm>
#include <cassert>
#include <iostream>
#include <tuple>
#include <vector>

namespace Opm {
template <class Scalar>
class Tabulated1DFunction
{
public:
    Tabulated1DFunction()
    {}

    template <class ScalarArrayX, class ScalarArrayY>
    Tabulated1DFunction(size_t nSamples,
                        const ScalarArrayX &x,
                        const ScalarArrayY &y,
                        bool sortInputs = false)
    { this->setXYArrays(nSamples, x, y, sortInputs); }

    template <class ScalarContainer>
    Tabulated1DFunction(const ScalarContainer &x,
                        const ScalarContainer &y,
                        bool sortInputs = false)
    { this->setXYContainers(x, y, sortInputs); }

    template <class PointContainer>
    Tabulated1DFunction(const PointContainer &points,
                        bool sortInputs = false)
    { this->setContainerOfTuples(points, sortInputs); }

    template <class ScalarArrayX, class ScalarArrayY>
    void setXYArrays(size_t nSamples,
                     const ScalarArrayX &x,
                     const ScalarArrayY &y,
                     bool sortInputs = false)
    {
        assert(nSamples > 1);

        resizeArrays_(nSamples);
        for (size_t i = 0; i < nSamples; ++i) {
            xValues_[i] = x[i];
            yValues_[i] = y[i];
        }

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

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

        resizeArrays_(x.size());
        std::copy(x.begin(), x.end(), xValues_.begin());
        std::copy(y.begin(), y.end(), yValues_.begin());

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

    template <class PointArray>
    void setArrayOfPoints(size_t nSamples,
                          const PointArray &points,
                          bool sortInputs = false)
    {
        // a linear function with less than two sampling points? what an incredible
        // bad idea!
        assert(nSamples > 1);

        resizeArrays_(nSamples);
        for (size_t i = 0; i < nSamples; ++i) {
            xValues_[i] = points[i][0];
            yValues_[i] = points[i][1];
        }

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

    template <class XYContainer>
    void setContainerOfTuples(const XYContainer &points,
                              bool sortInputs = false)
    {
        // a linear function with less than two sampling points? what an incredible
        // bad idea!
        assert(points.size() > 1);

        resizeArrays_(points.size());
        typename XYContainer::const_iterator it = points.begin();
        typename XYContainer::const_iterator endIt = points.end();
        for (int i = 0; it != endIt; ++i, ++it) {
            xValues_[i] = std::get<0>(*it);
            yValues_[i] = std::get<1>(*it);
        }

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

    size_t numSamples() const
    { return xValues_.size(); }

    Scalar xMin() const
    { return xValues_[0]; }

    Scalar xMax() const
    { return xValues_[numSamples() - 1]; }

    Scalar xAt(size_t i) const
    { return xValues_[i]; }

    Scalar valueAt(size_t i) const
    { return yValues_[i]; }

    bool applies(Scalar x) const
    { return xValues_[0] <= x && x <= xValues_[numSamples() - 1]; }

    Scalar eval(Scalar x, bool extrapolate=false) const
    {
        size_t segIdx;
        if (extrapolate && x < xValues_.front())
            segIdx = 0;
        else if (extrapolate && x > xValues_.back())
            segIdx = numSamples() - 2;
        else {
            assert(xValues_.front() <= x && x <= xValues_.back());
            segIdx = findSegmentIndex_(x);
        }

        Scalar x0 = xValues_[segIdx];
        Scalar x1 = xValues_[segIdx + 1];

        Scalar y0 = yValues_[segIdx];
        Scalar y1 = yValues_[segIdx + 1];

        return y0 + (y1 - y0)*(x - x0)/(x1 - x0);
    }

    template <class Evaluation>
    Evaluation eval(const Evaluation& x, bool extrapolate=false) const
    {
        size_t segIdx;
        if (extrapolate && x.value < xValues_.front())
            segIdx = 0;
        else if (extrapolate && x.value > xValues_.back())
            segIdx = numSamples() - 2;
        else {
            assert(xValues_.front() <= x.value && x.value <= xValues_.back());
            segIdx = findSegmentIndex_(x.value);
        }

        Scalar x0 = xValues_[segIdx];
        Scalar x1 = xValues_[segIdx + 1];

        Scalar y0 = yValues_[segIdx];
        Scalar y1 = yValues_[segIdx + 1];

        Scalar m = (y1 - y0)/(x1 - x0);

        Evaluation result;
        result.value = y0 + m*(x.value - x0);
        for (unsigned varIdx = 0; varIdx < result.derivatives.size(); ++varIdx)
            result.derivatives[varIdx] = m*x.derivatives[varIdx];

        return result;
    }

    Scalar evalDerivative(Scalar x, bool /*extrapolate*/=false) const
    {
        int segIdx = findSegmentIndex_(x);

        return evalDerivative(x, segIdx);
    }

    Scalar evalSecondDerivative(OPM_OPTIM_UNUSED Scalar x, OPM_OPTIM_UNUSED bool extrapolate=false) const
    {
        assert(extrapolate || applies(x));
        return 0.0;
    }

    Scalar evalThirdDerivative(OPM_OPTIM_UNUSED Scalar x, OPM_OPTIM_UNUSED bool extrapolate=false) const
    {
        assert(extrapolate || applies(x));
        return 0.0;
    }

    int monotonic(Scalar x0, Scalar x1, bool extrapolate OPM_OPTIM_UNUSED = false) const
    {
        assert(x0 != x1);

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

        assert(x0 < x1);

        int r = 3;
        if (x0 < xMin()) {
            assert(extrapolate);

            x0 = xMin();
        };

        size_t i = findSegmentIndex_(x0);
        if (xValues_[i + 1] >= x1) {
            // interval is fully contained within a single function
            // segment
            updateMonotonicity_(i, r);
            return r;
        }

        // the first segment overlaps with the specified interval
        // partially
        updateMonotonicity_(i, r);
        ++ i;

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

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

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

            return r;
        }

        // check for the last segment
        updateMonotonicity_(iEnd, r);

        return r;
    }

    int monotonic() const
    { return monotonic(xMin(), xMax()); }

    void printCSV(Scalar xi0, Scalar xi1, int k, std::ostream &os = std::cout) const
    {
        Scalar x0 = std::min(xi0, xi1);
        Scalar x1 = std::max(xi0, xi1);
        const int n = numSamples() - 1;
        for (int i = 0; i <= k; ++i) {
            double x = i*(x1 - x0)/k + x0;
            double y;
            double dy_dx;
            if (!applies(x)) {
                if (x < xValues_[0]) {
                    dy_dx = evalDerivative(xValues_[0]);
                    y = (x - xValues_[0])*dy_dx + yValues_[0];
                }
                else if (x > xValues_[n]) {
                    dy_dx = evalDerivative(xValues_[n]);
                    y = (x - xValues_[n])*dy_dx + yValues_[n];
                }
                else {
                    OPM_THROW(std::runtime_error,
                              "The sampling points given to a function must be sorted by their x value!");
                }
            }
            else {
                y = eval(x);
                dy_dx = evalDerivative(x);
            }

            os << x << " " << y << " " << dy_dx << "\n";
        }
    }

private:
    size_t findSegmentIndex_(Scalar x) const
    {
        // we need at least two sampling points!
        assert(xValues_.size() >= 2);

        if (x <= xValues_[1])
            return 0;
        else if (x >= xValues_[xValues_.size() - 2])
            return xValues_.size() - 2;
        else {
            // bisection
            size_t segmentIdx = 1;
            size_t upperIdx = xValues_.size() - 2;
            while (segmentIdx + 1 < upperIdx) {
                size_t pivotIdx = (segmentIdx + upperIdx) / 2;
                if (x < xValues_[pivotIdx])
                    upperIdx = pivotIdx;
                else
                    segmentIdx = pivotIdx;
            }

            assert(xValues_[segmentIdx] <= x);
            assert(x <= xValues_[segmentIdx + 1]);
            return segmentIdx;
        }
    }

    Scalar evalDerivative_(OPM_UNUSED Scalar x, int segIdx) const
    {
        Scalar x0 = xValues_[segIdx];
        Scalar x1 = xValues_[segIdx + 1];

        Scalar y0 = yValues_[segIdx];
        Scalar y1 = yValues_[segIdx + 1];

        return (y1 - y0)/(x1 - x0);
    }

    // returns the monotonicity of a segment
    //
    // The return value have the following meaning:
    //
    // 3: function is constant within interval [x0, x1]
    // 1: function is monotonously increasing in the specified interval
    // 9: function is not monotonic in the specified interval
    // -1: function is monotonously decreasing in the specified interval
    int updateMonotonicity_(int i, int &r) const
    {
        if (yValues_[i] < yValues_[i + 1]) {
            // monotonically increasing?
            if (r == 3 || r == 1)
                r = 1;
            else
                r = 0;
            return 1;
        }
        else if (yValues_[i] > yValues_[i + 1]) {
            // monotonically decreasing?
            if (r == 3 || r == -1)
                r = -1;
            else
                r = 0;
            return -1;
        }

        return 3;
    }

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

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

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

    void sortInput_()
    {
        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(xValues_);
        std::sort(idxVector.begin(), idxVector.end(), cmp);

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

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

    void resizeArrays_(size_t nSamples)
    {
        xValues_.resize(nSamples);
        yValues_.resize(nSamples);
    }

    std::vector<Scalar> xValues_;
    std::vector<Scalar> yValues_;
};
} // namespace Opm

#endif