us1zcljz
/
mathnet-numerics
forked from tsai/mathnet-numerics
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity
Browse Source

Interpolation: common api should not sort original data; alternative if data is already sorted #219

provider
Christoph Ruegg 13 years ago
parent
dd2df946b7
commit
8cf834e137
8 changed files with 474 additions and 204 deletions
Whitespace
Show all changes Ignore whitespace when comparing lines Ignore changes in amount of whitespace Ignore changes in whitespace at EOL
Split View
Diff Options
Show Stats Download Patch File Download Diff File
  1. 100
      src/Numerics/Interpolate.cs
  2. 151
      src/Numerics/Interpolation/Barycentric.cs
  3. 54
      src/Numerics/Interpolation/BulirschStoerRationalInterpolation.cs
  4. 239
      src/Numerics/Interpolation/CubicSpline.cs
  5. 47
      src/Numerics/Interpolation/LinearSpline.cs
  6. 60
      src/Numerics/Interpolation/NevillePolynomialInterpolation.cs
  7. 4
      src/UnitTests/InterpolationTests/BulirschStoerRationalTest.cs
  8. 23
      src/UnitTests/InterpolationTests/CubicSplineTest.cs

100
src/Numerics/Interpolate.cs
View File

@ -4,7 +4,7 @@
// http://github.com/mathnet/mathnet-numerics
// http://mathnetnumerics.codeplex.com
//
// Copyright (c) 2009-2013 Math.NET
// Copyright (c) 2009-2014 Math.NET
//
// Permission is hereby granted, free of charge, to any person
// obtaining a copy of this software and associated documentation
@ -41,13 +41,18 @@ namespace MathNet.Numerics
/// <summary>
/// Creates an interpolation based on arbitrary points.
/// </summary>
/// <param name="points">The sample points t. Supports both lists and arrays.</param>
/// <param name="values">The sample point values x(t). Supports both lists and arrays.</param>
/// <param name="points">The sample points t.</param>
/// <param name="values">The sample point values x(t).</param>
/// <returns>
/// An interpolation scheme optimized for the given sample points and values,
/// which can then be used to compute interpolations and extrapolations
/// on arbitrary points.
/// </returns>
/// <remarks>
/// if your data is already sorted in arrays, consider to use
/// MathNet.Numerics.Interpolation.Barycentric.InterpolateRationalFloaterHormannSorted
/// instead, which is more efficient.
/// </remarks>
public static IInterpolation Common(IEnumerable<double> points, IEnumerable<double> values)
{
return Barycentric.InterpolateRationalFloaterHormann(points, values);
@ -56,46 +61,57 @@ namespace MathNet.Numerics
/// <summary>
/// Create a floater hormann rational pole-free interpolation based on arbitrary points.
/// </summary>
/// <param name="points">The sample points t. Supports both lists and arrays.</param>
/// <param name="values">The sample point values x(t). Supports both lists and arrays.</param>
/// <param name="points">The sample points t.</param>
/// <param name="values">The sample point values x(t).</param>
/// <returns>
/// An interpolation scheme optimized for the given sample points and values,
/// which can then be used to compute interpolations and extrapolations
/// on arbitrary points.
/// </returns>
/// <remarks>
/// if your data is already sorted in arrays, consider to use
/// MathNet.Numerics.Interpolation.Barycentric.InterpolateRationalFloaterHormannSorted
/// instead, which is more efficient.
/// </remarks>
public static IInterpolation RationalWithoutPoles(IEnumerable<double> points, IEnumerable<double> values)
{
return Barycentric.InterpolateRationalFloaterHormann(points, values);
}
/// <summary>
/// Create a burlisch stoer rational interpolation based on arbitrary points.
/// Create a Bulirsch Stoer rational interpolation based on arbitrary points.
/// </summary>
/// <param name="points">The sample points t. Optimized for arrays..</param>
/// <param name="values">The sample point values x(t). Optimized for arrays.</param>
/// <param name="points">The sample points t.</param>
/// <param name="values">The sample point values x(t).</param>
/// <returns>
/// An interpolation scheme optimized for the given sample points and values,
/// which can then be used to compute interpolations and extrapolations
/// on arbitrary points.
/// </returns>
/// <remarks>
/// if your data is already sorted in arrays, consider to use
/// MathNet.Numerics.Interpolation.BulirschStoerRationalInterpolation.InterpolateSorted
/// instead, which is more efficient.
/// </remarks>
public static IInterpolation RationalWithPoles(IEnumerable<double> points, IEnumerable<double> values)
{
return new BulirschStoerRationalInterpolation(points, values);
return BulirschStoerRationalInterpolation.Interpolate(points, values);
}
/// <summary>
/// Create a barycentric polynomial interpolation where the given sample points are equidistant.
/// </summary>
/// <param name="points">The sample points t, must be equidistant. Optimized for arrays.</param>
/// <param name="values">The sample point values x(t). Optimized for arrays.</param>
/// <param name="points">The sample points t, must be equidistant.</param>
/// <param name="values">The sample point values x(t).</param>
/// <returns>
/// An interpolation scheme optimized for the given sample points and values,
/// which can then be used to compute interpolations and extrapolations
/// on arbitrary points.
/// </returns>
/// <remarks>
/// The value pairs do not have to be sorted, but if they are not sorted ascendingly
/// and the passed x and y arguments are arrays, they will be sorted inplace and thus modified.
/// if your data is already sorted in arrays, consider to use
/// MathNet.Numerics.Interpolation.Barycentric.InterpolatePolynomialEquidistantSorted
/// instead, which is more efficient.
/// </remarks>
public static IInterpolation PolynomialEquidistant(IEnumerable<double> points, IEnumerable<double> values)
{
@ -103,35 +119,41 @@ namespace MathNet.Numerics
}
/// <summary>
/// Create a neville polynomial interpolation based on arbitrary points.
/// Create a Neville polynomial interpolation based on arbitrary points.
/// If the points happen to be equidistant, consider to use the much more robust PolynomialEquidistant instead.
/// Otherwise, consider whether RationalWithoutPoles would not be a more robust alternative.
/// </summary>
/// <param name="points">The sample points t. Optimized for arrays.</param>
/// <param name="values">The sample point values x(t). Optimized for arrays.</param>
/// <param name="points">The sample points t.</param>
/// <param name="values">The sample point values x(t).</param>
/// <returns>
/// An interpolation scheme optimized for the given sample points and values,
/// which can then be used to compute interpolations and extrapolations
/// on arbitrary points.
/// </returns>
/// <remarks>
/// if your data is already sorted in arrays, consider to use
/// MathNet.Numerics.Interpolation.NevillePolynomialInterpolation.InterpolateSorted
/// instead, which is more efficient.
/// </remarks>
public static IInterpolation Polynomial(IEnumerable<double> points, IEnumerable<double> values)
{
return new NevillePolynomialInterpolation(points, values);
return NevillePolynomialInterpolation.Interpolate(points, values);
}
/// <summary>
/// Create a piecewise linear spline interpolation based on arbitrary points.
/// </summary>
/// <param name="points">The sample points t. Optimized for arrays.</param>
/// <param name="values">The sample point values x(t). Optimized for arrays.</param>
/// <param name="points">The sample points t.</param>
/// <param name="values">The sample point values x(t).</param>
/// <returns>
/// An interpolation scheme optimized for the given sample points and values,
/// which can then be used to compute interpolations and extrapolations
/// on arbitrary points.
/// </returns>
/// <remarks>
/// The value pairs do not have to be sorted, but if they are not sorted ascendingly
/// and the passed x and y arguments are arrays, they will be sorted inplace and thus modified.
/// if your data is already sorted in arrays, consider to use
/// MathNet.Numerics.Interpolation.LinearSpline.InterpolateSorted
/// instead, which is more efficient.
/// </remarks>
public static IInterpolation LinearSpline(IEnumerable<double> points, IEnumerable<double> values)
{
@ -139,18 +161,20 @@ namespace MathNet.Numerics
}
/// <summary>
/// Create an piecewise natural cubic spline interpolation based on arbitrary points, with zero secondary derivatives at the boundaries.
/// Create an piecewise natural cubic spline interpolation based on arbitrary points,
/// with zero secondary derivatives at the boundaries.
/// </summary>
/// <param name="points">The sample points t. Optimized for arrays.</param>
/// <param name="values">The sample point values x(t). Optimized for arrays.</param>
/// <param name="points">The sample points t.</param>
/// <param name="values">The sample point values x(t).</param>
/// <returns>
/// An interpolation scheme optimized for the given sample points and values,
/// which can then be used to compute interpolations and extrapolations
/// on arbitrary points.
/// </returns>
/// <remarks>
/// The value pairs do not have to be sorted, but if they are not sorted ascendingly
/// and the passed x and y arguments are arrays, they will be sorted inplace and thus modified.
/// if your data is already sorted in arrays, consider to use
/// MathNet.Numerics.Interpolation.CubicSpline.InterpolateNaturalSorted
/// instead, which is more efficient.
/// </remarks>
public static IInterpolation CubicSpline(IEnumerable<double> points, IEnumerable<double> values)
{
@ -158,18 +182,20 @@ namespace MathNet.Numerics
}
/// <summary>
/// Create an piecewise cubic Akima spline interpolation based on arbitrary points. Akima splines are robust to outliers.
/// Create an piecewise cubic Akima spline interpolation based on arbitrary points.
/// Akima splines are robust to outliers.
/// </summary>
/// <param name="points">The sample points t. Optimized for arrays.</param>
/// <param name="values">The sample point values x(t). Optimized for arrays.</param>
/// <param name="points">The sample points t.</param>
/// <param name="values">The sample point values x(t).</param>
/// <returns>
/// An interpolation scheme optimized for the given sample points and values,
/// which can then be used to compute interpolations and extrapolations
/// on arbitrary points.
/// </returns>
/// <remarks>
/// The value pairs do not have to be sorted, but if they are not sorted ascendingly
/// and the passed x and y arguments are arrays, they will be sorted inplace and thus modified.
/// if your data is already sorted in arrays, consider to use
/// MathNet.Numerics.Interpolation.CubicSpline.InterpolateAkimaSorted
/// instead, which is more efficient.
/// </remarks>
public static IInterpolation CubicSplineRobust(IEnumerable<double> points, IEnumerable<double> values)
{
@ -177,10 +203,11 @@ namespace MathNet.Numerics
}
/// <summary>
/// Create a piecewise cubic Hermite spline interpolation based on arbitrary points and their slopes/first derivative.
/// Create a piecewise cubic Hermite spline interpolation based on arbitrary points
/// and their slopes/first derivative.
/// </summary>
/// <param name="points">The sample points t. Optimized for arrays.</param>
/// <param name="values">The sample point values x(t). Optimized for arrays.</param>
/// <param name="points">The sample points t.</param>
/// <param name="values">The sample point values x(t).</param>
/// <param name="firstDerivatives">The slope at the sample points. Optimized for arrays.</param>
/// <returns>
/// An interpolation scheme optimized for the given sample points and values,
@ -188,8 +215,9 @@ namespace MathNet.Numerics
/// on arbitrary points.
/// </returns>
/// <remarks>
/// The value pairs do not have to be sorted, but if they are not sorted ascendingly
/// and the passed x and y arguments are arrays, they will be sorted inplace and thus modified.
/// if your data is already sorted in arrays, consider to use
/// MathNet.Numerics.Interpolation.CubicSpline.InterpolateHermiteSorted
/// instead, which is more efficient.
/// </remarks>
public static IInterpolation CubicSplineWithDerivatives(IEnumerable<double> points, IEnumerable<double> values, IEnumerable<double> firstDerivatives)
{

151
src/Numerics/Interpolation/Barycentric.cs
View File

@ -4,7 +4,7 @@
// http://github.com/mathnet/mathnet-numerics
// http://mathnetnumerics.codeplex.com
//
// Copyright (c) 2009-2013 Math.NET
// Copyright (c) 2009-2014 Math.NET
//
// Permission is hereby granted, free of charge, to any person
// obtaining a copy of this software and associated documentation
@ -45,9 +45,9 @@ namespace MathNet.Numerics.Interpolation
readonly double[] _y;
readonly double[] _w;
/// <param name="x">Sample points (N), no sorting assumed.</param>
/// <param name="y">Sample values (N).</param>
/// <param name="w">Barycentric weights (N).</param>
/// <param name="x">Sample points (N), sorted ascendingly.</param>
/// <param name="y">Sample values (N), sorted ascendingly by x.</param>
/// <param name="w">Barycentric weights (N), sorted ascendingly by x.</param>
public Barycentric(double[] x, double[] y, double[] w)
{
if (x.Length != y.Length || x.Length != w.Length)
@ -66,67 +66,80 @@ namespace MathNet.Numerics.Interpolation
}
/// <summary>
/// Create a barycentric polynomial interpolation from a set of (x,y) value pairs with equidistant x. No sorting is assumed.
/// Create a barycentric polynomial interpolation from a set of (x,y) value pairs with equidistant x, sorted ascendingly by x.
/// </summary>
/// <remarks>
/// The value pairs do not have to be sorted, but if they are not sorted ascendingly
/// and the passed x and y arguments are arrays, they will be sorted inplace and thus modified.
/// </remarks>
public static Barycentric InterpolatePolynomialEquidistant(IEnumerable<double> x, IEnumerable<double> y)
public static Barycentric InterpolatePolynomialEquidistantSorted(double[] x, double[] y)
{
var xx = (x as double[]) ?? x.ToArray();
var yy = (y as double[]) ?? y.ToArray();
if (xx.Length != yy.Length)
if (x.Length != y.Length)
{
throw new ArgumentException(Resources.ArgumentVectorsSameLength);
}
if (xx.Length < 1)
if (x.Length < 1)
{
throw new ArgumentOutOfRangeException("x");
}
Sorting.Sort(xx, yy);
var weights = new double[xx.Length];
var weights = new double[x.Length];
weights[0] = 1.0;
for (int i = 1; i < weights.Length; i++)
{
weights[i] = -(weights[i - 1]*(weights.Length - i))/i;
}
return new Barycentric(xx, yy, weights);
return new Barycentric(x, y, weights);
}
/// <summary>
/// Create a barycentric polynomial interpolation from an unordered set of (x,y) value pairs with equidistant x.
/// WARNING: Works in-place and can thus causes the data array to be reordered.
/// </summary>
public static Barycentric InterpolatePolynomialEquidistantInplace(double[] x, double[] y)
{
if (x.Length != y.Length)
{
throw new ArgumentException(Resources.ArgumentVectorsSameLength);
}
if (x.Length < 1)
{
throw new ArgumentOutOfRangeException("x");
}
Sorting.Sort(x, y);
return InterpolatePolynomialEquidistantSorted(x, y);
}
/// <summary>
/// Create a barycentric polynomial interpolation from an unsorted set of (x,y) value pairs with equidistant x.
/// </summary>
public static Barycentric InterpolatePolynomialEquidistant(IEnumerable<double> x, IEnumerable<double> y)
{
// note: we must make a copy, even if the input was arrays already
return InterpolatePolynomialEquidistantInplace(x.ToArray(), y.ToArray());
}
/// <summary>
/// Create a barycentric polynomial interpolation from a set of values related to linearly/equidistant spaced points within an interval.
/// </summary>
/// <remarks>
/// The value pairs do not have to be sorted, but if they are not sorted ascendingly
/// and the passed x and y arguments are arrays, they will be sorted inplace and thus modified.
/// </remarks>
public static Barycentric InterpolatePolynomialEquidistant(double leftBound, double rightBound, IEnumerable<double> y)
{
var yy = (y as double[]) ?? y.ToArray();
var xx = Generate.LinearSpaced(yy.Length, leftBound, rightBound);
return InterpolatePolynomialEquidistant(xx, yy);
return InterpolatePolynomialEquidistantSorted(xx, yy);
}
/// <summary>
/// Create a barycentric rational interpolation without poles, using Mike Floater and Kai Hormann's Algorithm.
/// The values are assumed to be sorted ascendingly by x.
/// </summary>
/// <param name="x">Sample points (N), no sorting assumed. Optimized for arrays.</param>
/// <param name="y">Sample values (N). Optimized for arrays.</param>
/// <param name="x">Sample points (N), sorted ascendingly.</param>
/// <param name="y">Sample values (N), sorted ascendingly by x.</param>
/// <param name="order">
/// Order of the interpolation scheme, 0 &lt;= order &lt;= N.
/// In most cases a value between 3 and 8 gives good results.
/// </param>
/// <remarks>
/// The value pairs do not have to be sorted, but if they are not sorted ascendingly
/// and the passed x and y arguments are arrays, they will be sorted inplace and thus modified.
/// </remarks>
public static Barycentric InterpolateRationalFloaterHormann(IEnumerable<double> x, IEnumerable<double> y, int order)
public static Barycentric InterpolateRationalFloaterHormannSorted(double[] x, double[] y, int order)
{
var xx = (x as double[]) ?? x.ToArray();
var yy = (y as double[]) ?? y.ToArray();
@ -175,18 +188,78 @@ namespace MathNet.Numerics.Interpolation
/// <summary>
/// Create a barycentric rational interpolation without poles, using Mike Floater and Kai Hormann's Algorithm.
/// WARNING: Works in-place and can thus causes the data array to be reordered.
/// </summary>
/// <param name="x">Sample points (N), no sorting assumed. Optimized for arrays.</param>
/// <param name="y">Sample values (N). Optimized for arrays.</param>
/// <remarks>
/// The value pairs do not have to be sorted, but if they are not sorted ascendingly
/// and the passed x and y arguments are arrays, they will be sorted inplace and thus modified.
/// </remarks>
/// <param name="x">Sample points (N), no sorting assumed.</param>
/// <param name="y">Sample values (N).</param>
/// <param name="order">
/// Order of the interpolation scheme, 0 &lt;= order &lt;= N.
/// In most cases a value between 3 and 8 gives good results.
/// </param>
public static Barycentric InterpolateRationalFloaterHormannInplace(double[] x, double[] y, int order)
{
if (x.Length != y.Length)
{
throw new ArgumentException(Resources.ArgumentVectorsSameLength);
}
if (0 > order || x.Length <= order)
{
throw new ArgumentOutOfRangeException("order");
}
Sorting.Sort(x, y);
return InterpolateRationalFloaterHormannSorted(x, y, order);
}
/// <summary>
/// Create a barycentric rational interpolation without poles, using Mike Floater and Kai Hormann's Algorithm.
/// </summary>
/// <param name="x">Sample points (N), no sorting assumed.</param>
/// <param name="y">Sample values (N).</param>
/// <param name="order">
/// Order of the interpolation scheme, 0 &lt;= order &lt;= N.
/// In most cases a value between 3 and 8 gives good results.
/// </param>
public static Barycentric InterpolateRationalFloaterHormann(IEnumerable<double> x, IEnumerable<double> y, int order)
{
// note: we must make a copy, even if the input was arrays already
return InterpolateRationalFloaterHormannInplace(x.ToArray(), y.ToArray(), order);
}
/// <summary>
/// Create a barycentric rational interpolation without poles, using Mike Floater and Kai Hormann's Algorithm.
/// The values are assumed to be sorted ascendingly by x.
/// </summary>
/// <param name="x">Sample points (N), sorted ascendingly.</param>
/// <param name="y">Sample values (N), sorted ascendingly by x.</param>
public static Barycentric InterpolateRationalFloaterHormannSorted(double[] x, double[] y)
{
return InterpolateRationalFloaterHormannSorted(x, y, Math.Min(3, x.Length - 1));
}
/// <summary>
/// Create a barycentric rational interpolation without poles, using Mike Floater and Kai Hormann's Algorithm.
/// WARNING: Works in-place and can thus causes the data array to be reordered.
/// </summary>
/// <param name="x">Sample points (N), no sorting assumed.</param>
/// <param name="y">Sample values (N).</param>
public static Barycentric InterpolateRationalFloaterHormannInplace(double[] x, double[] y)
{
return InterpolateRationalFloaterHormannInplace(x, y, Math.Min(3, x.Length - 1));
}
/// <summary>
/// Create a barycentric rational interpolation without poles, using Mike Floater and Kai Hormann's Algorithm.
/// </summary>
/// <param name="x">Sample points (N), no sorting assumed.</param>
/// <param name="y">Sample values (N).</param>
public static Barycentric InterpolateRationalFloaterHormann(IEnumerable<double> x, IEnumerable<double> y)
{
var xx = (x as double[]) ?? x.ToArray();
// note: we must make a copy, even if the input was arrays already
var xx = x.ToArray();
var order = Math.Min(3, xx.Length - 1);
return InterpolateRationalFloaterHormann(xx, y, order);
return InterpolateRationalFloaterHormannInplace(xx, y.ToArray(), order);
}
/// <summary>

54
src/Numerics/Interpolation/BulirschStoerRationalInterpolation.cs
View File

@ -4,7 +4,7 @@
// http://github.com/mathnet/mathnet-numerics
// http://mathnetnumerics.codeplex.com
//
// Copyright (c) 2009-2013 Math.NET
// Copyright (c) 2009-2014 Math.NET
//
// Permission is hereby granted, free of charge, to any person
// obtaining a copy of this software and associated documentation
@ -48,30 +48,54 @@ namespace MathNet.Numerics.Interpolation
readonly double[] _x;
readonly double[] _y;
/// <summary>
/// Initializes a new instance of the BulirschStoerRationalInterpolation class.
/// </summary>
/// <param name="x">Sample Points t</param>
/// <param name="y">Sample Values x(t)</param>
public BulirschStoerRationalInterpolation(IEnumerable<double> x, IEnumerable<double> y)
/// <param name="x">Sample Points t, sorted ascendingly.</param>
/// <param name="y">Sample Values x(t), sorted ascendingly by x.</param>
public BulirschStoerRationalInterpolation(double[] x, double[] y)
{
var xx = (x as double[]) ?? x.ToArray();
var yy = (y as double[]) ?? y.ToArray();
if (xx.Length != yy.Length)
if (x.Length != y.Length)
{
throw new ArgumentException(Resources.ArgumentVectorsSameLength);
}
if (xx.Length < 1)
if (x.Length < 1)
{
throw new ArgumentOutOfRangeException("x");
}
Sorting.Sort(xx, yy);
_x = x;
_y = y;
}
_x = xx;
_y = yy;
/// <summary>
/// Create a Bulirsch-Stoer rational interpolation from a set of (x,y) value pairs, sorted ascendingly by x.
/// </summary>
public static BulirschStoerRationalInterpolation InterpolateSorted(double[] x, double[] y)
{
return new BulirschStoerRationalInterpolation(x, y);
}
/// <summary>
/// Create a Bulirsch-Stoer rational interpolation from an unsorted set of (x,y) value pairs.
/// WARNING: Works in-place and can thus causes the data array to be reordered.
/// </summary>
public static BulirschStoerRationalInterpolation InterpolateInplace(double[] x, double[] y)
{
if (x.Length != y.Length)
{
throw new ArgumentException(Resources.ArgumentVectorsSameLength);
}
Sorting.Sort(x, y);
return InterpolateSorted(x, y);
}
/// <summary>
/// Create a Bulirsch-Stoer rational interpolation from an unsorted set of (x,y) value pairs.
/// </summary>
public static BulirschStoerRationalInterpolation Interpolate(IEnumerable<double> x, IEnumerable<double> y)
{
// note: we must make a copy, even if the input was arrays already
return InterpolateInplace(x.ToArray(), y.ToArray());
}
/// <summary>

239
src/Numerics/Interpolation/CubicSpline.cs
View File

@ -4,7 +4,7 @@
// http://github.com/mathnet/mathnet-numerics
// http://mathnetnumerics.codeplex.com
//
// Copyright (c) 2009-2013 Math.NET
// Copyright (c) 2009-2014 Math.NET
//
// Permission is hereby granted, free of charge, to any person
// obtaining a copy of this software and associated documentation
@ -69,79 +69,90 @@ namespace MathNet.Numerics.Interpolation
}
/// <summary>
/// Create a hermite cubic spline interpolation from a set of (x,y) value pairs and their slope (first derivative).
/// Create a hermite cubic spline interpolation from a set of (x,y) value pairs and their slope (first derivative), sorted ascendingly by x.
/// </summary>
/// <remarks>
/// The value pairs do not have to be sorted, but if they are not sorted ascendingly
/// and the passed x and y arguments are arrays, they will be sorted inplace and thus modified.
/// </remarks>
public static CubicSpline InterpolateHermite(IEnumerable<double> x, IEnumerable<double> y, IEnumerable<double> firstDerivatives)
public static CubicSpline InterpolateHermiteSorted(double[] x, double[] y, double[] firstDerivatives)
{
var xx = (x as double[]) ?? x.ToArray();
var yy = (y as double[]) ?? y.ToArray();
var dd = (firstDerivatives as double[]) ?? firstDerivatives.ToArray();
if (xx.Length != yy.Length || xx.Length != dd.Length)
if (x.Length != y.Length || x.Length != firstDerivatives.Length)
{
throw new ArgumentException(Resources.ArgumentVectorsSameLength);
}
if (xx.Length < 2)
if (x.Length < 2)
{
throw new ArgumentOutOfRangeException("x");
}
Sorting.Sort(xx, yy, dd);
var c0 = new double[xx.Length - 1];
var c1 = new double[xx.Length - 1];
var c2 = new double[xx.Length - 1];
var c3 = new double[xx.Length - 1];
var c0 = new double[x.Length - 1];
var c1 = new double[x.Length - 1];
var c2 = new double[x.Length - 1];
var c3 = new double[x.Length - 1];
for (int i = 0; i < c1.Length; i++)
{
double w = xx[i + 1] - xx[i];
double w = x[i + 1] - x[i];
double w2 = w*w;
c0[i] = yy[i];
c1[i] = dd[i];
c2[i] = (3*(yy[i + 1] - yy[i])/w - 2*dd[i] - dd[i + 1])/w;
c3[i] = (2*(yy[i] - yy[i + 1])/w + dd[i] + dd[i + 1])/w2;
c0[i] = y[i];
c1[i] = firstDerivatives[i];
c2[i] = (3*(y[i + 1] - y[i])/w - 2*firstDerivatives[i] - firstDerivatives[i + 1])/w;
c3[i] = (2*(y[i] - y[i + 1])/w + firstDerivatives[i] + firstDerivatives[i + 1])/w2;
}
return new CubicSpline(xx, c0, c1, c2, c3);
return new CubicSpline(x, c0, c1, c2, c3);
}
/// <summary>
/// Create an Akima cubic spline interpolation from a set of (x,y) value pairs. Akima splines are robust to outliers.
/// Create a hermite cubic spline interpolation from an unsorted set of (x,y) value pairs and their slope (first derivative).
/// WARNING: Works in-place and can thus causes the data array to be reordered.
/// </summary>
/// <remarks>
/// The value pairs do not have to be sorted, but if they are not sorted ascendingly
/// and the passed x and y arguments are arrays, they will be sorted inplace and thus modified.
/// </remarks>
public static CubicSpline InterpolateAkima(IEnumerable<double> x, IEnumerable<double> y)
public static CubicSpline InterpolateHermiteInplace(double[] x, double[] y, double[] firstDerivatives)
{
var xx = (x as double[]) ?? x.ToArray();
var yy = (y as double[]) ?? y.ToArray();
if (xx.Length != yy.Length)
if (x.Length != y.Length || x.Length != firstDerivatives.Length)
{
throw new ArgumentException(Resources.ArgumentVectorsSameLength);
}
if (xx.Length < 5)
if (x.Length < 2)
{
throw new ArgumentOutOfRangeException("x");
}
Sorting.Sort(xx, yy);
Sorting.Sort(x, y, firstDerivatives);
return InterpolateHermiteSorted(x, y, firstDerivatives);
}
/// <summary>
/// Create a hermite cubic spline interpolation from an unsorted set of (x,y) value pairs and their slope (first derivative).
/// </summary>
public static CubicSpline InterpolateHermite(IEnumerable<double> x, IEnumerable<double> y, IEnumerable<double> firstDerivatives)
{
// note: we must make a copy, even if the input was arrays already
return InterpolateHermiteInplace(x.ToArray(), y.ToArray(), firstDerivatives.ToArray());
}
/// <summary>
/// Create an Akima cubic spline interpolation from a set of (x,y) value pairs, sorted ascendingly by x.
/// Akima splines are robust to outliers.
/// </summary>
public static CubicSpline InterpolateAkimaSorted(double[] x, double[] y)
{
if (x.Length != y.Length)
{
throw new ArgumentException(Resources.ArgumentVectorsSameLength);
}
if (x.Length < 5)
{
throw new ArgumentOutOfRangeException("x");
}
/* Prepare divided differences (diff) and weights (w) */
var diff = new double[xx.Length - 1];
var weights = new double[xx.Length - 1];
var diff = new double[x.Length - 1];
var weights = new double[x.Length - 1];
for (int i = 0; i < diff.Length; i++)
{
diff[i] = (yy[i + 1] - yy[i])/(xx[i + 1] - xx[i]);
diff[i] = (y[i + 1] - y[i])/(x[i + 1] - x[i]);
}
for (int i = 1; i < weights.Length; i++)
@ -151,64 +162,75 @@ namespace MathNet.Numerics.Interpolation
/* Prepare Hermite interpolation scheme */
var dd = new double[xx.Length];
var dd = new double[x.Length];
for (int i = 2; i < dd.Length - 2; i++)
{
dd[i] = weights[i - 1].AlmostEqual(0.0) && weights[i + 1].AlmostEqual(0.0)
? (((xx[i + 1] - xx[i])*diff[i - 1]) + ((xx[i] - xx[i - 1])*diff[i]))/(xx[i + 1] - xx[i - 1])
? (((x[i + 1] - x[i])*diff[i - 1]) + ((x[i] - x[i - 1])*diff[i]))/(x[i + 1] - x[i - 1])
: ((weights[i + 1]*diff[i - 1]) + (weights[i - 1]*diff[i]))/(weights[i + 1] + weights[i - 1]);
}
dd[0] = DifferentiateThreePoint(xx, yy, 0, 0, 1, 2);
dd[1] = DifferentiateThreePoint(xx, yy, 1, 0, 1, 2);
dd[xx.Length - 2] = DifferentiateThreePoint(xx, yy, xx.Length - 2, xx.Length - 3, xx.Length - 2, xx.Length - 1);
dd[xx.Length - 1] = DifferentiateThreePoint(xx, yy, xx.Length - 1, xx.Length - 3, xx.Length - 2, xx.Length - 1);
dd[0] = DifferentiateThreePoint(x, y, 0, 0, 1, 2);
dd[1] = DifferentiateThreePoint(x, y, 1, 0, 1, 2);
dd[x.Length - 2] = DifferentiateThreePoint(x, y, x.Length - 2, x.Length - 3, x.Length - 2, x.Length - 1);
dd[x.Length - 1] = DifferentiateThreePoint(x, y, x.Length - 1, x.Length - 3, x.Length - 2, x.Length - 1);
/* Build Akima spline using Hermite interpolation scheme */
return InterpolateHermite(xx, yy, dd);
return InterpolateHermiteSorted(x, y, dd);
}
/// <summary>
/// Create a natural cubic spline interpolation from a set of (x,y) value pairs and zero second derivatives at the two boundaries.
/// Create an Akima cubic spline interpolation from an unsorted set of (x,y) value pairs.
/// Akima splines are robust to outliers.
/// WARNING: Works in-place and can thus causes the data array to be reordered.
/// </summary>
/// <remarks>
/// The value pairs do not have to be sorted, but if they are not sorted ascendingly
/// and the passed x and y arguments are arrays, they will be sorted inplace and thus modified.
/// </remarks>
public static CubicSpline InterpolateNatural(IEnumerable<double> x, IEnumerable<double> y)
public static CubicSpline InterpolateAkimaInplace(double[] x, double[] y)
{
return InterpolateBoundaries(x, y, SplineBoundaryCondition.SecondDerivative, 0.0, SplineBoundaryCondition.SecondDerivative, 0.0);
if (x.Length != y.Length)
{
throw new ArgumentException(Resources.ArgumentVectorsSameLength);
}
if (x.Length < 5)
{
throw new ArgumentOutOfRangeException("x");
}
Sorting.Sort(x, y);
return InterpolateAkimaSorted(x, y);
}
/// <summary>
/// Create a cubic spline interpolation from a set of (x,y) value pairs and custom boundary/termination conditions.
/// Create an Akima cubic spline interpolation from an unsorted set of (x,y) value pairs.
/// Akima splines are robust to outliers.
/// </summary>
/// <remarks>
/// The value pairs do not have to be sorted, but if they are not sorted ascendingly
/// and the passed x and y arguments are arrays, they will be sorted inplace and thus modified.
/// </remarks>
public static CubicSpline InterpolateBoundaries(IEnumerable<double> x, IEnumerable<double> y,
public static CubicSpline InterpolateAkima(IEnumerable<double> x, IEnumerable<double> y)
{
// note: we must make a copy, even if the input was arrays already
return InterpolateAkimaInplace(x.ToArray(), y.ToArray());
}
/// <summary>
/// Create a cubic spline interpolation from a set of (x,y) value pairs, sorted ascendingly by x,
/// and custom boundary/termination conditions.
/// </summary>
public static CubicSpline InterpolateBoundariesSorted(double[] x, double[] 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)
if (x.Length != y.Length)
{
throw new ArgumentException(Resources.ArgumentVectorsSameLength);
}
if (xx.Length < 2)
if (x.Length < 2)
{
throw new ArgumentOutOfRangeException("x");
}
Sorting.Sort(xx, yy);
int n = xx.Length;
int n = x.Length;
// normalize special cases
if ((n == 2)
@ -245,7 +267,7 @@ namespace MathNet.Numerics.Interpolation
a1[0] = 0;
a2[0] = 1;
a3[0] = 1;
b[0] = 2*(yy[1] - yy[0])/(xx[1] - xx[0]);
b[0] = 2*(y[1] - y[0])/(x[1] - x[0]);
break;
case SplineBoundaryCondition.FirstDerivative:
a1[0] = 0;
@ -257,19 +279,19 @@ namespace MathNet.Numerics.Interpolation
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]));
b[0] = (3*((y[1] - y[0])/(x[1] - x[0]))) - (0.5*leftBoundary*(x[1] - x[0]));
break;
default:
throw new NotSupportedException(Resources.InvalidLeftBoundaryCondition);
}
// Central Conditions
for (int i = 1; i < xx.Length - 1; i++)
for (int i = 1; i < x.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]));
a1[i] = x[i + 1] - x[i];
a2[i] = 2*(x[i + 1] - x[i - 1]);
a3[i] = x[i] - x[i - 1];
b[i] = (3*(y[i] - y[i - 1])/(x[i] - x[i - 1])*(x[i + 1] - x[i])) + (3*(y[i + 1] - y[i])/(x[i + 1] - x[i])*(x[i] - x[i - 1]));
}
// Right Boundary
@ -279,7 +301,7 @@ namespace MathNet.Numerics.Interpolation
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]);
b[n - 1] = 2*(y[n - 1] - y[n - 2])/(x[n - 1] - x[n - 2]);
break;
case SplineBoundaryCondition.FirstDerivative:
a1[n - 1] = 0;
@ -291,7 +313,7 @@ namespace MathNet.Numerics.Interpolation
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]));
b[n - 1] = (3*(y[n - 1] - y[n - 2])/(x[n - 1] - x[n - 2])) + (0.5*rightBoundary*(x[n - 1] - x[n - 2]));
break;
default:
throw new NotSupportedException(Resources.InvalidRightBoundaryCondition);
@ -299,7 +321,68 @@ namespace MathNet.Numerics.Interpolation
// Build Spline
double[] dd = SolveTridiagonal(a1, a2, a3, b);
return InterpolateHermite(xx, yy, dd);
return InterpolateHermiteSorted(x, y, dd);
}
/// <summary>
/// Create a cubic spline interpolation from an unsorted set of (x,y) value pairs and custom boundary/termination conditions.
/// WARNING: Works in-place and can thus causes the data array to be reordered.
/// </summary>
public static CubicSpline InterpolateBoundariesInplace(double[] x, double[] y,
SplineBoundaryCondition leftBoundaryCondition, double leftBoundary,
SplineBoundaryCondition rightBoundaryCondition, double rightBoundary)
{
if (x.Length != y.Length)
{
throw new ArgumentException(Resources.ArgumentVectorsSameLength);
}
if (x.Length < 2)
{
throw new ArgumentOutOfRangeException("x");
}
Sorting.Sort(x, y);
return InterpolateBoundariesSorted(x, y, leftBoundaryCondition, leftBoundary, rightBoundaryCondition, rightBoundary);
}
/// <summary>
/// Create a cubic spline interpolation from an unsorted set of (x,y) value pairs and custom boundary/termination conditions.
/// </summary>
public static CubicSpline InterpolateBoundaries(IEnumerable<double> x, IEnumerable<double> y,
SplineBoundaryCondition leftBoundaryCondition, double leftBoundary,
SplineBoundaryCondition rightBoundaryCondition, double rightBoundary)
{
// note: we must make a copy, even if the input was arrays already
return InterpolateBoundariesInplace(x.ToArray(), y.ToArray(), leftBoundaryCondition, leftBoundary, rightBoundaryCondition, rightBoundary);
}
/// <summary>
/// Create a natural cubic spline interpolation from a set of (x,y) value pairs
/// and zero second derivatives at the two boundaries, sorted ascendingly by x.
/// </summary>
public static CubicSpline InterpolateNaturalSorted(double[] x, double[] y)
{
return InterpolateBoundariesSorted(x, y, SplineBoundaryCondition.SecondDerivative, 0.0, SplineBoundaryCondition.SecondDerivative, 0.0);
}
/// <summary>
/// Create a natural cubic spline interpolation from an unsorted set of (x,y) value pairs
/// and zero second derivatives at the two boundaries.
/// WARNING: Works in-place and can thus causes the data array to be reordered.
/// </summary>
public static CubicSpline InterpolateNaturalInplace(double[] x, double[] y)
{
return InterpolateBoundariesInplace(x, y, SplineBoundaryCondition.SecondDerivative, 0.0, SplineBoundaryCondition.SecondDerivative, 0.0);
}
/// <summary>
/// Create a natural cubic spline interpolation from an unsorted set of (x,y) value pairs
/// and zero second derivatives at the two boundaries.
/// </summary>
public static CubicSpline InterpolateNatural(IEnumerable<double> x, IEnumerable<double> y)
{
return InterpolateBoundaries(x, y, SplineBoundaryCondition.SecondDerivative, 0.0, SplineBoundaryCondition.SecondDerivative, 0.0);
}
/// <summary>

47
src/Numerics/Interpolation/LinearSpline.cs
View File

@ -4,7 +4,7 @@
// http://github.com/mathnet/mathnet-numerics
// http://mathnetnumerics.codeplex.com
//
// Copyright (c) 2009-2013 Math.NET
// Copyright (c) 2009-2014 Math.NET
//
// Permission is hereby granted, free of charge, to any person
// obtaining a copy of this software and associated documentation
@ -63,31 +63,46 @@ namespace MathNet.Numerics.Interpolation
}
/// <summary>
/// Create a linear spline interpolation from a set of (x,y) value pairs.
/// Create a linear spline interpolation from a set of (x,y) value pairs, sorted ascendingly by x.
/// </summary>
/// <remarks>
/// The value pairs do not have to be sorted, but if they are not sorted ascendingly
/// and the passed x and y arguments are arrays, they will be sorted inplace and thus modified.
/// </remarks>
public static LinearSpline Interpolate(IEnumerable<double> x, IEnumerable<double> y)
public static LinearSpline InterpolateSorted(double[] x, double[] y)
{
var xx = (x as double[]) ?? x.ToArray();
var yy = (y as double[]) ?? y.ToArray();
if (xx.Length != yy.Length)
if (x.Length != y.Length)
{
throw new ArgumentException(Resources.ArgumentVectorsSameLength);
}
Sorting.Sort(xx, yy);
var c1 = new double[xx.Length - 1];
var c1 = new double[x.Length - 1];
for (int i = 0; i < c1.Length; i++)
{
c1[i] = (yy[i + 1] - yy[i])/(xx[i + 1] - xx[i]);
c1[i] = (y[i + 1] - y[i])/(x[i + 1] - x[i]);
}
return new LinearSpline(x, y, c1);
}
/// <summary>
/// Create a linear spline interpolation from an unsorted set of (x,y) value pairs.
/// WARNING: Works in-place and can thus causes the data array to be reordered.
/// </summary>
public static LinearSpline InterpolateInplace(double[] x, double[] y)
{
if (x.Length != y.Length)
{
throw new ArgumentException(Resources.ArgumentVectorsSameLength);
}
return new LinearSpline(xx, yy, c1);
Sorting.Sort(x, y);
return InterpolateSorted(x, y);
}
/// <summary>
/// Create a linear spline interpolation from an unsorted set of (x,y) value pairs.
/// </summary>
public static LinearSpline Interpolate(IEnumerable<double> x, IEnumerable<double> y)
{
// note: we must make a copy, even if the input was arrays already
return InterpolateInplace(x.ToArray(), y.ToArray());
}
/// <summary>

60
src/Numerics/Interpolation/NevillePolynomialInterpolation.cs
View File

@ -4,7 +4,7 @@
// http://github.com/mathnet/mathnet-numerics
// http://mathnetnumerics.codeplex.com
//
// Copyright (c) 2009-2013 Math.NET
// Copyright (c) 2009-2014 Math.NET
//
// Permission is hereby granted, free of charge, to any person
// obtaining a copy of this software and associated documentation
@ -53,38 +53,62 @@ namespace MathNet.Numerics.Interpolation
readonly double[] _x;
readonly double[] _y;
/// <summary>
/// Initializes a new instance of the NevillePolynomialInterpolation class.
/// </summary>
/// <param name="x">Sample Points t. Optimized for arrays.</param>
/// <param name="y">Sample Values x(t). Optimized for arrays.</param>
public NevillePolynomialInterpolation(IEnumerable<double> x, IEnumerable<double> y)
/// <param name="x">Sample Points t, sorted ascendingly.</param>
/// <param name="y">Sample Values x(t), sorted ascendingly by x.</param>
public NevillePolynomialInterpolation(double[] x, double[] y)
{
var xx = (x as double[]) ?? x.ToArray();
var yy = (y as double[]) ?? y.ToArray();
if (xx.Length != yy.Length)
if (x.Length != y.Length)
{
throw new ArgumentException(Resources.ArgumentVectorsSameLength);
}
if (xx.Length < 1)
if (x.Length < 1)
{
throw new ArgumentOutOfRangeException("x");
}
Sorting.Sort(xx, yy);
for (var i = 1; i < xx.Length; ++i)
for (var i = 1; i < x.Length; ++i)
{
if (xx[i] == xx[i - 1])
if (x[i] == x[i - 1])
{
throw new ArgumentException(Resources.Interpolation_Initialize_SamplePointsNotUnique, "x");
}
}
_x = xx;
_y = yy;
_x = x;
_y = y;
}
/// <summary>
/// Create a Neville polynomial interpolation from a set of (x,y) value pairs, sorted ascendingly by x.
/// </summary>
public static NevillePolynomialInterpolation InterpolateSorted(double[] x, double[] y)
{
return new NevillePolynomialInterpolation(x, y);
}
/// <summary>
/// Create a Neville polynomial interpolation from an unsorted set of (x,y) value pairs.
/// WARNING: Works in-place and can thus causes the data array to be reordered.
/// </summary>
public static NevillePolynomialInterpolation InterpolateInplace(double[] x, double[] y)
{
if (x.Length != y.Length)
{
throw new ArgumentException(Resources.ArgumentVectorsSameLength);
}
Sorting.Sort(x, y);
return InterpolateSorted(x, y);
}
/// <summary>
/// Create a Neville polynomial interpolation from an unsorted set of (x,y) value pairs.
/// </summary>
public static NevillePolynomialInterpolation Interpolate(IEnumerable<double> x, IEnumerable<double> y)
{
// note: we must make a copy, even if the input was arrays already
return InterpolateInplace(x.ToArray(), y.ToArray());
}
/// <summary>

4
src/UnitTests/InterpolationTests/BulirschStoerRationalTest.cs
View File

@ -55,7 +55,7 @@ namespace MathNet.Numerics.UnitTests.InterpolationTests
[Test]
public void FitsAtSamplePoints()
{
IInterpolation interpolation = new BulirschStoerRationalInterpolation(_t, _x);
IInterpolation interpolation = BulirschStoerRationalInterpolation.Interpolate(_t, _x);
for (int i = 0; i < _x.Length; i++)
{
@ -87,7 +87,7 @@ namespace MathNet.Numerics.UnitTests.InterpolationTests
[TestCase(-10.0, -3.6017584308603731307, 1e-13)]
public void FitsAtArbitraryPointsWithMaple(double t, double x, double maxAbsoluteError)
{
IInterpolation interpolation = new BulirschStoerRationalInterpolation(_t, _x);
IInterpolation interpolation = BulirschStoerRationalInterpolation.Interpolate(_t, _x);
Assert.AreEqual(x, interpolation.Interpolate(t), maxAbsoluteError, "Interpolation at {0}", t);
}

23
src/UnitTests/InterpolationTests/CubicSplineTest.cs
View File

@ -28,6 +28,9 @@
// OTHER DEALINGS IN THE SOFTWARE.
// </copyright>
using System;
using System.Linq;
using System.Threading.Tasks;
using MathNet.Numerics.Interpolation;
using NUnit.Framework;
@ -173,5 +176,25 @@ namespace MathNet.Numerics.UnitTests.InterpolationTests
Assert.AreEqual(ytest[i], it.Interpolate(xtest[i]), 1e-15, "Linear with {0} samples, sample {1}", samples, i);
}
}
[Test]
public void InterpolateAkimaSorted_MustBeThreadSafe_GitHub219([Values(8, 32, 256, 1024)] int samples)
{
var x = Generate.LinearSpaced(samples + 1, 0.0, 2.0*Math.PI);
var y = new double[samples][];
for (var i = 0; i < samples; ++i)
{
y[i] = x.Select(xx => Math.Sin(xx)/(i + 1)).ToArray();
}
var yipol = new double[samples];
Parallel.For(0, samples, i =>
{
var spline = CubicSpline.InterpolateAkimaSorted(x, y[i]);
yipol[i] = spline.Interpolate(1.0);
});
CollectionAssert.DoesNotContain(yipol, Double.NaN);
}
}
}

Write Preview
Loading…
Cancel
Save