diff --git a/src/Managed.UnitTests/InterpolationTests/InterpolationTest.cs b/src/Managed.UnitTests/InterpolationTests/InterpolationTest.cs index 686a50f7..84e2bf64 100644 --- a/src/Managed.UnitTests/InterpolationTests/InterpolationTest.cs +++ b/src/Managed.UnitTests/InterpolationTests/InterpolationTest.cs @@ -89,6 +89,14 @@ namespace MathNet.Numerics.UnitTests.InterpolationTests NonStandardParameters = true, }; + [VerifyContract] + public readonly IContract CubicHermiteSplineContractTests = new InterpolationContract() + { + Factory = (t, x) => new CubicHermiteSplineInterpolation(t, x, CubicSplineInterpolation.EvaluateSplineDerivatives(t, x, SplineBoundaryCondition.Natural, 0.0, SplineBoundaryCondition.Natural, 0.0)), + Order = new[] { 1, 2, 6 }, + NonStandardParameters = true, + }; + [VerifyContract] public readonly IContract LinearSplineContractTests = new InterpolationContract() { @@ -98,5 +106,15 @@ namespace MathNet.Numerics.UnitTests.InterpolationTests PolynomialBehavior = false, RationalBehavior = false }; + + [VerifyContract] + public readonly IContract CubicSplineContractTests = new InterpolationContract() + { + Factory = (t, x) => new CubicSplineInterpolation(t, x), + Order = new[] { 2, 3, 6 }, + LinearBehavior = true, + PolynomialBehavior = false, + RationalBehavior = false + }; } } diff --git a/src/Managed/Interpolation/Algorithms/CubicSplineInterpolation.cs b/src/Managed/Interpolation/Algorithms/CubicSplineInterpolation.cs new file mode 100644 index 00000000..7ec87e22 --- /dev/null +++ b/src/Managed/Interpolation/Algorithms/CubicSplineInterpolation.cs @@ -0,0 +1,383 @@ +// +// Math.NET Numerics, part of the Math.NET Project +// http://mathnet.opensourcedotnet.info +// +// Copyright (c) 2009 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. +// + +namespace MathNet.Numerics.Interpolation.Algorithms +{ + using System; + using System.Collections.Generic; + using Properties; + + /// + /// Cubic Spline Interpolation Algorithm with continuous first and second derivatives. + /// + /// + /// This algorithm supports both differentiation and integration. + /// + public class CubicSplineInterpolation : IInterpolation + { + /// + /// Internal Spline Interpolation + /// + private readonly CubicHermiteSplineInterpolation _spline; + + /// + /// Initializes a new instance of the CubicSplineInterpolation class. + /// + public CubicSplineInterpolation() + { + _spline = new CubicHermiteSplineInterpolation(); + } + + /// + /// Initializes a new instance of the CubicSplineInterpolation class. + /// + /// Sample Points t, sorted ascending. + /// Sample Values x(t) + public CubicSplineInterpolation( + IList samplePoints, + IList sampleValues) + { + _spline = new CubicHermiteSplineInterpolation(); + 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[] derivatives = EvaluateSplineDerivatives( + samplePoints, + sampleValues, + SplineBoundaryCondition.SecondDerivative, + 0.0, + SplineBoundaryCondition.SecondDerivative, + 0.0); + + _spline.Initialize(samplePoints, sampleValues, derivatives); + } + + /// + /// Initialize the interpolation method with the given spline coefficients (sorted by the sample points t). + /// + /// Sample Points t, sorted ascending. + /// Sample Values x(t) + /// Condition of the left boundary. + /// Left boundary value. Ignored in the parabolic case. + /// Condition of the right boundary. + /// Right boundary value. Ignored in the parabolic case. + public void Initialize( + IList samplePoints, + IList sampleValues, + SplineBoundaryCondition leftBoundaryCondition, + double leftBoundary, + SplineBoundaryCondition rightBoundaryCondition, + double rightBoundary) + { + double[] derivatives = EvaluateSplineDerivatives( + samplePoints, + sampleValues, + leftBoundaryCondition, + leftBoundary, + rightBoundaryCondition, + rightBoundary); + + _spline.Initialize(samplePoints, sampleValues, derivatives); + } + + /// + /// Evaluate the spline derivatives as used + /// internally by this interpolation algorithm. + /// + /// Sample Points t, sorted ascending. + /// Sample Values x(t) + /// Condition of the left boundary. + /// Left boundary value. Ignored in the parabolic case. + /// Condition of the right boundary. + /// Right boundary value. Ignored in the parabolic case. + /// Spline Derivative Vector + public static double[] EvaluateSplineDerivatives( + IList samplePoints, + IList sampleValues, + SplineBoundaryCondition leftBoundaryCondition, + double leftBoundary, + SplineBoundaryCondition rightBoundaryCondition, + double rightBoundary) + { + 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.ArgumentVectorsSameLengths); + } + + int n = samplePoints.Count; + + // 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; + } + + double[] a1 = new double[n]; + double[] a2 = new double[n]; + double[] a3 = new double[n]; + double[] b = new double[n]; + + // Left Boundary + switch (leftBoundaryCondition) + { + case SplineBoundaryCondition.ParabolicallyTerminated: + a1[0] = 0; + a2[0] = 1; + a3[0] = 1; + b[0] = 2 * (sampleValues[1] - sampleValues[0]) / (samplePoints[1] - samplePoints[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 * ((sampleValues[1] - sampleValues[0]) / (samplePoints[1] - samplePoints[0]))) - (0.5 * leftBoundary * (samplePoints[1] - samplePoints[0])); + break; + default: + throw new NotSupportedException(Resources.InvalidLeftBoundaryCondition); + } + + // Central Conditions + for (int i = 1; i < samplePoints.Count - 1; i++) + { + a1[i] = samplePoints[i + 1] - samplePoints[i]; + a2[i] = 2 * (samplePoints[i + 1] - samplePoints[i - 1]); + a3[i] = samplePoints[i] - samplePoints[i - 1]; + b[i] = (3 * (sampleValues[i] - sampleValues[i - 1]) / (samplePoints[i] - samplePoints[i - 1]) * (samplePoints[i + 1] - samplePoints[i])) + (3 * (sampleValues[i + 1] - sampleValues[i]) / (samplePoints[i + 1] - samplePoints[i]) * (samplePoints[i] - samplePoints[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 * (sampleValues[n - 1] - sampleValues[n - 2]) / (samplePoints[n - 1] - samplePoints[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 * (sampleValues[n - 1] - sampleValues[n - 2]) / (samplePoints[n - 1] - samplePoints[n - 2])) + (0.5 * rightBoundary * (samplePoints[n - 1] - samplePoints[n - 2])); + break; + default: + throw new NotSupportedException(Resources.InvalidRightBoundaryCondition); + } + + // Build Spline + return SolveTridiagonal(a1, a2, a3, b); + } + + /// + /// Evaluate the spline coefficients as used + /// internally by this interpolation algorithm. + /// + /// Sample Points t, sorted ascending. + /// Sample Values x(t) + /// Condition of the left boundary. + /// Left boundary value. Ignored in the parabolic case. + /// Condition of the right boundary. + /// Right boundary value. Ignored in the parabolic case. + /// Spline Coefficient Vector + public static double[] EvaluateSplineCoefficients( + IList samplePoints, + IList sampleValues, + SplineBoundaryCondition leftBoundaryCondition, + double leftBoundary, + SplineBoundaryCondition rightBoundaryCondition, + double rightBoundary) + { + double[] derivatives = EvaluateSplineDerivatives( + samplePoints, + sampleValues, + leftBoundaryCondition, + leftBoundary, + rightBoundaryCondition, + rightBoundary); + + return CubicHermiteSplineInterpolation.EvaluateSplineCoefficients( + samplePoints, + sampleValues, + derivatives); + } + + /// + /// 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] + private static double[] SolveTridiagonal( + double[] a, + double[] b, + double[] c, + double[] d) + { + double[] x = new double[a.Length]; + + 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]); + } + + 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; + } + + /// + /// 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); + } + + /// + /// Differentiate at point t. + /// + /// Point t to interpolate at. + /// Interpolated value x(t) + /// Interpolated second derivative at point t. + /// Interpolated first derivative at point t. + /// + /// + public double Differentiate( + double t, + out double interpolatedValue, + out double secondDerivative) + { + return _spline.Differentiate(t, out interpolatedValue, out secondDerivative); + } + + /// + /// 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); + } + } +} \ No newline at end of file diff --git a/src/Managed/Interpolation/SplineBoundaryCondition.cs b/src/Managed/Interpolation/SplineBoundaryCondition.cs new file mode 100644 index 00000000..45456d61 --- /dev/null +++ b/src/Managed/Interpolation/SplineBoundaryCondition.cs @@ -0,0 +1,56 @@ +// +// Math.NET Numerics, part of the Math.NET Project +// http://mathnet.opensourcedotnet.info +// +// Copyright (c) 2009 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. +// + +namespace MathNet.Numerics.Interpolation +{ + /// + /// Left and right boundary conditions. + /// + public enum SplineBoundaryCondition + { + /// + /// Natural Boundary (Zero second derivative). + /// + Natural = 0, + + /// + /// Parabolically Terminated boundary. + /// + ParabolicallyTerminated, + + /// + /// Fixed first derivative at the boundary. + /// + FirstDerivative, + + /// + /// Fixed second derivative at the boundary. + /// + SecondDerivative + } +} \ No newline at end of file diff --git a/src/Managed/Managed.csproj b/src/Managed/Managed.csproj index fc479e53..0a75cf59 100644 --- a/src/Managed/Managed.csproj +++ b/src/Managed/Managed.csproj @@ -56,6 +56,7 @@ + @@ -63,6 +64,7 @@ + diff --git a/src/Native/Native.csproj b/src/Native/Native.csproj index b3f9dba6..5f591f69 100644 --- a/src/Native/Native.csproj +++ b/src/Native/Native.csproj @@ -80,6 +80,9 @@ Interpolation\Algorithms\CubicHermiteSplineInterpolation.cs + + Interpolation\Algorithms\CubicSplineInterpolation.cs + Interpolation\Algorithms\FloaterHormannRationalInterpolation.cs @@ -98,6 +101,9 @@ Interpolation\Interpolate.cs + + Interpolation\SplineBoundaryCondition.cs + Precision.cs