TabFunction: Added optional interpolation type argument

lib/TabFunction.cpp

@@ -27,14 +27,16 @@ using namespace Latan;
  *                       TabFunction implementation                           *
  ******************************************************************************/
 // constructors ////////////////////////////////////////////////////////////////
-TabFunction::TabFunction(const DVec &x, const DVec &y)
-: TabFunction()
+TabFunction::TabFunction(const DVec &x, const DVec &y,
+                         const InterpType interpType)
+: interpType_(interpType)
 {
     setData(x, y);
 }
 
-TabFunction::TabFunction(const XYStatData &data, const Index i, const Index j)
-: TabFunction()
+TabFunction::TabFunction(const XYStatData &data, const Index i, const Index j,
+                         const InterpType interpType)
+: interpType_(interpType)
 {
     setData(data, i, j);
 }
@@ -60,26 +62,52 @@ void TabFunction::setData(const XYStatData &data, const Index i, const Index j)
 // function call ///////////////////////////////////////////////////////////////
 double TabFunction::operator()(const double *arg) const
 {
-    double x = arg[0], x_a, x_b, y_a, y_b;
-    
-    if ((x < value_.begin()->first)||(x >= value_.rbegin()->first))
-    {
-        LATAN_ERROR(Range, "tabulated function variable out of range (x= "
-                    + strFrom(x) + " not in ["
-                    + strFrom(value_.begin()->first) + ", "
-                    + strFrom(value_.rbegin()->first) + "])");
+    double result = 0.0, x = arg[0];
+
+
+    if ((x < value_.begin()->first) || (x >= value_.rbegin()->first)) {
+        LATAN_ERROR(Range, "tabulated function variable out of range "
+                               "(x= " + strFrom(x) + " not in ["
+                           + strFrom(value_.begin()->first) + ", "
+                           + strFrom(value_.rbegin()->first) + "])");
+    }
+
+    switch (interpType_) {
+        case InterpType::LINEAR: {
+            double x_a, x_b, y_a, y_b;
+
+            auto i = value_.equal_range(x);
+            auto low = (x == i.first->first) ? i.first : prev(i.first);
+            auto high = i.second;
+
+            x_a = low->first;
+            x_b = high->first;
+            y_a = low->second;
+            y_b = high->second;
+            result = y_a + (x - x_a) * (y_b - y_a) / (x_b - x_a);
+            break;
+        }
+        case InterpType::NEAREST: {
+            auto it = value_.upper_bound(x);
+            auto upper_pair = *it;
+            auto lower_pair = *(--it);
+            if (fabs(upper_pair.first - x) < fabs(lower_pair.first - x)) {
+                result = upper_pair.second;
+            }
+            else {
+                result = lower_pair.second;
+            }
+            break;
+        }
+        case InterpType::QUADRATIC:
+        default:
+            int intType = static_cast<int>(interpType_);
+            LATAN_ERROR(Implementation, "unsupported interpolation type in "
+                                        "tabulated function: "
+                                        + strFrom(intType));
     }
     
-    auto i    = value_.equal_range(x);
-    auto low  = (x == i.first->first) ? i.first : prev(i.first);
-    auto high = i.second;
-    
-    x_a = low->first;
-    x_b = high->first;
-    y_a = low->second;
-    y_b = high->second;
-    
-    return y_a + (x - x_a)*(y_b - y_a)/(x_b - x_a);
+    return result;
 }
 
 // DoubleFunction factory //////////////////////////////////////////////////////
@@ -101,13 +129,14 @@ DoubleFunction TabFunction::makeFunction(const bool makeHardCopy) const
     return res;
 }
 
-DoubleFunction Latan::interpolate(const DVec &x, const DVec &y)
+DoubleFunction Latan::interpolate(const DVec &x, const DVec &y,
+                                  const InterpType interpType)
 {
-    return TabFunction(x,y).makeFunction();
+    return TabFunction(x, y, interpType).makeFunction();
 }
 
 DoubleFunction Latan::interpolate(const XYStatData &data, const Index i,
-                                  const Index j)
+                                  const Index j, const InterpType interpType)
 {
-    return TabFunction(data, i, j).makeFunction();
+    return TabFunction(data, i, j, interpType).makeFunction();
 }

lib/TabFunction.hpp

@@ -20,9 +20,11 @@
 #ifndef Latan_TabFunction_hpp_
 #define Latan_TabFunction_hpp_
 
+#include <algorithm>
 #include <map>
 #include <LatAnalyze/Global.hpp>
 #include <LatAnalyze/Function.hpp>
+#include <LatAnalyze/Math.hpp>
 #include <LatAnalyze/XYStatData.hpp>
 
 BEGIN_LATAN_NAMESPACE
@@ -31,13 +33,22 @@ BEGIN_LATAN_NAMESPACE
  *                      tabulated function: 1D only                           *
  ******************************************************************************/
 
+enum class InterpType
+{
+  NEAREST,
+  LINEAR,
+  QUADRATIC
+};
+
 class TabFunction: public DoubleFunctionFactory
 {
 public:
     // constructors
     TabFunction(void) = default;
-    TabFunction(const DVec &x, const DVec &y);
-    TabFunction(const XYStatData &data, const Index i = 0, const Index j = 0);
+    TabFunction(const DVec &x, const DVec &y,
+                const InterpType interpType = InterpType::LINEAR);
+    TabFunction(const XYStatData &data, const Index i = 0, const Index j = 0,
+                const InterpType interpType = InterpType::LINEAR);
     // destructor
     virtual ~TabFunction(void) = default;
     // access
@@ -49,11 +60,14 @@ public:
     virtual DoubleFunction makeFunction(const bool makeHardCopy = true) const;
 private:
     std::map<double, double> value_;
+    InterpType interpType_;
 };
 
-DoubleFunction interpolate(const DVec &x, const DVec &y);
+DoubleFunction interpolate(const DVec &x, const DVec &y,
+                           const InterpType interpType = InterpType::LINEAR);
 DoubleFunction interpolate(const XYStatData &data, const Index i = 0,
-                           const Index j = 0);
+                           const Index j = 0,
+                           const InterpType interpType = InterpType::LINEAR);
 
 END_LATAN_NAMESPACE