Browse Source

Drop old redundant spline interpolation classes

optimization-3
Christoph Ruegg 13 years ago
parent
commit
27866408af
  1. 261
      src/Numerics/Interpolation/AkimaSplineInterpolation.cs
  2. 203
      src/Numerics/Interpolation/CubicHermiteSplineInterpolation.cs
  3. 373
      src/Numerics/Interpolation/CubicSplineInterpolation.cs
  4. 242
      src/Numerics/Interpolation/SplineInterpolation.cs
  5. 4
      src/Numerics/Numerics.csproj

261
src/Numerics/Interpolation/AkimaSplineInterpolation.cs

@ -1,261 +0,0 @@
// <copyright file="AkimaSplineInterpolation.cs" company="Math.NET">
// Math.NET Numerics, part of the Math.NET Project
// http://numerics.mathdotnet.com
// http://github.com/mathnet/mathnet-numerics
// http://mathnetnumerics.codeplex.com
//
// Copyright (c) 2009-2013 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.
// </copyright>
using System;
using System.Collections.Generic;
using MathNet.Numerics.Properties;
namespace MathNet.Numerics.Interpolation
{
/// <summary>
/// Akima Spline Interpolation Algorithm.
/// </summary>
/// <remarks>
/// This algorithm supports both differentiation and integration.
/// </remarks>
public class AkimaSplineInterpolation : IInterpolation
{
/// <summary>
/// Internal Spline Interpolation
/// </summary>
readonly CubicHermiteSplineInterpolation _spline;
/// <summary>
/// Initializes a new instance of the AkimaSplineInterpolation class.
/// </summary>
public AkimaSplineInterpolation()
{
_spline = new CubicHermiteSplineInterpolation();
}
/// <summary>
/// Initializes a new instance of the AkimaSplineInterpolation class.
/// </summary>
/// <param name="samplePoints">Sample Points t, sorted ascending.</param>
/// <param name="sampleValues">Sample Values x(t)</param>
public AkimaSplineInterpolation(IList<double> samplePoints, IList<double> sampleValues)
{
_spline = new CubicHermiteSplineInterpolation();
Initialize(samplePoints, sampleValues);
}
/// <summary>
/// Gets a value indicating whether the algorithm supports differentiation (interpolated derivative).
/// </summary>
bool IInterpolation.SupportsDifferentiation
{
get { return true; }
}
/// <summary>
/// Gets a value indicating whether the algorithm supports integration (interpolated quadrature).
/// </summary>
bool IInterpolation.SupportsIntegration
{
get { return true; }
}
/// <summary>
/// Initialize the interpolation method with the given spline coefficients (sorted by the sample points t).
/// </summary>
/// <param name="samplePoints">Sample Points t, sorted ascending.</param>
/// <param name="sampleValues">Sample Values x(t)</param>
public void Initialize(IList<double> samplePoints, IList<double> sampleValues)
{
double[] derivatives = EvaluateSplineDerivatives(samplePoints, sampleValues);
_spline.Initialize(samplePoints, sampleValues, derivatives);
}
/// <summary>
/// Evaluate the spline derivatives as used
/// internally by this interpolation algorithm.
/// </summary>
/// <param name="samplePoints">Sample Points t, sorted ascending.</param>
/// <param name="sampleValues">Sample Values x(t)</param>
/// <returns>Spline Derivative Vector</returns>
public static double[] EvaluateSplineDerivatives(IList<double> samplePoints, IList<double> 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.ArgumentVectorsSameLength);
}
for (var i = 1; i < samplePoints.Count; ++i)
if (samplePoints[i] <= samplePoints[i - 1])
throw new ArgumentException(Resources.Interpolation_Initialize_SamplePointsNotStrictlyAscendingOrder, "samplePoints");
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].AlmostEqual(0.0) && weights[i + 1].AlmostEqual(0.0)
? (((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;
}
/// <summary>
/// Evaluate the spline coefficients as used
/// internally by this interpolation algorithm.
/// </summary>
/// <param name="samplePoints">Sample Points t, sorted ascending.</param>
/// <param name="sampleValues">Sample Values x(t)</param>
/// <returns>Spline Coefficient Vector</returns>
public static double[] EvaluateSplineCoefficients(IList<double> samplePoints, IList<double> sampleValues)
{
double[] derivatives = EvaluateSplineDerivatives(samplePoints, sampleValues);
return CubicHermiteSplineInterpolation.EvaluateSplineCoefficients(samplePoints, sampleValues, derivatives);
}
/// <summary>
/// Three-Point Differentiation Helper.
/// </summary>
/// <param name="samplePoints">Sample Points t.</param>
/// <param name="sampleValues">Sample Values x(t).</param>
/// <param name="indexT">Index of the point of the differentiation.</param>
/// <param name="index0">Index of the first sample.</param>
/// <param name="index1">Index of the second sample.</param>
/// <param name="index2">Index of the third sample.</param>
/// <returns>The derivative approximation.</returns>
static double DifferentiateThreePoint(
IList<double> samplePoints, IList<double> 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;
}
/// <summary>
/// Interpolate at point t.
/// </summary>
/// <param name="t">Point t to interpolate at.</param>
/// <returns>Interpolated value x(t).</returns>
public double Interpolate(double t)
{
return _spline.Interpolate(t);
}
/// <summary>
/// Differentiate at point t.
/// </summary>
/// <param name="t">Point t to interpolate at.</param>
/// <returns>Interpolated first derivative at point t.</returns>
public double Differentiate(double t)
{
return _spline.Differentiate(t);
}
/// <summary>
/// Differentiate twice at point t.
/// </summary>
/// <param name="t">Point t to interpolate at.</param>
/// <returns>Interpolated second derivative at point t.</returns>
public double Differentiate2(double t)
{
return _spline.Differentiate2(t);
}
/// <summary>
/// Indefinite integral at point t.
/// </summary>
/// <param name="t">Point t to integrate at.</param>
public double Integrate(double t)
{
return _spline.Integrate(t);
}
/// <summary>
/// Definite integral between points a and b.
/// </summary>
/// <param name="a">Left bound of the integration interval [a,b].</param>
/// <param name="b">Right bound of the integration interval [a,b].</param>
public double Integrate(double a, double b)
{
return _spline.Integrate(a, b);
}
}
}

203
src/Numerics/Interpolation/CubicHermiteSplineInterpolation.cs

@ -1,203 +0,0 @@
// <copyright file="CubicHermiteSplineInterpolation.cs" company="Math.NET">
// Math.NET Numerics, part of the Math.NET Project
// http://numerics.mathdotnet.com
// http://github.com/mathnet/mathnet-numerics
// http://mathnetnumerics.codeplex.com
//
// Copyright (c) 2009-2013 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.
// </copyright>
using System;
using System.Collections.Generic;
using MathNet.Numerics.Properties;
namespace MathNet.Numerics.Interpolation
{
/// <summary>
/// Cubic Hermite Spline Interpolation Algorithm.
/// </summary>
/// <remarks>
/// This algorithm supports both differentiation and integration.
/// </remarks>
public class CubicHermiteSplineInterpolation : IInterpolation
{
/// <summary>
/// Internal Spline Interpolation
/// </summary>
readonly SplineInterpolation _spline;
/// <summary>
/// Initializes a new instance of the CubicHermiteSplineInterpolation class.
/// </summary>
public CubicHermiteSplineInterpolation()
{
_spline = new SplineInterpolation();
}
/// <summary>
/// Initializes a new instance of the CubicHermiteSplineInterpolation class.
/// </summary>
/// <param name="samplePoints">Sample Points t, sorted ascending.</param>
/// <param name="sampleValues">Sample Values x(t)</param>
/// <param name="sampleDerivatives">Sample Derivatives x'(t)</param>
public CubicHermiteSplineInterpolation(IList<double> samplePoints, IList<double> sampleValues, IList<double> sampleDerivatives)
{
_spline = new SplineInterpolation();
Initialize(samplePoints, sampleValues, sampleDerivatives);
}
/// <summary>
/// Gets a value indicating whether the algorithm supports differentiation (interpolated derivative).
/// </summary>
bool IInterpolation.SupportsDifferentiation
{
get { return true; }
}
/// <summary>
/// Gets a value indicating whether the algorithm supports integration (interpolated quadrature).
/// </summary>
bool IInterpolation.SupportsIntegration
{
get { return true; }
}
/// <summary>
/// Initialize the interpolation method with the given spline coefficients (sorted by the sample points t).
/// </summary>
/// <param name="samplePoints">Sample Points t, sorted ascending.</param>
/// <param name="sampleValues">Sample Values x(t)</param>
/// <param name="sampleDerivatives">Sample Derivatives x'(t)</param>
public void Initialize(IList<double> samplePoints, IList<double> sampleValues, IList<double> sampleDerivatives)
{
double[] coefficients = EvaluateSplineCoefficients(samplePoints, sampleValues, sampleDerivatives);
_spline.Initialize(samplePoints, coefficients);
}
/// <summary>
/// Evaluate the spline coefficients as used
/// internally by this interpolation algorithm.
/// </summary>
/// <param name="samplePoints">Sample Points t, sorted ascending.</param>
/// <param name="sampleValues">Sample Values x(t)</param>
/// <param name="sampleDerivatives">Sample Derivatives x'(t)</param>
/// <returns>Spline Coefficient Vector</returns>
public static double[] EvaluateSplineCoefficients(IList<double> samplePoints, IList<double> sampleValues, IList<double> sampleDerivatives)
{
if (null == samplePoints)
{
throw new ArgumentNullException("samplePoints");
}
if (null == sampleValues)
{
throw new ArgumentNullException("sampleValues");
}
if (null == sampleDerivatives)
{
throw new ArgumentNullException("sampleDerivatives");
}
if (samplePoints.Count < 2)
{
throw new ArgumentOutOfRangeException("samplePoints");
}
if (samplePoints.Count != sampleValues.Count
|| samplePoints.Count != sampleDerivatives.Count)
{
throw new ArgumentException(Resources.ArgumentVectorsSameLength);
}
for (var i = 1; i < samplePoints.Count; ++i)
if (samplePoints[i] <= samplePoints[i - 1])
throw new ArgumentException(Resources.Interpolation_Initialize_SamplePointsNotStrictlyAscendingOrder, "samplePoints");
var coefficients = new double[4*(samplePoints.Count - 1)];
for (int i = 0, j = 0; i < samplePoints.Count - 1; i++, j += 4)
{
double delta = samplePoints[i + 1] - samplePoints[i];
double delta2 = delta*delta;
double delta3 = delta*delta2;
coefficients[j] = sampleValues[i];
coefficients[j + 1] = sampleDerivatives[i];
coefficients[j + 2] = ((3*(sampleValues[i + 1] - sampleValues[i])) - (2*sampleDerivatives[i]*delta) - (sampleDerivatives[i + 1]*delta))/delta2;
coefficients[j + 3] = ((2*(sampleValues[i] - sampleValues[i + 1])) + (sampleDerivatives[i]*delta) + (sampleDerivatives[i + 1]*delta))/delta3;
}
return coefficients;
}
/// <summary>
/// Interpolate at point t.
/// </summary>
/// <param name="t">Point t to interpolate at.</param>
/// <returns>Interpolated value x(t).</returns>
public double Interpolate(double t)
{
return _spline.Interpolate(t);
}
/// <summary>
/// Differentiate at point t.
/// </summary>
/// <param name="t">Point t to interpolate at.</param>
/// <returns>Interpolated first derivative at point t.</returns>
public double Differentiate(double t)
{
return _spline.Differentiate(t);
}
/// <summary>
/// Differentiate twice at point t.
/// </summary>
/// <param name="t">Point t to interpolate at.</param>
/// <returns>Interpolated second derivative at point t.</returns>
public double Differentiate2(double t)
{
return _spline.Differentiate2(t);
}
/// <summary>
/// Indefinite integral at point t.
/// </summary>
/// <param name="t">Point t to integrate at.</param>
public double Integrate(double t)
{
return _spline.Integrate(t);
}
/// <summary>
/// Definite integral between points a and b.
/// </summary>
/// <param name="a">Left bound of the integration interval [a,b].</param>
/// <param name="b">Right bound of the integration interval [a,b].</param>
public double Integrate(double a, double b)
{
return _spline.Integrate(a, b);
}
}
}

373
src/Numerics/Interpolation/CubicSplineInterpolation.cs

@ -1,373 +0,0 @@
// <copyright file="CubicSplineInterpolation.cs" company="Math.NET">
// Math.NET Numerics, part of the Math.NET Project
// http://numerics.mathdotnet.com
// http://github.com/mathnet/mathnet-numerics
// http://mathnetnumerics.codeplex.com
//
// Copyright (c) 2009-2013 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.
// </copyright>
using System;
using System.Collections.Generic;
using MathNet.Numerics.Properties;
namespace MathNet.Numerics.Interpolation
{
/// <summary>
/// Cubic Spline Interpolation Algorithm with continuous first and second derivatives.
/// </summary>
/// <remarks>
/// This algorithm supports both differentiation and integration.
/// </remarks>
public class CubicSplineInterpolation : IInterpolation
{
/// <summary>
/// Internal Spline Interpolation
/// </summary>
readonly CubicHermiteSplineInterpolation _spline;
/// <summary>
/// Initializes a new instance of the CubicSplineInterpolation class.
/// </summary>
public CubicSplineInterpolation()
{
_spline = new CubicHermiteSplineInterpolation();
}
/// <summary>
/// Initializes a new instance of the CubicSplineInterpolation class.
/// </summary>
/// <param name="samplePoints">Sample Points t, sorted ascending.</param>
/// <param name="sampleValues">Sample Values x(t)</param>
public CubicSplineInterpolation(IList<double> samplePoints, IList<double> sampleValues)
{
_spline = new CubicHermiteSplineInterpolation();
Initialize(samplePoints, sampleValues);
}
/// <summary>
/// Initializes a new instance of the CubicSplineInterpolation class.
/// </summary>
/// <param name="samplePoints">Sample Points t, sorted ascending.</param>
/// <param name="sampleValues">Sample Values x(t)</param>
/// <param name="leftBoundaryCondition">Condition of the left boundary.</param>
/// <param name="leftBoundary">Left boundary value. Ignored in the parabolic case.</param>
/// <param name="rightBoundaryCondition">Condition of the right boundary.</param>
/// <param name="rightBoundary">Right boundary value. Ignored in the parabolic case.</param>
public CubicSplineInterpolation(
IList<double> samplePoints, IList<double> sampleValues,
SplineBoundaryCondition leftBoundaryCondition, double leftBoundary,
SplineBoundaryCondition rightBoundaryCondition, double rightBoundary)
{
_spline = new CubicHermiteSplineInterpolation();
Initialize(
samplePoints, sampleValues,
leftBoundaryCondition, leftBoundary,
rightBoundaryCondition, rightBoundary);
}
/// <summary>
/// Gets a value indicating whether the algorithm supports differentiation (interpolated derivative).
/// </summary>
bool IInterpolation.SupportsDifferentiation
{
get { return true; }
}
/// <summary>
/// Gets a value indicating whether the algorithm supports integration (interpolated quadrature).
/// </summary>
bool IInterpolation.SupportsIntegration
{
get { return true; }
}
/// <summary>
/// Initialize the interpolation method with the given spline coefficients (sorted by the sample points t).
/// </summary>
/// <param name="samplePoints">Sample Points t, sorted ascending.</param>
/// <param name="sampleValues">Sample Values x(t)</param>
public void Initialize(IList<double> samplePoints, IList<double> sampleValues)
{
double[] derivatives = EvaluateSplineDerivatives(samplePoints, sampleValues,
SplineBoundaryCondition.SecondDerivative, 0.0,
SplineBoundaryCondition.SecondDerivative, 0.0);
_spline.Initialize(samplePoints, sampleValues, derivatives);
}
/// <summary>
/// Initialize the interpolation method with the given spline coefficients (sorted by the sample points t).
/// </summary>
/// <param name="samplePoints">Sample Points t, sorted ascending.</param>
/// <param name="sampleValues">Sample Values x(t)</param>
/// <param name="leftBoundaryCondition">Condition of the left boundary.</param>
/// <param name="leftBoundary">Left boundary value. Ignored in the parabolic case.</param>
/// <param name="rightBoundaryCondition">Condition of the right boundary.</param>
/// <param name="rightBoundary">Right boundary value. Ignored in the parabolic case.</param>
public void Initialize(
IList<double> samplePoints, IList<double> sampleValues,
SplineBoundaryCondition leftBoundaryCondition, double leftBoundary,
SplineBoundaryCondition rightBoundaryCondition, double rightBoundary)
{
double[] derivatives = EvaluateSplineDerivatives(
samplePoints, sampleValues,
leftBoundaryCondition, leftBoundary,
rightBoundaryCondition, rightBoundary);
_spline.Initialize(samplePoints, sampleValues, derivatives);
}
/// <summary>
/// Evaluate the spline derivatives as used
/// internally by this interpolation algorithm.
/// </summary>
/// <param name="samplePoints">Sample Points t, sorted ascending.</param>
/// <param name="sampleValues">Sample Values x(t)</param>
/// <param name="leftBoundaryCondition">Condition of the left boundary.</param>
/// <param name="leftBoundary">Left boundary value. Ignored in the parabolic case.</param>
/// <param name="rightBoundaryCondition">Condition of the right boundary.</param>
/// <param name="rightBoundary">Right boundary value. Ignored in the parabolic case.</param>
/// <returns>Spline Derivative Vector</returns>
public static double[] EvaluateSplineDerivatives(
IList<double> samplePoints, IList<double> 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.ArgumentVectorsSameLength);
}
for (var i = 1; i < samplePoints.Count; ++i)
if (samplePoints[i] <= samplePoints[i - 1])
throw new ArgumentException(Resources.Interpolation_Initialize_SamplePointsNotStrictlyAscendingOrder, "samplePoints");
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;
}
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*(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);
}
/// <summary>
/// Evaluate the spline coefficients as used
/// internally by this interpolation algorithm.
/// </summary>
/// <param name="samplePoints">Sample Points t, sorted ascending.</param>
/// <param name="sampleValues">Sample Values x(t)</param>
/// <param name="leftBoundaryCondition">Condition of the left boundary.</param>
/// <param name="leftBoundary">Left boundary value. Ignored in the parabolic case.</param>
/// <param name="rightBoundaryCondition">Condition of the right boundary.</param>
/// <param name="rightBoundary">Right boundary value. Ignored in the parabolic case.</param>
/// <returns>Spline Coefficient Vector</returns>
public static double[] EvaluateSplineCoefficients(
IList<double> samplePoints, IList<double> sampleValues,
SplineBoundaryCondition leftBoundaryCondition, double leftBoundary,
SplineBoundaryCondition rightBoundaryCondition, double rightBoundary)
{
double[] derivatives = EvaluateSplineDerivatives(
samplePoints, sampleValues,
leftBoundaryCondition, leftBoundary,
rightBoundaryCondition, rightBoundary);
return CubicHermiteSplineInterpolation.EvaluateSplineCoefficients(samplePoints, sampleValues, derivatives);
}
/// <summary>
/// Tridiagonal Solve Helper.
/// </summary>
/// <param name="a">The a-vector[n].</param>
/// <param name="b">The b-vector[n], will be modified by this function.</param>
/// <param name="c">The c-vector[n].</param>
/// <param name="d">The d-vector[n], will be modified by this function.</param>
/// <returns>The x-vector[n]</returns>
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;
}
/// <summary>
/// Interpolate at point t.
/// </summary>
/// <param name="t">Point t to interpolate at.</param>
/// <returns>Interpolated value x(t).</returns>
public double Interpolate(double t)
{
return _spline.Interpolate(t);
}
/// <summary>
/// Differentiate at point t.
/// </summary>
/// <param name="t">Point t to interpolate at.</param>
/// <returns>Interpolated first derivative at point t.</returns>
public double Differentiate(double t)
{
return _spline.Differentiate(t);
}
/// <summary>
/// Differentiate twice at point t.
/// </summary>
/// <param name="t">Point t to interpolate at.</param>
/// <returns>Interpolated second derivative at point t.</returns>
public double Differentiate2(double t)
{
return _spline.Differentiate2(t);
}
/// <summary>
/// Indefinite integral at point t.
/// </summary>
/// <param name="t">Point t to integrate at.</param>
public double Integrate(double t)
{
return _spline.Integrate(t);
}
/// <summary>
/// Definite integral between points a and b.
/// </summary>
/// <param name="a">Left bound of the integration interval [a,b].</param>
/// <param name="b">Right bound of the integration interval [a,b].</param>
public double Integrate(double a, double b)
{
return _spline.Integrate(a, b);
}
}
}

242
src/Numerics/Interpolation/SplineInterpolation.cs

@ -1,242 +0,0 @@
// <copyright file="SplineInterpolation.cs" company="Math.NET">
// Math.NET Numerics, part of the Math.NET Project
// http://numerics.mathdotnet.com
// http://github.com/mathnet/mathnet-numerics
// http://mathnetnumerics.codeplex.com
//
// Copyright (c) 2009-2013 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.
// </copyright>
using System;
using System.Collections.Generic;
using MathNet.Numerics.Properties;
namespace MathNet.Numerics.Interpolation
{
/// <summary>
/// Third-Degree Spline Interpolation Algorithm.
/// </summary>
/// <remarks>
/// This algorithm supports both differentiation and integration.
/// </remarks>
public class SplineInterpolation : IInterpolation
{
/// <summary>
/// Sample Points t.
/// </summary>
IList<double> _points;
/// <summary>
/// Spline Coefficients c(t).
/// </summary>
IList<double> _coefficients;
/// <summary>
/// Number of samples.
/// </summary>
int _sampleCount;
/// <summary>
/// Initializes a new instance of the SplineInterpolation class.
/// </summary>
public SplineInterpolation()
{
}
/// <summary>
/// Initializes a new instance of the SplineInterpolation class.
/// </summary>
/// <param name="samplePoints">Sample Points t (length: N), sorted ascending.</param>
/// <param name="splineCoefficients">Spline Coefficients (length: 4*(N-1)).</param>
public SplineInterpolation(IList<double> samplePoints, IList<double> splineCoefficients)
{
Initialize(samplePoints, splineCoefficients);
}
/// <summary>
/// Gets a value indicating whether the algorithm supports differentiation (interpolated derivative).
/// </summary>
bool IInterpolation.SupportsDifferentiation
{
get { return true; }
}
/// <summary>
/// Gets a value indicating whether the algorithm supports integration (interpolated quadrature).
/// </summary>
bool IInterpolation.SupportsIntegration
{
get { return true; }
}
/// <summary>
/// Initialize the interpolation method with the given spline coefficients (sorted by the sample points t).
/// </summary>
/// <param name="samplePoints">Sample Points t (length: N), sorted ascending.</param>
/// <param name="splineCoefficients">Spline Coefficients (length: 4*(N-1)).</param>
public void Initialize(IList<double> samplePoints, IList<double> splineCoefficients)
{
if (null == samplePoints)
{
throw new ArgumentNullException("samplePoints");
}
if (null == splineCoefficients)
{
throw new ArgumentNullException("splineCoefficients");
}
if (samplePoints.Count < 2)
{
throw new ArgumentOutOfRangeException("samplePoints");
}
if (splineCoefficients.Count != 4*(samplePoints.Count - 1))
{
throw new ArgumentOutOfRangeException("splineCoefficients");
}
for (var i = 1; i < samplePoints.Count; ++i)
if (samplePoints[i] <= samplePoints[i - 1])
throw new ArgumentException(Resources.Interpolation_Initialize_SamplePointsNotStrictlyAscendingOrder, "samplePoints");
_points = samplePoints;
_coefficients = splineCoefficients;
_sampleCount = samplePoints.Count;
}
/// <summary>
/// Interpolate at point t.
/// </summary>
/// <param name="t">Point t to interpolate at.</param>
/// <returns>Interpolated value x(t).</returns>
public double Interpolate(double t)
{
int closestLeftIndex = LeftBracketIndex(t);
// Interpolation
double offset = t - _points[closestLeftIndex];
int k = closestLeftIndex << 2;
return _coefficients[k]
+ (offset*(_coefficients[k + 1]
+ (offset*(_coefficients[k + 2]
+ (offset*_coefficients[k + 3])))));
}
/// <summary>
/// Differentiate at point t.
/// </summary>
/// <param name="t">Point t to interpolate at.</param>
/// <returns>Interpolated first derivative at point t.</returns>
public double Differentiate(double t)
{
int closestLeftIndex = LeftBracketIndex(t);
double offset = t - _points[closestLeftIndex];
int k = closestLeftIndex << 2;
return _coefficients[k + 1]
+ (2*offset*_coefficients[k + 2])
+ (3*offset*offset*_coefficients[k + 3]);
}
/// <summary>
/// Differentiate twice at point t.
/// </summary>
/// <param name="t">Point t to interpolate at.</param>
/// <returns>Interpolated second derivative at point t.</returns>
public double Differentiate2(double t)
{
int closestLeftIndex = LeftBracketIndex(t);
double offset = t - _points[closestLeftIndex];
int k = closestLeftIndex << 2;
return (2*_coefficients[k + 2]) + (6*offset*_coefficients[k + 3]);
}
/// <summary>
/// Indefinite integral at point t.
/// </summary>
/// <param name="t">Point t to integrate at.</param>
public double Integrate(double t)
{
int closestLeftIndex = LeftBracketIndex(t);
// Integration
double result = 0;
for (int i = 0, j = 0; i < closestLeftIndex; i++, j += 4)
{
double w = _points[i + 1] - _points[i];
result += w*(_coefficients[j]
+ ((w*_coefficients[j + 1]*0.5)
+ (w*((_coefficients[j + 2]/3)
+ (w*_coefficients[j + 3]*0.25)))));
}
double offset = t - _points[closestLeftIndex];
int k = closestLeftIndex << 2;
return result + (offset*(_coefficients[k]
+ (offset*_coefficients[k + 1]*0.5)
+ (offset*_coefficients[k + 2]/3)
+ (offset*_coefficients[k + 3]*0.25)));
}
/// <summary>
/// Definite integral between points a and b.
/// </summary>
/// <param name="a">Left bound of the integration interval [a,b].</param>
/// <param name="b">Right bound of the integration interval [a,b].</param>
public double Integrate(double a, double b)
{
return Integrate(b) - Integrate(a);
}
/// <summary>
/// Find the index of the greatest sample point smaller than t.
/// </summary>
/// <param name="t">The value to look for.</param>
/// <returns>The sample point index.</returns>
int LeftBracketIndex(double t)
{
// Binary search in the [ t[0], ..., t[n-2] ] (t[n-1] is not included)
int low = 0;
int high = _sampleCount - 1;
while (low != high - 1)
{
int middle = (low + high)/2;
if (_points[middle] > t)
{
high = middle;
}
else
{
low = middle;
}
}
return low;
}
}
}

4
src/Numerics/Numerics.csproj

@ -379,15 +379,11 @@
<Compile Include="Integration\SimpsonRule.cs" />
<Compile Include="Integration\NewtonCotesTrapeziumRule.cs" />
<Compile Include="Integrate.cs" />
<Compile Include="Interpolation\AkimaSplineInterpolation.cs" />
<Compile Include="Interpolation\BarycentricInterpolation.cs" />
<Compile Include="Interpolation\BulirschStoerRationalInterpolation.cs" />
<Compile Include="Interpolation\CubicSplineInterpolation.cs" />
<Compile Include="Interpolation\CubicHermiteSplineInterpolation.cs" />
<Compile Include="Interpolation\FloaterHormannRationalInterpolation.cs" />
<Compile Include="Interpolation\LinearSpline.cs" />
<Compile Include="Interpolation\NevillePolynomialInterpolation.cs" />
<Compile Include="Interpolation\SplineInterpolation.cs" />
<Compile Include="Interpolation\IInterpolation.cs" />
<Compile Include="Interpolate.cs" />
<Compile Include="Interpolation\SplineBoundaryCondition.cs" />

Loading…
Cancel
Save