From 94b5c448e69bf1bd0313ff7e64ab6cc0561f6b5d Mon Sep 17 00:00:00 2001 From: Christoph Ruegg Date: Sat, 21 Dec 2013 01:30:41 +0100 Subject: [PATCH] Interpolation: migrate cubic spline to CubicSpline class --- src/Numerics/Interpolation/CubicSpline.cs | 155 ++++++++++++++++++ .../InterpolationTests/CubicSplineTest.cs | 52 ++---- 2 files changed, 173 insertions(+), 34 deletions(-) diff --git a/src/Numerics/Interpolation/CubicSpline.cs b/src/Numerics/Interpolation/CubicSpline.cs index 9dba29a9..cd4183cb 100644 --- a/src/Numerics/Interpolation/CubicSpline.cs +++ b/src/Numerics/Interpolation/CubicSpline.cs @@ -86,6 +86,11 @@ namespace MathNet.Numerics.Interpolation throw new ArgumentException(Resources.ArgumentVectorsSameLength); } + if (xx.Length < 2) + { + throw new ArgumentOutOfRangeException("x"); + } + Sorting.Sort(xx, yy, dd); var c0 = new double[xx.Length - 1]; @@ -122,6 +127,11 @@ namespace MathNet.Numerics.Interpolation throw new ArgumentException(Resources.ArgumentVectorsSameLength); } + if (xx.Length < 5) + { + throw new ArgumentOutOfRangeException("x"); + } + Sorting.Sort(xx, yy); /* Prepare divided differences (diff) and weights (w) */ @@ -160,6 +170,124 @@ namespace MathNet.Numerics.Interpolation return InterpolateHermite(xx, yy, dd); } + public static CubicSpline Interpolate(IEnumerable x, IEnumerable y) + { + return InterpolateBoundaries(x, y, SplineBoundaryCondition.SecondDerivative, 0.0, SplineBoundaryCondition.SecondDerivative, 0.0); + } + + public static CubicSpline InterpolateBoundaries(IEnumerable x, IEnumerable y, + SplineBoundaryCondition leftBoundaryCondition, double leftBoundary, + SplineBoundaryCondition rightBoundaryCondition, double rightBoundary) + { + var xx = (x as double[]) ?? x.ToArray(); + var yy = (y as double[]) ?? y.ToArray(); + + if (xx.Length != yy.Length) + { + throw new ArgumentException(Resources.ArgumentVectorsSameLength); + } + + if (xx.Length < 2) + { + throw new ArgumentOutOfRangeException("x"); + } + + Sorting.Sort(xx, yy); + + int n = xx.Length; + + // normalize special cases + if ((n == 2) + && (leftBoundaryCondition == SplineBoundaryCondition.ParabolicallyTerminated) + && (rightBoundaryCondition == SplineBoundaryCondition.ParabolicallyTerminated)) + { + leftBoundaryCondition = SplineBoundaryCondition.SecondDerivative; + leftBoundary = 0d; + rightBoundaryCondition = SplineBoundaryCondition.SecondDerivative; + rightBoundary = 0d; + } + + if (leftBoundaryCondition == SplineBoundaryCondition.Natural) + { + leftBoundaryCondition = SplineBoundaryCondition.SecondDerivative; + leftBoundary = 0d; + } + + if (rightBoundaryCondition == SplineBoundaryCondition.Natural) + { + rightBoundaryCondition = SplineBoundaryCondition.SecondDerivative; + rightBoundary = 0d; + } + + var a1 = new double[n]; + var a2 = new double[n]; + var a3 = new double[n]; + var b = new double[n]; + + // Left Boundary + switch (leftBoundaryCondition) + { + case SplineBoundaryCondition.ParabolicallyTerminated: + a1[0] = 0; + a2[0] = 1; + a3[0] = 1; + b[0] = 2*(yy[1] - yy[0])/(xx[1] - xx[0]); + break; + case SplineBoundaryCondition.FirstDerivative: + a1[0] = 0; + a2[0] = 1; + a3[0] = 0; + b[0] = leftBoundary; + break; + case SplineBoundaryCondition.SecondDerivative: + a1[0] = 0; + a2[0] = 2; + a3[0] = 1; + b[0] = (3*((yy[1] - yy[0])/(xx[1] - xx[0]))) - (0.5*leftBoundary*(xx[1] - xx[0])); + break; + default: + throw new NotSupportedException(Resources.InvalidLeftBoundaryCondition); + } + + // Central Conditions + for (int i = 1; i < xx.Length - 1; i++) + { + a1[i] = xx[i + 1] - xx[i]; + a2[i] = 2*(xx[i + 1] - xx[i - 1]); + a3[i] = xx[i] - xx[i - 1]; + b[i] = (3*(yy[i] - yy[i - 1])/(xx[i] - xx[i - 1])*(xx[i + 1] - xx[i])) + (3*(yy[i + 1] - yy[i])/(xx[i + 1] - xx[i])*(xx[i] - xx[i - 1])); + } + + // Right Boundary + switch (rightBoundaryCondition) + { + case SplineBoundaryCondition.ParabolicallyTerminated: + a1[n - 1] = 1; + a2[n - 1] = 1; + a3[n - 1] = 0; + b[n - 1] = 2*(yy[n - 1] - yy[n - 2])/(xx[n - 1] - xx[n - 2]); + break; + case SplineBoundaryCondition.FirstDerivative: + a1[n - 1] = 0; + a2[n - 1] = 1; + a3[n - 1] = 0; + b[n - 1] = rightBoundary; + break; + case SplineBoundaryCondition.SecondDerivative: + a1[n - 1] = 1; + a2[n - 1] = 2; + a3[n - 1] = 0; + b[n - 1] = (3*(yy[n - 1] - yy[n - 2])/(xx[n - 1] - xx[n - 2])) + (0.5*rightBoundary*(xx[n - 1] - xx[n - 2])); + break; + default: + throw new NotSupportedException(Resources.InvalidRightBoundaryCondition); + } + + // Build Spline + double[] dd = SolveTridiagonal(a1, a2, a3, b); + return InterpolateHermite(xx, yy, dd); + } + /// /// Three-Point Differentiation Helper. /// @@ -185,6 +313,33 @@ namespace MathNet.Numerics.Interpolation return (2*a*t) + b; } + /// + /// Tridiagonal Solve Helper. + /// + /// The a-vector[n]. + /// The b-vector[n], will be modified by this function. + /// The c-vector[n]. + /// The d-vector[n], will be modified by this function. + /// The x-vector[n] + static double[] SolveTridiagonal(double[] a, double[] b, double[] c, double[] d) + { + for (int k = 1; k < a.Length; k++) + { + double t = a[k]/b[k - 1]; + b[k] = b[k] - (t*c[k - 1]); + d[k] = d[k] - (t*d[k - 1]); + } + + var x = new double[a.Length]; + x[x.Length - 1] = d[d.Length - 1]/b[b.Length - 1]; + for (int k = x.Length - 2; k >= 0; k--) + { + x[k] = (d[k] - (c[k]*x[k + 1]))/b[k]; + } + + return x; + } + /// /// Gets a value indicating whether the algorithm supports differentiation (interpolated derivative). /// diff --git a/src/UnitTests/InterpolationTests/CubicSplineTest.cs b/src/UnitTests/InterpolationTests/CubicSplineTest.cs index 4e1c27af..4185b5d3 100644 --- a/src/UnitTests/InterpolationTests/CubicSplineTest.cs +++ b/src/UnitTests/InterpolationTests/CubicSplineTest.cs @@ -33,21 +33,11 @@ using NUnit.Framework; namespace MathNet.Numerics.UnitTests.InterpolationTests { - /// - /// CubicSpline Test case. - /// [TestFixture, Category("Interpolation")] public class CubicSplineTest { - /// - /// Sample points. - /// readonly double[] _t = { -2.0, -1.0, 0.0, 1.0, 2.0 }; - - /// - /// Sample values. - /// - readonly double[] _x = { 1.0, 2.0, -1.0, 0.0, 1.0 }; + readonly double[] _y = { 1.0, 2.0, -1.0, 0.0, 1.0 }; /// /// Verifies that the interpolation matches the given value at all the provided sample points. @@ -55,11 +45,10 @@ namespace MathNet.Numerics.UnitTests.InterpolationTests [Test] public void NaturalFitsAtSamplePoints() { - IInterpolation interpolation = new CubicSplineInterpolation(_t, _x); - - for (int i = 0; i < _x.Length; i++) + IInterpolation it = CubicSpline.Interpolate(_t, _y); + for (int i = 0; i < _y.Length; i++) { - Assert.AreEqual(_x[i], interpolation.Interpolate(_t[i]), "A Exact Point " + i); + Assert.AreEqual(_y[i], it.Interpolate(_t[i]), "A Exact Point " + i); } } @@ -85,9 +74,8 @@ namespace MathNet.Numerics.UnitTests.InterpolationTests [TestCase(-10.0, 677, 1e-12)] public void NaturalFitsAtArbitraryPointsWithMaple(double t, double x, double maxAbsoluteError) { - IInterpolation interpolation = new CubicSplineInterpolation(_t, _x); - - Assert.AreEqual(x, interpolation.Interpolate(t), maxAbsoluteError, "Interpolation at {0}", t); + IInterpolation it = CubicSpline.Interpolate(_t, _y); + Assert.AreEqual(x, it.Interpolate(t), maxAbsoluteError, "Interpolation at {0}", t); } /// @@ -96,11 +84,10 @@ namespace MathNet.Numerics.UnitTests.InterpolationTests [Test] public void FixedFirstDerivativeFitsAtSamplePoints() { - IInterpolation interpolation = new CubicSplineInterpolation(_t, _x, SplineBoundaryCondition.FirstDerivative, 1.0, SplineBoundaryCondition.FirstDerivative, -1.0); - - for (int i = 0; i < _x.Length; i++) + IInterpolation it = CubicSpline.InterpolateBoundaries(_t, _y, SplineBoundaryCondition.FirstDerivative, 1.0, SplineBoundaryCondition.FirstDerivative, -1.0); + for (int i = 0; i < _y.Length; i++) { - Assert.AreEqual(_x[i], interpolation.Interpolate(_t[i]), "A Exact Point " + i); + Assert.AreEqual(_y[i], it.Interpolate(_t[i]), "A Exact Point " + i); } } @@ -126,9 +113,8 @@ namespace MathNet.Numerics.UnitTests.InterpolationTests [TestCase(-10.0, 1330.1428571428571429, 1e-12)] public void FixedFirstDerivativeFitsAtArbitraryPointsWithMaple(double t, double x, double maxAbsoluteError) { - IInterpolation interpolation = new CubicSplineInterpolation(_t, _x, SplineBoundaryCondition.FirstDerivative, 1.0, SplineBoundaryCondition.FirstDerivative, -1.0); - - Assert.AreEqual(x, interpolation.Interpolate(t), maxAbsoluteError, "Interpolation at {0}", t); + IInterpolation it = CubicSpline.InterpolateBoundaries(_t, _y, SplineBoundaryCondition.FirstDerivative, 1.0, SplineBoundaryCondition.FirstDerivative, -1.0); + Assert.AreEqual(x, it.Interpolate(t), maxAbsoluteError, "Interpolation at {0}", t); } /// @@ -137,11 +123,10 @@ namespace MathNet.Numerics.UnitTests.InterpolationTests [Test] public void FixedSecondDerivativeFitsAtSamplePoints() { - IInterpolation interpolation = new CubicSplineInterpolation(_t, _x, SplineBoundaryCondition.SecondDerivative, -5.0, SplineBoundaryCondition.SecondDerivative, -1.0); - - for (int i = 0; i < _x.Length; i++) + IInterpolation it = CubicSpline.InterpolateBoundaries(_t, _y, SplineBoundaryCondition.SecondDerivative, -5.0, SplineBoundaryCondition.SecondDerivative, -1.0); + for (int i = 0; i < _y.Length; i++) { - Assert.AreEqual(_x[i], interpolation.Interpolate(_t[i]), "A Exact Point " + i); + Assert.AreEqual(_y[i], it.Interpolate(_t[i]), "A Exact Point " + i); } } @@ -167,9 +152,8 @@ namespace MathNet.Numerics.UnitTests.InterpolationTests [TestCase(-10.0, -37, 1e-12)] public void FixedSecondDerivativeFitsAtArbitraryPointsWithMaple(double t, double x, double maxAbsoluteError) { - IInterpolation interpolation = new CubicSplineInterpolation(_t, _x, SplineBoundaryCondition.SecondDerivative, -5.0, SplineBoundaryCondition.SecondDerivative, -1.0); - - Assert.AreEqual(x, interpolation.Interpolate(t), maxAbsoluteError, "Interpolation at {0}", t); + IInterpolation it = CubicSpline.InterpolateBoundaries(_t, _y, SplineBoundaryCondition.SecondDerivative, -5.0, SplineBoundaryCondition.SecondDerivative, -1.0); + Assert.AreEqual(x, it.Interpolate(t), maxAbsoluteError, "Interpolation at {0}", t); } /// @@ -183,10 +167,10 @@ namespace MathNet.Numerics.UnitTests.InterpolationTests { double[] x, y, xtest, ytest; LinearInterpolationCase.Build(out x, out y, out xtest, out ytest, samples); - IInterpolation interpolation = new CubicSplineInterpolation(x, y); + IInterpolation it = CubicSpline.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); + Assert.AreEqual(ytest[i], it.Interpolate(xtest[i]), 1e-15, "Linear with {0} samples, sample {1}", samples, i); } } }