From ce9f979652762e124c5c942ca0b28cbdcbe4520b Mon Sep 17 00:00:00 2001
From: Michael Wedel <michael.wedel@psi.ch>
Date: Mon, 23 Jun 2014 16:53:36 +0200
Subject: [PATCH] Re #9682. Implemented binary search

Replaced the linear search with a binary search and added unit tests.
---
 .../Kernel/inc/MantidKernel/Interpolation.h   |  4 ++
 .../Framework/Kernel/src/Interpolation.cpp    | 38 ++++++++++++++-----
 .../Framework/Kernel/test/InterpolationTest.h | 38 +++++++++++++++++++
 3 files changed, 70 insertions(+), 10 deletions(-)

diff --git a/Code/Mantid/Framework/Kernel/inc/MantidKernel/Interpolation.h b/Code/Mantid/Framework/Kernel/inc/MantidKernel/Interpolation.h
index 7cf862d611b5..417c166cf542 100644
--- a/Code/Mantid/Framework/Kernel/inc/MantidKernel/Interpolation.h
+++ b/Code/Mantid/Framework/Kernel/inc/MantidKernel/Interpolation.h
@@ -62,10 +62,14 @@ class MANTID_KERNEL_DLL Interpolation
   /// unit of y-axis
   Unit_sptr m_yUnit;
 
+protected:
+  size_t findIndexOfNextLargerValue(const std::vector<double> &data, double key, size_t range_start, size_t range_end) const;
+
 public:
 
   /// Constructor default to linear interpolation and x-unit set to TOF
   Interpolation();
+  virtual ~Interpolation() { }
 
   /// add data point
   void addPoint(const double& xx, const double& yy);
diff --git a/Code/Mantid/Framework/Kernel/src/Interpolation.cpp b/Code/Mantid/Framework/Kernel/src/Interpolation.cpp
index 6d5bb1756325..346b13a01df1 100644
--- a/Code/Mantid/Framework/Kernel/src/Interpolation.cpp
+++ b/Code/Mantid/Framework/Kernel/src/Interpolation.cpp
@@ -21,6 +21,26 @@ namespace Kernel
     m_yUnit = UnitFactory::Instance().create("TOF");
   }
 
+  size_t Interpolation::findIndexOfNextLargerValue(const std::vector<double> &data, double key, size_t range_start, size_t range_end) const
+  {
+      if(range_end < range_start)
+      {
+          throw std::range_error("Value is outside array range.");
+      }
+
+      size_t center = range_start + (range_end - range_start) / 2;
+
+      if(data[center] > key && data[center - 1] <= key) {
+          return center;
+      }
+
+      if(data[center] <= key) {
+          return findIndexOfNextLargerValue(data, key, center + 1, range_end);
+      } else {
+          return findIndexOfNextLargerValue(data, key, range_start, center - 1);
+      }
+  }
+
   void Interpolation::setXUnit(const std::string& unit) 
   { 
     m_xUnit = UnitFactory::Instance().create(unit);
@@ -48,7 +68,7 @@ namespace Kernel
     
     // check first if at is within the limits of interpolation interval
 
-    if ( at <= m_x[0] )
+    if ( at < m_x[0] )
     {
       return m_y[0]-(m_x[0]-at)*(m_y[1]-m_y[0])/(m_x[1]-m_x[0]);
     }
@@ -58,16 +78,14 @@ namespace Kernel
       return m_y[N-1]+(at-m_x[N-1])*(m_y[N-1]-m_y[N-2])/(m_x[N-1]-m_x[N-2]);
     }
 
-    // otherwise
-
-    for (unsigned int i = 1; i < N; i++)
-    {
-      if ( m_x[i] > at )
-      {
-        return m_y[i-1] + (at-m_x[i-1])*(m_y[i]-m_y[i-1])/(m_x[i]-m_x[i-1]);
-      }
+    try {
+        // otherwise
+        // General case. Find index of next largest value by binary search.
+        size_t idx = findIndexOfNextLargerValue(m_x, at, 1, N - 1);
+        return m_y[idx-1] + (at - m_x[idx-1])*(m_y[idx]-m_y[idx-1])/(m_x[idx]-m_x[idx-1]);
+    } catch(std::range_error) {
+        return 0.0;
     }
-    return 0.0;
   }
 
   /** Add point in the interpolation.
diff --git a/Code/Mantid/Framework/Kernel/test/InterpolationTest.h b/Code/Mantid/Framework/Kernel/test/InterpolationTest.h
index 45909fa77a96..636b0fc468e9 100644
--- a/Code/Mantid/Framework/Kernel/test/InterpolationTest.h
+++ b/Code/Mantid/Framework/Kernel/test/InterpolationTest.h
@@ -163,6 +163,34 @@ class InterpolationTest : public CxxTest::TestSuite
         checkInterpolationResults(interpolation);
     }
 
+    void testFindIndexOfNextLargerValue()
+    {
+        TestableInterpolation interpolation;
+
+        size_t N = m_tableXValues.size();
+
+        // lower limit - can be treated like general case
+        TS_ASSERT_EQUALS(interpolation.findIndexOfNextLargerValue(m_tableXValues, 200.0, 1, N - 1), 1);
+
+        // Exact interpolation points
+        TS_ASSERT_EQUALS(interpolation.findIndexOfNextLargerValue(m_tableXValues, 201.0, 1, N - 1), 2);
+        TS_ASSERT_EQUALS(interpolation.findIndexOfNextLargerValue(m_tableXValues, 202.0, 1, N - 1), 3);
+        TS_ASSERT_EQUALS(interpolation.findIndexOfNextLargerValue(m_tableXValues, 203.0, 1, N - 1), 4);
+
+        // Arbitrary interpolation points
+        TS_ASSERT_EQUALS(interpolation.findIndexOfNextLargerValue(m_tableXValues, 200.5, 1, N - 1), 1);
+        TS_ASSERT_EQUALS(interpolation.findIndexOfNextLargerValue(m_tableXValues, 201.25, 1, N - 1), 2);
+        TS_ASSERT_EQUALS(interpolation.findIndexOfNextLargerValue(m_tableXValues, 203.5, 1, N - 1), 4);
+
+
+        // upper limit - must be covered as edge case before this can ever be called.
+        TS_ASSERT_THROWS(interpolation.findIndexOfNextLargerValue(m_tableXValues, 204.0, 1, N - 1), std::range_error);
+
+        // outside interpolation limits - edge cases as well
+        TS_ASSERT_THROWS(interpolation.findIndexOfNextLargerValue(m_tableXValues, 199, 1, N - 1), std::range_error)
+        TS_ASSERT_THROWS(interpolation.findIndexOfNextLargerValue(m_tableXValues, 2000.0, 1, N - 1), std::range_error)
+    }
+
 private:
     Interpolation getInitializedInterpolation(std::string xUnit, std::string yUnit)
     {
@@ -240,6 +268,16 @@ class InterpolationTest : public CxxTest::TestSuite
     // Values outside interpolation range
     std::vector<double> m_outsideXValues;
     std::vector<double> m_outsideYValues;
+
+    // For the test of findIndexOfNextLargerValue access to protected member is needed
+    class TestableInterpolation : public Interpolation {
+        friend class InterpolationTest;
+
+    public:
+        TestableInterpolation() : Interpolation()
+        { }
+        ~TestableInterpolation() { }
+    };
 };
 
 #endif /*INTERPOLATIONTEST_H_*/