diff --git a/QuantLib.vcxproj b/QuantLib.vcxproj
index a3d24659604..54d1b917cb6 100644
--- a/QuantLib.vcxproj
+++ b/QuantLib.vcxproj
@@ -2073,6 +2073,7 @@
+
diff --git a/QuantLib.vcxproj.filters b/QuantLib.vcxproj.filters
index cf3d3abb3ff..327ce92335c 100644
--- a/QuantLib.vcxproj.filters
+++ b/QuantLib.vcxproj.filters
@@ -6173,6 +6173,9 @@
legacy\libormarketmodels
+
+ experimental\math
+
experimental\volatility
diff --git a/ql/CMakeLists.txt b/ql/CMakeLists.txt
index 4647a233626..d378e6a5a4c 100644
--- a/ql/CMakeLists.txt
+++ b/ql/CMakeLists.txt
@@ -173,6 +173,7 @@ set(QL_SOURCES
experimental/math/fireflyalgorithm.cpp
experimental/math/gaussiancopulapolicy.cpp
experimental/math/gaussiannoncentralchisquaredpolynomial.cpp
+ experimental/math/laplaceinterpolation.cpp
experimental/math/multidimintegrator.cpp
experimental/math/multidimquadrature.cpp
experimental/math/particleswarmoptimization.cpp
diff --git a/ql/experimental/math/Makefile.am b/ql/experimental/math/Makefile.am
index 0c4827c8050..9a0473eb326 100644
--- a/ql/experimental/math/Makefile.am
+++ b/ql/experimental/math/Makefile.am
@@ -32,6 +32,7 @@ cpp_files = \
fireflyalgorithm.cpp \
gaussiancopulapolicy.cpp \
gaussiannoncentralchisquaredpolynomial.cpp \
+ laplaceinterpolation.cpp \
multidimintegrator.cpp \
multidimquadrature.cpp \
particleswarmoptimization.cpp \
diff --git a/ql/experimental/math/laplaceinterpolation.cpp b/ql/experimental/math/laplaceinterpolation.cpp
new file mode 100644
index 00000000000..93748bd9a65
--- /dev/null
+++ b/ql/experimental/math/laplaceinterpolation.cpp
@@ -0,0 +1,245 @@
+/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
+
+/*
+ Copyright (C) 2015, 2024 Peter Caspers
+
+ This file is part of QuantLib, a free-software/open-source library
+ for financial quantitative analysts and developers - http://quantlib.org/
+
+ QuantLib is free software: you can redistribute it and/or modify it
+ under the terms of the QuantLib license. You should have received a
+ copy of the license along with this program; if not, please email
+ . The license is also available online at
+ .
+
+ This program is distributed in the hope that it will be useful, but WITHOUT
+ ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+ FOR A PARTICULAR PURPOSE. See the license for more details.
+*/
+
+/*! \file laplaceinterpolation.hpp
+ \brief Laplace interpolation of missing values
+*/
+
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+
+namespace QuantLib {
+
+ LaplaceInterpolation::LaplaceInterpolation(std::function&)> y,
+ std::vector> x,
+ const Real relTol)
+ : y_(std::move(y)), x_(std::move(x)), relTol_(relTol) {
+
+ // set up the mesher
+
+ std::vector dim;
+ coordinateIncluded_.resize(x_.size());
+ for (Size i = 0; i < x_.size(); ++i) {
+ coordinateIncluded_[i] = x_[i].size() > 1;
+ if (coordinateIncluded_[i])
+ dim.push_back(x_[i].size());
+ }
+
+ numberOfCoordinatesIncluded_ = dim.size();
+
+ if (numberOfCoordinatesIncluded_ == 0) {
+ return;
+ }
+
+ QL_REQUIRE(!dim.empty(), "LaplaceInterpolation: singular point or no points given");
+
+ layout_ = ext::make_shared(dim);
+
+ std::vector> meshers;
+ for (Size i = 0; i < x_.size(); ++i) {
+ if (x_[i].size() > 1)
+ meshers.push_back(ext::make_shared(x_[i]));
+ }
+
+ auto mesher = ext::make_shared(layout_, meshers);
+
+ // set up the Laplace operator and convert it to matrix
+
+ struct LaplaceOp : public FdmLinearOpComposite {
+ explicit LaplaceOp(const ext::shared_ptr& mesher) {
+ for (Size direction = 0; direction < mesher->layout()->dim().size(); ++direction) {
+ if (mesher->layout()->dim()[direction] > 1)
+ map_.push_back(SecondDerivativeOp(direction, mesher));
+ }
+ }
+ std::vector map_;
+
+ Size size() const override { QL_FAIL("no impl"); }
+ void setTime(Time t1, Time t2) override { QL_FAIL("no impl"); }
+ Array apply(const array_type& r) const override { QL_FAIL("no impl"); }
+ Array apply_mixed(const Array& r) const override { QL_FAIL("no impl"); }
+ Array apply_direction(Size direction, const Array& r) const override {
+ QL_FAIL("no impl");
+ }
+ Array solve_splitting(Size direction, const Array& r, Real s) const override {
+ QL_FAIL("no impl");
+ }
+ Array preconditioner(const Array& r, Real s) const override { QL_FAIL("no impl"); }
+ std::vector toMatrixDecomp() const override {
+ std::vector decomp;
+ for (auto const& m : map_)
+ decomp.push_back(m.toMatrix());
+ return decomp;
+ }
+ };
+
+ SparseMatrix op = LaplaceOp(mesher).toMatrix();
+
+ // set up the linear system to solve
+
+ Size N = layout_->size();
+
+ SparseMatrix g(N, N, 5 * N);
+ Array rhs(N, 0.0), guess(N, 0.0);
+ Real guessTmp = 0.0;
+
+ struct f_A {
+ const SparseMatrix& g;
+ explicit f_A(const SparseMatrix& g) : g(g) {}
+ Array operator()(const Array& x) const { return prod(g, x); }
+ };
+
+ auto rowit = op.begin1();
+ Size count = 0;
+ std::vector corner_h(dim.size());
+ std::vector corner_neighbour_index(dim.size());
+ for (auto const& pos : *layout_) {
+ auto coord = pos.coordinates();
+ Real val =
+ y_(numberOfCoordinatesIncluded_ == x_.size() ? coord : fullCoordinates(coord));
+ QL_REQUIRE(rowit != op.end1() && rowit.index1() == count,
+ "LaplaceInterpolation: op matrix row iterator ("
+ << (rowit != op.end1() ? std::to_string(rowit.index1()) : "na")
+ << ") does not match expected row count (" << count << ")");
+ if (val == Null()) {
+ bool isCorner = true;
+ for (Size d = 0; d < dim.size() && isCorner; ++d) {
+ if (coord[d] == 0) {
+ corner_h[d] = meshers[d]->dplus(0);
+ corner_neighbour_index[d] = 1;
+ } else if (coord[d] == layout_->dim()[d] - 1) {
+ corner_h[d] = meshers[d]->dminus(dim[d] - 1);
+ corner_neighbour_index[d] = dim[d] - 2;
+ } else {
+ isCorner = false;
+ }
+ }
+ if (isCorner) {
+ // handling of the "corners", all second derivs are zero in the op
+ // this directly generalizes Numerical Recipes, 3rd ed, eq 3.8.6
+ Real sum_corner_h =
+ std::accumulate(corner_h.begin(), corner_h.end(), 0.0, std::plus());
+ for (Size j = 0; j < dim.size(); ++j) {
+ std::vector coord_j(coord);
+ coord_j[j] = corner_neighbour_index[j];
+ Real weight = 0.0;
+ for (Size i = 0; i < dim.size(); ++i) {
+ if (i != j)
+ weight += corner_h[i];
+ }
+ weight = dim.size() == 1 ? 1.0 : weight / sum_corner_h;
+ g(count, layout_->index(coord_j)) = -weight;
+ }
+ g(count, count) = 1.0;
+ } else {
+ // point with at least one dimension with non-trivial second derivative
+ for (auto colit = rowit.begin(); colit != rowit.end(); ++colit)
+ g(count, colit.index2()) = *colit;
+ }
+ rhs[count] = 0.0;
+ guess[count] = guessTmp;
+ } else {
+ g(count, count) = 1;
+ rhs[count] = val;
+ guess[count] = guessTmp = val;
+ }
+ ++count;
+ ++rowit;
+ }
+
+ interpolatedValues_ = BiCGstab(f_A(g), 10 * N, relTol_).solve(rhs, guess).x;
+ }
+
+ std::vector
+ LaplaceInterpolation::projectedCoordinates(const std::vector& coordinates) const {
+ std::vector tmp;
+ for (Size i = 0; i < coordinates.size(); ++i) {
+ if (coordinateIncluded_[i])
+ tmp.push_back(coordinates[i]);
+ }
+ return tmp;
+ }
+
+ std::vector
+ LaplaceInterpolation::fullCoordinates(const std::vector& projectedCoordinates) const {
+ std::vector tmp(coordinateIncluded_.size(), 0);
+ for (Size i = 0, count = 0; i < coordinateIncluded_.size(); ++i) {
+ if (coordinateIncluded_[i])
+ tmp[i] = projectedCoordinates[count++];
+ }
+ return tmp;
+ }
+
+ Real LaplaceInterpolation::operator()(const std::vector& coordinates) const {
+ QL_REQUIRE(coordinates.size() == x_.size(), "LaplaceInterpolation::operator(): expected "
+ << x_.size() << " coordinates, got "
+ << coordinates.size());
+ if (numberOfCoordinatesIncluded_ == 0) {
+ Real val = y_(coordinates);
+ return val == Null() ? 0.0 : val;
+ } else {
+ return interpolatedValues_[layout_->index(numberOfCoordinatesIncluded_ == x_.size() ?
+ coordinates :
+ projectedCoordinates(coordinates))];
+ }
+ }
+
+ void laplaceInterpolation(Matrix& A,
+ const std::vector& x,
+ const std::vector& y,
+ Real relTol) {
+
+ std::vector> tmp;
+ tmp.push_back(y);
+ tmp.push_back(x);
+
+ if (y.empty()) {
+ tmp[0].resize(A.rows());
+ std::iota(tmp[0].begin(), tmp[0].end(), 0.0);
+ }
+
+ if (x.empty()) {
+ tmp[1].resize(A.columns());
+ std::iota(tmp[1].begin(), tmp[1].end(), 0.0);
+ }
+
+ LaplaceInterpolation interpolation(
+ [&A](const std::vector& coordinates) {
+ return A(coordinates[0], coordinates[1]);
+ },
+ tmp, relTol);
+
+ for (Size i = 0; i < A.rows(); ++i) {
+ for (Size j = 0; j < A.columns(); ++j) {
+ if (A(i, j) == Null())
+ A(i, j) = interpolation({i, j});
+ }
+ }
+ }
+
+} // namespace QuantLib
diff --git a/ql/experimental/math/laplaceinterpolation.hpp b/ql/experimental/math/laplaceinterpolation.hpp
index 831488de5c2..ee0745802eb 100644
--- a/ql/experimental/math/laplaceinterpolation.hpp
+++ b/ql/experimental/math/laplaceinterpolation.hpp
@@ -1,7 +1,7 @@
/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
/*
- Copyright (C) 2015 Peter Caspers
+ Copyright (C) 2015, 2024 Peter Caspers
This file is part of QuantLib, a free-software/open-source library
for financial quantitative analysts and developers - http://quantlib.org/
@@ -24,140 +24,51 @@
#ifndef quantlib_laplace_interpolation
#define quantlib_laplace_interpolation
-#include
-#include
+#include
+#include
+#include
+#include
+
+#include
namespace QuantLib {
-/*! reference: Numerical Recipes, 3rd edition, ch. 3.8
- two dimensional reconstruction of missing (i.e. null)
- values using laplace interpolation assuming an
- equidistant grid */
+ class FdmLinearOpLayout;
+
+ /*! Reconstruction of missing values using Laplace interpolation. We support an arbitrary number
+ of dimensions n >= 1 and non-equidistant grids. For n = 1 the method is identical to linear
+ interpolation with flat extrapolation. Reference: Numerical Recipes, 3rd edition, ch. 3.8. */
+
+ class LaplaceInterpolation {
+ public:
+ /*! Missing values y should be encoded as Null(). */
+ LaplaceInterpolation(std::function&)> y,
+ std::vector> x,
+ const Real relTol = 1E-6);
+ Real operator()(const std::vector& coordinates) const;
+
+ private:
+ std::vector projectedCoordinates(const std::vector& coordinates) const;
+ std::vector fullCoordinates(const std::vector& projectedCoordinates) const;
-template void laplaceInterpolation(M &A, Real relTol = 1E-6) {
+ std::function&)> y_;
+ std::vector> x_;
+ Real relTol_;
- struct f_A {
- const SparseMatrix &g;
- explicit f_A(const SparseMatrix &g) : g(g) {}
- Array operator()(const Array &x) const {
- return prod(g, x);
- }
+ std::vector coordinateIncluded_;
+ Size numberOfCoordinatesIncluded_;
+
+ ext::shared_ptr layout_;
+ Array interpolatedValues_;
};
- Size m = A.rows();
- Size n = A.columns();
-
- QL_REQUIRE(n > 1 && m > 1, "matrix (" << m << "," << n
- << ") must at least be 2x2");
-
- SparseMatrix g(m * n, m * n, 5 * m * n);
- Array rhs(m * n, 0.0), guess(m * n, 0.0);
- Real guessTmp = 0.0;
- Size i1, i2, i3, i4, j1, j2, j3, j4;
- bool inner;
-
- for (Size l = 0, i = 0; i < m; ++i) {
- for (Size j = 0; j < n; ++j) {
-
- inner = false;
-
- // top
- if (i == 0) {
- if (j == 0) {
- i1 = 0;
- j1 = 1;
- i2 = 1;
- j2 = 0;
- } else {
- if (j == n - 1) {
- i1 = 0;
- j1 = n - 2;
- i2 = 1;
- j2 = n - 1;
- } else {
- i1 = i2 = 0;
- j1 = j - 1;
- j2 = j + 1;
- }
- }
- }
-
- // bottom
- if (i == m - 1) {
- if (j == 0) {
- i1 = m - 1;
- j1 = 1;
- i2 = m - 2;
- j2 = 0;
- } else {
- if (j == n - 1) {
- i1 = m - 1;
- j1 = n - 2;
- i2 = m - 2;
- j2 = n - 1;
- } else {
- i1 = i2 = m - 1;
- j1 = j - 1;
- j2 = j + 1;
- }
- }
- }
-
- // left / right
- if (i > 0 && i < m - 1) {
- if (j == 0 || j == n - 1) {
- j1 = j2 = j;
- i1 = i - 1;
- i2 = i + 1;
- } else {
- inner = true;
- i1 = i - 1;
- i2 = i - 1;
- i3 = i + 1;
- i4 = i + 1;
- j1 = j - 1;
- j2 = j + 1;
- j3 = j - 1;
- j4 = j + 1;
- }
- }
-
- g(l, i * n + j) = 1.0;
- if (A[i][j] == Null()) {
- if (inner) {
- g(l, i1 * n + j1) = -0.25;
- g(l, i2 * n + j2) = -0.25;
- g(l, i3 * n + j3) = -0.25;
- g(l, i4 * n + j4) = -0.25;
- } else {
- g(l, i1 * n + j1) = -0.5;
- g(l, i2 * n + j2) = -0.5;
- }
- rhs[l] = 0.0;
- guess[l] = guessTmp;
- } else {
- rhs[l] = A[i][j];
- guess[l] = guessTmp = A[i][j];
- }
- ++l;
- }
- }
-
- // solve the equation (preconditioner is identiy)
- Array s = BiCGstab(f_A(g), 10 * m * n, relTol)
- .solve(rhs, guess)
- .x;
-
- // replace missing values by solution
- for (Size i = 0; i < m; ++i) {
- for (Size j = 0; j < n; ++j) {
- if (A[i][j] == Null()) {
- A[i][j] = s[i * n + j];
- }
- }
- }
-};
-
-} // namespace QuantLib
-
-#endif // include guard
+ /*! Convenience function that Laplace-interpolates null values in a given matrix.
+ If the x or y grid or both are not given, an equidistant grid is assumed. */
+
+ void laplaceInterpolation(Matrix& A,
+ const std::vector& x = {},
+ const std::vector& y = {},
+ Real relTol = 1E-6);
+}
+
+#endif
diff --git a/test-suite/interpolations.cpp b/test-suite/interpolations.cpp
index acfc872114a..6f0bb135dc4 100644
--- a/test-suite/interpolations.cpp
+++ b/test-suite/interpolations.cpp
@@ -26,6 +26,7 @@
#include "toplevelfixture.hpp"
#include "utilities.hpp"
#include
+#include
#include
#include
#include
@@ -2501,6 +2502,266 @@ BOOST_AUTO_TEST_CASE(testChebyshevInterpolationUpdateY) {
}
}
+BOOST_AUTO_TEST_CASE(testLaplaceInterpolation) {
+ BOOST_TEST_MESSAGE("Testing Laplace interpolation...");
+
+ Real tol = 1E-12;
+ Real na = Null();
+
+ // full matrix
+
+ Matrix m1 = {
+ {1.0, 2.0, 4.0},
+ {6.0, 6.5, 7.0},
+ {5.0, 3.0, 2.0},
+ };
+
+ laplaceInterpolation(m1);
+
+ BOOST_CHECK_CLOSE(m1(0, 0), 1.0, tol);
+ BOOST_CHECK_CLOSE(m1(0, 2), 4.0, tol);
+ BOOST_CHECK_CLOSE(m1(2, 0), 5.0, tol);
+ BOOST_CHECK_CLOSE(m1(2, 1), 3.0, tol);
+
+ // inner point
+
+ Matrix m2 = {
+ {1.0, 2.0, 4.0},
+ {6.0, na, 7.0},
+ {5.0, 3.0, 2.0},
+ };
+
+ laplaceInterpolation(m2);
+
+ BOOST_CHECK_CLOSE(m2(1, 1), 4.5, tol);
+ BOOST_CHECK_CLOSE(m2(2, 1), 3.0, tol);
+
+ // boundaries
+
+ Matrix m3 = {
+ {1.0, na, 4.0},
+ {6.0, 6.5, 7.0},
+ {5.0, 3.0, 2.0},
+ };
+
+ laplaceInterpolation(m3);
+
+ BOOST_CHECK_CLOSE(m3(0, 1), 2.5, tol);
+
+ Matrix m4 = {
+ {1.0, 2.0, 4.0},
+ {na, 6.5, 7.0},
+ {5.0, 3.0, 2.0},
+ };
+
+ laplaceInterpolation(m4);
+
+ BOOST_CHECK_CLOSE(m4(1, 0), 3.0, tol);
+
+ Matrix m5 = {
+ {1.0, 2.0, 4.0},
+ {6.0, 6.5, 7.0},
+ {5.0, na, 2.0},
+ };
+
+ laplaceInterpolation(m5);
+
+ BOOST_CHECK_CLOSE(m5(2, 1), 3.5, tol);
+
+ Matrix m6 = {
+ {1.0, 2.0, 4.0},
+ {6.0, 6.5, na},
+ {5.0, 3.0, 2.0},
+ };
+
+ laplaceInterpolation(m6);
+
+ BOOST_CHECK_CLOSE(m6(1, 2), 3.0, tol);
+
+ // corners
+
+ Matrix m7 = {
+ {na, 2.0, 4.0},
+ {6.0, 6.5, 7.0},
+ {5.0, 3.0, 2.0},
+ };
+
+ laplaceInterpolation(m7);
+
+ BOOST_CHECK_CLOSE(m7(0, 0), 4.0, tol);
+
+ Matrix m8 = {
+ {1.0, 2.0, 4.0},
+ {6.0, 6.5, 7.0},
+ {na, 3.0, 2.0},
+ };
+
+ laplaceInterpolation(m8);
+
+ BOOST_CHECK_CLOSE(m8(2, 0), 4.5, tol);
+
+ Matrix m9 = {
+ {1.0, 2.0, 4.0},
+ {6.0, 6.5, 7.0},
+ {5.0, 3.0, na},
+ };
+
+ laplaceInterpolation(m9);
+
+ BOOST_CHECK_CLOSE(m9(2, 2), 5.0, tol);
+
+ Matrix m10 = {
+ {1.0, 2.0, na},
+ {6.0, 6.5, 7.0},
+ {5.0, 3.0, 2.0},
+ };
+
+ laplaceInterpolation(m10);
+
+ BOOST_CHECK_CLOSE(m10(0, 2), 4.5, tol);
+
+ // one dim (col vector)
+
+ Matrix m20 = { {na}, {na}, {3.0}, {5.0}, {7.0}, {na} };
+
+ laplaceInterpolation(m20);
+
+ BOOST_CHECK_CLOSE(m20(0, 0), 3.0, tol);
+ BOOST_CHECK_CLOSE(m20(1, 0), 3.0, tol);
+ BOOST_CHECK_CLOSE(m20(5, 0), 7.0, tol);
+
+ // one dim (row vector)
+
+ Matrix m21 = {{na, na, 3.0, 5.0, 7.0, na}};
+
+ laplaceInterpolation(m21);
+
+ BOOST_CHECK_CLOSE(m21(0, 0), 3.0, tol);
+ BOOST_CHECK_CLOSE(m21(0, 1), 3.0, tol);
+ BOOST_CHECK_CLOSE(m21(0, 5), 7.0, tol);
+
+ // non equidistant grid, inner point
+
+ Matrix m30 = {
+ {1.0, 2.0, 4.0},
+ {6.0, na, 7.0},
+ {5.0, 3.0, 2.0},
+ };
+
+ laplaceInterpolation(m30, {1.0, 2.0, 4.0}, {1.0, 2.0, 4.0});
+
+ BOOST_CHECK_CLOSE(m30(1, 1), 26.0 / 6.0, tol);
+
+ // non equidistant grid, boundaries
+
+ Matrix m31 = {
+ {1.0, na, 4.0},
+ {6.0, 6.5, 7.0},
+ {5.0, 3.0, 2.0},
+ };
+
+ laplaceInterpolation(m31, {1.0, 2.0, 4.0}, {1.0, 2.0, 4.0});
+
+ BOOST_CHECK_CLOSE(m31(0, 1), 6.0 / 3.0, tol);
+
+ Matrix m32 = {
+ {1.0, 2.0, 4.0},
+ {na, 6.5, 7.0},
+ {5.0, 3.0, 2.0},
+ };
+
+ laplaceInterpolation(m32, {1.0, 2.0, 4.0}, {1.0, 2.0, 4.0});
+
+ BOOST_CHECK_CLOSE(m32(1, 0), 7.0 / 3.0, tol);
+
+ Matrix m33 = {
+ {1.0, 2.0, 4.0},
+ {6.0, 6.5, 7.0},
+ {5.0, na, 2.0},
+ };
+
+ laplaceInterpolation(m33, {1.0, 2.0, 4.0}, {1.0, 2.0, 4.0});
+
+ BOOST_CHECK_CLOSE(m33(2, 1), 12.0 / 3.0, tol);
+
+ Matrix m34 = {
+ {1.0, 2.0, 4.0},
+ {6.0, 6.5, na},
+ {5.0, 3.0, 2.0},
+ };
+
+ laplaceInterpolation(m34, {1.0, 2.0, 4.0}, {1.0, 2.0, 4.0});
+
+ BOOST_CHECK_CLOSE(m34(1, 2), 10.0 / 3.0, tol);
+
+ // non equidistant grid, corners
+
+ Matrix m35 = {
+ {na, 2.0, 4.0},
+ {6.0, 6.5, 7.0},
+ {5.0, 3.0, 2.0},
+ };
+
+ laplaceInterpolation(m35, {1.0, 2.0, 4.0}, {1.0, 3.0, 7.0});
+
+ BOOST_CHECK_CLOSE(m35(0, 0), 10.0 / 3.0, tol);
+
+ Matrix m36 = {
+ {1.0, 2.0, 4.0},
+ {6.0, 6.5, 7.0},
+ {na, 3.0, 2.0},
+ };
+
+ laplaceInterpolation(m36, {1.0, 2.0, 4.0}, {1.0, 3.0, 7.0});
+
+ BOOST_CHECK_CLOSE(m36(2, 0), 18.0 / 5.0, tol);
+
+ Matrix m37 = {
+ {1.0, 2.0, 4.0},
+ {6.0, 6.5, 7.0},
+ {5.0, 3.0, na},
+ };
+
+ laplaceInterpolation(m37, {1.0, 2.0, 4.0}, {1.0, 3.0, 7.0});
+
+ BOOST_CHECK_CLOSE(m37(2, 2), 13.0 / 3.0, tol);
+
+ Matrix m38 = {
+ {1.0, 2.0, na},
+ {6.0, 6.5, 7.0},
+ {5.0, 3.0, 2.0},
+ };
+
+ laplaceInterpolation(m38, {1.0, 2.0, 4.0}, {1.0, 2.0, 3.0});
+
+ BOOST_CHECK_CLOSE(m38(0, 2), 16.0 / 3.0, tol);
+
+ // single point with given value
+
+ Matrix m50 = {
+ {1.0},
+ };
+
+ laplaceInterpolation(m50);
+
+ BOOST_CHECK_CLOSE(m50(0, 0), 1.0, tol);
+
+ // single point with missing value
+
+ Matrix m51 = {
+ {Null()},
+ };
+
+ laplaceInterpolation(m51);
+
+ BOOST_CHECK_CLOSE(m51(0, 0), 0.0, tol);
+
+ // no point
+
+ LaplaceInterpolation l0([](const std::vector& x) { return Null(); }, {});
+ BOOST_CHECK_CLOSE(l0({}), 0.0, tol);
+}
+
BOOST_AUTO_TEST_SUITE_END()
BOOST_AUTO_TEST_SUITE_END()