basis.hpp

Go to the documentation of this file.

// Copyright (C) 2021-2022 Petter Nilsson. MIT License.
 
#pragma once
 
#include <cassert>
#include <limits>
 
#include "static_matrix.hpp"
 
SMOOTH_BEGIN_NAMESPACE
 
template<std::size_t K, typename Scalar>
constexpr StaticMatrix<Scalar, 1, K + 1> monomial_derivative(Scalar u, std::size_t p = 0)
{
  StaticMatrix<Scalar, 1, K + 1> ret;
 
  if (p > K) { return ret; }
 
  for (auto i = 0u; i < p; ++i) { ret[0][i] = Scalar(0); }
  Scalar P1      = 1;
  std::size_t P2 = 1;
  for (auto j = 2u; j <= p; ++j) { P2 *= j; }
  ret[0][p] = P1 * static_cast<Scalar>(P2);
  for (auto i = p + 1; i <= K; ++i) {
    P1 *= u;
    P2 *= i;
    P2 /= i - p;
    ret[0][i] = P1 * static_cast<Scalar>(P2);
  }
 
  return ret;
}
 
template<std::size_t K, std::size_t P, typename Scalar>
constexpr StaticMatrix<Scalar, P + 1, K + 1> monomial_derivatives(Scalar u)
{
  StaticMatrix<Scalar, P + 1, K + 1> ret;
  for (auto p = 0u; p <= P; ++p) { ret[p] = monomial_derivative<K, Scalar>(u, p)[0]; }
  return ret;
}
 
namespace detail {
 
template<std::size_t K, typename Scalar = double>
constexpr StaticMatrix<Scalar, K + 1, K + 1> bspline_basis()
{
  StaticMatrix<Scalar, K + 1, K + 1> ret;
  if constexpr (K == 0) {
    ret[0][0] = 1;
    return ret;
  } else {
    const auto coeff_mat_km1 = bspline_basis<K - 1, Scalar>();
    StaticMatrix<Scalar, K + 1, K> low, high;
    StaticMatrix<Scalar, K, K + 1> left, right;
 
    for (auto i = 0u; i != K; ++i) {
      for (auto j = 0u; j != K; ++j) {
        low[i][j]      = coeff_mat_km1[i][j];
        high[i + 1][j] = coeff_mat_km1[i][j];
      }
    }
 
    for (auto k = 0u; k != K; ++k) {
      left[k][k + 1] = static_cast<Scalar>(K - (k + 1)) / static_cast<Scalar>(K);
      left[k][k]     = Scalar(1) - left[k][k + 1];
 
      right[k][k + 1] = Scalar(1) / static_cast<Scalar>(K);
      right[k][k]     = -right[k][k + 1];
    }
 
    return low * left + high * right;
  }
}
 
template<std::size_t K, typename Scalar = double>
constexpr StaticMatrix<Scalar, K + 1, K + 1> bernstein_basis()
{
  StaticMatrix<Scalar, K + 1, K + 1> ret;
  if constexpr (K == 0) {
    ret[0][0] = 1;
    return ret;
  } else {
    constexpr auto coeff_mat_km1 = bernstein_basis<K - 1, Scalar>();
    StaticMatrix<Scalar, K + 1, K> low, high;
    StaticMatrix<Scalar, K, K + 1> left, right;
 
    for (auto i = 0u; i != K; ++i) {
      for (auto j = 0u; j != K; ++j) {
        low[i][j]      = coeff_mat_km1[i][j];
        high[i + 1][j] = coeff_mat_km1[i][j];
      }
    }
 
    for (auto k = 0u; k != K; ++k) {
      left[k][k] = Scalar(1);
 
      right[k][k]     = Scalar(-1);
      right[k][k + 1] = Scalar(1);
    }
 
    return low * left + high * right;
  }
}
 
template<std::size_t K, typename Scalar = double>
constexpr StaticMatrix<Scalar, K + 1, K + 1> hermite_basis()
{
  StaticMatrix<Scalar, K + 1, K + 1> ret;
  ret[0][0] = 1;
 
  if constexpr (K > 0) { ret[1][1] = 2; }
 
  if constexpr (K > 1) {
    for (auto k = 2u; k < K + 1; ++k) {
      for (auto i = 0u; i < k; ++i) { ret[i + 1][k] += 2 * ret[i][k - 1]; }
      for (auto i = 0u; i + 1 < k; ++i) { ret[i][k] -= 2 * (k - 1) * ret[i][k - 2]; }
    }
  }
  return ret;
}
 
template<std::size_t K, typename Scalar = double>
constexpr StaticMatrix<Scalar, K + 1, K + 1> laguerre_basis()
{
  StaticMatrix<Scalar, K + 1, K + 1> ret;
  ret[0][0] = 1;
 
  if constexpr (K > 0) {
    ret[0][1] = 1.;
    ret[1][1] = -1;
  }
 
  if constexpr (K > 1) {
    for (auto k = 2u; k < K + 1; ++k) {
      for (auto i = 0u; i < k; ++i) {
        ret[i][k] += (2 * k - 1) * ret[i][k - 1] / k;
        ret[i + 1][k] -= ret[i][k - 1] / k;
      }
      for (auto i = 0u; i + 1 < k; ++i) { ret[i][k] -= static_cast<double>(k - 1) * ret[i][k - 2] / k; }
    }
  }
  return ret;
}
 
template<std::size_t K, typename Scalar = double>
constexpr StaticMatrix<Scalar, K + 1, K + 1> jacobi_basis(double alpha, double beta)
{
  StaticMatrix<Scalar, K + 1, K + 1> ret;
  ret[0][0] = 1;
 
  if constexpr (K > 0) {
    ret[0][1] = alpha + 1. - (alpha + beta + 2) / 2;
    ret[1][1] = (alpha + beta + 2.) / 2;
  }
 
  if constexpr (K > 1) {
    for (auto k = 2u; k < K + 1; ++k) {
      const auto frac = 1. / ((2 * k) * (k + alpha + beta) * (2 * k + alpha + beta - 2));
      const auto c1   = (2 * k + alpha + beta - 1) * (alpha * alpha - beta * beta);
      const auto c2   = (2 * k + alpha + beta - 1) * (2 * k + alpha + beta) * (2 * k + alpha + beta - 2);
      const auto c3   = 2 * (k + alpha - 1) * (k + beta - 1) * (2 * k + alpha + beta);
      for (auto i = 0u; i < k; ++i) {
        ret[i][k] += c1 * ret[i][k - 1] * frac;
        ret[i + 1][k] += c2 * ret[i][k - 1] * frac;
      }
      for (auto i = 0u; i + 1 < k; ++i) { ret[i][k] -= c3 * ret[i][k - 2] * frac; }
    }
  }
 
  return ret;
}
 
}  // namespace detail
 
enum class PolynomialBasis {
  Bernstein,
  Bspline,
  Chebyshev1st,
  Chebyshev2nd,
  Hermite,
  Laguerre,
  Legendre,
  Monomial,
};
 
template<PolynomialBasis Basis, std::size_t K, typename Scalar = double>
constexpr StaticMatrix<Scalar, K + 1, K + 1> polynomial_basis()
{
  if constexpr (Basis == PolynomialBasis::Monomial) {
    StaticMatrix<Scalar, K + 1, K + 1> ret;
    for (auto k = 0u; k <= K; ++k) { ret[k][k] = Scalar(1); }
    return ret;
  }
  if constexpr (Basis == PolynomialBasis::Bernstein) { return detail::bernstein_basis<K, Scalar>(); }
  if constexpr (Basis == PolynomialBasis::Laguerre) { return detail::laguerre_basis<K, Scalar>(); }
  if constexpr (Basis == PolynomialBasis::Hermite) { return detail::hermite_basis<K, Scalar>(); }
  if constexpr (Basis == PolynomialBasis::Legendre) { return detail::jacobi_basis<K, Scalar>(0, 0); }
  if constexpr (Basis == PolynomialBasis::Chebyshev1st) {
    auto ret       = detail::jacobi_basis<K, Scalar>(-0.5, -0.5);
    const auto fac = monomial_derivative<K>(1.) * ret;
    for (auto k = 0u; k < K + 1; ++k) {
      for (auto r = 0u; r < K + 1; ++r) { ret[r][k] /= fac[0][k]; }
    }
    return ret;
  }
  if constexpr (Basis == PolynomialBasis::Chebyshev2nd) {
    auto ret       = detail::jacobi_basis<K, Scalar>(0.5, 0.5);
    const auto fac = monomial_derivative<K>(1.) * ret;
    for (auto k = 0u; k < K + 1; ++k) {
      for (auto r = 0u; r < K + 1; ++r) { ret[r][k] *= static_cast<Scalar>(k + 1) / fac[0][k]; }
    }
    return ret;
  }
  if constexpr (Basis == PolynomialBasis::Bspline) { return detail::bspline_basis<K, Scalar>(); }
}
 
template<std::size_t K, std::ranges::random_access_range R, typename Scalar = std::ranges::range_value_t<R>>
constexpr StaticMatrix<Scalar, K + 1, K + 1> lagrange_basis(const R & ts)
{
  StaticMatrix<Scalar, K + 1, K + 1> ret;
 
  for (auto row = 0u; row < K + 1; ++row) {
    ret[row][0] = 1;
    for (auto col = 0u; col < K + 1; ++col) {
      if (col != row) {
        const auto row_copy = ret[row];
        ret[row].fill(0);
        for (auto i = 0u; i <= col - (col > row ? 1 : 0); ++i) {
          ret[row][i + 1] += row_copy[i] / (*(ts.begin() + row) - *(ts.begin() + col));
          ret[row][i] -= *(ts.begin() + col) * row_copy[i] / (*(ts.begin() + row) - *(ts.begin() + col));
        }
      }
    }
  }
 
  return ret.transpose();
}
 
template<std::size_t K, std::size_t N, typename Scalar, std::ranges::sized_range R>
constexpr StaticMatrix<Scalar, K + 1, N>
polynomial_basis_derivatives(const StaticMatrix<Scalar, K + 1, K + 1> & B, const R & ts)
{
  assert(std::ranges::size(ts) == N);
 
  StaticMatrix<Scalar, K + 1, N> ret;
 
  for (auto j = 0u; const auto t : ts) {
    const auto U   = monomial_derivative<K>(t, 1);
    const auto U_B = U * B;
    for (auto i = 0u; i < K + 1; ++i) { ret[i][j] = U_B[0][i]; }
    ++j;
  }
 
  return ret;
}
 
template<PolynomialBasis Basis, std::size_t K, typename Scalar = double>
constexpr StaticMatrix<Scalar, K + 1, K + 1> polynomial_cumulative_basis()
{
  auto M = polynomial_basis<Basis, K, Scalar>();
  for (auto i = 0u; i != K + 1; ++i) {
    for (auto j = 0u; j != K; ++j) { M[i][K - 1 - j] += M[i][K - j]; }
  }
  return M;
}
 
template<std::size_t K, std::size_t P, typename Scalar = double>
constexpr StaticMatrix<Scalar, K + 1, K + 1> monomial_integral()
{
  StaticMatrix<Scalar, K + 1, K + 1> ret;
  for (auto i = 0u; i <= K; ++i) {
    for (auto j = i; j <= K; ++j) {
      if (i >= P && j >= P) {
        std::size_t c = 1;
        for (auto _i = i - P + 1; _i <= i; ++_i) { c *= _i; }
        for (auto _j = j - P + 1; _j <= j; ++_j) { c *= _j; }
        ret[i][j] = static_cast<Scalar>(c) / static_cast<Scalar>(i + j - 2 * P + 1);
      } else {
        ret[i][j] = 0;
      }
      ret[j][i] = ret[i][j];
    }
  }
  return ret;
}
 
inline constexpr double integrate_absolute_polynomial(double t0, double t1, double A, double B, double C)
{
  // location of first zero (if any)
  double mid1 = std::numeric_limits<double>::infinity();
  // location of second zero (if any)
  double mid2 = std::numeric_limits<double>::infinity();
 
  if (std::abs(A) < 1e-9 && std::abs(B) > 1e-9) {
    // linear non-constant function
    mid1 = std::clamp(-C / B, t0, t1);
  } else if (std::abs(A) > 1e-9) {
    // quadratic function
    const double res = B * B / (4 * A * A) - C / A;
 
    if (res > 0) {
      mid1 = -B / (2 * A) - std::sqrt(res);
      mid2 = -B / (2 * A) + std::sqrt(res);
    }
  }
 
  const auto integ = [&](double u) { return A * u * u * u / 3 + B * u * u / 2 + C * u; };
 
  const double mid1cl = std::clamp(mid1, t0, t1);
  const double mid2cl = std::clamp(mid2, t0, t1);
 
  return std::abs(integ(t1) - integ(t0) + 2 * integ(mid1cl) - 2 * integ(mid2cl));
}
 
SMOOTH_END_NAMESPACE