diff --git a/src/Numerics.Tests/InterpolationTests/PchipSplineTest.cs b/src/Numerics.Tests/InterpolationTests/PchipSplineTest.cs new file mode 100644 index 00000000..806c191f --- /dev/null +++ b/src/Numerics.Tests/InterpolationTests/PchipSplineTest.cs @@ -0,0 +1,130 @@ +// +// Math.NET Numerics, part of the Math.NET Project +// http://numerics.mathdotnet.com +// http://github.com/mathnet/mathnet-numerics +// +// Copyright (c) 2009-2016 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 MathNet.Numerics.Interpolation; +using NUnit.Framework; + +namespace MathNet.Numerics.UnitTests.InterpolationTests +{ + [TestFixture, Category("Interpolation")] + public class PchipSplineTest + { + readonly double[] _t = { -2.0, -1.0, 0.0, 1.0, 2.0 }; + readonly double[] _y = { 1.0, 2.0, -1.0, 0.0, 1.0 }; + + readonly double[] _tNag = new double[] { 7.99, 8.09, 8.19, 8.70, 9.20, 10.00, 12.00, 15.00, 20.00 }; + readonly double[] _yNag = new double[] { 0.00000E+0, 0.27643E-4, 0.43750E-1, 0.16918E+0, 0.46943E+0, 0.94374E+0, 0.99864E+0, 0.99992E+0, 0.99999E+0 }; + + /// + /// Verifies that the interpolation matches the given value at all the provided sample points. + /// + [Test] + public void FitsAtSamplePoints() + { + IInterpolation it = CubicSpline.InterpolatePchip(_t, _y); + for (int i = 0; i < _y.Length; i++) + { + Assert.AreEqual(_y[i], it.Interpolate(_t[i]), "A Exact Point " + i); + } + } + + /// + /// Verifies that at points other than the provided sample points, the interpolation matches the one computed by Octave as a reference. + /// + /// Sample point. + /// Sample value. + /// Maximum absolute error. + [TestCase(-2.4, -0.7440, 1e-15)] + [TestCase(-0.9, 1.9160, 1e-15)] + [TestCase(-0.5, 0.5000, 1e-15)] + [TestCase(-0.1, -0.9160, 1e-15)] + [TestCase(0.1, -0.9810, 1e-15)] + [TestCase(0.4, -0.7440, 1e-15)] + [TestCase(1.2, 0.2000, 1e-15)] + [TestCase(10.0, 9.0000, 1e-15)] + [TestCase(-10.0, -727.0000, 1e-15)] + public void FitsAtArbitraryPoints(double t, double x, double maxAbsoluteError) + { + IInterpolation it = CubicSpline.InterpolatePchip(_t, _y); + + Assert.AreEqual(x, it.Interpolate(t), maxAbsoluteError, "Interpolation at {0}", t); + } + + /// + /// Verifies that at points other than the provided sample points, the interpolation matches the one computed by NAG as a reference. + /// Reference: https://www.nag.com/numeric/cl/nagdoc_cl25/html/e01/e01bec.html + /// + /// Sample point. + /// Sample value. + /// Maximum absolute error. + [TestCase(7.9900, 0.0000, 5e-5)] + [TestCase(9.1910, 0.4640, 5e-5)] + [TestCase(10.3920, 0.9645, 5e-5)] + [TestCase(11.5930, 0.9965, 5e-5)] + [TestCase(12.7940, 0.9992, 5e-5)] + [TestCase(13.9950, 0.9998, 5e-5)] + [TestCase(15.1960, 0.9999, 5e-5)] + [TestCase(16.3970, 1.0000, 5e-5)] + [TestCase(17.5980, 1.0000, 5e-5)] + [TestCase(18.7990, 1.0000, 5e-5)] + [TestCase(20.0000, 1.0000, 5e-5)] + public void FitsAtExamplePoints(double t, double x, double maxAbsoluteError) + { + IInterpolation it = CubicSpline.InterpolatePchip(_tNag, _yNag); + + Assert.AreEqual(x, it.Interpolate(t), maxAbsoluteError, "Interpolation at {0}", t); + } + + /// + /// Verifies that the interpolation supports the linear case appropriately + /// + /// Samples array. + [TestCase(5)] + [TestCase(7)] + [TestCase(15)] + public void SupportsLinearCase(int samples) + { + double[] x, y, xtest, ytest; + LinearInterpolationCase.Build(out x, out y, out xtest, out ytest, samples); + IInterpolation it = CubicSpline.InterpolatePchip(x, y); + for (int i = 0; i < xtest.Length; i++) + { + 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.InterpolatePchip(new double[0], new double[0]), Throws.ArgumentException); + Assert.That(() => CubicSpline.InterpolatePchip(new double[2], new double[2]), Throws.ArgumentException); + Assert.That(CubicSpline.InterpolatePchip(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)); + } + } +} diff --git a/src/Numerics/Interpolation/CubicSpline.cs b/src/Numerics/Interpolation/CubicSpline.cs index 5da1a2e6..29159bd1 100644 --- a/src/Numerics/Interpolation/CubicSpline.cs +++ b/src/Numerics/Interpolation/CubicSpline.cs @@ -212,6 +212,7 @@ namespace MathNet.Numerics.Interpolation /// /// Create a piecewise cubic Hermite interpolating polynomial from an unsorted set of (x,y) value pairs. + /// Monotone-preserving interpolation with continuous first derivative. /// public static CubicSpline InterpolatePchipSorted(double[] x, double[] y) { @@ -219,54 +220,69 @@ namespace MathNet.Numerics.Interpolation { throw new ArgumentException("All vectors must have the same dimensionality."); } - // TODO: minsize? - if (x.Length < 5) + + if (x.Length < 3) { - throw new ArgumentException("The given array is too small. It must be at least 5 long.", nameof(x)); + throw new ArgumentException("The given array is too small. It must be at least 3 long.", nameof(x)); } - var h = new double[x.Length - 1]; - // The slopes between each x. var m = new double[x.Length - 1]; - // The slope of the interpolant at each x. - var d = new double[x.Length]; - for (var k = 0; k < x.Length - 1; ++k) + for (int i = 0; i < m.Length; i++) + { + m[i] = (y[i + 1] - y[i])/(x[i + 1] - x[i]); + } + + var dd = new double[x.Length]; + var hPrev = x[1] - x[0]; + // This check is quite costly as it usually involves a Math.Pow(). + var mPrevIs0 = m[0].AlmostEqual(0.0); + + for (var i = 1; i < x.Length - 1; ++i) { - h[k] = x[k + 1] - x[k]; - m[k] = (y[k + 1] - y[k]) / h[k]; - if (k == 0) - continue; - if (m[k].AlmostEqual(0.0) || m[k - 1].AlmostEqual(0.0) || Math.Sign(m[k]) != Math.Sign(m[k-1])) - d[k] = 0; + var h = x[i + 1] - x[i]; + var mIs0 = m[i].AlmostEqual(0.0); + + if (mIs0 || mPrevIs0 || Math.Sign(m[i]) != Math.Sign(m[i - 1])) + dd[i] = 0; else { // Weighted harmonic mean of each slope. - var w1 = 2 * h[k] + h[k - 1]; - var w2 = h[k] + 2 * h[k - 1]; - d[k] = (w1 + w2) / (w1 / m[k - 1] + w2 / m[k]); + var w1 = 2 * h + hPrev; + var w2 = h + 2 * hPrev; + dd[i] = (w1 + w2) / (w1 / m[i - 1] + w2 / m[i]); } + + hPrev = h; + mPrevIs0 = mIs0; } // Special case end-points. - d[0] = PchipEndPoints(h[0], h[1], m[0], m[1]); - d[d.Length - 1] = PchipEndPoints(h[h.Length - 1], h[h.Length - 2], m[m.Length - 1], m[m.Length - 2]); + dd[0] = PchipEndPoints(x[1] - x[0], x[2] - x[1], m[0], m[1]); + dd[dd.Length - 1] = PchipEndPoints(x[x.Length - 1] - x[x.Length - 2], x[x.Length - 2] - x[x.Length - 3], + m[m.Length - 1], m[m.Length - 2]); - return InterpolateHermiteSorted(x, y, d); + return InterpolateHermiteSorted(x, y, dd); } - private static double PchipEndPoints(double h0, double h1, double m0, double m1) + static double PchipEndPoints(double h0, double h1, double m0, double m1) { + // One-sided three-point estimate for the derivative. var d = ((2 * h0 + h1) * m0 - h0 * m1) / (h0 + h1); + if (Math.Sign(d) != Math.Sign(m0)) return 0.0; + if (Math.Sign(m0) != Math.Sign(m1) && (Math.Abs(d) > 3 * Math.Abs(m0))) return 3 * m0; + return d; } + /// /// Create a piecewise cubic Hermite interpolating polynomial from an unsorted set of (x,y) value pairs. + /// Monotone-preserving interpolation with continuous first derivative. /// WARNING: Works in-place and can thus causes the data array to be reordered. /// public static CubicSpline InterpolatePchipInplace(double[] x, double[] y) @@ -282,6 +298,7 @@ namespace MathNet.Numerics.Interpolation /// /// Create a piecewise cubic Hermite interpolating polynomial from an unsorted set of (x,y) value pairs. + /// Monotone-preserving interpolation with continuous first derivative. /// public static CubicSpline InterpolatePchip(IEnumerable x, IEnumerable y) { @@ -477,7 +494,7 @@ namespace MathNet.Numerics.Interpolation double t1 = xx[index1] - xx[index0]; double t2 = xx[index2] - xx[index0]; - double a = (x2 - x0 - (t2/t1*(x1 - x0)))/(t2*t2 - t1*t2); + double a = (x2 - x0 - (t2/t1*(x1 - x0)))/(t2*(t2 - t1)); double b = (x1 - x0 - a*t1*t1)/t1; return (2*a*t) + b; }