expose max iterations multiplier, add unit test case

pcaspers · May 19, 2024 · 541f74a · 541f74a
1 parent 14377c2
commit 541f74a
Show file tree

Hide file tree

Showing 3 changed files with 38 additions and 7 deletions.
diff --git a/ql/experimental/math/laplaceinterpolation.cpp b/ql/experimental/math/laplaceinterpolation.cpp
@@ -37,8 +37,9 @@ namespace QuantLib {
 
     LaplaceInterpolation::LaplaceInterpolation(std::function<Real(const std::vector<Size>&)> y,
                                                std::vector<std::vector<Real>> x,
-                                               const Real relTol)
-    : y_(std::move(y)), x_(std::move(x)), relTol_(relTol) {
+                                               Real relTol,
+                                               Size maxIterMultiplier)
+    : y_(std::move(y)), x_(std::move(x)), relTol_(relTol), maxIterMultiplier_(maxIterMultiplier) {
 
         // set up the mesher
 
@@ -173,7 +174,7 @@ for (auto const& m : map_)
             ++rowit;
         }
 
-        interpolatedValues_ = BiCGstab(f_A(g), 10 * N, relTol_).solve(rhs, guess).x;
+        interpolatedValues_ = BiCGstab(f_A(g), maxIterMultiplier_ * N, relTol_).solve(rhs, guess).x;
     }
 
     std::vector<Size>
@@ -213,7 +214,8 @@ for (auto const& m : map_)
     void laplaceInterpolation(Matrix& A,
                               const std::vector<Real>& x,
                               const std::vector<Real>& y,
-                              Real relTol) {
+                              Real relTol,
+                              Size maxIterMultiplier) {
 
         std::vector<std::vector<Real>> tmp;
         tmp.push_back(y);
@@ -233,7 +235,7 @@ for (auto const& m : map_)
             [&A](const std::vector<Size>& coordinates) {
                 return A(coordinates[0], coordinates[1]);
             },
-            tmp, relTol);
+            tmp, relTol, maxIterMultiplier);
 
         for (Size i = 0; i < A.rows(); ++i) {
             for (Size j = 0; j < A.columns(); ++j) {

diff --git a/ql/experimental/math/laplaceinterpolation.hpp b/ql/experimental/math/laplaceinterpolation.hpp
@@ -44,7 +44,8 @@ namespace QuantLib {
         /*! Missing values y should be encoded as Null<Real>(). */
         LaplaceInterpolation(std::function<Real(const std::vector<Size>&)> y,
                              std::vector<std::vector<Real>> x,
-                             Real relTol = 1E-6);
+                             Real relTol = 1E-6,
+                             Size maxIterMultiplier = 10);
         Real operator()(const std::vector<Size>& coordinates) const;
 
       private:
@@ -54,6 +55,7 @@ namespace QuantLib {
         std::function<Real(const std::vector<Size>&)> y_;
         std::vector<std::vector<Real>> x_;
         Real relTol_;
+        Size maxIterMultiplier_;
 
         std::vector<bool> coordinateIncluded_;
         Size numberOfCoordinatesIncluded_;
@@ -68,7 +70,8 @@ namespace QuantLib {
     void laplaceInterpolation(Matrix& A,
                               const std::vector<Real>& x = {},
                               const std::vector<Real>& y = {},
-                              Real relTol = 1E-6);
+                              Real relTol = 1E-6,
+                              Size maxIterMultiplier = 10);
 }
 
 #endif
diff --git a/test-suite/interpolations.cpp b/test-suite/interpolations.cpp
@@ -2760,6 +2760,32 @@ BOOST_AUTO_TEST_CASE(testLaplaceInterpolation) {
 
     LaplaceInterpolation l0([](const std::vector<Size>& x) { return Null<Real>(); }, {});
     BOOST_CHECK_CLOSE(l0({}), 0.0, tol);
+
+    // single test cases from actual issues observed in the field
+
+    std::vector<Real> tx = {0.0849315, 0.257534, 0.509589, 1.00548, 2.00274, 3.00274, 4.00274,
+                            5.00548,   7.00822,  10.0082,  15.011,  20.0137, 30.0219, 70.0493};
+    std::vector<Real> ty = {0.25, 1, 2, 3, 4, 5, 7, 10, 15, 20, 30, 100};
+    Matrix m52 = {{na, na, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, na},
+                  {na, na, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, na},
+                  {na, na, na, na, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, na},
+                  {na, na, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, na},
+                  {na, na, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, na},
+                  {na, na, na, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, na},
+                  {na, na, na, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, na},
+                  {na, na, na, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, na},
+                  {na, na, na, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, na, na},
+                  {na, na, na, na, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, na, na},
+                  {na, na, na, na, na, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, na, na},
+                  {na, na, na, na, na, na, na, na, na, na, na, na, na, na}};
+
+    // we need to allow for more iterations to achieve the desired accuracy
+    laplaceInterpolation(m52, tx, ty, 1E-6, 100);
+
+    for (auto const& v : m52) {
+        BOOST_CHECK_CLOSE(v, 1.0, 0.1);
+    }
+
 }
 
 BOOST_AUTO_TEST_SUITE_END()