Spline.h

/*
 * template<typename T, size_t order>
 * class Spline
 *
 * Represents a spline of dataytype T and order order. The datatype has to
 * fulfill the following requirements:
 *   - comparisons <, <=, >, >=, == and != have to be implemented.
 *   - arithmetic operators + - * /  += -= *= /= have to be implemented
 *   - A pathway must exist, such that integer values of type T can be
 * constructed via static_cast<T>(int) (e. g. via a constructor taking an int).
 * BSplines can be generated via the method generateBspline(...).
 *
 * All methods accessing two splines assume that these splines are defined on
 * the same grid (i.e. that both splines have the same interval boundaries
 * within the intersection of their respective supports). This may also cause
 * problems when splines are constructed by adding up multiple splines. To be
 * safe, make sure that the supports of two splines being added overlap at least
 * in one grid point.
 *
 *
 * ########################################################################
 * The contents of this file is free and unencumbered software released into the
 * public domain. For more information, please refer to <http://unlicense.org/>
 * ########################################################################
 */
 
#ifndef BSPLINE_SPLINE_H
#define BSPLINE_SPLINE_H
 
#include <bspline/exceptions/BSplineException.h>
#include <bspline/internal/misc.h>
#include <bspline/support/Support.h>
 
#include <algorithm>
#include <array>
#include <optional>
#include <type_traits>
#include <vector>
 
namespace bspline {
 
using namespace support;
using namespace bspline::exceptions;
 
template <typename T, size_t order>
class Spline final {
 private:
  static constexpr size_t ARRAY_SIZE = order + 1;
  Support<T> _support;
  std::vector<std::array<T, ARRAY_SIZE>> _coefficients;
 
  std::optional<size_t> findInterval(const T &x) const {
    if (_support.size() < 2 || x > _support.back() || x < _support.front())
      return std::nullopt;  // x is not part of the spline's support
 
    const auto begin = _support.begin();
    const auto it = std::lower_bound(begin, _support.end(), x);
 
    // Type of intervalEndIndex is signed.
    const auto intervalEndIndex = std::distance(begin, it);
    return std::max<decltype(intervalEndIndex)>(0, intervalEndIndex - 1);
  };
 
  void checkValidity(
      const Support<T> &support,
      const std::vector<std::array<T, ARRAY_SIZE>> &coefficients) const {
    const bool isValid =
        (!support.containsIntervals() && coefficients.size() == 0) ||
        (support.size() >= 2 &&
         coefficients.size() == support.numberOfIntervals());
    if (!isValid) throw BSplineException(ErrorCode::INCONSISTENT_DATA);
  };
 
  void checkValidity() const { checkValidity(_support, _coefficients); };
 
  void setData(Support<T> support,
               std::vector<std::array<T, ARRAY_SIZE>> coefficients) {
    checkValidity(support, coefficients);
    _support = std::move(support);
    _coefficients = std::move(coefficients);
  };
 
 public:
  using data_type = T;
 
  static constexpr size_t spline_order = order;
 
  Spline(Support<T> support,
         std::vector<std::array<T, ARRAY_SIZE>> coefficients)
      : _support(std::move(support)), _coefficients(std::move(coefficients)) {
    checkValidity(_support, _coefficients);
  };
 
  explicit Spline(Grid<T> grid)
      : Spline(Support<T>::createEmpty(std::move(grid)), {}){};
 
  const Support<T> &getSupport() const noexcept {
    DURING_TEST_CHECK_VALIDITY();
    return _support;
  };
 
  const std::vector<std::array<T, ARRAY_SIZE>> &getCoefficients()
      const noexcept {
    DURING_TEST_CHECK_VALIDITY();
    return _coefficients;
  };
 
  T operator()(const T &x) const {
    DURING_TEST_CHECK_VALIDITY();
    const auto intervalIndex = findInterval(x);
 
    if (!intervalIndex) return static_cast<T>(0);
 
    const T xm = (_support[*intervalIndex + 1] + _support[*intervalIndex]) /
                 static_cast<T>(2);
 
    return internal::evaluateInterval(x, _coefficients[*intervalIndex], xm);
  };
 
  const T &front() const {
    DURING_TEST_CHECK_VALIDITY();
    return _support.front();
  };
 
  const T &back() const {
    DURING_TEST_CHECK_VALIDITY();
    return _support.back();
  };
 
  template <size_t order2>
  bool checkOverlap(const Spline<T, order2> &m2) const {
    DURING_TEST_CHECK_VALIDITY();
    if (!_support.containsIntervals() || !m2.getSupport().containsIntervals())
      return false;
    const bool isNotOverlapping = m2.getSupport().back() <= _support.front() ||
                                  m2.getSupport().front() >= _support.back();
    return !isNotOverlapping;
  }
 
  bool isZero() const {
    DURING_TEST_CHECK_VALIDITY();
    static const T ZERO = static_cast<T>(0);
    if (!_support.containsIntervals()) return true;
    for (const auto &cs : _coefficients) {
      for (const auto &c : cs) {
        if (c != ZERO) return false;
      }
    }
    return true;
  };
 
  // ######################## Operator definitions ########################
  // ######################################################################
 
  Spline<T, order> operator/(const T &d) const {
    DURING_TEST_CHECK_VALIDITY();
    return (*this) * (static_cast<T>(1) / d);
  };
 
  Spline<T, order> operator*(const T &d) const {
    DURING_TEST_CHECK_VALIDITY();
    Spline<T, order> ret(*this);
    for (auto &cs : ret._coefficients) {
      for (auto &c : cs) {
        c *= d;
      }
    }
    return ret;
  };
 
  Spline<T, order> &operator*=(const T &d) {
    DURING_TEST_CHECK_VALIDITY();
    for (auto &cs : _coefficients) {
      for (auto &c : cs) {
        c *= d;
      }
    }
    return *this;
  };
 
  Spline<T, order> &operator/=(const T &d) {
    DURING_TEST_CHECK_VALIDITY();
    (*this) *= (static_cast<T>(1) / d);
    return *this;
  };
 
  Spline<T, order> operator-() const {
    DURING_TEST_CHECK_VALIDITY();
    return (*this) * static_cast<T>(-1);
  };
 
  template <size_t ordera>
  Spline<T, order> &operator=(const Spline<T, ordera> &a) {
    DURING_TEST_CHECK_VALIDITY();
    // The case ordera == order should be handled by the default assignment
    // operator which is automatically generated.
    static_assert(
        ordera < order,
        "The assignment operator is only defined if the order of the rhs "
        "spline is lower than or equal to that of the lhs spline.");
 
    std::vector<std::array<T, ARRAY_SIZE>> ncoefficients(
        a.getCoefficients().size(),
        internal::make_array<T, ARRAY_SIZE>(static_cast<T>(0)));
    for (size_t i = 0; i < a.getCoefficients().size(); i++) {
      const auto &coeffsi = a.getCoefficients()[i];
      auto &ncoeffsi = ncoefficients[i];
      for (size_t j = 0; j < coeffsi.size(); j++) ncoeffsi[j] = coeffsi[j];
    }
    setData(a.getSupport(), std::move(ncoefficients));
    return *this;
  }
 
  template <size_t ordera>
  Spline<T, order + ordera> operator*(const Spline<T, ordera> &a) const {
    DURING_TEST_CHECK_VALIDITY();
    static constexpr size_t NEW_ORDER = order + ordera;
    static constexpr size_t NEW_ARRAY_SIZE = NEW_ORDER + 1;
 
    // Will also check whether the two grids are equivalent.
    Support newSupport = _support.calcIntersection(a.getSupport());
    const size_t nintervals = newSupport.numberOfIntervals();
 
    if (nintervals == 0)
      return Spline<T, NEW_ORDER>(std::move(newSupport), {});  // No overlap
 
    std::vector<std::array<T, NEW_ARRAY_SIZE>> newCoefficients(
        nintervals, internal::make_array<T, NEW_ARRAY_SIZE>(static_cast<T>(0)));
 
    for (size_t i = 0; i < nintervals; i++) {
      const auto ai = newSupport.absoluteFromRelative(i);
 
      auto thisRelIndex = _support.intervalIndexFromAbsolute(ai).value();
      auto aRelIndex = a.getSupport().intervalIndexFromAbsolute(ai).value();
 
      const auto &thiscoeffs = _coefficients[thisRelIndex];
      const auto &acoeffs = a.getCoefficients()[aRelIndex];
      auto &coeffsi = newCoefficients[i];
 
      for (size_t j = 0; j < order + 1; j++) {
        for (size_t k = 0; k < ordera + 1; k++) {
          coeffsi[j + k] += thiscoeffs[j] * acoeffs[k];
        }
      }
    }
    return Spline<T, NEW_ORDER>(std::move(newSupport),
                                std::move(newCoefficients));
  }
 
  template <size_t ordera>
  Spline<T, std::max(order, ordera)> operator+(
      const Spline<T, ordera> &a) const {
    DURING_TEST_CHECK_VALIDITY();
    static constexpr size_t NEW_ORDER = std::max(order, ordera);
    static constexpr size_t NEW_ARRAY_SIZE = NEW_ORDER + 1;
 
    // Will also check whether the two grids are equivalent.
    Support newSupport = _support.calcUnion(a.getSupport());
    const size_t nintervals = newSupport.numberOfIntervals();
 
    std::vector<std::array<T, NEW_ARRAY_SIZE>> ncoefficients(nintervals);
 
    for (size_t i = 0; i < nintervals; i++) {
      const auto absIndex = newSupport.absoluteFromRelative(i);
      const auto thisRelIndex = _support.intervalIndexFromAbsolute(absIndex);
      const auto aRelIndex = a.getSupport().intervalIndexFromAbsolute(absIndex);
 
      if (thisRelIndex && !aRelIndex) {
        ncoefficients[i] =
            internal::changearraysize<T, order + 1, NEW_ARRAY_SIZE>(
                _coefficients[*thisRelIndex]);
      } else if (aRelIndex && !thisRelIndex) {
        ncoefficients[i] =
            internal::changearraysize<T, ordera + 1, NEW_ARRAY_SIZE>(
                a.getCoefficients()[*aRelIndex]);
      } else if (thisRelIndex && aRelIndex) {
        ncoefficients[i] = internal::add<T, ordera + 1, order + 1>(
            a.getCoefficients()[*aRelIndex], _coefficients[*thisRelIndex]);
      } else {
        ncoefficients[i] =
            internal::make_array<T, NEW_ARRAY_SIZE>(static_cast<T>(0));
      }
    }
    return Spline<T, NEW_ORDER>(std::move(newSupport),
                                std::move(ncoefficients));
  }
 
  template <size_t ordera>
  Spline<T, order> &operator+=(const Spline<T, ordera> &a) {
    DURING_TEST_CHECK_VALIDITY();
    static_assert(
        ordera <= order,
        "The operators += and -= are only defined if the order of the rhs "
        "spline is lower than or equal to that of the lhs spline.");
    (*this) = (*this) + a;
    return *this;
  }
 
  template <size_t ordera>
  Spline<T, order> &operator-=(const Spline<T, ordera> &a) {
    (*this) += (static_cast<T>(-1) * a);
    return *this;
  }
 
  template <size_t ordera>
  Spline<T, std::max(order, ordera)> operator-(
      const Spline<T, ordera> &a) const {
    return (*this) + (static_cast<T>(-1) * a);
  }
 
  bool operator==(const Spline &other) const {
    return _support == other._support && _coefficients == other._coefficients;
  }
 
  bool operator!=(const Spline &other) const { return !(*this == other); }
};  // class Spline
 
template <typename T, size_t ARRAY_SIZE>
Spline(Support<T> support, std::vector<std::array<T, ARRAY_SIZE>> coefficients)
    -> Spline<T, ARRAY_SIZE - 1>;
 
// ################### End of defintion of Spline class ###################
// ########################################################################
 
template <typename T, size_t order>
inline Spline<T, order> operator*(const T &d, const Spline<T, order> &b) {
  return b * d;
}
 
template <typename CoeffIter, typename SplineIter>
auto linearCombination(CoeffIter coeffsBegin, CoeffIter coeffsEnd,
                       SplineIter splinesBegin, SplineIter splinesEnd) {
  // The spline data type.
  using Spline = typename std::remove_cv_t<
      typename std::iterator_traits<SplineIter>::value_type>;
 
  // The coefficient data type.
  using T = typename std::remove_cv_t<
      typename std::iterator_traits<CoeffIter>::value_type>;
 
  // The order of the spline.
  constexpr size_t order = Spline::spline_order;
 
  // Check, the data type of the spline and the coefficients are consistent.
  static_assert(
      std::is_same<typename Spline::data_type, T>::value,
      "Coefficients must be of the same type as the data type of the spline.");
 
  {
    // The number of coefficients and splines.
    const auto coeffsSize = std::distance(coeffsBegin, coeffsEnd);
    const auto splinesSize = std::distance(splinesBegin, splinesEnd);
 
    // Check, the data is consistent.
    if (coeffsSize != splinesSize) {
      throw BSplineException(
          ErrorCode::INCONSISTENT_DATA,
          "The number of coefficients and splines must coincide.");
    }
 
    // std::distance() may return negative values.
    if (coeffsSize <= 0) {
      throw BSplineException(
          ErrorCode::MISSING_DATA,
          "The number of coefficients and splines may not be zero.");
    }
  }
 
  const auto &support0 = splinesBegin->getSupport();
  std::optional<size_t> startIndex;
  std::optional<size_t> endIndex;
 
  // Determine the union of all the splines' supports.
  for (auto it = splinesBegin; it < splinesEnd; it++) {
    if (!it->getSupport().hasSameGrid(support0)) {
      throw BSplineException(ErrorCode::DIFFERING_GRIDS);
    }
 
    if (!it->getSupport().empty()) {
      const size_t si = it->getSupport().getStartIndex();
      if (!startIndex || si < *startIndex) {
        startIndex = si;
      }
 
      const size_t ei = it->getSupport().getEndIndex();
      if (!endIndex || ei > *endIndex) {
        endIndex = ei;
      }
    }
  }
 
  // Set up support and coefficients vector for the returned spline.
  Support newSupport(support0.getGrid(), startIndex.value_or(0),
                     endIndex.value_or(0));
  std::vector<std::array<T, order + 1>> newCoefficients(
      newSupport.numberOfIntervals(),
      internal::make_array<T, order + 1>(static_cast<T>(0)));
 
  auto coeffIt = coeffsBegin;
  auto splineIt = splinesBegin;
  while (splineIt < splinesEnd) {
    // Get spline and coefficient.
    const auto &spline = *splineIt;
    const T coeff = *coeffIt;
 
    for (size_t j = 0; j < spline.getSupport().numberOfIntervals(); j++) {
      // Index of the interval relative to the global grid.
      const size_t absoluteIndex = spline.getSupport().absoluteFromRelative(j);
 
      // Index of the interval relative to the newSupport.
      const size_t newSupportIndex =
          newSupport.intervalIndexFromAbsolute(absoluteIndex).value();
 
      const std::array<T, order + 1> &splineCoeffs =
          spline.getCoefficients().at(j);
      std::array<T, order + 1> &newCoeffs = newCoefficients.at(newSupportIndex);
 
      for (size_t k = 0; k < order + 1; k++) {
        newCoeffs[k] += coeff * splineCoeffs[k];
      }
    }
 
    splineIt++;
    coeffIt++;
  }
 
  return Spline(std::move(newSupport), std::move(newCoefficients));
}
 
template <typename CoeffCollection, typename SplineCollection>
auto linearCombination(const CoeffCollection &coeffs,
                       const SplineCollection &splines) {
  return linearCombination(coeffs.begin(), coeffs.end(), splines.begin(),
                           splines.end());
}
 
#ifndef BSPLINE_DOXYGEN_IGNORE
template <typename S>
struct is_spline : std::false_type {};
 
template <typename T, size_t order>
struct is_spline<Spline<T, order>> : std::true_type {};
#endif  // BSPLINE_DOXYGEN_IGNORE
 
template <typename S>
inline constexpr bool is_spline_v = is_spline<S>::value;
 
}  // namespace bspline
#endif  // BSPLINE_SPLINE_H