MultivariateHermiteInterpolation.h

Go to the documentation of this file.

 
#pragma once
 
#ifndef __UTILS_MULTIVARIATE_HERMITE_INTERPOLATION_H__
#define __UTILS_MULTIVARIATE_HERMITE_INTERPOLATION_H__
 
#include <cassert>
#include <memory>
#include <vector>
 
#include "BivariateHermiteInterpolation.h"
 
namespace Utils
{
 
class MultivariateHermiteInterpolation
{
  public:
    // Samples are assumed to be sorted lexicographically by the coordinates in x.
    // All x[i] must have the same dimension (>= 1).
    void SetSamples(const std::vector<std::vector<double>> &x, const std::vector<double> &y)
    {
        assert(x.size() == y.size());
        assert(!x.empty());
 
        if (x.empty() || y.empty())
            return;
 
        assert(x[0].size() >= 2);
 
        leafInterpolator.reset();
        children.clear();
        xValues.clear();
 
        dimension = x[0].size();
 
        if (dimension == 2)
        {
            // degenerate to the bivariate case
            std::vector<std::vector<double>> xv;
            xv.reserve(x.size());
            for (size_t i = 0; i < x.size(); ++i)
                xv.push_back(x[i]);
 
            leafInterpolator = std::make_unique<BivariateHermiteInterpolation>();
            leafInterpolator->SetTrueInterpolation(trueInterpolation);
            leafInterpolator->SetSamples(xv, y);
            return;
        }
 
        // Group samples by the first coordinate, then recurse on the rest.
        std::vector<std::vector<double>> subX;
        std::vector<double> subY;
 
        auto flushGroup = [&]() {
            if (subX.empty())
                return;
            children.emplace_back(std::make_unique<MultivariateHermiteInterpolation>());
            children.back()->SetTrueInterpolation(trueInterpolation);
            children.back()->SetSamples(subX, subY);
            subX.clear();
            subY.clear();
        };
 
        xValues.push_back(x[0][0]);
        for (size_t i = 0; i < x.size(); ++i)
        {
            assert(x[i].size() == dimension);
 
            if (x[i][0] != xValues.back())
            {
                flushGroup();
                xValues.push_back(x[i][0]);
            }
 
            subX.emplace_back(x[i].begin() + 1, x[i].end());
            subY.push_back(y[i]);
        }
        flushGroup();
    }
 
    double Predict(const std::vector<double> &x) const
    {
        assert(x.size() == dimension);
 
        if (dimension < 2 || x.size() != dimension)
            return 0;
 
        if (leafInterpolator)
            return leafInterpolator->Predict(x);
 
        const std::vector<double> tail(x.begin() + 1, x.end());
 
        std::vector<double> vals;
        vals.reserve(children.size());
        for (size_t i = 0; i < children.size(); ++i)
            vals.push_back(children[i]->Predict(tail));
 
        HermiteInterpolation firstInterpolator;
        firstInterpolator.SetTrueInterpolation(trueInterpolation);
        firstInterpolator.SetSamples(xValues, vals);
 
        const double val = firstInterpolator.Predict(x[0]);
 
        if (trueInterpolation)
            return val;
 
        const auto m = 1E-12;
        if (val < m)
            return m;
 
        return val;
    }
 
    void SetTrueInterpolation(bool reg)
    {
        trueInterpolation = reg;
 
        if (leafInterpolator)
            leafInterpolator->SetTrueInterpolation(reg);
 
        for (auto &child : children)
            if (child)
                child->SetTrueInterpolation(reg);
    }
 
  private:
    // Used when dimension == 2.
    std::unique_ptr<BivariateHermiteInterpolation> leafInterpolator;
 
    // Used when dimension > 2: one child per unique value of the first coordinate.
    std::vector<std::unique_ptr<MultivariateHermiteInterpolation>> children;
    std::vector<double> xValues;
 
    size_t dimension = 0;
 
    bool trueInterpolation = false;
};
 
} // namespace Utils
 
#endif // __UTILS_MULTIVARIATE_HERMITE_INTERPOLATION_H__