From 5d9b0143e353fc4904ae8cbe3ceba55b76d26079 Mon Sep 17 00:00:00 2001 From: Christoph Ruegg Date: Fri, 20 Dec 2013 17:22:57 +0100 Subject: [PATCH] Interpolation: more reasonable LinearSpline implementation --- src/Numerics/Interpolate.cs | 6 +- src/Numerics/Interpolation/LinearSpline.cs | 182 +++++++++++++++++ .../LinearSplineInterpolation.cs | 188 ------------------ src/Numerics/Numerics.csproj | 2 +- .../InterpolationTests/LinearSplineTest.cs | 29 +-- 5 files changed, 193 insertions(+), 214 deletions(-) create mode 100644 src/Numerics/Interpolation/LinearSpline.cs delete mode 100644 src/Numerics/Interpolation/LinearSplineInterpolation.cs diff --git a/src/Numerics/Interpolate.cs b/src/Numerics/Interpolate.cs index 475aa423..72e63eba 100644 --- a/src/Numerics/Interpolate.cs +++ b/src/Numerics/Interpolate.cs @@ -63,11 +63,9 @@ namespace MathNet.Numerics /// which can then be used to compute interpolations and extrapolations /// on arbitrary points. /// - public static IInterpolation LinearBetweenPoints(IList points, IList values) + public static IInterpolation LinearBetweenPoints(IEnumerable points, IEnumerable values) { - var method = new LinearSplineInterpolation(); - method.Initialize(points, values); - return method; + return LinearSpline.Interpolate(points, values); } /// diff --git a/src/Numerics/Interpolation/LinearSpline.cs b/src/Numerics/Interpolation/LinearSpline.cs new file mode 100644 index 00000000..8996429f --- /dev/null +++ b/src/Numerics/Interpolation/LinearSpline.cs @@ -0,0 +1,182 @@ +// +// Math.NET Numerics, part of the Math.NET Project +// http://numerics.mathdotnet.com +// http://github.com/mathnet/mathnet-numerics +// http://mathnetnumerics.codeplex.com +// +// Copyright (c) 2009-2013 Math.NET +// +// Permission is hereby granted, free of charge, to any person +// obtaining a copy of this software and associated documentation +// files (the "Software"), to deal in the Software without +// restriction, including without limitation the rights to use, +// copy, modify, merge, publish, distribute, sublicense, and/or sell +// copies of the Software, and to permit persons to whom the +// Software is furnished to do so, subject to the following +// conditions: +// +// The above copyright notice and this permission notice shall be +// included in all copies or substantial portions of the Software. +// +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, +// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES +// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND +// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT +// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, +// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING +// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR +// OTHER DEALINGS IN THE SOFTWARE. +// + +using System; +using System.Collections.Generic; +using System.Linq; +using MathNet.Numerics.Properties; + +namespace MathNet.Numerics.Interpolation +{ + /// + /// Linear Spline Interpolation. + /// + /// Supports both differentiation and integration. + public class LinearSpline : IInterpolation + { + readonly double[] _x; + readonly double[] _c0; + readonly double[] _c1; + + /// Sample points (N+1), sorted ascending + /// Sample values (N or N+1) at the corresponding points; intercept, zero order coefficients + /// Slopes (N) at the sample points (first order coefficients): N + public LinearSpline(double[] x, double[] c0, double[] c1) + { + _x = x; + _c0 = c0; + _c1 = c1; + } + + /// + /// Create a linear spline interpolation from a set of (x,y) value pairs. + /// + /// + /// The value pairs do not have to be sorted, but if they are not sorted ascendingly + /// and the passed x and y arguments are arrays, they will be sorted inplace and thus modified. + /// + public static LinearSpline Interpolate(IEnumerable x, IEnumerable y) + { + var xx = (x as double[]) ?? x.ToArray(); + var yy = (y as double[]) ?? y.ToArray(); + + if (xx.Length != yy.Length) + { + throw new ArgumentException(Resources.ArgumentVectorsSameLength); + } + + Sorting.Sort(xx, yy); + + var c1 = new double[xx.Length - 1]; + for (int i = 0; i < c1.Length; i++) + { + c1[i] = (yy[i + 1] - yy[i])/(xx[i + 1] - xx[i]); + } + + return new LinearSpline(xx, yy, c1); + } + + /// + /// Gets a value indicating whether the algorithm supports differentiation (interpolated derivative). + /// + public bool SupportsDifferentiation + { + get { return true; } + } + + /// + /// Gets a value indicating whether the algorithm supports integration (interpolated quadrature). + /// + public bool SupportsIntegration + { + get { return true; } + } + + /// + /// Interpolate at point t. + /// + /// Point t to interpolate at. + /// Interpolated value x(t). + public double Interpolate(double t) + { + int k = LeftBracketIndex(t); + return _c0[k] + (t - _x[k])*_c1[k]; + } + + /// + /// Differentiate at point t. + /// + /// Point t to interpolate at. + /// Interpolated first derivative at point t. + public double Differentiate(double t) + { + int k = LeftBracketIndex(t); + return _c1[k]; + } + + /// + /// Differentiate twice at point t. + /// + /// Point t to interpolate at. + /// Interpolated second derivative at point t. + public double Differentiate2(double t) + { + return 0d; + } + + /// + /// Integrate up to point t. + /// + /// Right bound of the integration interval [a,t]. + /// Interpolated definite integral over the interval [a,t]. + public double Integrate(double t) + { + int k = LeftBracketIndex(t); + var x = (t - _x[k]); + return x*(_c0[k] + x*0.5*_c1[k]); + } + + /// + /// Interpolate, differentiate and 2nd differentiate at point t. + /// + /// Point t to interpolate at. + /// Interpolated first derivative at point t. + public Tuple DifferentiateAll(double t) + { + int k = LeftBracketIndex(t); + var x = (t - _x[k]); + return new Tuple(_c0[k] + x*_c1[k], _c1[k], 0d); + } + + /// + /// Find the index of the greatest sample point smaller than t. + /// + int LeftBracketIndex(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) + { + int middle = (low + high)/2; + if (_x[middle] > t) + { + high = middle; + } + else + { + low = middle; + } + } + + return low; + } + } +} diff --git a/src/Numerics/Interpolation/LinearSplineInterpolation.cs b/src/Numerics/Interpolation/LinearSplineInterpolation.cs deleted file mode 100644 index 7b413d80..00000000 --- a/src/Numerics/Interpolation/LinearSplineInterpolation.cs +++ /dev/null @@ -1,188 +0,0 @@ -// -// Math.NET Numerics, part of the Math.NET Project -// http://numerics.mathdotnet.com -// http://github.com/mathnet/mathnet-numerics -// http://mathnetnumerics.codeplex.com -// -// Copyright (c) 2009-2010 Math.NET -// -// Permission is hereby granted, free of charge, to any person -// obtaining a copy of this software and associated documentation -// files (the "Software"), to deal in the Software without -// restriction, including without limitation the rights to use, -// copy, modify, merge, publish, distribute, sublicense, and/or sell -// copies of the Software, and to permit persons to whom the -// Software is furnished to do so, subject to the following -// conditions: -// -// The above copyright notice and this permission notice shall be -// included in all copies or substantial portions of the Software. -// -// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, -// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES -// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND -// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT -// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, -// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING -// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR -// OTHER DEALINGS IN THE SOFTWARE. -// - -using System; -using System.Collections.Generic; -using MathNet.Numerics.Properties; - -namespace MathNet.Numerics.Interpolation -{ - /// - /// Linear Spline Interpolation Algorithm. - /// - /// - /// This algorithm supports both differentiation and integration. - /// - public class LinearSplineInterpolation : IInterpolation - { - /// - /// Internal Spline Interpolation - /// - readonly SplineInterpolation _spline; - - /// - /// Initializes a new instance of the LinearSplineInterpolation class. - /// - public LinearSplineInterpolation() - { - _spline = new SplineInterpolation(); - } - - /// - /// Initializes a new instance of the LinearSplineInterpolation class. - /// - /// Sample Points t, sorted ascending. - /// Sample Values x(t) - public LinearSplineInterpolation(IList samplePoints, IList sampleValues) - { - _spline = new SplineInterpolation(); - Initialize(samplePoints, sampleValues); - } - - /// - /// Gets a value indicating whether the algorithm supports differentiation (interpolated derivative). - /// - /// - /// - bool IInterpolation.SupportsDifferentiation - { - get { return true; } - } - - /// - /// Gets a value indicating whether the algorithm supports integration (interpolated quadrature). - /// - /// - bool IInterpolation.SupportsIntegration - { - get { return true; } - } - - /// - /// Initialize the interpolation method with the given spline coefficients (sorted by the sample points t). - /// - /// Sample Points t, sorted ascending. - /// Sample Values x(t) - public void Initialize(IList samplePoints, IList sampleValues) - { - double[] coefficients = EvaluateSplineCoefficients(samplePoints, sampleValues); - _spline.Initialize(samplePoints, coefficients); - } - - /// - /// Evaluate the spline coefficients as used - /// internally by this interpolation algorithm. - /// - /// Sample Points t, sorted ascending. - /// Sample Values x(t) - /// Spline Coefficient Vector - public static double[] EvaluateSplineCoefficients(IList samplePoints, IList sampleValues) - { - if (null == samplePoints) - { - throw new ArgumentNullException("samplePoints"); - } - - if (null == sampleValues) - { - throw new ArgumentNullException("sampleValues"); - } - - if (samplePoints.Count < 2) - { - throw new ArgumentOutOfRangeException("samplePoints"); - } - - if (samplePoints.Count != sampleValues.Count) - { - 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"); - - var coefficients = new double[4*(samplePoints.Count - 1)]; - for (int i = 0, j = 0; i < samplePoints.Count - 1; i++, j += 4) - { - coefficients[j] = sampleValues[i]; - coefficients[j + 1] = (sampleValues[i + 1] - sampleValues[i])/(samplePoints[i + 1] - samplePoints[i]); - coefficients[j + 2] = 0; - coefficients[j + 3] = 0; - } - return coefficients; - } - - /// - /// Interpolate at point t. - /// - /// Point t to interpolate at. - /// Interpolated value x(t). - public double Interpolate(double t) - { - return _spline.Interpolate(t); - } - - /// - /// Differentiate at point t. - /// - /// Point t to interpolate at. - /// Interpolated first derivative at point t. - /// - /// - public double Differentiate(double t) - { - return _spline.Differentiate(t); - } - - /// - /// Interpolate, differentiate and 2nd differentiate at point t. - /// - /// Point t to interpolate at. - /// Interpolated first derivative at point t. - /// - /// - public Tuple DifferentiateAll(double t) - { - return _spline.DifferentiateAll(t); - } - - /// - /// Integrate up to point t. - /// - /// Right bound of the integration interval [a,t]. - /// Interpolated definite integral over the interval [a,t]. - /// - public double Integrate(double t) - { - return _spline.Integrate(t); - } - } -} diff --git a/src/Numerics/Numerics.csproj b/src/Numerics/Numerics.csproj index b338bde0..37cbe59c 100644 --- a/src/Numerics/Numerics.csproj +++ b/src/Numerics/Numerics.csproj @@ -383,7 +383,7 @@ - + diff --git a/src/UnitTests/InterpolationTests/LinearSplineTest.cs b/src/UnitTests/InterpolationTests/LinearSplineTest.cs index 192e0c9e..fea52514 100644 --- a/src/UnitTests/InterpolationTests/LinearSplineTest.cs +++ b/src/UnitTests/InterpolationTests/LinearSplineTest.cs @@ -28,12 +28,11 @@ // OTHER DEALINGS IN THE SOFTWARE. // +using MathNet.Numerics.Interpolation; +using NUnit.Framework; + namespace MathNet.Numerics.UnitTests.InterpolationTests { - using System; - using Interpolation; - using NUnit.Framework; - /// /// LinearSpline test case /// @@ -43,12 +42,12 @@ namespace MathNet.Numerics.UnitTests.InterpolationTests /// /// Sample points. /// - private readonly double[] _t = new[] { -2.0, -1.0, 0.0, 1.0, 2.0 }; + readonly double[] _t = { -2.0, -1.0, 0.0, 1.0, 2.0 }; /// /// Sample values. /// - private readonly double[] _x = new[] { 1.0, 2.0, -1.0, 0.0, 1.0 }; + readonly double[] _x = { 1.0, 2.0, -1.0, 0.0, 1.0 }; /// /// Verifies that the interpolation matches the given value at all the provided sample points. @@ -56,7 +55,7 @@ namespace MathNet.Numerics.UnitTests.InterpolationTests [Test] public void FitsAtSamplePoints() { - IInterpolation interpolation = new LinearSplineInterpolation(_t, _x); + IInterpolation interpolation = LinearSpline.Interpolate(_t, _x); for (int i = 0; i < _x.Length; i++) { @@ -89,7 +88,7 @@ namespace MathNet.Numerics.UnitTests.InterpolationTests [TestCase(-10.0, -7.0, 1e-15)] public void FitsAtArbitraryPointsWithMaple(double t, double x, double maxAbsoluteError) { - IInterpolation interpolation = new LinearSplineInterpolation(_t, _x); + IInterpolation interpolation = LinearSpline.Interpolate(_t, _x); Assert.AreEqual(x, interpolation.Interpolate(t), maxAbsoluteError, "Interpolation at {0}", t); @@ -108,23 +107,11 @@ namespace MathNet.Numerics.UnitTests.InterpolationTests { double[] x, y, xtest, ytest; LinearInterpolationCase.Build(out x, out y, out xtest, out ytest, samples); - IInterpolation interpolation = new LinearSplineInterpolation(x, y); + IInterpolation interpolation = LinearSpline.Interpolate(x, y); for (int i = 0; i < xtest.Length; i++) { 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); - } } }