Updated interpolation classes to prevent non-unique and/or non-monotonically ascending sample points. Fixes issue #1.

data/Github-Cureos-1.csv

@@ -0,0 +1,49 @@
+0.001,3719,-6.907755279,8.221210094
+0.00103542,3393,-6.872948137,8.129469765
+0.0010721,3094,-6.838135937,8.037220031
+0.0010721,3100,-6.838135937,8.03915739
+0.0015,1251,-6.502290171,7.13169851
+0.002,559.4,-6.214608098,6.326864781
+0.0021455,458.1,-6.144382654,6.127087501
+0.0021455,462.6,-6.144382654,6.13686275
+0.00230297,377.6,-6.073555685,5.933835434
+0.002472,308.5,-6.002727739,5.731721843
+0.002472,314,-6.002727739,5.749392986
+0.0026414,259.7,-5.936446199,5.559527119
+0.0028224,214.5,-5.870167692,5.368309738
+0.0028224,216,-5.870167692,5.375278408
+0.003,181.2,-5.80914299,5.199601394
+0.0036074,105.7,-5.624767988,4.660604893
+0.0036074,110,-5.624767988,4.700480366
+0.004,81.27,-5.521460918,4.397776945
+0.005,42.06,-5.298317367,3.73909717
+0.006,24.46,-5.11599581,3.19703913
+0.008,10.37,-4.828313737,2.338917022
+0.01,5.356,-4.605170186,1.678217428
+0.015,1.693,-4.199705078,0.526502103
+0.02,0.8205,-3.912023005,-0.197841368
+0.03,0.3783,-3.506557897,-0.972067747
+0.04,0.2685,-3.218875825,-1.314904365
+0.05,0.2262,-2.995732274,-1.486335715
+0.06,0.2048,-2.813410717,-1.585721386
+0.08,0.1823,-2.525728644,-1.702101597
+0.1,0.1693,-2.302585093,-1.77608299
+0.15,0.1492,-1.897119985,-1.902467591
+0.2,0.1358,-1.609437912,-1.996572064
+0.3,0.1176,-1.203972804,-2.140466244
+0.4,0.1052,-0.916290732,-2.251891979
+0.5,0.09598,-0.693147181,-2.343615443
+0.6,0.08874,-0.510825624,-2.422044533
+0.8,0.07793,-0.223143551,-2.551944291
+1,0.07007,0,-2.658260537
+1.25,0.06265,0.223143551,-2.770191598
+1.5,0.05701,0.405465108,-2.864528588
+2,0.04896,0.693147181,-3.016751641
+3,0.03931,1.098612289,-3.23627634
+4,0.03369,1.386294361,-3.390554222
+5,0.03,1.609437912,-3.506557897
+6,0.02741,1.791759469,-3.596847369
+8,0.02401,2.079441542,-3.729284869
+10,0.02192,2.302585093,-3.820355817
+15,0.01915,2.708050201,-3.955452563
+20,0.01786,2.995732274,-4.025191704

src/Numerics/Interpolation/Algorithms/AkimaSplineInterpolation.cs

@@ -134,6 +134,10 @@ namespace MathNet.Numerics.Interpolation.Algorithms
                 throw new ArgumentException(Resources.ArgumentVectorsSameLength);
             }
 
+            for (var i = 1; i < samplePoints.Count; ++i)
+                if (samplePoints[i] <= samplePoints[i - 1])
+                    throw new ArgumentException(Resources.Interpolation_Initialize_SamplePointsNotStrictlyAscendingOrder, "samplePoints");
+
             int n = samplePoints.Count;
 
             /* Prepare divided differences (diff) and weights (w) */

src/Numerics/Interpolation/Algorithms/CubicHermiteSplineInterpolation.cs

@@ -32,6 +32,7 @@ namespace MathNet.Numerics.Interpolation.Algorithms
 {
     using System;
     using System.Collections.Generic;
+    using Properties;
 
     /// <summary>
     /// Cubic Hermite Spline Interpolation Algorithm.
@@ -146,6 +147,10 @@ namespace MathNet.Numerics.Interpolation.Algorithms
                 throw new ArgumentException(Properties.Resources.ArgumentVectorsSameLength);
             }
 
+            for (var i = 1; i < samplePoints.Count; ++i)
+                if (samplePoints[i] <= samplePoints[i - 1])
+                    throw new ArgumentException(Resources.Interpolation_Initialize_SamplePointsNotStrictlyAscendingOrder, "samplePoints");
+
             double[] coefficients = new double[4 * (samplePoints.Count - 1)];
 
             for (int i = 0, j = 0; i < samplePoints.Count - 1; i++, j += 4)

src/Numerics/Interpolation/Algorithms/CubicSplineInterpolation.cs

@@ -205,6 +205,10 @@ namespace MathNet.Numerics.Interpolation.Algorithms
                 throw new ArgumentException(Resources.ArgumentVectorsSameLength);
             }
 
+            for (var i = 1; i < samplePoints.Count; ++i)
+                if (samplePoints[i] <= samplePoints[i - 1])
+                    throw new ArgumentException(Resources.Interpolation_Initialize_SamplePointsNotStrictlyAscendingOrder, "samplePoints");
+
             int n = samplePoints.Count;
 
             // normalize special cases

src/Numerics/Interpolation/Algorithms/LinearSplineInterpolation.cs

@@ -32,6 +32,7 @@ namespace MathNet.Numerics.Interpolation.Algorithms
 {
     using System;
     using System.Collections.Generic;
+    using Properties;
 
     /// <summary>
     /// Linear Spline Interpolation Algorithm.
@@ -133,6 +134,10 @@ namespace MathNet.Numerics.Interpolation.Algorithms
                 throw new ArgumentException(Properties.Resources.ArgumentVectorsSameLength);
             }
 
+            for (var i = 1; i < samplePoints.Count; ++i)
+                if (samplePoints[i] <= samplePoints[i - 1])
+                    throw new ArgumentException(Resources.Interpolation_Initialize_SamplePointsNotStrictlyAscendingOrder, "samplePoints");
+
             double[] coefficients = new double[4 * (samplePoints.Count - 1)];
 
             for (int i = 0, j = 0; i < samplePoints.Count - 1; i++, j += 4)

src/Numerics/Interpolation/Algorithms/NevillePolynomialInterpolation.cs

@@ -32,6 +32,7 @@ namespace MathNet.Numerics.Interpolation.Algorithms
 {
     using System;
     using System.Collections.Generic;
+    using Properties;
 
     /// <summary>
     /// Lagrange Polynomial Interpolation using Neville's Algorithm.
@@ -125,6 +126,10 @@ namespace MathNet.Numerics.Interpolation.Algorithms
                 throw new ArgumentException(Properties.Resources.ArgumentVectorsSameLength);
             }
 
+            for (var i = 1; i < samplePoints.Count; ++i)
+                if (samplePoints[i] <= samplePoints[i - 1])
+                    throw new ArgumentException(Resources.Interpolation_Initialize_SamplePointsNotUnique, "samplePoints");
+
             _points = samplePoints;
             _values = sampleValues;
         }

src/Numerics/Interpolation/Algorithms/SplineInterpolation.cs

@@ -32,6 +32,7 @@ namespace MathNet.Numerics.Interpolation.Algorithms
 {
     using System;
     using System.Collections.Generic;
+    using Properties;
 
     /// <summary>
     /// Third-Degree Spline Interpolation Algorithm.
@@ -123,6 +124,10 @@ namespace MathNet.Numerics.Interpolation.Algorithms
                 throw new ArgumentOutOfRangeException("splineCoefficients");
             }
 
+            for (var i = 1; i < samplePoints.Count; ++i)
+                if (samplePoints[i] <= samplePoints[i - 1])
+                    throw new ArgumentException(Resources.Interpolation_Initialize_SamplePointsNotStrictlyAscendingOrder, "samplePoints");
+
             _points = samplePoints;
             _coefficients = splineCoefficients;
             _sampleCount = samplePoints.Count;

src/Numerics/Properties/Resources.Designer.cs

@@ -1,7 +1,7 @@
 //------------------------------------------------------------------------------
 // <auto-generated>
 //     This code was generated by a tool.
-//     Runtime Version:4.0.30319.239
+//     Runtime Version:4.0.30319.488
 //
 //     Changes to this file may cause incorrect behavior and will be lost if
 //     the code is regenerated.
@@ -564,6 +564,24 @@ namespace MathNet.Numerics.Properties {
             }
         }
         
+        /// <summary>
+        ///   Looks up a localized string similar to Sample points should be sorted in strictly ascending order.
+        /// </summary>
+        internal static string Interpolation_Initialize_SamplePointsNotStrictlyAscendingOrder {
+            get {
+                return ResourceManager.GetString("Interpolation_Initialize_SamplePointsNotStrictlyAscendingOrder", resourceCulture);
+            }
+        }
+        
+        /// <summary>
+        ///   Looks up a localized string similar to All sample points should be unique..
+        /// </summary>
+        internal static string Interpolation_Initialize_SamplePointsNotUnique {
+            get {
+                return ResourceManager.GetString("Interpolation_Initialize_SamplePointsNotUnique", resourceCulture);
+            }
+        }
+        
         /// <summary>
         ///   Looks up a localized string similar to Invalid parameterization for the distribution..
         /// </summary>

src/Numerics/Properties/Resources.resx

@@ -369,4 +369,10 @@
   <data name="TooManyElements" xml:space="preserve">
     <value>We only support sparse matrix with less than int.MaxValue elements.</value>
   </data>
+  <data name="Interpolation_Initialize_SamplePointsNotStrictlyAscendingOrder">
+    <value xml:space="preserve">Sample points should be sorted in strictly ascending order</value>
+  </data>
+  <data name="Interpolation_Initialize_SamplePointsNotUnique">
+    <value xml:space="preserve">All sample points should be unique.</value>
+  </data>
 </root>

src/Numerics/Properties/Resources1.Designer.cs

@@ -1,7 +1,7 @@
 //------------------------------------------------------------------------------
 // <auto-generated>
 //     This code was generated by a tool.
-//     Runtime Version:4.0.30319.239
+//     Runtime Version:4.0.30319.488
 //
 //     Changes to this file may cause incorrect behavior and will be lost if
 //     the code is regenerated.
@@ -564,6 +564,24 @@ namespace MathNet.Numerics.Properties {
             }
         }
         
+        /// <summary>
+        ///   Looks up a localized string similar to Sample points should be sorted in strictly ascending order.
+        /// </summary>
+        internal static string Interpolation_Initialize_SamplePointsNotStrictlyAscendingOrder {
+            get {
+                return ResourceManager.GetString("Interpolation_Initialize_SamplePointsNotStrictlyAscendingOrder", resourceCulture);
+            }
+        }
+        
+        /// <summary>
+        ///   Looks up a localized string similar to All sample points should be unique..
+        /// </summary>
+        internal static string Interpolation_Initialize_SamplePointsNotUnique {
+            get {
+                return ResourceManager.GetString("Interpolation_Initialize_SamplePointsNotUnique", resourceCulture);
+            }
+        }
+        
         /// <summary>
         ///   Looks up a localized string similar to Invalid parameterization for the distribution..
         /// </summary>

src/UnitTests/InterpolationTests/CubicSplineTest.cs

@@ -30,6 +30,10 @@
 
 namespace MathNet.Numerics.UnitTests.InterpolationTests
 {
+    using System;
+    using System.Globalization;
+    using System.IO;
+    using System.Linq;
     using Interpolation;
     using Interpolation.Algorithms;
     using NUnit.Framework;
@@ -220,5 +224,17 @@ namespace MathNet.Numerics.UnitTests.InterpolationTests
                 Assert.AreEqual(ytest[i], interpolation.Interpolate(xtest[i]), 1e-15, "Linear with {0} samples, sample {1}", samples, i);
             }
         }
+
+        /// <summary>
+        /// Verifies that sample points are required to be sorted in strictly monotonically ascending order.
+        /// </summary>
+        [Test]
+        [ExpectedException(typeof(ArgumentException))]
+        public void Constructor_SamplePointsNotStrictlyAscending_Throws()
+        {
+            var x = new[] { -1.0, 0.0, 1.5, 1.5, 2.5, 4.0 };
+            var y = new[] { 1.0, 0.3, -0.7, -0.6, -0.1, 0.4 };
+            var interpolation = new CubicSplineInterpolation(x, y);
+        }
     }
 }

src/UnitTests/InterpolationTests/LinearSplineTest.cs

@@ -30,6 +30,10 @@
 
 namespace MathNet.Numerics.UnitTests.InterpolationTests
 {
+    using System;
+    using System.Globalization;
+    using System.IO;
+    using System.Linq;
     using Interpolation;
     using Interpolation.Algorithms;
     using NUnit.Framework;
@@ -118,5 +122,17 @@ namespace MathNet.Numerics.UnitTests.InterpolationTests
                 Assert.AreEqual(ytest[i], interpolation.Interpolate(xtest[i]), 1e-15, "Linear with {0} samples, sample {1}", samples, i);
             }
         }
+
+        /// <summary>
+        /// Verifies that sample points are required to be sorted in strictly monotonically ascending order.
+        /// </summary>
+        [Test]
+        [ExpectedException(typeof(ArgumentException))]
+        public void Constructor_SamplePointsNotStrictlyAscending_Throws()
+        {
+            var x = new[] { -1.0, 0.0, 1.5, 1.5, 2.5, 4.0 };
+            var y = new[] { 1.0, 0.3, -0.7, -0.6, -0.1, 0.4 };
+            var interpolation = new LinearSplineInterpolation(x, y);
+        }
     }
 }

src/UnitTests/InterpolationTests/NevillePolynomialTest.cs

@@ -30,6 +30,10 @@
 
 namespace MathNet.Numerics.UnitTests.InterpolationTests
 {
+    using System;
+    using System.Globalization;
+    using System.IO;
+    using System.Linq;
     using Interpolation;
     using Interpolation.Algorithms;
     using NUnit.Framework;
@@ -116,5 +120,42 @@ namespace MathNet.Numerics.UnitTests.InterpolationTests
                 Assert.AreEqual(ytest[i], interpolation.Interpolate(xtest[i]), 1e-13, "Linear with {0} samples, sample {1}", samples, i);
             }
         }
+
+        /// <summary>
+        /// Verifies that all sample points must be unique.
+        /// </summary>
+        [Test]
+        [ExpectedException(typeof(ArgumentException))]
+        public void Constructor_SamplePointsNotUnique_Throws()
+        {
+            var x = new[] { -1.0, 0.0, 1.5, 1.5, 2.5, 4.0 };
+            var y = new[] { 1.0, 0.3, -0.7, -0.6, -0.1, 0.4 };
+            var interpolation = new NevillePolynomialInterpolation(x, y);
+        }
+
+        /// <summary>
+        /// Verifies that interpolation does not yield NaN values for a case where some sample points
+        /// has an upper and a lower sample value.
+        /// </summary>
+        /// <param name="value">Value for which interpolation is requested.</param>
+        /// <remarks>Columns 3 and 4 of the .csv file contain the logarithms of energy and mass attenuation
+        /// coefficients, respectively.</remarks>
+        [Test, Sequential, Ignore]
+        public void Interpolate_LogLogAttenuationData_InterpolationShouldNotYieldNaN(
+            [Values(0.0025, 0.035, 0.45, 5.5, 18.5, 35.0)] double value)
+        {
+            var data = File.ReadLines(@"./data/Github-Cureos-1.csv").
+                Select(line =>
+                {
+                    var vals = line.Split(new[] { ',' }, StringSplitOptions.RemoveEmptyEntries);
+                    return Tuple.Create(vals[2], vals[3]);
+                }).ToArray();
+            var x = data.Select(tuple => Double.Parse(tuple.Item1, CultureInfo.InvariantCulture)).ToArray();
+            var y = data.Select(tuple => Double.Parse(tuple.Item2, CultureInfo.InvariantCulture)).ToArray();
+            IInterpolation interpolation = new NevillePolynomialInterpolation(x, y);
+
+            var actual = interpolation.Interpolate(Math.Log(value));
+            Assert.That(actual, Is.Not.NaN);
+        }
     }
 }

src/UnitTests/UnitTests.csproj

@@ -769,6 +769,10 @@
       <Link>data\Codeplex-5667.csv</Link>
       <CopyToOutputDirectory>Always</CopyToOutputDirectory>
     </None>
+    <None Include="..\..\data\Github-Cureos-1.csv">
+      <Link>data\Github-Cureos-1.csv</Link>
+      <CopyToOutputDirectory>Always</CopyToOutputDirectory>
+    </None>
     <None Include="..\..\data\Matlab\A.mat">
       <Link>data\Matlab\A.mat</Link>
       <CopyToOutputDirectory>Always</CopyToOutputDirectory>