From 118da972b47ada15301f18b1957014fbf73d11f0 Mon Sep 17 00:00:00 2001
From: Christoph Ruegg <git@cdrnet.ch>
Date: Sun, 20 Jul 2014 13:27:51 +0200
Subject: [PATCH] Interpolation: modify step interpolation to allow non-zero
 right-open segment value.

---
 .../Interpolation/StepInterpolation.cs        | 31 +++++--------------
 .../StepInterpolationTest.cs                  |  4 ++-
 2 files changed, 11 insertions(+), 24 deletions(-)

diff --git a/src/Numerics/Interpolation/StepInterpolation.cs b/src/Numerics/Interpolation/StepInterpolation.cs
index 7f90d7e4..b851eb8f 100644
--- a/src/Numerics/Interpolation/StepInterpolation.cs
+++ b/src/Numerics/Interpolation/StepInterpolation.cs
@@ -36,9 +36,9 @@ using MathNet.Numerics.Properties;
 namespace MathNet.Numerics.Interpolation
 {
     /// <summary>
-    /// A step function where the start of each segment is included, and last is excluded.
-    /// Segment i is [x_i, x_i+1).
-    /// The domain of the function is all real numbers, such that y = 0 where x &lt; x_0 or x gt; x_n
+    /// A step function where the start of each segment is included, and the last segment is open-ended.
+    /// Segment i is [x_i, x_i+1) for i &lt; N, or [x_i, infinity] for i = N.
+    /// The domain of the function is all real numbers, such that y = 0 where x &lt;.
     /// </summary>
     /// <remarks>Supports both differentiation and integration.</remarks>
     public class StepInterpolation : IInterpolation
@@ -47,11 +47,11 @@ namespace MathNet.Numerics.Interpolation
         readonly double[] _y;
         readonly Lazy<double[]> _indefiniteIntegral;
 
-        /// <param name="x">Sample points (N+1), sorted ascending</param>
-        /// <param name="sy">Functional values (N) of each segment</param>
+        /// <param name="x">Sample points (N), sorted ascending</param>
+        /// <param name="sy">Samples values (N) of each segment starting at the corresponding sample point.</param>
         public StepInterpolation(double[] x, double[] sy)
         {
-            if (x.Length != sy.Length + 1)
+            if (x.Length != sy.Length)
             {
                 throw new ArgumentException(Resources.ArgumentVectorsSameLength);
             }
@@ -63,25 +63,14 @@ namespace MathNet.Numerics.Interpolation
 
         /// <summary>
         /// Create a linear spline interpolation from a set of (x,y) value pairs, sorted ascendingly by x.
-        /// The y-value corresponding to the largest x-sample is ignored.
         /// </summary>
         public static StepInterpolation InterpolateSorted(double[] x, double[] y)
         {
-            if (x.Length != y.Length)
-            {
-                throw new ArgumentException(Resources.ArgumentVectorsSameLength);
-            }
-
-            // drop the last value which is not part of any segment.
-            var segmentValues = new double[x.Length - 1];
-            Array.Copy(y, 0, segmentValues, 0, segmentValues.Length);
-
-            return new StepInterpolation(x, segmentValues);
+            return new StepInterpolation(x, y);
         }
 
         /// <summary>
         /// Create a linear spline interpolation from an unsorted set of (x,y) value pairs.
-        /// The y-value corresponding to the largest x-sample is ignored.
         /// WARNING: Works in-place and can thus causes the data array to be reordered.
         /// </summary>
         public static StepInterpolation InterpolateInplace(double[] x, double[] y)
@@ -97,7 +86,6 @@ namespace MathNet.Numerics.Interpolation
 
         /// <summary>
         /// Create a linear spline interpolation from an unsorted set of (x,y) value pairs.
-        /// The y-value corresponding to the largest x-sample is ignored.
         /// </summary>
         public static StepInterpolation Interpolate(IEnumerable<double> x, IEnumerable<double> y)
         {
@@ -122,7 +110,7 @@ namespace MathNet.Numerics.Interpolation
         /// <returns>Interpolated value x(t).</returns>
         public double Interpolate(double t)
         {
-            if (t < _x[0] || t >= _x[_x.Length - 1])
+            if (t < _x[0])
                 return 0.0;
 
             int k = LeftBracketIndex(t);
@@ -160,9 +148,6 @@ namespace MathNet.Numerics.Interpolation
         {
             if (t <= _x[0])
                 return 0.0;
-            int last = _x.Length - 1;
-            if (t >= _x[last])
-                return _indefiniteIntegral.Value[last];
 
             int k = LeftBracketIndex(t);
             var x = (t - _x[k]);
diff --git a/src/UnitTests/InterpolationTests/StepInterpolationTest.cs b/src/UnitTests/InterpolationTests/StepInterpolationTest.cs
index d84ff19f..15a4c930 100644
--- a/src/UnitTests/InterpolationTests/StepInterpolationTest.cs
+++ b/src/UnitTests/InterpolationTests/StepInterpolationTest.cs
@@ -37,7 +37,7 @@ namespace MathNet.Numerics.UnitTests.InterpolationTests
     public class StepInterpolationTest
     {
         readonly double[] _t = { -2.0, -1.0, 0.0, 1.0, 2.0, 3.0 };
-        readonly double[] _y = { 1.0, 2.0, -1.0, 0.0, 1.0 };
+        readonly double[] _y = { 1.0, 2.0, -1.0, 0.0, 1.0, 0.0 };
 
         [Test]
         public void FirstDerivative()
@@ -53,6 +53,7 @@ namespace MathNet.Numerics.UnitTests.InterpolationTests
             Assert.That(ip.Differentiate(1.0), Is.EqualTo(double.NaN));
             Assert.That(ip.Differentiate(2.0), Is.EqualTo(double.NaN));
             Assert.That(ip.Differentiate(3.0), Is.EqualTo(double.NaN));
+            Assert.That(ip.Differentiate(4.0), Is.EqualTo(0.0));
         }
 
         [Test]
@@ -71,6 +72,7 @@ namespace MathNet.Numerics.UnitTests.InterpolationTests
             Assert.That(ip.Integrate(-3.0, 4.0), Is.EqualTo(3.0));
             Assert.That(ip.Integrate(0.5, 1.5), Is.EqualTo(-0.5));
             Assert.That(ip.Integrate(-1.5, -0.5), Is.EqualTo(1.5));
+            Assert.That(ip.Integrate(3.0, 4.0), Is.EqualTo(0.0));
         }
 
         /// <summary>