From 7e17434889bbba3e699d98a7be1938a120de0d05 Mon Sep 17 00:00:00 2001 From: cureos Date: Tue, 7 Feb 2012 09:54:27 +0100 Subject: [PATCH] Updated interpolation classes to prevent non-unique and/or non-monotonically ascending sample points. Fixes issue #1. --- data/Github-Cureos-1.csv | 49 +++++++++++++++++++ .../Algorithms/AkimaSplineInterpolation.cs | 4 ++ .../CubicHermiteSplineInterpolation.cs | 5 ++ .../Algorithms/CubicSplineInterpolation.cs | 4 ++ .../Algorithms/LinearSplineInterpolation.cs | 5 ++ .../NevillePolynomialInterpolation.cs | 5 ++ .../Algorithms/SplineInterpolation.cs | 5 ++ src/Numerics/Properties/Resources.Designer.cs | 20 +++++++- src/Numerics/Properties/Resources.resx | 6 +++ .../Properties/Resources1.Designer.cs | 20 +++++++- .../InterpolationTests/CubicSplineTest.cs | 16 ++++++ .../InterpolationTests/LinearSplineTest.cs | 16 ++++++ .../NevillePolynomialTest.cs | 41 ++++++++++++++++ src/UnitTests/UnitTests.csproj | 4 ++ 14 files changed, 198 insertions(+), 2 deletions(-) create mode 100644 data/Github-Cureos-1.csv diff --git a/data/Github-Cureos-1.csv b/data/Github-Cureos-1.csv new file mode 100644 index 00000000..62bde724 --- /dev/null +++ b/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 diff --git a/src/Numerics/Interpolation/Algorithms/AkimaSplineInterpolation.cs b/src/Numerics/Interpolation/Algorithms/AkimaSplineInterpolation.cs index 1a1161f6..98b3016a 100644 --- a/src/Numerics/Interpolation/Algorithms/AkimaSplineInterpolation.cs +++ b/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) */ diff --git a/src/Numerics/Interpolation/Algorithms/CubicHermiteSplineInterpolation.cs b/src/Numerics/Interpolation/Algorithms/CubicHermiteSplineInterpolation.cs index 0e4d257d..8e22b6a0 100644 --- a/src/Numerics/Interpolation/Algorithms/CubicHermiteSplineInterpolation.cs +++ b/src/Numerics/Interpolation/Algorithms/CubicHermiteSplineInterpolation.cs @@ -32,6 +32,7 @@ namespace MathNet.Numerics.Interpolation.Algorithms { using System; using System.Collections.Generic; + using Properties; /// /// 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) diff --git a/src/Numerics/Interpolation/Algorithms/CubicSplineInterpolation.cs b/src/Numerics/Interpolation/Algorithms/CubicSplineInterpolation.cs index 33f66163..85b5ecf8 100644 --- a/src/Numerics/Interpolation/Algorithms/CubicSplineInterpolation.cs +++ b/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 diff --git a/src/Numerics/Interpolation/Algorithms/LinearSplineInterpolation.cs b/src/Numerics/Interpolation/Algorithms/LinearSplineInterpolation.cs index e403bf40..54146c3a 100644 --- a/src/Numerics/Interpolation/Algorithms/LinearSplineInterpolation.cs +++ b/src/Numerics/Interpolation/Algorithms/LinearSplineInterpolation.cs @@ -32,6 +32,7 @@ namespace MathNet.Numerics.Interpolation.Algorithms { using System; using System.Collections.Generic; + using Properties; /// /// 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) diff --git a/src/Numerics/Interpolation/Algorithms/NevillePolynomialInterpolation.cs b/src/Numerics/Interpolation/Algorithms/NevillePolynomialInterpolation.cs index a0596bdb..755b525a 100644 --- a/src/Numerics/Interpolation/Algorithms/NevillePolynomialInterpolation.cs +++ b/src/Numerics/Interpolation/Algorithms/NevillePolynomialInterpolation.cs @@ -32,6 +32,7 @@ namespace MathNet.Numerics.Interpolation.Algorithms { using System; using System.Collections.Generic; + using Properties; /// /// 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; } diff --git a/src/Numerics/Interpolation/Algorithms/SplineInterpolation.cs b/src/Numerics/Interpolation/Algorithms/SplineInterpolation.cs index f7546aa8..c654b21f 100644 --- a/src/Numerics/Interpolation/Algorithms/SplineInterpolation.cs +++ b/src/Numerics/Interpolation/Algorithms/SplineInterpolation.cs @@ -32,6 +32,7 @@ namespace MathNet.Numerics.Interpolation.Algorithms { using System; using System.Collections.Generic; + using Properties; /// /// 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; diff --git a/src/Numerics/Properties/Resources.Designer.cs b/src/Numerics/Properties/Resources.Designer.cs index 37163f68..5ba402b7 100644 --- a/src/Numerics/Properties/Resources.Designer.cs +++ b/src/Numerics/Properties/Resources.Designer.cs @@ -1,7 +1,7 @@ //------------------------------------------------------------------------------ // // 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 { } } + /// + /// Looks up a localized string similar to Sample points should be sorted in strictly ascending order. + /// + internal static string Interpolation_Initialize_SamplePointsNotStrictlyAscendingOrder { + get { + return ResourceManager.GetString("Interpolation_Initialize_SamplePointsNotStrictlyAscendingOrder", resourceCulture); + } + } + + /// + /// Looks up a localized string similar to All sample points should be unique.. + /// + internal static string Interpolation_Initialize_SamplePointsNotUnique { + get { + return ResourceManager.GetString("Interpolation_Initialize_SamplePointsNotUnique", resourceCulture); + } + } + /// /// Looks up a localized string similar to Invalid parameterization for the distribution.. /// diff --git a/src/Numerics/Properties/Resources.resx b/src/Numerics/Properties/Resources.resx index b6925415..2d7c5c10 100644 --- a/src/Numerics/Properties/Resources.resx +++ b/src/Numerics/Properties/Resources.resx @@ -369,4 +369,10 @@ We only support sparse matrix with less than int.MaxValue elements. + + Sample points should be sorted in strictly ascending order + + + All sample points should be unique. + \ No newline at end of file diff --git a/src/Numerics/Properties/Resources1.Designer.cs b/src/Numerics/Properties/Resources1.Designer.cs index 37163f68..5ba402b7 100644 --- a/src/Numerics/Properties/Resources1.Designer.cs +++ b/src/Numerics/Properties/Resources1.Designer.cs @@ -1,7 +1,7 @@ //------------------------------------------------------------------------------ // // 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 { } } + /// + /// Looks up a localized string similar to Sample points should be sorted in strictly ascending order. + /// + internal static string Interpolation_Initialize_SamplePointsNotStrictlyAscendingOrder { + get { + return ResourceManager.GetString("Interpolation_Initialize_SamplePointsNotStrictlyAscendingOrder", resourceCulture); + } + } + + /// + /// Looks up a localized string similar to All sample points should be unique.. + /// + internal static string Interpolation_Initialize_SamplePointsNotUnique { + get { + return ResourceManager.GetString("Interpolation_Initialize_SamplePointsNotUnique", resourceCulture); + } + } + /// /// Looks up a localized string similar to Invalid parameterization for the distribution.. /// diff --git a/src/UnitTests/InterpolationTests/CubicSplineTest.cs b/src/UnitTests/InterpolationTests/CubicSplineTest.cs index 4b249274..ccaa1894 100644 --- a/src/UnitTests/InterpolationTests/CubicSplineTest.cs +++ b/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); } } + + /// + /// Verifies that sample points are required to be sorted in strictly monotonically ascending order. + /// + [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); + } } } diff --git a/src/UnitTests/InterpolationTests/LinearSplineTest.cs b/src/UnitTests/InterpolationTests/LinearSplineTest.cs index 5fee9710..dba3da9b 100644 --- a/src/UnitTests/InterpolationTests/LinearSplineTest.cs +++ b/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); } } + + /// + /// Verifies that sample points are required to be sorted in strictly monotonically ascending order. + /// + [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); + } } } diff --git a/src/UnitTests/InterpolationTests/NevillePolynomialTest.cs b/src/UnitTests/InterpolationTests/NevillePolynomialTest.cs index 1c5abfcb..533ef58a 100644 --- a/src/UnitTests/InterpolationTests/NevillePolynomialTest.cs +++ b/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); } } + + /// + /// Verifies that all sample points must be unique. + /// + [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); + } + + /// + /// Verifies that interpolation does not yield NaN values for a case where some sample points + /// has an upper and a lower sample value. + /// + /// Value for which interpolation is requested. + /// Columns 3 and 4 of the .csv file contain the logarithms of energy and mass attenuation + /// coefficients, respectively. + [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); + } } } diff --git a/src/UnitTests/UnitTests.csproj b/src/UnitTests/UnitTests.csproj index 8408f2b2..d9c822ce 100644 --- a/src/UnitTests/UnitTests.csproj +++ b/src/UnitTests/UnitTests.csproj @@ -769,6 +769,10 @@ data\Codeplex-5667.csv Always + + data\Github-Cureos-1.csv + Always + data\Matlab\A.mat Always