From 628ebc5810accc29eab68de65fb78689682e0af0 Mon Sep 17 00:00:00 2001 From: Peter Caspers Date: Tue, 19 Mar 2024 14:54:08 +0100 Subject: [PATCH 1/5] handle non-equidistant grids and arbitrary dimensions --- ql/CMakeLists.txt | 1 + ql/experimental/math/laplaceinterpolation.cpp | 253 +++++++++++++++++ ql/experimental/math/laplaceinterpolation.hpp | 173 +++--------- test-suite/interpolations.cpp | 261 ++++++++++++++++++ 4 files changed, 557 insertions(+), 131 deletions(-) create mode 100644 ql/experimental/math/laplaceinterpolation.cpp 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/laplaceinterpolation.cpp b/ql/experimental/math/laplaceinterpolation.cpp new file mode 100644 index 00000000000..9f5f2a9b6ab --- /dev/null +++ b/ql/experimental/math/laplaceinterpolation.cpp @@ -0,0 +1,253 @@ +/* -*- 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 { + 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) { + // shandling 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))]; + } + } + + template + void laplaceInterpolation(M& 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}); + } + } + } + + // template instantiations for matrix classes we want to support + + template void laplaceInterpolation(Matrix& A, + const std::vector& x, + const std::vector& y, + Real relTol); + +} // namespace QuantLib diff --git a/ql/experimental/math/laplaceinterpolation.hpp b/ql/experimental/math/laplaceinterpolation.hpp index 831488de5c2..6ce1d8cc8d3 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 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. */ + + template + void laplaceInterpolation(M& 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() From 3e1e10ffea614dd3da1f4ce1b86db83bcbd27858 Mon Sep 17 00:00:00 2001 From: Peter Caspers Date: Tue, 19 Mar 2024 15:02:48 +0100 Subject: [PATCH 2/5] update project files, fix typo --- QuantLib.vcxproj | 1 + QuantLib.vcxproj.filters | 3 +++ ql/experimental/math/laplaceinterpolation.cpp | 2 +- 3 files changed, 5 insertions(+), 1 deletion(-) 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/experimental/math/laplaceinterpolation.cpp b/ql/experimental/math/laplaceinterpolation.cpp index 9f5f2a9b6ab..38a74064283 100644 --- a/ql/experimental/math/laplaceinterpolation.cpp +++ b/ql/experimental/math/laplaceinterpolation.cpp @@ -140,7 +140,7 @@ namespace QuantLib { } } if (isCorner) { - // shandling of the "corners", all second derivs are zero in the op + // 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()); From 93586e09011074f8d01976315ea067e5e3abe3d8 Mon Sep 17 00:00:00 2001 From: Peter Caspers Date: Tue, 19 Mar 2024 15:10:48 +0100 Subject: [PATCH 3/5] update make file --- ql/experimental/math/Makefile.am | 1 + 1 file changed, 1 insertion(+) 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 \ From 85448fc9950e3bf71e891804d37ff4440a2565ab Mon Sep 17 00:00:00 2001 From: Peter Caspers Date: Tue, 19 Mar 2024 15:27:48 +0100 Subject: [PATCH 4/5] make ctor explicit --- ql/experimental/math/laplaceinterpolation.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ql/experimental/math/laplaceinterpolation.cpp b/ql/experimental/math/laplaceinterpolation.cpp index 38a74064283..278ee534757 100644 --- a/ql/experimental/math/laplaceinterpolation.cpp +++ b/ql/experimental/math/laplaceinterpolation.cpp @@ -71,7 +71,7 @@ namespace QuantLib { // set up the Laplace operator and convert it to matrix struct LaplaceOp : public FdmLinearOpComposite { - LaplaceOp(const ext::shared_ptr& mesher) { + 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)); From 5727400e1738fdf0935c0e600eaad01c77ab8922 Mon Sep 17 00:00:00 2001 From: Peter Caspers Date: Thu, 21 Mar 2024 15:07:03 +0100 Subject: [PATCH 5/5] remove template version, just support QuantLib::Matrix --- ql/experimental/math/laplaceinterpolation.cpp | 10 +--------- ql/experimental/math/laplaceinterpolation.hpp | 4 ++-- 2 files changed, 3 insertions(+), 11 deletions(-) diff --git a/ql/experimental/math/laplaceinterpolation.cpp b/ql/experimental/math/laplaceinterpolation.cpp index 278ee534757..93748bd9a65 100644 --- a/ql/experimental/math/laplaceinterpolation.cpp +++ b/ql/experimental/math/laplaceinterpolation.cpp @@ -209,8 +209,7 @@ namespace QuantLib { } } - template - void laplaceInterpolation(M& A, + void laplaceInterpolation(Matrix& A, const std::vector& x, const std::vector& y, Real relTol) { @@ -243,11 +242,4 @@ namespace QuantLib { } } - // template instantiations for matrix classes we want to support - - template void laplaceInterpolation(Matrix& A, - const std::vector& x, - const std::vector& y, - Real relTol); - } // namespace QuantLib diff --git a/ql/experimental/math/laplaceinterpolation.hpp b/ql/experimental/math/laplaceinterpolation.hpp index 6ce1d8cc8d3..ee0745802eb 100644 --- a/ql/experimental/math/laplaceinterpolation.hpp +++ b/ql/experimental/math/laplaceinterpolation.hpp @@ -25,6 +25,7 @@ #define quantlib_laplace_interpolation #include +#include #include #include @@ -64,8 +65,7 @@ namespace QuantLib { /*! 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. */ - template - void laplaceInterpolation(M& A, + void laplaceInterpolation(Matrix& A, const std::vector& x = {}, const std::vector& y = {}, Real relTol = 1E-6);