diff --git a/src/pyinterp/core/include/pyinterp/detail/interpolation/akima.hpp b/src/pyinterp/core/include/pyinterp/detail/interpolation/akima.hpp new file mode 100644 index 00000000..6e027d01 --- /dev/null +++ b/src/pyinterp/core/include/pyinterp/detail/interpolation/akima.hpp @@ -0,0 +1,135 @@ +// Copyright (c) 2024 CNES +// +// All rights reserved. Use of this source code is governed by a +// BSD-style license that can be found in the LICENSE file. +#pragma once +#include "pyinterp/detail/interpolation/interpolator_1d.hpp" + +namespace pyinterp::detail::interpolation { + +/// Akima interpolation +template +class Akima : public Interpolator1D { + public: + using Interpolator1D::Interpolator1D; + using Interpolator1D::operator(); + using Interpolator1D::derivative; + + /// Returns the minimum number of points required for the interpolation. + auto min_size() const -> Eigen::Index override { return 5; } + + private: + Vector m_{}; + Vector s_{}; + + /// Compute the boundary conditions. + virtual auto boundary_condition(T* m, const size_t size) -> void { + m[-2] = 3 * m[0] - 2 * m[1]; + m[-1] = 2 * m[0] - m[1]; + m[size - 1] = 2 * m[size - 2] - m[size - 3]; + m[size] = 3 * m[size - 2] - 2 * m[size - 3]; + } + + /// @brief Compute the coefficients of the interpolation + /// @param xa X-coordinates of the data points. + /// @param ya Y-coordinates of the data points. + auto compute_coefficients(const Eigen::Ref>& xa, + const Eigen::Ref>& ya) + -> void override; + + /// Interpolation + /// @param xa X-coordinates of the data points. + /// @param ya Y-coordinates of the data points. + /// @param x The point where the interpolation must be calculated. + /// @param index The index of the last point found in the search. + /// @return The interpolated value at the point x. + auto operator()(const Eigen::Ref>& xa, + const Eigen::Ref>& ya, const T& x, + Eigen::Index* index) const -> T override; + + /// @brief Returns the derivative of the interpolation function at the point + /// x. + /// @param xa X-coordinates of the data points. + /// @param ya Y-coordinates of the data points. + /// @param x The point where the derivative must be calculated. + /// @param index The index of the last point found in the search. + /// @return The derivative of the interpolation function at the point x. + auto derivative(const Eigen::Ref>& xa, + const Eigen::Ref>& ya, const T& x, + Eigen::Index* index) const -> T override; +}; + +template +auto Akima::compute_coefficients(const Eigen::Ref>& xa, + const Eigen::Ref>& ya) + -> void { + Interpolator1D::compute_coefficients(xa, ya); + auto size = xa.size(); + if (m_.size() < size + 4) { + m_.resize(size + 4); + s_.resize(size); + } + + // m contains the slopes of the lines between the points. Two extra points + // are added at the beginning and end to handle the boundary conditions. + auto* m = m_.data() + 2; + for (Eigen::Index ix = 0; ix < size - 1; ++ix) { + m[ix] = (ya[ix + 1] - ya[ix]) / (xa[ix + 1] - xa[ix]); + } + + boundary_condition(m, size); + + // Compute the spline slopes of the lines between the points. + for (Eigen::Index ix = 2; ix < size - 2; ++ix) { + auto denominator = + std::abs(m[ix + 1] - m[ix]) + std::abs(m[ix - 1] - m[ix - 2]); + if (denominator != 0) { + s_(ix) = std::abs(m[ix + 1] - m[ix]) * m[ix - 1] + + std::abs(m[ix - 1] - m[ix - 2]) * m[ix]; + s_(ix) /= denominator; + } else { + s_(ix) = (m[ix - 1] + m[ix]) * 0.5; + } + } + s_(0) = m[0]; + s_(1) = (m[0] + m[2]) * 0.5; + s_(size - 2) = (m[size - 3] + m[size - 1]) * 0.5; + s_(size - 1) = m[size - 1]; +} + +template +auto Akima::operator()(const Eigen::Ref>& xa, + const Eigen::Ref>& ya, const T& x, + Eigen::Index* index) const -> T { + auto search = this->search(xa, x, index); + if (!search) { + throw std::numeric_limits::quiet_NaN(); + } + auto [i0, i1] = *search; + const auto dx = xa(i1) - xa(i0); + const auto h = x - xa(i0); + const auto ai = ya(i0); + const auto bi = s_[i0]; + const auto ci = (3 * m_[i0] - 2 * s_[i0] - s_[i1]) / dx; + const auto di = (s_[i0] + s_[i1] - 2 * m_[i0]) / (dx * dx); + return ai + h * (bi + h * (ci + h * di)); +} + +template +auto Akima::derivative(const Eigen::Ref>& xa, + const Eigen::Ref>& ya, const T& x, + Eigen::Index* index) const -> T { + auto search = this->search(xa, x, index); + if (!search) { + throw std::numeric_limits::quiet_NaN(); + } + auto [i0, i1] = *search; + const auto dx = xa(i1) - xa(i0); + const auto h = x - xa(i0); + const auto bi = s_[i0]; + const auto ci = (3 * m_[i0] - 2 * s_[i0] - s_[i1]) / dx; + const auto di = (s_[i0] + s_[i1] - 2 * m_[i0]) / (dx * dx); + return bi + h * (2 * ci + h * 3 * di); +} + +} // namespace pyinterp::detail::interpolation diff --git a/src/pyinterp/core/include/pyinterp/detail/interpolation/akima_periodic.hpp b/src/pyinterp/core/include/pyinterp/detail/interpolation/akima_periodic.hpp new file mode 100644 index 00000000..bdcbe723 --- /dev/null +++ b/src/pyinterp/core/include/pyinterp/detail/interpolation/akima_periodic.hpp @@ -0,0 +1,26 @@ +// Copyright (c) 2024 CNES +// +// All rights reserved. Use of this source code is governed by a +// BSD-style license that can be found in the LICENSE file. +#pragma once +#include "pyinterp/detail/interpolation/akima.hpp" + +namespace pyinterp::detail::interpolation { + +/// Akima periodic interpolation +template +class AkimaPeriodic : public Akima { + public: + using Akima::Akima; + + private: + /// Compute the boundary conditions. + auto boundary_condition(T* m, const size_t size) -> void override { + m[-2] = m[size - 3]; + m[-1] = m[size - 2]; + m[size - 1] = m[0]; + m[size] = m[1]; + } +}; + +} // namespace pyinterp::detail::interpolation diff --git a/src/pyinterp/core/include/pyinterp/detail/interpolation/bicubic.hpp b/src/pyinterp/core/include/pyinterp/detail/interpolation/bicubic.hpp new file mode 100644 index 00000000..693a4a45 --- /dev/null +++ b/src/pyinterp/core/include/pyinterp/detail/interpolation/bicubic.hpp @@ -0,0 +1,181 @@ +// Copyright (c) 2024 CNES +// +// All rights reserved. Use of this source code is governed by a +// BSD-style license that can be found in the LICENSE file. +#pragma once + +#include "pyinterp/detail/interpolation/cspline.hpp" +#include "pyinterp/detail/interpolation/interpolator_2d.hpp" + +namespace pyinterp::detail::interpolation { + +/// Bicubic interpolation +/// @tparam T type of the data +template +class Bicubic : public Interpolator2D { + public: + using Interpolator2D::Interpolator2D; + using Interpolator2D::operator(); + + /// Returns the minimum number of points required for the interpolation. + auto min_size() const -> Eigen::Index override { return 2; } + + private: + /// Interpolation + /// @param xa X-coordinates of the data points. + /// @param ya Y-coordinates of the data points. + /// @param za Z-values of the data points. + /// @param x The point where the interpolation must be calculated. + /// @param y The point where the interpolation must be calculated. + /// @param ix The index of the last point found in the search. + /// @param jx The index of the last point found in the search. + /// @return The interpolated value at the coordinates x, y. + auto operator()(const Eigen::Ref> &xa, + const Eigen::Ref> &ya, + const Eigen::Ref> &za, const T &x, const T &y, + Eigen::Index *ix, Eigen::Index *jx) const -> T override; + + /// Compute the coefficients of the bicubic interpolation + /// @param xa X-coordinates of the data points. + /// @param ya Y-coordinates of the data points. + /// @param za Z-values of the data points. + auto compute_coefficients(const Eigen::Ref> &xa, + const Eigen::Ref> &ya, + const Eigen::Ref> &za) + -> void override; + + Matrix zx_{}; + Matrix zy_{}; + Matrix zxy_{}; +}; + +template +auto Bicubic::compute_coefficients(const Eigen::Ref> &xa, + const Eigen::Ref> &ya, + const Eigen::Ref> &za) + -> void { + Interpolator2D::compute_coefficients(xa, ya, za); + auto xsize = xa.size(); + auto ysize = ya.size(); + + if (zx_.rows() != xsize || zx_.cols() != ysize) { + zx_ = Matrix(xsize, ysize); + zy_ = Matrix(xsize, ysize); + zxy_ = Matrix(xsize, ysize); + } + + auto spline = CSpline(); + for (Eigen::Index j = 0; j < ysize; ++j) { + zx_.col(j) = spline.derivative(xa, za.col(j), xa); + } + for (Eigen::Index i = 0; i < xsize; ++i) { + zy_.row(i) = spline.derivative(ya, za.row(i), ya); + } + for (Eigen::Index j = 0; j < ysize; ++j) { + zxy_.col(j) = spline.derivative(xa, zy_.col(j), xa); + } +} + +template +auto Bicubic::operator()(const Eigen::Ref> &xa, + const Eigen::Ref> &ya, + const Eigen::Ref> &za, const T &x, + const T &y, Eigen::Index *ix, + Eigen::Index *jx) const -> T { + auto search_x = this->search(xa, x, ix); + auto search_y = this->search(ya, y, jx); + if (!search_x || !search_y) { + throw std::numeric_limits::quiet_NaN(); + } + const auto [i0, i1] = *search_x; + const auto [j0, j1] = *search_y; + const auto x0 = xa(i0); + const auto x1 = xa(i1); + const auto y0 = ya(j0); + const auto y1 = ya(j1); + const auto z00 = za(i0, j0); + const auto z01 = za(i0, j1); + const auto z10 = za(i1, j0); + const auto z11 = za(i1, j1); + const auto dx = x1 - x0; + const auto dy = y1 - y0; + const auto t = (x - x0) / dx; + const auto u = (y - y0) / dy; + const auto zx00 = zx_(i0, j0) * dx; + const auto zx01 = zx_(i0, j1) * dx; + const auto zx10 = zx_(i1, j0) * dx; + const auto zx11 = zx_(i1, j1) * dx; + const auto zy00 = zy_(i0, j0) * dy; + const auto zy01 = zy_(i0, j1) * dy; + const auto zy10 = zy_(i1, j0) * dy; + const auto zy11 = zy_(i1, j1) * dy; + const auto zxy00 = zxy_(i0, j0) * (dx * dy); + const auto zxy01 = zxy_(i0, j1) * (dx * dy); + const auto zxy10 = zxy_(i1, j0) * (dx * dy); + const auto zxy11 = zxy_(i1, j1) * (dx * dy); + const auto t0 = 1; + const auto t1 = t; + const auto t2 = t * t; + const auto t3 = t * t2; + const auto u0 = 1; + const auto u1 = u; + const auto u2 = u * u; + const auto u3 = u * u2; + + auto v = z00; + auto z = v * t0 * u0; + + v = zy00; + z += v * t0 * u1; + + v = 3 * (-z00 + z01) - 2 * zy00 - zy01; + z += v * t0 * u2; + + v = 2 * (z00 - z01) + zy00 + zy01; + z += v * t0 * u3; + + v = zx00; + z += v * t1 * u0; + + v = zxy00; + z += v * t1 * u1; + + v = 3 * (-zx00 + zx01) - 2 * zxy00 - zxy01; + z += v * t1 * u2; + + v = 2 * (zx00 - zx01) + zxy00 + zxy01; + z += v * t1 * u3; + + v = 3 * (-z00 + z10) - 2 * zx00 - zx10; + z += v * t2 * u0; + + v = 3 * (-zy00 + zy10) - 2 * zxy00 - zxy10; + z += v * t2 * u1; + + v = 9 * (z00 - z10 + z11 - z01) + 6 * (zx00 - zx01 + zy00 - zy10) + + 3 * (zx10 - zx11 - zy11 + zy01) + 4 * zxy00 + 2 * (zxy10 + zxy01) + zxy11; + z += v * t2 * u2; + + v = 6 * (-z00 + z10 - z11 + z01) + 4 * (-zx00 + zx01) + + 3 * (-zy00 + zy10 + zy11 - zy01) + 2 * (-zx10 + zx11 - zxy00 - zxy01) - + zxy10 - zxy11; + z += v * t2 * u3; + + v = 2 * (z00 - z10) + zx00 + zx10; + z += v * t3 * u0; + + v = 2 * (zy00 - zy10) + zxy00 + zxy10; + z += v * t3 * u1; + + v = 6 * (-z00 + z10 - z11 + z01) + 3 * (-zx00 - zx10 + zx11 + zx01) + + 4 * (-zy00 + zy10) + 2 * (zy11 - zy01 - zxy00 - zxy10) - zxy11 - zxy01; + z += v * t3 * u2; + + v = 4 * (z00 - z10 + z11 - z01) + + 2 * (zx00 + zx10 - zx11 - zx01 + zy00 - zy10 - zy11 + zy01) + zxy00 + + zxy10 + zxy11 + zxy01; + z += v * t3 * u3; + return z; +} + +} // namespace pyinterp::detail::interpolation diff --git a/src/pyinterp/core/include/pyinterp/detail/interpolation/bilinear.hpp b/src/pyinterp/core/include/pyinterp/detail/interpolation/bilinear.hpp new file mode 100644 index 00000000..0df05c8a --- /dev/null +++ b/src/pyinterp/core/include/pyinterp/detail/interpolation/bilinear.hpp @@ -0,0 +1,60 @@ +// Copyright (c) 2024 CNES +// +// All rights reserved. Use of this source code is governed by a +// BSD-style license that can be found in the LICENSE file. +#pragma once +#include "pyinterp/detail/interpolation/interpolator_2d.hpp" + +namespace pyinterp::detail::interpolation { + +/// Bilinear interpolation +/// @tparam T type of the data +template +class Bilinear : public Interpolator2D { + public: + using Interpolator2D::Interpolator2D; + using Interpolator2D::operator(); + + /// Returns the minimum number of points required for the interpolation. + auto min_size() const -> Eigen::Index override { return 2; } + + private: + /// Interpolation + /// @param xa X-coordinates of the data points. + /// @param ya Y-coordinates of the data points. + /// @param za Z-values of the data points. + /// @param x The point where the interpolation must be calculated. + /// @param y The point where the interpolation must be calculated. + /// @param ix The index of the last point found in the search. + auto operator()(const Eigen::Ref> &xa, + const Eigen::Ref> &ya, + const Eigen::Ref> &za, const T &x, const T &y, + Eigen::Index *ix, Eigen::Index *jx) const -> T override; +}; + +template +auto Bilinear::operator()(const Eigen::Ref> &xa, + const Eigen::Ref> &ya, + const Eigen::Ref> &za, const T &x, + const T &y, Eigen::Index *ix, + Eigen::Index *jx) const -> T { + auto search_x = this->search(xa, x, ix); + auto search_y = this->search(ya, y, jx); + if (!search_x || !search_y) { + throw std::numeric_limits::quiet_NaN(); + } + auto [i0, i1] = *search_x; + auto [j0, j1] = *search_y; + auto x0 = xa[i0]; + auto x1 = xa[i1]; + auto y0 = ya[j0]; + auto y1 = ya[j1]; + auto dx = x1 - x0; + auto dy = y1 - y0; + auto t = (x - x0) / dx; + auto u = (y - y0) / dy; + return (T(1) - t) * (T(1) - u) * za(i0, j0) + t * (T(1) - u) * za(i1, j0) + + (T(1) - t) * u * za(i0, j1) + t * u * za(i1, j1); +} + +} // namespace pyinterp::detail::interpolation diff --git a/src/pyinterp/core/include/pyinterp/detail/interpolation/cspline.hpp b/src/pyinterp/core/include/pyinterp/detail/interpolation/cspline.hpp new file mode 100644 index 00000000..6620cc36 --- /dev/null +++ b/src/pyinterp/core/include/pyinterp/detail/interpolation/cspline.hpp @@ -0,0 +1,156 @@ +// Copyright (c) 2024 CNES +// +// All rights reserved. Use of this source code is governed by a +// BSD-style license that can be found in the LICENSE file. +#pragma once +#include + +#include "pyinterp/detail/interpolation/interpolator_1d.hpp" + +namespace pyinterp::detail::interpolation { + +/// Coefficients of the cubic spline interpolation +template +class CSplineCoefficients { + public: + /// Constructor + /// @param c0 The first derivative at the first point + /// @param c1 The first derivative at the last point + constexpr CSplineCoefficients(const T &c0, const T &c1) : c0_(c0), c1_(c1){}; + + /// Compute the coefficients of the cubic spline interpolation + /// @param dx The distance between the two points + /// @param dy The difference between the two points + /// @return The coefficients of the cubic spline interpolation + constexpr auto operator()(const T &dx, const T &dy) const + -> std::tuple { + return {(dy / dx) - dx * (c1_ + 2 * c0_) / 3, c0_, (c1_ - c0_) / (3 * dx)}; + } + + private: + T c0_; + T c1_; +}; + +/// Cubic spline interpolation +template +class CSpline : public Interpolator1D { + public: + using Interpolator1D::Interpolator1D; + using Interpolator1D::operator(); + using Interpolator1D::derivative; + + /// Returns the minimum number of points required for the interpolation. + auto min_size() const -> Eigen::Index override { return 3; } + + private: + /// Interpolation + /// @param xa X-coordinates of the data points. + /// @param ya Y-coordinates of the data points. + /// @param x The point where the interpolation must be calculated. + /// @param i The index of the last point found in the search. + /// @return The interpolated value at the point x. + auto operator()(const Eigen::Ref> &xa, + const Eigen::Ref> &ya, const T &x, + Eigen::Index *i) const -> T override; + + /// @brief Returns the derivative of the interpolation function at the point + /// x. + /// @param xa X-coordinates of the data points. + /// @param ya Y-coordinates of the data points. + /// @param x The point where the derivative must be calculated. + /// @param i The index of the last point found in the search. + /// @return The derivative of the interpolation function at the point x. + auto derivative(const Eigen::Ref> &xa, + const Eigen::Ref> &ya, const T &x, + Eigen::Index *i) const -> T override; + + protected: + /// @brief Compute the coefficients of the interpolation + /// @param xa X-coordinates of the data points. + /// @param ya Y-coordinates of the data points. + auto compute_coefficients(const Eigen::Ref> &xa, + const Eigen::Ref> &ya) + -> void override; + + Matrix A_; + Vector b_; + Vector x_; +}; + +template +auto CSpline::compute_coefficients(const Eigen::Ref> &xa, + const Eigen::Ref> &ya) + -> void { + Interpolator1D::compute_coefficients(xa, ya); + auto size = xa.size(); + if (x_.size() != size) { + A_.resize(size - 2, size - 2); + b_.resize(size - 2); + x_.resize(size); + A_.setZero(); + } + + for (auto i = 0; i < size - 2; i++) { + const auto h_i0 = xa[i + 1] - xa[i]; + const auto h_i1 = xa[i + 2] - xa[i + 1]; + const auto y_i0 = ya[i + 1] - ya[i]; + const auto y_i1 = ya[i + 2] - ya[i + 1]; + const auto g_i0 = (h_i0 != 0) ? 1 / h_i0 : 0; + const auto g_i1 = (h_i1 != 0) ? 1 / h_i1 : 0; + if (i > 0) { + A_(i, i - 1) = h_i0; + } + A_(i, i) = 2 * (h_i0 + h_i1); + if (i < size - 3) { + A_(i, i + 1) = h_i1; + } + b_(i) = 3 * (y_i1 * g_i1 - y_i0 * g_i0); + } + x_.segment(1, size - 2) = A_.fullPivLu().solve(b_); + x_(0) = x_(size - 1) = 0; +} + +template +auto CSpline::operator()(const Eigen::Ref> &xa, + const Eigen::Ref> &ya, const T &x, + Eigen::Index *i) const -> T { + auto where = this->search(xa, x, i); + if (!where) { + return std::numeric_limits::quiet_NaN(); + } + auto x_lo = xa(where->first); + auto y_lo = ya(where->first); + auto x_hi = xa(where->second); + auto y_hi = ya(where->second); + auto dx = x_hi - x_lo; + auto dy = y_hi - y_lo; + auto h = x - x_lo; + + auto [b, c, d] = + CSplineCoefficients(x_(where->first), x_(where->second))(dx, dy); + return y_lo + h * (b + h * (c + h * d)); +} + +template +auto CSpline::derivative(const Eigen::Ref> &xa, + const Eigen::Ref> &ya, const T &x, + Eigen::Index *i) const -> T { + auto where = this->search(xa, x, i); + if (!where) { + return std::numeric_limits::quiet_NaN(); + } + auto x_lo = xa(where->first); + auto y_lo = ya(where->first); + auto x_hi = xa(where->second); + auto y_hi = ya(where->second); + auto dx = x_hi - x_lo; + auto dy = y_hi - y_lo; + auto h = x - x_lo; + + auto [b, c, d] = + CSplineCoefficients(x_(where->first), x_(where->second))(dx, dy); + return b + h * (2 * c + h * 3 * d); +} + +} // namespace pyinterp::detail::interpolation diff --git a/src/pyinterp/core/include/pyinterp/detail/interpolation/cspline_periodic.hpp b/src/pyinterp/core/include/pyinterp/detail/interpolation/cspline_periodic.hpp new file mode 100644 index 00000000..bdd570aa --- /dev/null +++ b/src/pyinterp/core/include/pyinterp/detail/interpolation/cspline_periodic.hpp @@ -0,0 +1,71 @@ +// Copyright (c) 2024 CNES +// +// All rights reserved. Use of this source code is governed by a +// BSD-style license that can be found in the LICENSE file. +#pragma once +#include + +#include "pyinterp/detail/interpolation/cspline.hpp" + +namespace pyinterp::detail::interpolation { + +/// Periodic cubic spline interpolation +template +class CSplinePeriodic : public CSpline { + public: + using CSpline::CSpline; + + /// Returns the minimum number of points required for the interpolation. + auto min_size() const -> Eigen::Index override { return 2; } + + private: + /// Compute the coefficients of the interpolation + /// @param xa X-coordinates of the data points. + /// @param ya Y-coordinates of the data points. + auto compute_coefficients(const Eigen::Ref> &xa, + const Eigen::Ref> &ya) + -> void override; +}; + +template +auto CSplinePeriodic::compute_coefficients( + const Eigen::Ref> &xa, + const Eigen::Ref> &ya) -> void { + Interpolator1D::compute_coefficients(xa, ya); + auto size = xa.size(); + if (this->x_.size() < size) { + this->A_.resize(size - 1, size - 1); + this->b_.resize(size - 1); + this->x_.resize(size); + this->A_.setZero(); + } + + for (auto i = 0; i < size - 2; i++) { + const auto h_i0 = xa[i + 1] - xa[i]; + const auto h_i1 = xa[i + 2] - xa[i + 1]; + const auto y_i0 = ya[i + 1] - ya[i]; + const auto y_i1 = ya[i + 2] - ya[i + 1]; + const auto g_i0 = (h_i0 != 0) ? 1 / h_i0 : 0; + const auto g_i1 = (h_i1 != 0) ? 1 / h_i1 : 0; + this->A_(i + 1, i) = h_i1; + this->A_(i, i) = 2 * (h_i0 + h_i1); + this->A_(i, i + 1) = h_i1; + this->b_(i) = 3 * (y_i1 * g_i1 - y_i0 * g_i0); + } + auto i = size - 2; + const auto h_i0 = xa[i + 1] - xa[i]; + const auto h_i1 = xa[1] - xa[0]; + const auto y_i0 = ya[i + 1] - ya[i]; + const auto y_i1 = ya[1] - ya[0]; + const auto g_i0 = (h_i0 != 0) ? 1 / h_i0 : 0; + const auto g_i1 = (h_i1 != 0) ? 1 / h_i1 : 0; + + this->A_(i, 0) = h_i1; + this->A_(i, i) = 2 * (h_i0 + h_i1); + this->A_(0, i) = h_i1; + this->b_(i) = 3 * (y_i1 * g_i1 - y_i0 * g_i0); + this->x_.segment(1, size - 1) = this->A_.fullPivLu().solve(this->b_); + this->x_(0) = this->x_(size - 1); +} + +} // namespace pyinterp::detail::interpolation diff --git a/src/pyinterp/core/include/pyinterp/detail/interpolation/interpolator.hpp b/src/pyinterp/core/include/pyinterp/detail/interpolation/interpolator.hpp new file mode 100644 index 00000000..ef2c4e3a --- /dev/null +++ b/src/pyinterp/core/include/pyinterp/detail/interpolation/interpolator.hpp @@ -0,0 +1,79 @@ +// Copyright (c) 2024 CNES +// +// All rights reserved. Use of this source code is governed by a +// BSD-style license that can be found in the LICENSE file. +#pragma once +#include +#include + +#include "pyinterp/eigen.hpp" + +namespace pyinterp::detail::interpolation { + +/// @brief Base class for all interpolators +template +class Interpolator { + public: + /// Constructor. + Interpolator() = default; + + /// Destructor. + virtual ~Interpolator() = default; + + /// @brief Search for the index of the first element in xa that is greater + /// than x. + /// @param xa The array to search. + /// @param x The value to search for. + /// @return An optional pair of indices (i, j) such that xa[i] < x <= xa[j]. + constexpr auto search(const Vector &xa, const T &x) const + -> std::optional>; + + /// @brief Search for the index of the first element in xa that is greater + /// than x. + /// @param xa The array to search. + /// @param x The value to search for. + /// @param i The index to start the search from. + /// @return An optional pair of indices (i, j) such that xa[i] < x <= xa[j]. + constexpr auto search(const Vector &xa, const T &x, Eigen::Index *i) const + -> std::optional>; +}; + +template +constexpr auto Interpolator::search(const Vector &xa, const T &x) const + -> std::optional> { + auto index = Eigen::Index{}; + return search(xa, x, &index); +} + +template +constexpr auto Interpolator::search(const Vector &xa, const T &x, + Eigen::Index *index) const + -> std::optional> { + // If the index is within the bounds of the array, check if x is within the + // last interval. + if (*index >= 0 && *index < xa.size() - 1 && xa(*index) < x && + x <= xa(*index + 1)) { + return std::make_pair(*index, *index + 1); + } + const auto begin = xa.array().begin(); + const auto end = xa.array().end(); + for (auto i : {std::max(Eigen::Index(0), (*index) - 1), Eigen::Index(0)}) { + const auto it = std::lower_bound(begin + i, end, x); + if (it == end) { + return std::nullopt; + } + if (it == begin) { + return *it == x + ? std::optional>({0, 1}) + : std::nullopt; + } + i = static_cast(std::distance(begin, it)); + if (xa(i - 1) < x && x <= xa(i)) { + *index = i - 1; + return std::make_pair(*index, i); + } + } + return {}; +} + +} // namespace pyinterp::detail::interpolation diff --git a/src/pyinterp/core/include/pyinterp/detail/interpolation/interpolator_1d.hpp b/src/pyinterp/core/include/pyinterp/detail/interpolation/interpolator_1d.hpp new file mode 100644 index 00000000..6dce39c5 --- /dev/null +++ b/src/pyinterp/core/include/pyinterp/detail/interpolation/interpolator_1d.hpp @@ -0,0 +1,101 @@ +// Copyright (c) 2024 CNES +// +// All rights reserved. Use of this source code is governed by a +// BSD-style license that can be found in the LICENSE file. +#pragma once +#include "pyinterp/detail/interpolation/interpolator.hpp" + +namespace pyinterp::detail::interpolation { + +/// @brief One-dimensional interpolation (interpolation along the x-axis). +/// @tparam T type of data +template +class Interpolator1D : public Interpolator { + public: + /// The minimum size of the arrays to be interpolated. + virtual auto min_size() const -> Eigen::Index = 0; + + /// Interpolate the value of y at x. + /// @param xa X-coordinates of the data points. + /// @param ya Y-coordinates of the data points. + /// @param x The point where the interpolation must be calculated. + /// @return The interpolated value at the point x. + auto operator()(const Eigen::Ref> &xa, + const Eigen::Ref> &ya, const T &x) -> T { + compute_coefficients(xa, ya); + auto index = Eigen::Index{}; + return (*this)(xa, ya, x, &index); + } + + /// Interpolate the values of y at x. + /// @param xa X-coordinates of the data points. + /// @param ya Y-coordinates of the data points. + /// @param x The points where the interpolation must be calculated. + /// @return The interpolated values at the points x. + auto operator()(const Eigen::Ref> &xa, + const Eigen::Ref> &ya, + const Eigen::Ref> &x) -> Vector { + compute_coefficients(xa, ya); + auto ix = Eigen::Index{}; + auto y = Vector(x.size()); + for (Eigen::Index i = 0; i < x.size(); ++i) { + y(i) = (*this)(xa, ya, x(i), &ix); + } + return y; + } + + /// Calculate the derivative of y at x. + /// @param xa X-coordinates of the data points. + /// @param ya Y-coordinates of the data points. + /// @param x The point where the derivative must be calculated. + /// @return The derivative of the interpolation function at the point x. + auto derivative(const Eigen::Ref> &xa, + const Eigen::Ref> &ya, const T &x) -> T { + compute_coefficients(xa, ya); + auto index = Eigen::Index{}; + return derivative(xa, ya, x, &index); + } + + /// Calculate the derivatives of y at x. + /// @param xa X-coordinates of the data points. + /// @param ya Y-coordinates of the data points. + /// @param x The points where the derivative must be calculated. + /// @return The derivatives of the interpolation function at the points x. + auto derivative(const Eigen::Ref> &xa, + const Eigen::Ref> &ya, + const Eigen::Ref> &x) -> Vector { + compute_coefficients(xa, ya); + auto ix = Eigen::Index{}; + auto y = Vector(x.size()); + for (Eigen::Index i = 0; i < x.size(); ++i) { + y(i) = derivative(xa, ya, x(i), &ix); + } + return y; + } + + protected: + /// Interpolate the value of y at x using the index of the last search. + virtual auto operator()(const Eigen::Ref> &xa, + const Eigen::Ref> &ya, const T &x, + Eigen::Index *index) const -> T = 0; + + /// Calculate the derivative of y at x using the index of the last search. + virtual auto derivative(const Eigen::Ref> &xa, + const Eigen::Ref> &ya, const T &x, + Eigen::Index *index) const -> T = 0; + + /// Check if the arrays are valid. + virtual auto compute_coefficients(const Eigen::Ref> &xa, + const Eigen::Ref> &ya) + -> void { + if (xa.size() != ya.size()) { + throw std::invalid_argument("xa and ya must have the same size"); + } + if (xa.size() < min_size()) { + throw std::invalid_argument("xa and ya must have at least " + + std::to_string(min_size()) + " elements"); + } + } +}; + +} // namespace pyinterp::detail::interpolation diff --git a/src/pyinterp/core/include/pyinterp/detail/interpolation/interpolator_2d.hpp b/src/pyinterp/core/include/pyinterp/detail/interpolation/interpolator_2d.hpp new file mode 100644 index 00000000..83a7b66b --- /dev/null +++ b/src/pyinterp/core/include/pyinterp/detail/interpolation/interpolator_2d.hpp @@ -0,0 +1,90 @@ +// Copyright (c) 2024 CNES +// +// All rights reserved. Use of this source code is governed by a +// BSD-style license that can be found in the LICENSE file. +#pragma once + +#include "pyinterp/detail/interpolation/interpolator.hpp" + +namespace pyinterp::detail::interpolation { + +/// Base class for all 2D interpolators +/// @tparam T type of the data +template +class Interpolator2D : public Interpolator { + public: + /// The minimum size of the arrays to be interpolated. + virtual auto min_size() const -> Eigen::Index = 0; + + /// Interpolate the value of y at x. + /// @param xa X-coordinates of the data points. + /// @param ya Y-coordinates of the data points. + /// @param za Z-values of the data points. + /// @param x The point where the interpolation must be calculated. + /// @param y The point where the interpolation must be calculated. + /// @return The interpolated value at the point x. + auto operator()(const Eigen::Ref> &xa, + const Eigen::Ref> &ya, + const Eigen::Ref> &za, const T &x, const T &y) + -> T { + compute_coefficients(xa, ya, za); + auto ix = Eigen::Index{}; + auto jx = Eigen::Index{}; + return (*this)(xa, ya, za, x, y, &ix, &jx); + } + + /// Interpolate the values of y at x. + /// @param xa X-coordinates of the data points. + /// @param ya Y-coordinates of the data points. + /// @param za Z-values of the data points. + /// @param x The point where the interpolation must be calculated. + /// @param y The point where the interpolation must be calculated. + /// @return The interpolated value at the point x. + auto operator()(const Eigen::Ref> &xa, + const Eigen::Ref> &ya, + const Eigen::Ref> &za, + const Eigen::Ref> &x, + const Eigen::Ref> &y) -> Vector { + compute_coefficients(xa, ya, za); + auto ix = Eigen::Index{}; + auto jx = Eigen::Index{}; + auto z = Vector(x.size()); + for (Eigen::Index i = 0; i < x.size(); ++i) { + z(i) = (*this)(xa, ya, za, x(i), y(i), &ix, &jx); + } + return z; + } + + protected: + /// Interpolate the value of y at x using the index of the last search. + virtual auto operator()(const Eigen::Ref> &xa, + const Eigen::Ref> &ya, + const Eigen::Ref> &za, const T &x, + const T &y, Eigen::Index *ix, Eigen::Index *jx) const + -> T = 0; + + /// Check if the arrays are valid. + virtual auto compute_coefficients(const Eigen::Ref> &xa, + const Eigen::Ref> &ya, + const Eigen::Ref> &za) + -> void { + if (xa.size() != za.rows()) { + throw std::invalid_argument( + "xa and za must have the same number of rows"); + } + if (ya.size() != za.cols()) { + throw std::invalid_argument( + "ya and za must have the same number of columns"); + } + if (xa.size() < min_size()) { + throw std::invalid_argument("xa must have at least " + + std::to_string(min_size()) + " elements"); + } + if (ya.size() < min_size()) { + throw std::invalid_argument("ya must have at least " + + std::to_string(min_size()) + " elements"); + } + } +}; + +} // namespace pyinterp::detail::interpolation diff --git a/src/pyinterp/core/include/pyinterp/detail/interpolation/linear.hpp b/src/pyinterp/core/include/pyinterp/detail/interpolation/linear.hpp new file mode 100644 index 00000000..5fe4fdad --- /dev/null +++ b/src/pyinterp/core/include/pyinterp/detail/interpolation/linear.hpp @@ -0,0 +1,58 @@ +// Copyright (c) 2024 CNES +// +// All rights reserved. Use of this source code is governed by a +// BSD-style license that can be found in the LICENSE file. +#pragma once +#include "pyinterp/detail/interpolation/interpolator_1d.hpp" + +namespace pyinterp::detail::interpolation { + +/// Linear interpolation +template +class Linear : public Interpolator1D { + public: + using Interpolator1D::Interpolator1D; + using Interpolator1D::operator(); + using Interpolator1D::derivative; + + /// Returns the minimum number of points required for the interpolation. + auto min_size() const -> Eigen::Index override { return 2; } + + private: + /// Interpolation + /// @param xa X-coordinates of the data points. + /// @param ya Y-coordinates of the data points. + /// @param x The point where the interpolation must be calculated. + /// @param i The index of the last point found in the search. + /// @return The interpolated value at the point x. + auto operator()(const Eigen::Ref> &xa, + const Eigen::Ref> &ya, const T &x, + Eigen::Index *i) const -> T override { + auto where = this->search(xa, x, i); + if (!where) { + return std::numeric_limits::quiet_NaN(); + } + auto [i0, i1] = *where; + return ya(i0) + (ya(i1) - ya(i0)) / (xa(i1) - xa(i0)) * (x - xa(i0)); + } + + /// @brief Returns the derivative of the interpolation function at the point + /// x. + /// @param xa X-coordinates of the data points. + /// @param ya Y-coordinates of the data points. + /// @param x The point where the derivative must be calculated. + /// @param i The index of the last point found in the search. + /// @return The derivative of the interpolation function at the point x. + auto derivative(const Eigen::Ref> &xa, + const Eigen::Ref> &ya, const T &x, + Eigen::Index *i) const -> T override { + auto where = this->search(xa, x, i); + if (!where) { + return std::numeric_limits::quiet_NaN(); + } + auto [i0, i1] = *where; + return (ya(i1) - ya(i0)) / (xa(i1) - xa(i0)); + } +}; + +} // namespace pyinterp::detail::interpolation diff --git a/src/pyinterp/core/include/pyinterp/detail/interpolation/polynomial.hpp b/src/pyinterp/core/include/pyinterp/detail/interpolation/polynomial.hpp new file mode 100644 index 00000000..54ccdbb2 --- /dev/null +++ b/src/pyinterp/core/include/pyinterp/detail/interpolation/polynomial.hpp @@ -0,0 +1,129 @@ +// Copyright (c) 2024 CNES +// +// All rights reserved. Use of this source code is governed by a +// BSD-style license that can be found in the LICENSE file. +#pragma once +#include "pyinterp/detail/interpolation/interpolator_1d.hpp" + +namespace pyinterp::detail::interpolation { + +/// Polynomial interpolation +template +class Polynomial : public Interpolator1D { + public: + using Interpolator1D::Interpolator1D; + using Interpolator1D::operator(); + using Interpolator1D::derivative; + + /// The minimum size of the arrays to be interpolated. + auto min_size() const -> Eigen::Index override { return 3; } + + private: + /// Coefficients of the interpolation + Vector work_{}; + + /// Compute the coefficients of the interpolation + /// @param xa X-coordinates of the data points. + /// @param ya Y-coordinates of the data points. + auto compute_coefficients(const Eigen::Ref> &xa, + const Eigen::Ref> &ya) + -> void override; + + /// Compute the coefficients of the interpolation + /// @param xa X-coordinates of the data points. + /// @param x The point where the interpolation must be calculated. + auto taylor(const Eigen::Ref> &xa, const T &x) const + -> Vector; + + /// Interpolation + /// @param xa X-coordinates of the data points. + /// @param ya Y-coordinates of the data points. + /// @param x The point where the interpolation must be calculated. + /// @param index The index of the last point found in the search. + /// @return The interpolated value at the point x. + auto operator()(const Eigen::Ref> &xa, + const Eigen::Ref> &ya, const T &x, + Eigen::Index *index) const -> T override; + + /// @brief Returns the derivative of the interpolation function at the point + /// x. + /// @param xa X-coordinates of the data points. + /// @param ya Y-coordinates of the data points. + /// @param x The point where the derivative must be calculated. + /// @param index The index of the last point found in the search. + /// @return The derivative of the interpolation function at the point x. + auto derivative(const Eigen::Ref> &xa, + const Eigen::Ref> &ya, const T &x, + Eigen::Index *index) const -> T override; +}; + +template +auto Polynomial::compute_coefficients(const Eigen::Ref> &xa, + const Eigen::Ref> &ya) + -> void { + Interpolator1D::compute_coefficients(xa, ya); + auto size = xa.size(); + if (work_.size() < size) { + work_.resize(size); + } + work_[0] = ya[0]; + work_.segment(1, size - 1) = + (ya.segment(1, size - 1) - ya.segment(0, size - 1)).array() / + (xa.segment(1, size - 1) - xa.segment(0, size - 1)).array(); + for (Eigen::Index ix = 2; ix < size; ix++) { + work_.segment(ix, size - ix) = + (work_.segment(ix, size - ix) - work_.segment(ix - 1, size - ix)) + .array() / + (xa.segment(ix, size - ix) - xa.segment(ix - 1, size - ix)).array(); + } +} + +template +auto Polynomial::taylor(const Eigen::Ref> &xa, + const T &x) const -> Vector { + auto size = xa.size(); + auto c = Vector(size); + auto w = Vector(size); + w(size - 1) = T(1); + c(0) = work_(0); + for (Eigen::Index ix = size - 1; ix--;) { + w(ix) = -w(ix + 1) * (xa(size - 2 - ix) - x); + for (Eigen::Index jx = ix + 1; jx < size - 1; ++jx) { + w(jx) -= w(jx + 1) * (xa(size - 2 - ix) - x); + } + for (Eigen::Index jx = ix; jx < size; ++jx) { + c(jx - ix) += w(jx) * work_(size - ix - 1); + } + } + return c; +} + +template +auto Polynomial::operator()(const Eigen::Ref> &xa, + const Eigen::Ref> &ya, + const T &x, Eigen::Index *index) const -> T { + auto search = this->search(xa, x, index); + if (!search) { + throw std::numeric_limits::quiet_NaN(); + } + auto size = xa.size(); + auto y = work_[size - 1]; + for (Eigen::Index ix = size - 1; ix--;) { + y = work_[ix] + (x - xa(ix)) * y; + } + return y; +} + +template +auto Polynomial::derivative(const Eigen::Ref> &xa, + const Eigen::Ref> &ya, + const T &x, Eigen::Index *index) const -> T { + auto search = this->search(xa, x, index); + if (!search) { + throw std::numeric_limits::quiet_NaN(); + } + auto coefficients = taylor(xa, x); + return coefficients(1); +} + +} // namespace pyinterp::detail::interpolation diff --git a/src/pyinterp/core/include/pyinterp/detail/interpolation/steffen.hpp b/src/pyinterp/core/include/pyinterp/detail/interpolation/steffen.hpp new file mode 100644 index 00000000..6a44ca4b --- /dev/null +++ b/src/pyinterp/core/include/pyinterp/detail/interpolation/steffen.hpp @@ -0,0 +1,138 @@ +// Copyright (c) 2024 CNES +// +// All rights reserved. Use of this source code is governed by a +// BSD-style license that can be found in the LICENSE file. +#pragma once +#include "pyinterp/detail/interpolation/interpolator_1d.hpp" + +namespace pyinterp::detail::interpolation { + +/// Steffen interpolation +template +class Steffen : public Interpolator1D { + public: + using Interpolator1D::Interpolator1D; + using Interpolator1D::operator(); + using Interpolator1D::derivative; + + /// Returns the minimum number of points required for the interpolation. + auto min_size() const -> Eigen::Index override { return 3; } + + private: + /// Interpolation + /// @param xa X-coordinates of the data points. + /// @param ya Y-coordinates of the data points. + /// @param x The point where the interpolation must be calculated. + /// @param i The index of the last point found in the search. + /// @return The interpolated value at the point x. + auto operator()(const Eigen::Ref> &xa, + const Eigen::Ref> &ya, const T &x, + Eigen::Index *i) const -> T override; + + /// @brief Returns the derivative of the interpolation function at the point + /// x. + /// @param xa X-coordinates of the data points. + /// @param ya Y-coordinates of the data points. + /// @param x The point where the derivative must be calculated. + /// @param i The index of the last point found in the search. + auto derivative(const Eigen::Ref> &xa, + const Eigen::Ref> &ya, const T &x, + Eigen::Index *i) const -> T override; + + private: + /// Compute the coefficients of the interpolation + /// @param xa X-coordinates of the data points. + /// @param ya Y-coordinates of the data points. + auto compute_coefficients(const Eigen::Ref> &xa, + const Eigen::Ref> &ya) + -> void override; + + /// Return the sign of x multiplied by the sign of y + static constexpr auto copysign(const T &x, const T &y) -> T { + if ((x < T(0) && y > T(0)) || (x > T(0) && y < T(0))) { + return -x; + } + return x; + } + /// The slopes of the interpolation + Vector y_prime_; +}; + +template +auto Steffen::compute_coefficients(const Eigen::Ref> &xa, + const Eigen::Ref> &ya) + -> void { + Interpolator1D::compute_coefficients(xa, ya); + auto size = xa.size(); + if (y_prime_.size() < size) { + y_prime_.resize(size); + } + // First assign the interval and slopes for the left boundary. + // We use the "simplest possibility" method described in the paper + // in section 2.2 + auto h0 = (xa[1] - xa[0]); + auto s0 = (ya[1] - ya[0]) / h0; + + y_prime_[0] = s0; + + // Now we calculate all the necessary s, h, p, and y' variables from 1 to + // size-2 (0 to size - 2 inclusive) + for (Eigen::Index i = 1; i < size - 1; ++i) { + // Eq. 6 + auto hi = (xa[i + 1] - xa[i]); + auto him1 = (xa[i] - xa[i - 1]); + // Eq. 7 + auto si = (ya[i + 1] - ya[i]) / hi; + auto sim1 = (ya[i] - ya[i - 1]) / him1; + // Eq. 8 + auto pi = (sim1 * hi + si * him1) / (him1 + hi); + + y_prime_[i] = + (Steffen::copysign(T(1), sim1) + Steffen::copysign(T(1), si)) * + std::min(std::fabs(sim1), std::min(std::fabs(si), 0.5 * std::fabs(pi))); + } + // We also need y' for the rightmost boundary; we use the + // "simplest possibility" method described in the paper in + // section 2.2 + y_prime_[size - 1] = + (ya[size - 1] - ya[size - 2]) / (xa[size - 1] - xa[size - 2]); +} + +template +auto Steffen::operator()(const Eigen::Ref> &xa, + const Eigen::Ref> &ya, const T &x, + Eigen::Index *i) const -> T { + auto where = this->search(xa, x, i); + if (!where) { + return std::numeric_limits::quiet_NaN(); + } + auto [i0, i1] = *where; + const auto h = x - xa[i0]; + const auto hi = (xa[i1] - xa[i0]); + const auto si = (ya[i1] - ya[i0]) / hi; + const auto a = (y_prime_[i0] + y_prime_[i1] - 2 * si) / hi / hi; + const auto b = (3 * si - 2 * y_prime_[i0] - y_prime_[i1]) / hi; + const auto c = y_prime_[i0]; + const auto d = ya[i0]; + return d + h * (c + h * (b + h * a)); +} + +template +auto Steffen::derivative(const Eigen::Ref> &xa, + const Eigen::Ref> &ya, const T &x, + Eigen::Index *i) const -> T { + auto where = this->search(xa, x, i); + if (!where) { + return std::numeric_limits::quiet_NaN(); + } + auto [i0, i1] = *where; + const auto h = x - xa[i0]; + const auto hi = (xa[i1] - xa[i0]); + const auto si = (ya[i1] - ya[i0]) / hi; + const auto a = (y_prime_[i0] + y_prime_[i1] - 2 * si) / hi / hi; + const auto b = (3 * si - 2 * y_prime_[i0] - y_prime_[i1]) / hi; + const auto c = y_prime_[i0]; + return c + h * (2 * b + h * 3 * a); +} + +} // namespace pyinterp::detail::interpolation diff --git a/src/pyinterp/core/tests/CMakeLists.txt b/src/pyinterp/core/tests/CMakeLists.txt index 7baa57a9..29cc8dc3 100755 --- a/src/pyinterp/core/tests/CMakeLists.txt +++ b/src/pyinterp/core/tests/CMakeLists.txt @@ -20,12 +20,20 @@ add_testcase(geodetic_coordinates) add_testcase(geodetic_system) add_testcase(geometry_rtree) add_testcase(gsl GSL::gsl GSL::gslcblas) +add_testcase(interpolation_akima) +add_testcase(interpolation_bicubic ${BLAS_LIBRARIES}) +add_testcase(interpolation_bilinear) +add_testcase(interpolation_cspline ${BLAS_LIBRARIES}) +add_testcase(interpolation_linear) +add_testcase(interpolation_polynomial) +add_testcase(interpolation_search) +add_testcase(interpolation_steffen) add_testcase(math_bicubic GSL::gsl GSL::gslcblas) add_testcase(math_binning) add_testcase(math_bivariate) add_testcase(math_descriptive_statistics) -add_testcase(math_linear) add_testcase(math_kriging) +add_testcase(math_linear) add_testcase(math_rbf) add_testcase(math_spline GSL::gsl GSL::gslcblas) add_testcase(math_streaming_histogram) diff --git a/src/pyinterp/core/tests/interpolation_akima.cpp b/src/pyinterp/core/tests/interpolation_akima.cpp new file mode 100644 index 00000000..95771a42 --- /dev/null +++ b/src/pyinterp/core/tests/interpolation_akima.cpp @@ -0,0 +1,65 @@ +// Copyright (c) 2024 CNES +// +// All rights reserved. Use of this source code is governed by a +// BSD-style license that can be found in the LICENSE file. +#include + +#include "pyinterp/detail/interpolation/akima.hpp" +#include "pyinterp/detail/interpolation/akima_periodic.hpp" + +TEST(Akima, interpolate) { + Eigen::Matrix xa(5); + Eigen::Matrix ya(5); + Eigen::Matrix xp(4); + Eigen::Matrix yp(4); + + xa << 0.0, 1.0, 2.0, 3.0, 4.0; + ya << 0.0, 1.0, 2.0, 3.0, 4.0; + + xp << 0.0, 0.5, 1.0, 2.0; + yp << 0.0, 0.5, 1.0, 2.0; + + auto interpolator = pyinterp::detail::interpolation::Akima(); + auto y = interpolator(xa, ya, xp); + for (auto ix = 0; ix < xp.size(); ix++) { + EXPECT_DOUBLE_EQ(y(ix), yp(ix)); + } +} + +TEST(Akima, derivative) { + Eigen::Matrix xa(5); + Eigen::Matrix ya(5); + Eigen::Matrix xp(4); + Eigen::Matrix dyp(4); + + xa << 0.0, 1.0, 2.0, 3.0, 4.0; + ya << 0.0, 1.0, 2.0, 3.0, 4.0; + + xp << 0.0, 0.5, 1.0, 2.0; + dyp << 1.0, 1.0, 1.0, 1.0; + + auto interpolator = pyinterp::detail::interpolation::Akima(); + auto dy = interpolator.derivative(xa, ya, xp); + for (auto ix = 0; ix < xp.size(); ix++) { + EXPECT_DOUBLE_EQ(dy(ix), dyp(ix)); + } +} + +TEST(AkimaPeriodic, interpolate) { + Eigen::Matrix xa(5); + Eigen::Matrix ya(5); + Eigen::Matrix xp(4); + Eigen::Matrix yp(4); + + xa << 0.0, 1.0, 2.0, 3.0, 4.0; + ya << 0.0, 1.0, 2.0, 3.0, 4.0; + + xp << 0.0, 0.5, 1.0, 2.0; + yp << 0.0, 0.5, 1.0, 2.0; + + auto interpolator = pyinterp::detail::interpolation::AkimaPeriodic(); + auto y = interpolator(xa, ya, xp); + for (auto ix = 0; ix < xp.size(); ix++) { + EXPECT_DOUBLE_EQ(y(ix), yp(ix)); + } +} diff --git a/src/pyinterp/core/tests/interpolation_bicubic.cpp b/src/pyinterp/core/tests/interpolation_bicubic.cpp new file mode 100644 index 00000000..f72666a8 --- /dev/null +++ b/src/pyinterp/core/tests/interpolation_bicubic.cpp @@ -0,0 +1,80 @@ +// Copyright (c) 2024 CNES +// +// All rights reserved. Use of this source code is governed by a +// BSD-style license that can be found in the LICENSE file. +#include + +#include "pyinterp/detail/interpolation/bicubic.hpp" + +TEST(Bicubic, case_one) { + Eigen::Matrix za(4, 4); + Eigen::Matrix xa(4); + Eigen::Matrix ya(4); + Eigen::Matrix xp(3); + Eigen::Matrix yp(3); + Eigen::Matrix zp(3); + + xa << 0.0, 1.0, 2.0, 3.0; + ya << 0.0, 1.0, 2.0, 3.0; + za << 1.0, 1.1, 1.2, 1.3, 1.1, 1.2, 1.3, 1.4, 1.2, 1.3, 1.4, 1.5, 1.3, 1.4, + 1.5, 1.6; + xp << 1.0, 1.5, 2.0; + yp << 1.0, 1.5, 2.0; + zp << 1.2, 1.3, 1.4; + auto bicubic = pyinterp::detail::interpolation::Bicubic(); + auto z = bicubic(xa, ya, za, xp, yp); + for (Eigen::Index i = 0; i < z.size(); ++i) { + EXPECT_NEAR(z(i), zp(i), 1.0e-12); + } +} + +TEST(Bicubic, non_linear) { + Eigen::Matrix za(8, 8); + Eigen::Matrix xa(8); + Eigen::Matrix ya(8); + Eigen::Matrix xp(7); + Eigen::Matrix yp(7); + Eigen::Matrix zp(7); + + xa << 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0; + ya << 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0; + za << 1, 2, 3, 4, 5, 6, 7, 8, 2, 2, 6, 4, 10, 6, 14, 8, 3, 6, 3, 12, 15, 6, + 21, 24, 4, 4, 12, 4, 20, 12, 28, 8, 5, 10, 15, 20, 5, 30, 35, 40, 6, 6, 6, + 12, 30, 6, 42, 24, 7, 14, 21, 28, 35, 42, 7, 56, 8, 8, 24, 8, 40, 24, 56, + 8; + xp << 1.4, 2.3, 4.7, 3.3, 7.5, 6.6, 5.1; + yp << 1.0, 1.8, 1.9, 2.5, 2.7, 4.1, 3.3; + zp << 1.4, 3.11183531264736, 8.27114315792559, 5.03218982537718, + 22.13230634702637, 23.63206834997871, 17.28553080971182; + auto bicubic = pyinterp::detail::interpolation::Bicubic(); + auto z = bicubic(xa, ya, za, xp, yp); + for (Eigen::Index i = 0; i < z.size(); ++i) { + EXPECT_NEAR(z(i), zp(i), 1.0e-12); + } +} + +TEST(Bicubic, non_square) { + Eigen::Matrix za(10, 8); + Eigen::Matrix xa(10); + Eigen::Matrix ya(8); + Eigen::Matrix xp(7); + Eigen::Matrix yp(7); + Eigen::Matrix zp(7); + + xa << 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0; + ya << 1.0, 4.0, 6.0, 8.0, 10.0, 12.0, 14.0, 16.0; + za << 1, 2, 3, 4, 5, 6, 7, 8, 2, 2, 6, 4, 10, 6, 14, 8, 3, 6, 3, 12, 15, 6, + 21, 24, 4, 4, 12, 4, 20, 12, 28, 8, 5, 10, 15, 20, 5, 30, 35, 40, 6, 6, 6, + 12, 30, 6, 42, 24, 7, 14, 21, 28, 35, 42, 7, 56, 8, 8, 24, 8, 40, 24, 56, + 8, 9, 11, 13, 15, 17, 19, 21, 23, 10, 12, 14, 16, 18, 20, 22, 24; + xp << 1.4, 2.3, 9.7, 3.3, 9.5, 6.6, 5.1; + yp << 1.0, 1.8, 1.9, 2.5, 2.7, 4.1, 3.3; + zp << 1.4, 2.46782030941187003, 10.7717721621846465, 4.80725067958096375, + 11.6747032398627297, 11.2619968682970111, 9.00168877916872567; + + auto bicubic = pyinterp::detail::interpolation::Bicubic(); + auto z = bicubic(xa, ya, za, xp, yp); + for (Eigen::Index i = 0; i < z.size(); ++i) { + EXPECT_NEAR(z(i), zp(i), 1.0e-12); + } +} diff --git a/src/pyinterp/core/tests/interpolation_bilinear.cpp b/src/pyinterp/core/tests/interpolation_bilinear.cpp new file mode 100644 index 00000000..6eeb8088 --- /dev/null +++ b/src/pyinterp/core/tests/interpolation_bilinear.cpp @@ -0,0 +1,56 @@ +// Copyright (c) 2024 CNES +// +// All rights reserved. Use of this source code is governed by a +// BSD-style license that can be found in the LICENSE file. +#include + +#include "pyinterp/detail/interpolation/bilinear.hpp" + +TEST(Bilinear, symmetric) { + Eigen::Matrix za(4, 4); + Eigen::Matrix xa(4); + Eigen::Matrix ya(4); + Eigen::Matrix xp(6); + Eigen::Matrix yp(6); + Eigen::Matrix zp(6); + + xa << 0.0, 1.0, 2.0, 3.0; + ya << 0.0, 1.0, 2.0, 3.0; + za << 1.0, 1.1, 1.2, 1.3, 1.1, 1.2, 1.3, 1.4, 1.2, 1.3, 1.4, 1.5, 1.3, 1.4, + 1.5, 1.6; + xp << 0.0, 0.5, 1.0, 1.5, 2.5, 3.0; + yp << 0.0, 0.5, 1.0, 1.5, 2.5, 3.0; + zp << 1.0, 1.1, 1.2, 1.3, 1.5, 1.6; + + auto bilinear = pyinterp::detail::interpolation::Bilinear(); + auto z = bilinear(xa, ya, za, xp, yp); + for (Eigen::Index i = 0; i < z.size(); ++i) { + EXPECT_NEAR(z(i), zp(i), 1.0e-12); + } +} + +TEST(Bilinear, asymmetric) { + Eigen::Matrix za(4, 4); + Eigen::Matrix xa(4); + Eigen::Matrix ya(4); + Eigen::Matrix xp(12); + Eigen::Matrix yp(12); + Eigen::Matrix zp(12); + + xa << 0.0, 1.0, 2.0, 3.0; + ya << 0.0, 1.0, 2.0, 3.0; + za << 1.0, 1.3, 1.5, 1.6, 1.1, 1.4, 1.6, 1.9, 1.2, 1.5, 1.7, 2.2, 1.4, 1.7, + 1.9, 2.3; + xp << 0.0, 0.5, 1.0, 1.5, 2.5, 3.0, 1.3954, 1.6476, 0.824957, 2.41108, + 2.98619, 1.36485; + yp << 0.0, 0.5, 1.0, 1.5, 2.5, 3.0, 0.265371, 2.13849, 1.62114, 1.22198, + 0.724681, 0.0596087; + zp << 1.0, 1.2, 1.4, 1.55, 2.025, 2.3, 1.2191513, 1.7242442248, 1.5067237, + 1.626612, 1.6146423, 1.15436761; + + auto bilinear = pyinterp::detail::interpolation::Bilinear(); + auto z = bilinear(xa, ya, za, xp, yp); + for (Eigen::Index i = 0; i < z.size(); ++i) { + EXPECT_NEAR(z(i), zp(i), 1.0e-12); + } +} diff --git a/src/pyinterp/core/tests/interpolation_cspline.cpp b/src/pyinterp/core/tests/interpolation_cspline.cpp new file mode 100644 index 00000000..2a54ab75 --- /dev/null +++ b/src/pyinterp/core/tests/interpolation_cspline.cpp @@ -0,0 +1,166 @@ +// Copyright (c) 2024 CNES +// +// All rights reserved. Use of this source code is governed by a +// BSD-style license that can be found in the LICENSE file. +#include + +#include "pyinterp/detail/interpolation/cspline.hpp" +#include "pyinterp/detail/interpolation/cspline_periodic.hpp" + +TEST(CSpline, case_one) { + Eigen::Matrix xa(6); + Eigen::Matrix ya(6); + Eigen::Matrix xp(50); + Eigen::Matrix yp(50); + Eigen::Matrix dyp(50); + + xa << 0.0, 0.2, 0.4, 0.6, 0.8, 1.0; + ya << 1.0, 0.961538461538461, 0.862068965517241, 0.735294117647059, + 0.609756097560976, 0.500000000000000; + + xp << 0.00, 0.02, 0.04, 0.06, 0.08, 0.10, 0.12, 0.14, 0.16, 0.18, 0.20, 0.22, + 0.24, 0.26, 0.28, 0.30, 0.32, 0.34, 0.36, 0.38, 0.40, 0.42, 0.44, 0.46, + 0.48, 0.50, 0.52, 0.54, 0.56, 0.58, 0.60, 0.62, 0.64, 0.66, 0.68, 0.70, + 0.72, 0.74, 0.76, 0.78, 0.80, 0.82, 0.84, 0.86, 0.88, 0.90, 0.92, 0.94, + 0.96, 0.98; + yp << 1.000000000000000, 0.997583282975581, 0.995079933416512, + 0.992403318788142, 0.989466806555819, 0.986183764184894, + 0.982467559140716, 0.978231558888635, 0.973389130893999, + 0.967853642622158, 0.961538461538461, 0.954382579685350, + 0.946427487413627, 0.937740299651188, 0.928388131325928, + 0.918438097365742, 0.907957312698524, 0.897012892252170, + 0.885671950954575, 0.874001603733634, 0.862068965517241, + 0.849933363488199, 0.837622973848936, 0.825158185056786, + 0.812559385569085, 0.799846963843167, 0.787041308336369, + 0.774162807506023, 0.761231849809467, 0.748268823704033, + 0.735294117647059, 0.722328486073082, 0.709394147325463, + 0.696513685724764, 0.683709685591549, 0.671004731246381, + 0.658421407009825, 0.645982297202442, 0.633709986144797, + 0.621627058157454, 0.609756097560976, 0.598112015427308, + 0.586679029833925, 0.575433685609685, 0.564352527583445, + 0.553412100584061, 0.542588949440392, 0.531859618981294, + 0.521200654035625, 0.510588599432241; + + dyp << -0.120113913432180, -0.122279726798445, -0.128777166897241, + -0.139606233728568, -0.154766927292426, -0.174259247588814, + -0.198083194617734, -0.226238768379184, -0.258725968873165, + -0.295544796099676, -0.336695250058719, -0.378333644186652, + -0.416616291919835, -0.451543193258270, -0.483114348201955, + -0.511329756750890, -0.536189418905076, -0.557693334664512, + -0.575841504029200, -0.590633926999137, -0.602070603574326, + -0.611319695518765, -0.619549364596455, -0.626759610807396, + -0.632950434151589, -0.638121834629033, -0.642273812239728, + -0.645406366983674, -0.647519498860871, -0.648613207871319, + -0.648687494015019, -0.647687460711257, -0.645558211379322, + -0.642299746019212, -0.637912064630930, -0.632395167214473, + -0.625749053769843, -0.617973724297039, -0.609069178796061, + -0.599035417266910, -0.587872439709585, -0.576731233416743, + -0.566762785681043, -0.557967096502484, -0.550344165881066, + -0.543893993816790, -0.538616580309654, -0.534511925359660, + -0.531580028966807, -0.529820891131095; + + auto interpolator = pyinterp::detail::interpolation::CSpline(); + auto y = interpolator(xa, ya, xp); + for (auto ix = 0; ix < xp.size(); ix++) { + EXPECT_NEAR(y(ix), yp(ix), 1e-6); + } + auto dy = interpolator.derivative(xa, ya, xp); + for (auto ix = 0; ix < xp.size(); ix++) { + EXPECT_NEAR(dy(ix), dyp(ix), 1e-6); + } +} + +TEST(CSpline, case_two) { + Eigen::Matrix xa(7); + Eigen::Matrix ya(7); + Eigen::Matrix xp(19); + Eigen::Matrix yp(19); + Eigen::Matrix dyp(19); + + xa << -1.2139767065644265, -0.792590494453907, -0.250954683125019, + 0.665867809951305, 0.735655088722706, 0.827622053027153, + 1.426592227816582; + ya << -0.00453877449035645, 0.49763182550668716, 0.17805472016334534, + 0.40514493733644485, -0.21595209836959839, 0.47405586764216423, + 0.46561462432146072; + + xp << -1.2139767065644265, -1.0735146358609200, -0.9330525651574135, + -0.7925904944539071, -0.6120452240109444, -0.4314999535679818, + -0.2509546831250191, 0.0546528145670890, 0.3602603122591972, + 0.6658678099513053, 0.6891302362084388, 0.7123926624655723, + 0.7356550887227058, 0.7663107434908548, 0.7969663982590039, + 0.8276220530271530, 1.0272787779569625, 1.2269355028867721, + 1.4265922278165817; + yp << -0.00453877449035645, 0.25816917628390590, 0.44938881397673230, + 0.49763182550668716, 0.31389980410075147, 0.09948951681196887, + 0.17805472016334534, 1.27633142487980233, 2.04936553432792001, + 0.40514493733644485, 0.13322324792901385, -0.09656315924697809, + -0.21595209836959839, -0.13551147728045118, 0.13466779030061801, + 0.47405586764216423, 1.68064089899304370, 1.43594739539458649, + 0.46561462432146072; + dyp << 1.955137555965937, 1.700662049790549, 0.937235531264386, + -0.335141999612553, -1.401385073563169, -0.674982149482761, + 1.844066772628670, 4.202528085784793, -0.284432022227558, + -11.616813551408383, -11.272731243226174, -7.994209291156876, + -1.781247695200491, 6.373970868827501, 10.597456848997197, + 10.889210245308570, 1.803124267866902, -3.648527318598099, + -5.465744514086432; + + auto interpolator = pyinterp::detail::interpolation::CSpline(); + auto y = interpolator(xa, ya, xp); + for (auto ix = 0; ix < xp.size(); ix++) { + EXPECT_NEAR(y(ix), yp(ix), 1e-6); + } + auto dy = interpolator.derivative(xa, ya, xp); + for (auto ix = 0; ix < xp.size(); ix++) { + EXPECT_NEAR(dy(ix), dyp(ix), 1e-6); + } +} + +TEST(CSplinePeriodic, case_one) { + Eigen::Matrix xa(11); + Eigen::Matrix ya(11); + Eigen::Matrix xp(31); + Eigen::Matrix yp(31); + + xa << 0.000000000000000, 0.130153674349869, 0.164545962312740, + 0.227375461261537, 0.256465324353657, 0.372545206874658, + 0.520820016781720, 0.647654717733075, 0.753429306654340, + 0.900873984827658, 1.000000000000000; + ; + ya << 0.000000000000000, 0.729629261832041, 0.859286331568207, + 0.989913099419008, 0.999175006262120, 0.717928599519215, + -0.130443237213363, -0.800267961158980, -0.999767873040527, + -0.583333769240853, -0.000000000000000; + ; + + xp << 0.000000000000000, 0.043384558116623, 0.086769116233246, + 0.130153674349869, 0.141617770337492, 0.153081866325116, + 0.164545962312740, 0.185489128629005, 0.206432294945271, + 0.227375461261537, 0.237072082292243, 0.246768703322951, + 0.256465324353657, 0.295158618527324, 0.333851912700991, + 0.372545206874658, 0.421970143510346, 0.471395080146033, + 0.520820016781720, 0.563098250432172, 0.605376484082623, + 0.647654717733075, 0.682912914040164, 0.718171110347252, + 0.753429306654340, 0.802577532712113, 0.851725758769885, + 0.900873984827658, 0.933915989885105, 0.966957994942553, + 1.000000000000000; + yp << 0.000000000000000, 0.268657574670719, 0.517940878523929, + 0.729629261832041, 0.777012551497867, 0.820298314554859, + 0.859286331568207, 0.918833991960315, 0.962624749226346, + 0.989913099419008, 0.996756196601349, 0.999858105635752, + 0.999175006262120, 0.959248551766306, 0.863713527741856, + 0.717928599519215, 0.470065187871106, 0.177694938589523, + -0.130443237213363, -0.385093922365765, -0.613840011545983, + -0.800267961158980, -0.912498361131651, -0.980219217412290, + -0.999767873040528, -0.943635958253643, -0.800314354800596, + -0.583333769240853, -0.403689914131666, -0.206151346799382, + -0.000000000000000; + + auto interpolator = + pyinterp::detail::interpolation::CSplinePeriodic(); + auto y = interpolator(xa, ya, xp); + for (auto ix = 0; ix < xp.size(); ix++) { + EXPECT_NEAR(y(ix), yp(ix), 1e-6); + } +} diff --git a/src/pyinterp/core/tests/interpolation_linear.cpp b/src/pyinterp/core/tests/interpolation_linear.cpp new file mode 100644 index 00000000..4415fcbd --- /dev/null +++ b/src/pyinterp/core/tests/interpolation_linear.cpp @@ -0,0 +1,45 @@ +// Copyright (c) 2024 CNES +// +// All rights reserved. Use of this source code is governed by a +// BSD-style license that can be found in the LICENSE file. +#include + +#include "pyinterp/detail/interpolation/linear.hpp" + +TEST(Linear, interpolate) { + Eigen::Matrix xa(4); + Eigen::Matrix ya(4); + Eigen::Matrix xp(6); + Eigen::Matrix yp(6); + + xa << 0.0, 1.0, 2.0, 3.0; + ya << 0.0, 1.0, 2.0, 3.0; + + xp << 0.0, 0.5, 1.0, 1.5, 2.5, 3.0; + yp << 0.0, 0.5, 1.0, 1.5, 2.5, 3.0; + + auto interpolator = pyinterp::detail::interpolation::Linear(); + auto y = interpolator(xa, ya, xp); + for (Eigen::Index i = 0; i < xp.size(); ++i) { + EXPECT_DOUBLE_EQ(y(i), yp(i)); + } +} + +TEST(Linear, derivative) { + Eigen::Matrix xa(4); + Eigen::Matrix ya(4); + Eigen::Matrix xp(6); + Eigen::Matrix dyp(6); + + xa << 0.0, 1.0, 2.0, 3.0; + ya << 0.0, 1.0, 2.0, 3.0; + + xp << 0.0, 0.5, 1.0, 1.5, 2.5, 3.0; + dyp << 1.0, 1.0, 1.0, 1.0, 1.0, 1.0; + + auto interpolator = pyinterp::detail::interpolation::Linear(); + auto dy = interpolator.derivative(xa, ya, xp); + for (Eigen::Index i = 0; i < xp.size(); ++i) { + EXPECT_DOUBLE_EQ(dy(i), dyp(i)); + } +} diff --git a/src/pyinterp/core/tests/interpolation_polynomial.cpp b/src/pyinterp/core/tests/interpolation_polynomial.cpp new file mode 100644 index 00000000..745e8828 --- /dev/null +++ b/src/pyinterp/core/tests/interpolation_polynomial.cpp @@ -0,0 +1,45 @@ +// Copyright (c) 2024 CNES +// +// All rights reserved. Use of this source code is governed by a +// BSD-style license that can be found in the LICENSE file. +#include + +#include "pyinterp/detail/interpolation/polynomial.hpp" + +TEST(Polynomial, interpolate) { + Eigen::Matrix xa(4); + Eigen::Matrix ya(4); + Eigen::Matrix xp(6); + Eigen::Matrix yp(6); + + xa << 0.0, 1.0, 2.0, 3.0; + ya << 0.0, 1.0, 2.0, 3.0; + + xp << 0.0, 0.5, 1.0, 1.5, 2.5, 3.0; + yp << 0.0, 0.5, 1.0, 1.5, 2.5, 3.0; + + auto interpolator = pyinterp::detail::interpolation::Polynomial(); + auto y = interpolator(xa, ya, xp); + for (auto ix = 0; ix < xp.size(); ix++) { + EXPECT_DOUBLE_EQ(y(ix), yp(ix)); + } +} + +TEST(Polynomial, derivative) { + Eigen::Matrix xa(4); + Eigen::Matrix ya(4); + Eigen::Matrix xp(6); + Eigen::Matrix dyp(6); + + xa << 0.0, 1.0, 2.0, 3.0; + ya << 0.0, 1.0, 2.0, 3.0; + + xp << 0.0, 0.5, 1.0, 1.5, 2.5, 3.0; + dyp << 1.0, 1.0, 1.0, 1.0, 1.0, 1.0; + + auto interpolator = pyinterp::detail::interpolation::Polynomial(); + auto dy = interpolator.derivative(xa, ya, xp); + for (auto ix = 0; ix < xp.size(); ix++) { + EXPECT_DOUBLE_EQ(dy(ix), dyp(ix)); + } +} diff --git a/src/pyinterp/core/tests/interpolation_search.cpp b/src/pyinterp/core/tests/interpolation_search.cpp new file mode 100644 index 00000000..6448b3a9 --- /dev/null +++ b/src/pyinterp/core/tests/interpolation_search.cpp @@ -0,0 +1,61 @@ +// Copyright (c) 2024 CNES +// +// All rights reserved. Use of this source code is governed by a +// BSD-style license that can be found in the LICENSE file. +#include + +#include "pyinterp/detail/interpolation/interpolator.hpp" + +template +using Vector = Eigen::Array; + +TEST(Interpolator, search) { + Vector x(10); + x << 0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9; + + auto interpolator = pyinterp::detail::interpolation::Interpolator(); + auto where = interpolator.search(x, 0.5); + ASSERT_TRUE(where.has_value()); + auto [i0, i1] = *where; + EXPECT_EQ(i0, 4); + EXPECT_EQ(i1, 5); + auto index = static_cast(5); + where = interpolator.search(x, 0.5, &index); + ASSERT_TRUE(where.has_value()); + auto [i02, i12] = *where; + EXPECT_EQ(i02, 4); + EXPECT_EQ(i12, 5); + index = 6; + where = interpolator.search(x, 0.5, &index); + ASSERT_TRUE(where.has_value()); + index = 7; + where = interpolator.search(x, 0.5, &index); + ASSERT_TRUE(where.has_value()); + where = interpolator.search(x, 0.0); + ASSERT_TRUE(where.has_value()); + EXPECT_EQ(where->first, 0); + EXPECT_EQ(where->second, 1); + where = interpolator.search(x, 0.9); + ASSERT_TRUE(where.has_value()); + EXPECT_EQ(where->first, 8); + EXPECT_EQ(where->second, 9); + index = 8; + where = interpolator.search(x, 0.9, &index); + ASSERT_TRUE(where.has_value()); + EXPECT_EQ(where->first, 8); + EXPECT_EQ(where->second, 9); + index = 9; + where = interpolator.search(x, 0.9, &index); + ASSERT_TRUE(where.has_value()); + EXPECT_EQ(where->first, 8); + EXPECT_EQ(where->second, 9); + index = 10; + where = interpolator.search(x, 0.9, &index); + ASSERT_TRUE(where.has_value()); + index = 11; + where = interpolator.search(x, 0.9, &index); + ASSERT_FALSE(where.has_value()); + where = interpolator.search(x, -0.1); + ASSERT_FALSE(where.has_value()); + where = interpolator.search(x, 1.0); +} diff --git a/src/pyinterp/core/tests/interpolation_steffen.cpp b/src/pyinterp/core/tests/interpolation_steffen.cpp new file mode 100644 index 00000000..dcd18234 --- /dev/null +++ b/src/pyinterp/core/tests/interpolation_steffen.cpp @@ -0,0 +1,197 @@ +// Copyright (c) 2024 CNES +// +// All rights reserved. Use of this source code is governed by a +// BSD-style license that can be found in the LICENSE file. +#include + +#include "pyinterp/detail/interpolation/steffen.hpp" + +TEST(Steffen, case_one) { + Eigen::Matrix xa(5); + Eigen::Matrix ya(5); + Eigen::Matrix xp(6); + Eigen::Matrix yp(6); + Eigen::Matrix dyp(6); + + xa << 0.0, 1.0, 2.0, 3.0, 4.0; + ya << 0.0, 1.0, 2.0, 3.0, 4.0; + xp << 0.0, 0.5, 1.0, 2.0, 2.5, 3.95; + yp << 0.0, 0.5, 1.0, 2.0, 2.5, 3.95; + dyp << 1.0, 1.0, 1.0, 1.0, 1.0, 1.0; + + auto interpolator = pyinterp::detail::interpolation::Steffen(); + auto y = interpolator(xa, ya, xp); + for (auto ix = 0; ix < xp.size(); ix++) { + EXPECT_NEAR(y(ix), yp(ix), 1e-6); + } + + auto dy = interpolator.derivative(xa, ya, xp); + for (auto ix = 0; ix < xp.size(); ix++) { + EXPECT_NEAR(dy(ix), dyp(ix), 1e-6); + } +} + +TEST(Steffen, case_two) { + Eigen::Matrix xa(30); + Eigen::Matrix ya(30); + Eigen::Matrix xp(129); + Eigen::Matrix yp(129); + Eigen::Matrix dyp(129); + + xa << 4.673405471947611, 4.851675778029557, 6.185473620119991, + 7.066003430727031, 7.236222118389267, 7.82067161881444, 8.054504031621493, + 8.615033792800752, 9.16666529764867, 9.442092828126395, + 10.719577968900381, 11.741967577609238, 12.564881194673035, + 12.938424721515084, 12.985892202096888, 13.995771653255744, + 14.459157783380315, 15.834857667614848, 15.887779086499942, + 18.64645279534696, 19.256361741244287, 19.85100457105047, + 19.943627313550987, 20.5033910115468, 20.895380493234843, + 21.759942483053962, 21.948984866760803, 22.2504057728461, + 22.741860290952104, 23.088720630939164; + + ya << -6.131584089652438, 4.997132870704837, 9.50359639889188, + -3.788474042493764, 1.98825324341826, -6.553460294945401, + 8.44565760012345, -3.643687601520438, -4.720262269613782, + -5.31956608504474, 4.7484698519996655, 1.8191111043001218, + 2.7234424260251124, -7.282776099988966, -6.966038807035955, + 2.6047260012112545, -3.6363485275348033, -4.124714595822157, + -4.820986899771757, -2.8550035950284602, 6.017204757558943, + 6.7275899382734075, -9.423359755713639, -0.3645325139870703, + 6.925739188986089, -0.07843256857855074, 1.6744315603866244, + -8.332083151442944, 6.059246349756446, -7.360159198725659; + + xp << 4.6734054719476106, 4.8516757780295574, 4.8575586235375265, + 5.0417117751274407, 5.2258649267173567, 5.4100180783072727, + 5.5941712298971886, 5.7783243814871028, 5.9624775330770188, + 6.1466306846669347, 6.1854736201199909, 6.3307838362568507, + 6.5149369878467667, 6.6990901394366809, 6.8832432910265968, + 7.0660034307270312, 7.0673964426165128, 7.236222118389267, + 7.2515495942064288, 7.4357027457963438, 7.6198558973862589, + 7.8040090489761749, 7.8206716188144396, 7.9881622005660908, + 8.0545040316214926, 8.172315352156005, 8.356468503745921, + 8.5406216553358369, 8.6150337928007517, 8.7247748069257529, + 8.9089279585156689, 9.0930811101055831, 9.1666652976486702, + 9.277234261695499, 9.4420928281263947, 9.461387413285415, + 9.645540564875331, 9.8296937164652469, 10.013846868055161, + 10.198000019645077, 10.382153171234993, 10.566306322824907, + 10.719577968900381, 10.750459474414823, 10.934612626004739, + 11.118765777594655, 11.302918929184571, 11.487072080774485, + 11.671225232364399, 11.741967577609238, 11.855378383954315, + 12.039531535544231, 12.223684687134147, 12.407837838724063, + 12.564881194673035, 12.591990990313979, 12.776144141903895, + 12.938424721515084, 12.960297293493811, 12.985892202096888, + 13.144450445083727, 13.328603596673641, 13.512756748263556, + 13.696909899853472, 13.881063051443387, 13.995771653255744, + 14.065216203033303, 14.249369354623219, 14.433522506213134, + 14.459157783380315, 14.617675657803051, 14.801828809392966, + 14.985981960982881, 15.170135112572796, 15.354288264162712, + 15.538441415752628, 15.722594567342544, 15.834857667614848, + 15.90674771893246, 15.887779086499942, 16.090900870522375, + 16.27505402211229, 16.459207173702204, 16.643360325292122, + 16.827513476882036, 17.01166662847195, 17.195819780061868, + 17.379972931651782, 17.5641260832417, 17.748279234831614, + 17.932432386421532, 18.116585538011446, 18.300738689601364, + 18.484891841191278, 18.64645279534696, 18.669044992781188, + 18.853198144371106, 19.03735129596102, 19.221504447550938, + 19.256361741244287, 19.405657599140852, 19.58981075073077, + 19.773963902320684, 19.851004571050471, 19.958117053910598, + 19.943627313550987, 20.142270205500516, 20.32642335709043, + 20.5033910115468, 20.510576508680348, 20.694729660270262, + 20.87888281186018, 20.895380493234843, 21.063035963450094, + 21.247189115040012, 21.431342266629926, 21.61549541821984, + 21.759942483053962, 21.799648569809754, 21.948984866760803, + 21.983801721399669, 22.167954872989586, 22.250405772846101, + 22.352108024579501, 22.536261176169418, 22.720414327759332, + 22.741860290952104, 22.90456747934925, 23.088720630939164; + yp << -6.1315840896524376, 4.9971328707048368, 5.0367975988827167, + 6.1897906340782765, 7.170975366812117, 7.9803517970842286, + 8.6179199248946077, 9.083679750243256, 9.3776312731301772, + 9.4997744935553676, 9.5035963988918795, 8.5371013292095554, + 5.3135045284256019, 1.2120062346303317, -2.3083133782071208, + -3.788474042493764, -3.7873197326712797, 1.98825324341826, + 1.9709370138809248, -0.31768851396215103, -4.2211558425051745, + -6.5330277558649303, -6.5534602949454008, 5.5087083141507147, + 8.4456576001234502, 7.1443430435268214, 1.9932653799771245, + -2.8426372908118083, -3.6436876015204378, -3.9926784426485034, + -4.3315242293057175, -4.5849882013553547, -4.7202622696137819, + -5.0157009305506106, -5.3195660850447402, -5.3127453670702938, + -4.6348437675141847, -3.1014810884327102, -1.0745634258401102, + 1.0840031242494419, 3.0123124658217071, 4.3484585028624627, + 4.7484698519996655, 4.740613456600772, 4.4142234654988952, + 3.7574719541934956, 2.9757785771330334, 2.274562988765974, + 1.8592448435407634, 1.8191111043001218, 1.8659054809992697, + 2.0883298868527498, 2.3859686670095699, 2.6372078307738107, + 2.7234424260251124, 2.5729816013335247, -3.258105131647655, + -7.2827760999889657, -7.179945071809863, -6.9660388070359547, + -5.5662805880783406, -3.3905902235190721, -1.0497645733987406, + 1.0095869363327665, 2.3408548797255095, 2.6047260012112545, + 2.2325150500282476, -0.91241316251151927, -3.5649171021138555, + -3.6363485275348033, -3.7309363021244821, -3.8038342533632736, + -3.8503835649797429, -3.8846412514674595, -3.920664327319995, + -3.9725098070309213, -4.0542347050938083, -4.1247145958221569, + -4.8208939493861829, -4.8209868997717571, -4.8103284962540727, + -4.7822416973267625, -4.7366335526042516, -4.6735040620865389, + -4.5928532257736272, -4.4946810436655156, -4.3789875157622005, + -4.2457726420636881, -4.0950364225699714, -3.9267788572810582, + -3.7409999461969403, -3.537699689317626, -3.3168780866431069, + -3.0785351381733914, -2.8550035950284602, -2.7914506045621672, + -0.46968123806956008, 3.263658433044764, 5.8622198554846658, + 6.0172047575589431, 6.3291356295622618, 6.5905310151592165, + 6.7156659481449141, 6.7275899382734075, -9.4118959280855066, + -9.4233597557136388, -7.6111870705157543, -3.9651685725226056, + -0.36453251398707032, -0.23537781181431175, 4.0332797539781309, + 6.899794364641159, 6.9257391889860891, 6.2377210430025594, + 4.3902665602806792, 2.1878722453596344, 0.44278317319805183, + -0.078432568578550743, 0.12107113682074738, 1.6744315603866244, + 1.3047436323966959, -6.4955024539509765, -8.3320831514429443, + -6.7382472104931805, 0.61052993339902706, 5.9794239504090552, + 6.0592463497564459, 1.5387265272596853, -7.3601591987256612; + dyp << 62.426083204466224, 6.757341159173202, 6.727537246743899, + 5.7945730211610407, 4.8616087955781726, 3.9286445699953054, + 2.9956803444124378, 2.0627161188295791, 1.1297518932467119, + 0.1967876676638447, 0, -12.480296842307126, -21.209126982426163, + -22.014768399615317, -14.897221093874615, 0, 1.65274069980919, 0, + -2.2393981239975669, -19.714311699236049, -19.77742868994725, + -2.4287490961310456, 0, 78.213312807624717, 0, -20.358420305108574, + -31.350476855816961, -16.93545425712459, -3.9032385156832157, + -2.5401833171916435, -1.3740294011308083, -1.6128919418291603, + -2.1012124848241789, -2.8800570317456784, 0, 0.70341285683675259, + 6.3314138906742139, 9.9941697310971538, 11.69168037810557, + 11.42394583169949, 9.1909660918788916, 4.992741158643824, 0, + -0.50358091648162351, -2.8552721239834269, -4.0914806331173992, + -4.2122064438835425, -3.2174495562818706, -1.1072099703123885, 0, + 0.78347408424733667, 1.5221058615312819, 1.6003417148468082, + 1.018181644193916, 0, -10.817925943273289, -39.490030870551372, 0, + 8.0126635338650072, 6.7986203849007927, 10.557804177353285, + 12.667135433334296, 12.351260158262027, 9.6101783521364421, + 4.44389001495753, 0, -10.1306309231183, -19.882866043398995, + -4.7826445467139278, -0.70998925548225966, -0.49283947499262398, + -0.31159279753453034, -0.20667940230488324, -0.17809928930368468, + -0.22585245853093405, -0.3499389099866323, -0.55035864367077902, + -0.70998925548225966, 0.0098004308855327432, 0, 0.10494594234665755, + 0.20009145380778143, 0.2952369652689053, 0.39038247673003107, + 0.48552798819115495, 0.58067349965227877, 0.67581901111340459, + 0.77096452257452841, 0.86611003403565412, 0.96125554549677805, + 1.0564010569579039, 1.1515465684190276, 1.2466920798801533, + 1.3418375913412772, 1.4253105022449177, 4.1661049863071913, + 18.74497348837053, 19.496500284294363, 6.4206853740787171, + 2.3892836005304425, 1.7894106566060179, 1.0494805960296585, + 0.30955053545331346, 0, 1.5724451183167569, 0, 16.386389851530527, + 21.613477536080744, 17.603560096681914, 18.338573561497682, + 23.697110333020628, 3.1105642331329868, 0, -7.5982202651552164, + -11.73098890423519, -11.453053366134201, -6.7644136508523696, 0, + 9.2309078320594402, 0, -20.350273706124355, -39.581660196515507, 0, + 28.835095151183879, 42.753317235684989, 7.3325239511975031, 0, + -47.053315974301256, -38.688209637869583; + + auto interpolator = pyinterp::detail::interpolation::Steffen(); + auto y = interpolator(xa, ya, xp); + for (auto ix = 0; ix < xp.size(); ix++) { + EXPECT_NEAR(y(ix), yp(ix), 1e-6); + } + + auto dy = interpolator.derivative(xa, ya, xp); + for (auto ix = 0; ix < xp.size(); ix++) { + EXPECT_NEAR(dy(ix), dyp(ix), 1e-6); + } +}