Interpolation: explicitly check for min required number of samples #232

src/Numerics/Interpolation/Barycentric.cs

@@ -57,7 +57,7 @@ namespace MathNet.Numerics.Interpolation
 
             if (x.Length < 1)
             {
-                throw new ArgumentOutOfRangeException("x");
+                throw new ArgumentException(string.Format(Resources.ArrayTooSmall, 1), "x");
             }
 
             _x = x;
@@ -77,7 +77,7 @@ namespace MathNet.Numerics.Interpolation
 
             if (x.Length < 1)
             {
-                throw new ArgumentOutOfRangeException("x");
+                throw new ArgumentException(string.Format(Resources.ArrayTooSmall, 1), "x");
             }
 
             var weights = new double[x.Length];
@@ -101,11 +101,6 @@ namespace MathNet.Numerics.Interpolation
                 throw new ArgumentException(Resources.ArgumentVectorsSameLength);
             }
 
-            if (x.Length < 1)
-            {
-                throw new ArgumentOutOfRangeException("x");
-            }
-
             Sorting.Sort(x, y);
             return InterpolatePolynomialEquidistantSorted(x, y);
         }
@@ -141,28 +136,28 @@ namespace MathNet.Numerics.Interpolation
         /// </param>
         public static Barycentric InterpolateRationalFloaterHormannSorted(double[] x, double[] y, int order)
         {
-            var xx = (x as double[]) ?? x.ToArray();
+            if (x.Length != y.Length)
-            var yy = (y as double[]) ?? y.ToArray();
-
-            if (xx.Length != yy.Length)
             {
                 throw new ArgumentException(Resources.ArgumentVectorsSameLength);
             }
 
-            if (0 > order || xx.Length <= order)
+            if (x.Length < 1)
             {
-                throw new ArgumentOutOfRangeException("order");
+                throw new ArgumentException(string.Format(Resources.ArrayTooSmall, 1), "x");
             }
 
-            Sorting.Sort(xx, yy);
+            if (0 > order || x.Length <= order)
+            {
+                throw new ArgumentOutOfRangeException("order");
+            }
 
-            var weights = new double[xx.Length];
+            var weights = new double[x.Length];
 
             // order: odd -> negative, even -> positive
             double sign = ((order & 0x1) == 0x1) ? -1.0 : 1.0;
 
             // compute barycentric weights
-            for (int k = 0; k < xx.Length; k++)
+            for (int k = 0; k < x.Length; k++)
             {
                 double s = 0;
                 for (int i = Math.Max(k - order, 0); i <= Math.Min(k, weights.Length - 1 - order); i++)
@@ -172,7 +167,7 @@ namespace MathNet.Numerics.Interpolation
                     {
                         if (j != k)
                         {
-                            v = v/Math.Abs(xx[k] - xx[j]);
+                            v = v/Math.Abs(x[k] - x[j]);
                         }
                     }
 
@@ -183,7 +178,7 @@ namespace MathNet.Numerics.Interpolation
                 sign = -sign;
             }
 
-            return new Barycentric(xx, yy, weights);
+            return new Barycentric(x, y, weights);
         }
 
         /// <summary>
@@ -203,11 +198,6 @@ namespace MathNet.Numerics.Interpolation
                 throw new ArgumentException(Resources.ArgumentVectorsSameLength);
             }
 
-            if (0 > order || x.Length <= order)
-            {
-                throw new ArgumentOutOfRangeException("order");
-            }
-
             Sorting.Sort(x, y);
             return InterpolateRationalFloaterHormannSorted(x, y, order);
         }

src/Numerics/Interpolation/BulirschStoerRationalInterpolation.cs

@@ -59,7 +59,7 @@ namespace MathNet.Numerics.Interpolation
 
             if (x.Length < 1)
             {
-                throw new ArgumentOutOfRangeException("x");
+                throw new ArgumentException(string.Format(Resources.ArrayTooSmall, 1), "x");
             }
 
             _x = x;

src/Numerics/Interpolation/CubicSpline.cs

@@ -60,6 +60,11 @@ namespace MathNet.Numerics.Interpolation
                 throw new ArgumentException(Resources.ArgumentVectorsSameLength);
             }
 
+            if (x.Length < 2)
+            {
+                throw new ArgumentException(string.Format(Resources.ArrayTooSmall, 2), "x");
+            }
+
             _x = x;
             _c0 = c0;
             _c1 = c1;
@@ -80,7 +85,7 @@ namespace MathNet.Numerics.Interpolation
 
             if (x.Length < 2)
             {
-                throw new ArgumentOutOfRangeException("x");
+                throw new ArgumentException(string.Format(Resources.ArrayTooSmall, 2), "x");
             }
 
             var c0 = new double[x.Length - 1];
@@ -113,7 +118,7 @@ namespace MathNet.Numerics.Interpolation
 
             if (x.Length < 2)
             {
-                throw new ArgumentOutOfRangeException("x");
+                throw new ArgumentException(string.Format(Resources.ArrayTooSmall, 2), "x");
             }
 
             Sorting.Sort(x, y, firstDerivatives);
@@ -142,7 +147,7 @@ namespace MathNet.Numerics.Interpolation
 
             if (x.Length < 5)
             {
-                throw new ArgumentOutOfRangeException("x");
+                throw new ArgumentException(string.Format(Resources.ArrayTooSmall, 5), "x");
             }
 
             /* Prepare divided differences (diff) and weights (w) */
@@ -193,11 +198,6 @@ namespace MathNet.Numerics.Interpolation
                 throw new ArgumentException(Resources.ArgumentVectorsSameLength);
             }
 
-            if (x.Length < 5)
-            {
-                throw new ArgumentOutOfRangeException("x");
-            }
-
             Sorting.Sort(x, y);
             return InterpolateAkimaSorted(x, y);
         }
@@ -227,7 +227,7 @@ namespace MathNet.Numerics.Interpolation
 
             if (x.Length < 2)
             {
-                throw new ArgumentOutOfRangeException("x");
+                throw new ArgumentException(string.Format(Resources.ArrayTooSmall, 2), "x");
             }
 
             int n = x.Length;
@@ -337,11 +337,6 @@ namespace MathNet.Numerics.Interpolation
                 throw new ArgumentException(Resources.ArgumentVectorsSameLength);
             }
 
-            if (x.Length < 2)
-            {
-                throw new ArgumentOutOfRangeException("x");
-            }
-
             Sorting.Sort(x, y);
             return InterpolateBoundariesSorted(x, y, leftBoundaryCondition, leftBoundary, rightBoundaryCondition, rightBoundary);
         }
@@ -460,7 +455,7 @@ namespace MathNet.Numerics.Interpolation
         /// <returns>Interpolated value x(t).</returns>
         public double Interpolate(double t)
         {
-            int k = LeftBracketIndex(t);
+            int k = LeftSegmentIndex(t);
             var x = t - _x[k];
             return _c0[k] + x*(_c1[k] + x*(_c2[k] + x*_c3[k]));
         }
@@ -472,7 +467,7 @@ namespace MathNet.Numerics.Interpolation
         /// <returns>Interpolated first derivative at point t.</returns>
         public double Differentiate(double t)
         {
-            int k = LeftBracketIndex(t);
+            int k = LeftSegmentIndex(t);
             var x = t - _x[k];
             return _c1[k] + x*(2*_c2[k] + x*3*_c3[k]);
         }
@@ -484,7 +479,7 @@ namespace MathNet.Numerics.Interpolation
         /// <returns>Interpolated second derivative at point t.</returns>
         public double Differentiate2(double t)
         {
-            int k = LeftBracketIndex(t);
+            int k = LeftSegmentIndex(t);
             var x = t - _x[k];
             return 2*_c2[k] + x*6*_c3[k];
         }
@@ -495,7 +490,7 @@ namespace MathNet.Numerics.Interpolation
         /// <param name="t">Point t to integrate at.</param>
         public double Integrate(double t)
         {
-            int k = LeftBracketIndex(t);
+            int k = LeftSegmentIndex(t);
             var x = t - _x[k];
             return _indefiniteIntegral.Value[k] + x*(_c0[k] + x*(_c1[k]/2 + x*(_c2[k]/3 + x*_c3[k]/4)));
         }
@@ -523,11 +518,11 @@ namespace MathNet.Numerics.Interpolation
         }
 
         /// <summary>
-        /// Find the index of the greatest sample point smaller than t.
+        /// Find the index of the greatest sample point smaller than t,
+        /// or the left index of the closest segment for extrapolation.
         /// </summary>
-        int LeftBracketIndex(double t)
+        int LeftSegmentIndex(double t)
         {
-            // Binary search in the [ t[0], ..., t[n-2] ] (t[n-1] is not included)
             int low = 0;
             int high = _x.Length - 1;
             while (low != high - 1)

src/Numerics/Interpolation/LinearSpline.cs

@@ -56,6 +56,11 @@ namespace MathNet.Numerics.Interpolation
                 throw new ArgumentException(Resources.ArgumentVectorsSameLength);
             }
 
+            if (x.Length < 2)
+            {
+                throw new ArgumentException(string.Format(Resources.ArrayTooSmall, 2), "x");
+            }
+
             _x = x;
             _c0 = c0;
             _c1 = c1;
@@ -72,6 +77,11 @@ namespace MathNet.Numerics.Interpolation
                 throw new ArgumentException(Resources.ArgumentVectorsSameLength);
             }
 
+            if (x.Length < 2)
+            {
+                throw new ArgumentException(string.Format(Resources.ArrayTooSmall, 2), "x");
+            }
+
             var c1 = new double[x.Length - 1];
             for (int i = 0; i < c1.Length; i++)
             {
@@ -128,7 +138,7 @@ namespace MathNet.Numerics.Interpolation
         /// <returns>Interpolated value x(t).</returns>
         public double Interpolate(double t)
         {
-            int k = LeftBracketIndex(t);
+            int k = LeftSegmentIndex(t);
             return _c0[k] + (t - _x[k])*_c1[k];
         }
 
@@ -139,7 +149,7 @@ namespace MathNet.Numerics.Interpolation
         /// <returns>Interpolated first derivative at point t.</returns>
         public double Differentiate(double t)
         {
-            int k = LeftBracketIndex(t);
+            int k = LeftSegmentIndex(t);
             return _c1[k];
         }
 
@@ -159,7 +169,7 @@ namespace MathNet.Numerics.Interpolation
         /// <param name="t">Point t to integrate at.</param>
         public double Integrate(double t)
         {
-            int k = LeftBracketIndex(t);
+            int k = LeftSegmentIndex(t);
             var x = t - _x[k];
             return _indefiniteIntegral.Value[k] + x*(_c0[k] + x*_c1[k]/2);
         }
@@ -187,11 +197,11 @@ namespace MathNet.Numerics.Interpolation
         }
 
         /// <summary>
-        /// Find the index of the greatest sample point smaller than t.
+        /// Find the index of the greatest sample point smaller than t,
+        /// or the left index of the closest segment for extrapolation.
         /// </summary>
-        int LeftBracketIndex(double t)
+        int LeftSegmentIndex(double t)
         {
-            // Binary search in the [ t[0], ..., t[n-2] ] (t[n-1] is not included)
             int low = 0;
             int high = _x.Length - 1;
             while (low != high - 1)

src/Numerics/Interpolation/NevillePolynomialInterpolation.cs

@@ -64,7 +64,7 @@ namespace MathNet.Numerics.Interpolation
 
             if (x.Length < 1)
             {
-                throw new ArgumentOutOfRangeException("x");
+                throw new ArgumentException(string.Format(Resources.ArrayTooSmall, 1), "x");
             }
 
             for (var i = 1; i < x.Length; ++i)

src/Numerics/Interpolation/QuadraticSpline.cs

@@ -56,6 +56,11 @@ namespace MathNet.Numerics.Interpolation
                 throw new ArgumentException(Resources.ArgumentVectorsSameLength);
             }
 
+            if (x.Length < 2)
+            {
+                throw new ArgumentException(string.Format(Resources.ArrayTooSmall, 2), "x");
+            }
+
             _x = x;
             _c0 = c0;
             _c1 = c1;
@@ -86,7 +91,7 @@ namespace MathNet.Numerics.Interpolation
         /// <returns>Interpolated value x(t).</returns>
         public double Interpolate(double t)
         {
-            int k = LeftBracketIndex(t);
+            int k = LeftSegmentIndex(t);
             var x = t - _x[k];
             return _c0[k] + x*(_c1[k] + x*_c2[k]);
         }
@@ -98,7 +103,7 @@ namespace MathNet.Numerics.Interpolation
         /// <returns>Interpolated first derivative at point t.</returns>
         public double Differentiate(double t)
         {
-            int k = LeftBracketIndex(t);
+            int k = LeftSegmentIndex(t);
             return _c1[k] + (t - _x[k])*2*_c2[k];
         }
 
@@ -109,7 +114,7 @@ namespace MathNet.Numerics.Interpolation
         /// <returns>Interpolated second derivative at point t.</returns>
         public double Differentiate2(double t)
         {
-            int k = LeftBracketIndex(t);
+            int k = LeftSegmentIndex(t);
             return 2*_c2[k];
         }
 
@@ -119,7 +124,7 @@ namespace MathNet.Numerics.Interpolation
         /// <param name="t">Point t to integrate at.</param>
         public double Integrate(double t)
         {
-            int k = LeftBracketIndex(t);
+            int k = LeftSegmentIndex(t);
             var x = t - _x[k];
             return _indefiniteIntegral.Value[k] + x*(_c0[k] + x*(_c1[k]/2 + x*_c2[k]/3));
         }
@@ -147,11 +152,11 @@ namespace MathNet.Numerics.Interpolation
         }
 
         /// <summary>
-        /// Find the index of the greatest sample point smaller than t.
+        /// Find the index of the greatest sample point smaller than t,
+        /// or the left index of the closest segment for extrapolation.
         /// </summary>
-        int LeftBracketIndex(double t)
+        int LeftSegmentIndex(double t)
         {
-            // Binary search in the [ t[0], ..., t[n-2] ] (t[n-1] is not included)
             int low = 0;
             int high = _x.Length - 1;
             while (low != high - 1)

src/Numerics/Interpolation/StepInterpolation.cs

@@ -56,6 +56,11 @@ namespace MathNet.Numerics.Interpolation
                 throw new ArgumentException(Resources.ArgumentVectorsSameLength);
             }
 
+            if (x.Length < 1)
+            {
+                throw new ArgumentException(string.Format(Resources.ArrayTooSmall, 1), "x");
+            }
+
             _x = x;
             _y = sy;
             _indefiniteIntegral = new Lazy<double[]>(ComputeIndefiniteIntegral);

src/UnitTests/InterpolationTests/AkimaSplineTest.cs

@@ -92,5 +92,13 @@ namespace MathNet.Numerics.UnitTests.InterpolationTests
                 Assert.AreEqual(ytest[i], it.Interpolate(xtest[i]), 1e-15, "Linear with {0} samples, sample {1}", samples, i);
             }
         }
+
+        [Test]
+        public void FewSamples()
+        {
+            Assert.That(() => CubicSpline.InterpolateAkima(new double[0], new double[0]), Throws.ArgumentException);
+            Assert.That(() => CubicSpline.InterpolateAkima(new double[4], new double[4]), Throws.ArgumentException);
+            Assert.That(CubicSpline.InterpolateAkima(new[] { 1.0, 2.0, 3.0, 4.0, 5.0 }, new[] { 2.0, 2.0, 2.0, 2.0, 2.0 }).Interpolate(1.0), Is.EqualTo(2.0));
+        }
     }
 }

src/UnitTests/InterpolationTests/BulirschStoerRationalTest.cs

@@ -92,6 +92,13 @@ namespace MathNet.Numerics.UnitTests.InterpolationTests
             Assert.AreEqual(x, interpolation.Interpolate(t), maxAbsoluteError, "Interpolation at {0}", t);
         }
 
+        [Test]
+        public void FewSamples()
+        {
+            Assert.That(() => BulirschStoerRationalInterpolation.Interpolate(new double[0], new double[0]), Throws.ArgumentException);
+            Assert.That(BulirschStoerRationalInterpolation.Interpolate(new[] { 1.0 }, new[] { 2.0 }).Interpolate(1.0), Is.EqualTo(2.0));
+        }
+
         // NOTE: No test for the linear case because this algorithms is incredibly bad at this.
     }
 }

src/UnitTests/InterpolationTests/CubicSplineTest.cs

@@ -176,6 +176,14 @@ namespace MathNet.Numerics.UnitTests.InterpolationTests
             }
         }
 
+        [Test]
+        public void FewSamples()
+        {
+            Assert.That(() => CubicSpline.InterpolateNatural(new double[0], new double[0]), Throws.ArgumentException);
+            Assert.That(() => CubicSpline.InterpolateNatural(new double[1], new double[1]), Throws.ArgumentException);
+            Assert.That(CubicSpline.InterpolateNatural(new[] { 1.0, 2.0 }, new[] { 2.0, 2.0 }).Interpolate(1.0), Is.EqualTo(2.0));
+        }
+
 #if !NET35 && !PORTABLE
         [Test]
         public void InterpolateAkimaSorted_MustBeThreadSafe_GitHub219([Values(8, 32, 256, 1024)] int samples)

src/UnitTests/InterpolationTests/EquidistantPolynomialTest.cs

@@ -94,5 +94,12 @@ namespace MathNet.Numerics.UnitTests.InterpolationTests
                 Assert.AreEqual(ytest[i], it.Interpolate(xtest[i]), 1e-12, "Linear with {0} samples, sample {1}", samples, i);
             }
         }
+
+        [Test]
+        public void FewSamples()
+        {
+            Assert.That(() => Barycentric.InterpolatePolynomialEquidistant(new double[0], new double[0]), Throws.ArgumentException);
+            Assert.That(Barycentric.InterpolatePolynomialEquidistant(new[] { 1.0 }, new[] { 2.0 }).Interpolate(1.0), Is.EqualTo(2.0));
+        }
     }
 }

src/UnitTests/InterpolationTests/FloaterHormannRationalTest.cs

@@ -119,5 +119,12 @@ namespace MathNet.Numerics.UnitTests.InterpolationTests
                 Assert.AreEqual(ytest[i], it.Interpolate(xtest[i]), 1e-14, "Linear with {0} samples, sample {1}", samples, i);
             }
         }
+
+        [Test]
+        public void FewSamples()
+        {
+            Assert.That(() => Barycentric.InterpolateRationalFloaterHormann(new double[0], new double[0]), Throws.ArgumentException);
+            Assert.That(Barycentric.InterpolateRationalFloaterHormann(new[] { 1.0 }, new[] { 2.0 }).Interpolate(1.0), Is.EqualTo(2.0));
+        }
     }
 }

src/UnitTests/InterpolationTests/LinearSplineTest.cs

@@ -131,5 +131,13 @@ namespace MathNet.Numerics.UnitTests.InterpolationTests
                 Assert.AreEqual(ytest[i], ip.Interpolate(xtest[i]), 1e-15, "Linear with {0} samples, sample {1}", samples, i);
             }
         }
+
+        [Test]
+        public void FewSamples()
+        {
+            Assert.That(() => LinearSpline.Interpolate(new double[0], new double[0]), Throws.ArgumentException);
+            Assert.That(() => LinearSpline.Interpolate(new double[1], new double[1]), Throws.ArgumentException);
+            Assert.That(LinearSpline.Interpolate(new[] { 1.0, 2.0 }, new[] { 2.0, 2.0 }).Interpolate(1.0), Is.EqualTo(2.0));
+        }
     }
 }

src/UnitTests/InterpolationTests/NevillePolynomialTest.cs

@@ -146,5 +146,12 @@ namespace MathNet.Numerics.UnitTests.InterpolationTests
             var actual = interpolation.Interpolate(Math.Log(value));
             Assert.That(actual, Is.Not.NaN);
         }
+
+        [Test]
+        public void FewSamples()
+        {
+            Assert.That(() => NevillePolynomialInterpolation.Interpolate(new double[0], new double[0]), Throws.ArgumentException);
+            Assert.That(NevillePolynomialInterpolation.Interpolate(new[] { 1.0 }, new[] { 2.0 }).Interpolate(1.0), Is.EqualTo(2.0));
+        }
     }
 }

src/UnitTests/InterpolationTests/StepInterpolationTest.cs

@@ -28,6 +28,7 @@
 // OTHER DEALINGS IN THE SOFTWARE.
 // </copyright>
 
+using System;
 using MathNet.Numerics.Interpolation;
 using NUnit.Framework;
 
@@ -87,5 +88,12 @@ namespace MathNet.Numerics.UnitTests.InterpolationTests
                 Assert.AreEqual(_y[i], ip.Interpolate(_t[i]), "A Exact Point " + i);
             }
         }
+
+        [Test]
+        public void FewSamples()
+        {
+            Assert.That(() => StepInterpolation.Interpolate(new double[0], new double[0]), Throws.ArgumentException);
+            Assert.That(StepInterpolation.Interpolate(new[] { 1.0 }, new[] { 2.0 }).Interpolate(1.0), Is.EqualTo(2.0));
+        }
     }
 }