From 3e60249dfd5d8c22c869506c6e4e5e9a10daaff8 Mon Sep 17 00:00:00 2001 From: Christoph Ruegg Date: Sat, 1 Aug 2009 05:16:05 +0800 Subject: [PATCH] interpolation: ported akima spline Signed-off-by: Christoph Ruegg --- .../InterpolationTests/InterpolationTest.cs | 10 + .../Algorithms/AkimaSplineInterpolation.cs | 280 ++++++++++++++++++ src/Managed/Managed.csproj | 1 + src/Native/Native.csproj | 3 + 4 files changed, 294 insertions(+) create mode 100644 src/Managed/Interpolation/Algorithms/AkimaSplineInterpolation.cs diff --git a/src/Managed.UnitTests/InterpolationTests/InterpolationTest.cs b/src/Managed.UnitTests/InterpolationTests/InterpolationTest.cs index 28667e62..1c733696 100644 --- a/src/Managed.UnitTests/InterpolationTests/InterpolationTest.cs +++ b/src/Managed.UnitTests/InterpolationTests/InterpolationTest.cs @@ -116,5 +116,15 @@ namespace MathNet.Numerics.UnitTests.InterpolationTests PolynomialBehavior = false, RationalBehavior = false }; + + [VerifyContract] + public readonly IContract AkimaSplineContractTests = new InterpolationContract() + { + Factory = (t, x) => new AkimaSplineInterpolation(t, x), + MinimumSampleCount = 5, + LinearBehavior = true, + PolynomialBehavior = false, + RationalBehavior = false + }; } } diff --git a/src/Managed/Interpolation/Algorithms/AkimaSplineInterpolation.cs b/src/Managed/Interpolation/Algorithms/AkimaSplineInterpolation.cs new file mode 100644 index 00000000..65318bd6 --- /dev/null +++ b/src/Managed/Interpolation/Algorithms/AkimaSplineInterpolation.cs @@ -0,0 +1,280 @@ +// +// 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; + + /// + /// Akima Spline Interpolation Algorithm. + /// + /// + /// This algorithm supports both differentiation and integration. + /// + public class AkimaSplineInterpolation : IInterpolation + { + /// + /// Internal Spline Interpolation + /// + private readonly CubicHermiteSplineInterpolation _spline; + + /// + /// Initializes a new instance of the AkimaSplineInterpolation class. + /// + public AkimaSplineInterpolation() + { + _spline = new CubicHermiteSplineInterpolation(); + } + + /// + /// Initializes a new instance of the AkimaSplineInterpolation class. + /// + /// Sample Points t, sorted ascending. + /// Sample Values x(t) + public AkimaSplineInterpolation( + 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); + + _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) + /// Spline Derivative Vector + public static double[] EvaluateSplineDerivatives( + IList samplePoints, + IList sampleValues) + { + if (null == samplePoints) + { + throw new ArgumentNullException("samplePoints"); + } + + if (null == sampleValues) + { + throw new ArgumentNullException("sampleValues"); + } + + if (samplePoints.Count < 5) + { + throw new ArgumentOutOfRangeException("samplePoints"); + } + + if (samplePoints.Count != sampleValues.Count) + { + throw new ArgumentException(Resources.ArgumentVectorsSameLengths); + } + + int n = samplePoints.Count; + + /* Prepare divided differences (diff) and weights (w) */ + + var differences = new double[n - 1]; + var weights = new double[n - 1]; + + for (int i = 0; i < differences.Length; i++) + { + differences[i] = (sampleValues[i + 1] - sampleValues[i]) / (samplePoints[i + 1] - samplePoints[i]); + } + + for (int i = 1; i < weights.Length; i++) + { + weights[i] = Math.Abs(differences[i] - differences[i - 1]); + } + + /* Prepare Hermite interpolation scheme */ + + var derivatives = new double[n]; + + for (int i = 2; i < derivatives.Length - 2; i++) + { + derivatives[i] = + weights[i - 1].AlmostZero() && weights[i + 1].AlmostZero() + ? (((samplePoints[i + 1] - samplePoints[i]) * differences[i - 1]) + + ((samplePoints[i] - samplePoints[i - 1]) * differences[i])) + / (samplePoints[i + 1] - samplePoints[i - 1]) + : ((weights[i + 1] * differences[i - 1]) + + (weights[i - 1] * differences[i])) + / (weights[i + 1] + weights[i - 1]); + } + + derivatives[0] = DifferentiateThreePoint(samplePoints, sampleValues, 0, 0, 1, 2); + derivatives[1] = DifferentiateThreePoint(samplePoints, sampleValues, 1, 0, 1, 2); + derivatives[n - 2] = DifferentiateThreePoint(samplePoints, sampleValues, n - 2, n - 3, n - 2, n - 1); + derivatives[n - 1] = DifferentiateThreePoint(samplePoints, sampleValues, n - 1, n - 3, n - 2, n - 1); + + /* Build Akima spline using Hermite interpolation scheme */ + + return derivatives; + } + + /// + /// 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) + { + double[] derivatives = EvaluateSplineDerivatives( + samplePoints, + sampleValues); + + return CubicHermiteSplineInterpolation.EvaluateSplineCoefficients( + samplePoints, + sampleValues, + derivatives); + } + + /// + /// Three-Point Differentiation Helper. + /// + /// Sample Points t. + /// Sample Values x(t). + /// Index of the point of the differentiation. + /// Index of the first sample. + /// Index of the second sample. + /// Index of the third sample. + /// The derivative approximation. + private static double DifferentiateThreePoint( + IList samplePoints, + IList sampleValues, + int indexT, + int index0, + int index1, + int index2) + { + double x0 = sampleValues[index0]; + double x1 = sampleValues[index1]; + double x2 = sampleValues[index2]; + + double t = samplePoints[indexT] - samplePoints[index0]; + double t1 = samplePoints[index1] - samplePoints[index0]; + double t2 = samplePoints[index2] - samplePoints[index0]; + + double a = (x2 - x0 - (t2 / t1 * (x1 - x0))) / ((t2 * t2) - (t1 * t2)); + double b = (x1 - x0 - (a * t1 * t1)) / t1; + return (2 * a * t) + b; + } + + /// + /// 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/Managed.csproj b/src/Managed/Managed.csproj index ed96e700..d392c239 100644 --- a/src/Managed/Managed.csproj +++ b/src/Managed/Managed.csproj @@ -58,6 +58,7 @@ + diff --git a/src/Native/Native.csproj b/src/Native/Native.csproj index 10948882..81515625 100644 --- a/src/Native/Native.csproj +++ b/src/Native/Native.csproj @@ -83,6 +83,9 @@ Integration\Integrate.cs + + Interpolation\Algorithms\AkimaSplineInterpolation.cs + Interpolation\Algorithms\BarycentricInterpolation.cs