From 9e1734638c9a8afbd521e62d31a117e16193559e Mon Sep 17 00:00:00 2001 From: Febin <25066330+febkor@users.noreply.github.com> Date: Mon, 28 Dec 2020 17:13:40 +0200 Subject: [PATCH 1/5] Add piecewise cubic interpolating polynomial function. --- src/Numerics/Interpolation/CubicSpline.cs | 79 +++++++++++++++++++++++ 1 file changed, 79 insertions(+) diff --git a/src/Numerics/Interpolation/CubicSpline.cs b/src/Numerics/Interpolation/CubicSpline.cs index f2284506..5da1a2e6 100644 --- a/src/Numerics/Interpolation/CubicSpline.cs +++ b/src/Numerics/Interpolation/CubicSpline.cs @@ -210,6 +210,85 @@ namespace MathNet.Numerics.Interpolation return InterpolateAkimaInplace(x.ToArray(), y.ToArray()); } + /// + /// Create a piecewise cubic Hermite interpolating polynomial from an unsorted set of (x,y) value pairs. + /// + public static CubicSpline InterpolatePchipSorted(double[] x, double[] y) + { + if (x.Length != y.Length) + { + throw new ArgumentException("All vectors must have the same dimensionality."); + } + // TODO: minsize? + if (x.Length < 5) + { + throw new ArgumentException("The given array is too small. It must be at least 5 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) + { + 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; + 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]); + } + } + + // 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]); + + return InterpolateHermiteSorted(x, y, d); + } + + private static double PchipEndPoints(double h0, double h1, double m0, double m1) + { + 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. + /// WARNING: Works in-place and can thus causes the data array to be reordered. + /// + public static CubicSpline InterpolatePchipInplace(double[] x, double[] y) + { + if (x.Length != y.Length) + { + throw new ArgumentException("All vectors must have the same dimensionality."); + } + + Sorting.Sort(x, y); + return InterpolatePchipSorted(x, y); + } + + /// + /// Create a piecewise cubic Hermite interpolating polynomial from an unsorted set of (x,y) value pairs. + /// + public static CubicSpline InterpolatePchip(IEnumerable x, IEnumerable y) + { + // note: we must make a copy, even if the input was arrays already + return InterpolatePchipInplace(x.ToArray(), y.ToArray()); + } + /// /// Create a cubic spline interpolation from a set of (x,y) value pairs, sorted ascendingly by x, /// and custom boundary/termination conditions. From 8e5c6136b66fb40271cc9a94064a6c6c7798bc11 Mon Sep 17 00:00:00 2001 From: Febin <25066330+febkor@users.noreply.github.com> Date: Mon, 28 Dec 2020 23:01:01 +0200 Subject: [PATCH 2/5] Refactored to minimize allocations; added tests; make the code look more similar to Akima --- .../InterpolationTests/PchipSplineTest.cs | 130 ++++++++++++++++++ src/Numerics/Interpolation/CubicSpline.cs | 61 +++++--- 2 files changed, 169 insertions(+), 22 deletions(-) create mode 100644 src/Numerics.Tests/InterpolationTests/PchipSplineTest.cs 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; } From 054292fe3122c754b57a5398055c6e0c47c0abd7 Mon Sep 17 00:00:00 2001 From: Febin <25066330+febkor@users.noreply.github.com> Date: Tue, 29 Dec 2020 09:39:09 +0200 Subject: [PATCH 3/5] Minor clean-up. --- .../InterpolationTests/PchipSplineTest.cs | 10 +++++----- src/Numerics/Interpolation/CubicSpline.cs | 6 +++--- 2 files changed, 8 insertions(+), 8 deletions(-) diff --git a/src/Numerics.Tests/InterpolationTests/PchipSplineTest.cs b/src/Numerics.Tests/InterpolationTests/PchipSplineTest.cs index 806c191f..1f74751e 100644 --- a/src/Numerics.Tests/InterpolationTests/PchipSplineTest.cs +++ b/src/Numerics.Tests/InterpolationTests/PchipSplineTest.cs @@ -3,7 +3,7 @@ // http://numerics.mathdotnet.com // http://github.com/mathnet/mathnet-numerics // -// Copyright (c) 2009-2016 Math.NET +// Copyright (c) 2009-2021 Math.NET // // Permission is hereby granted, free of charge, to any person // obtaining a copy of this software and associated documentation @@ -38,8 +38,8 @@ namespace MathNet.Numerics.UnitTests.InterpolationTests 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 }; + readonly double[] _tNag = { 7.99, 8.09, 8.19, 8.70, 9.20, 10.00, 12.00, 15.00, 20.00 }; + readonly double[] _yNag = { 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. @@ -94,7 +94,7 @@ namespace MathNet.Numerics.UnitTests.InterpolationTests [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) + public void FitsAtNagExamplePoints(double t, double x, double maxAbsoluteError) { IInterpolation it = CubicSpline.InterpolatePchip(_tNag, _yNag); @@ -124,7 +124,7 @@ namespace MathNet.Numerics.UnitTests.InterpolationTests { 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)); + Assert.That(CubicSpline.InterpolatePchip(new[] { 1.0, 2.0, 3.0 }, new[] { 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 29159bd1..5dc52146 100644 --- a/src/Numerics/Interpolation/CubicSpline.cs +++ b/src/Numerics/Interpolation/CubicSpline.cs @@ -259,7 +259,8 @@ namespace MathNet.Numerics.Interpolation // Special case end-points. 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], + 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, dd); @@ -267,7 +268,7 @@ namespace MathNet.Numerics.Interpolation static double PchipEndPoints(double h0, double h1, double m0, double m1) { - // One-sided three-point estimate for the derivative. + // One-sided, shape-preserving, three-point estimate for the derivative. var d = ((2 * h0 + h1) * m0 - h0 * m1) / (h0 + h1); if (Math.Sign(d) != Math.Sign(m0)) @@ -279,7 +280,6 @@ namespace MathNet.Numerics.Interpolation 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. From e3cf8d12f5f06a05418d6b56561d95b0f6499639 Mon Sep 17 00:00:00 2001 From: Febin <25066330+febkor@users.noreply.github.com> Date: Thu, 31 Dec 2020 17:58:58 +0200 Subject: [PATCH 4/5] Add CubicSplineMonotone to Interpolate. --- src/Numerics/Interpolate.cs | 23 ++++++++++++++++++++++- src/Numerics/Interpolation/CubicSpline.cs | 6 ++++++ 2 files changed, 28 insertions(+), 1 deletion(-) diff --git a/src/Numerics/Interpolate.cs b/src/Numerics/Interpolate.cs index af551b96..7c81c638 100644 --- a/src/Numerics/Interpolate.cs +++ b/src/Numerics/Interpolate.cs @@ -201,7 +201,7 @@ namespace MathNet.Numerics } /// - /// Create an piecewise cubic Akima spline interpolation based on arbitrary points. + /// Create a piecewise cubic Akima spline interpolation based on arbitrary points. /// Akima splines are robust to outliers. /// /// The sample points t. @@ -221,6 +221,27 @@ namespace MathNet.Numerics return Interpolation.CubicSpline.InterpolateAkima(points, values); } + /// + /// Create a piecewise cubic monotone spline interpolation based on arbitrary points. + /// This is a shape-preserving spline with continuous first derivative. + /// + /// The sample points t. + /// The sample point values x(t). + /// + /// An interpolation scheme optimized for the given sample points and values, + /// which can then be used to compute interpolations and extrapolations + /// on arbitrary points. + /// + /// + /// if your data is already sorted in arrays, consider to use + /// MathNet.Numerics.Interpolation.CubicSpline.InterpolatePchipSorted + /// instead, which is more efficient. + /// + public static IInterpolation CubicSplineMonotone(IEnumerable points, IEnumerable values) + { + return Interpolation.CubicSpline.InterpolatePchip(points, values); + } + /// /// Create a piecewise cubic Hermite spline interpolation based on arbitrary points /// and their slopes/first derivative. diff --git a/src/Numerics/Interpolation/CubicSpline.cs b/src/Numerics/Interpolation/CubicSpline.cs index 5dc52146..cba2c54c 100644 --- a/src/Numerics/Interpolation/CubicSpline.cs +++ b/src/Numerics/Interpolation/CubicSpline.cs @@ -244,7 +244,9 @@ namespace MathNet.Numerics.Interpolation 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. @@ -272,10 +274,14 @@ namespace MathNet.Numerics.Interpolation 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; } From ea01dbbd296afbd212522dd691c9603442e6c8e6 Mon Sep 17 00:00:00 2001 From: Febin <25066330+febkor@users.noreply.github.com> Date: Thu, 31 Dec 2020 18:07:47 +0200 Subject: [PATCH 5/5] Add monotone spline to interpolation doc; update site menu to add links. --- docs/content/Interpolation.md | 3 ++- src/Numerics/Interpolation/CubicSpline.cs | 2 ++ 2 files changed, 4 insertions(+), 1 deletion(-) diff --git a/docs/content/Interpolation.md b/docs/content/Interpolation.md index 384e1d71..5f298153 100644 --- a/docs/content/Interpolation.md +++ b/docs/content/Interpolation.md @@ -37,12 +37,13 @@ Interpolation on arbitrary sample points ---------------------------------------- * *Rational pole-free*: Barycentric Floater-Hormann Algorithm -* **Rational with poles**: Bulirsch & Stoer Algorithm +* *Rational with poles*: Bulirsch & Stoer Algorithm * *Neville Polynomial*: Neville Algorithm. Note that the Neville algorithm performs very badly on equidistant points. If you need to interpolate a polynomial on equidistant points, we recommend to use the barycentric algorithm instead. * *Linear Spline* * *Cubic Spline* with boundary conditions * *Natural Cubic Spline* * *Akima Cubic Spline* +* *Monotone Cubic Spline*: Monotone-preserving piecewise cubic Hermite interpolating polynomial (PCHIP), based on Fritsch & Carlson (1980). Interpolation with additional data diff --git a/src/Numerics/Interpolation/CubicSpline.cs b/src/Numerics/Interpolation/CubicSpline.cs index cba2c54c..d2344901 100644 --- a/src/Numerics/Interpolation/CubicSpline.cs +++ b/src/Numerics/Interpolation/CubicSpline.cs @@ -216,6 +216,8 @@ namespace MathNet.Numerics.Interpolation /// public static CubicSpline InterpolatePchipSorted(double[] x, double[] y) { + // Implementation based on "Numerical Computing with Matlab" (Moler, 2004). + if (x.Length != y.Length) { throw new ArgumentException("All vectors must have the same dimensionality.");