diff --git a/src/Numerics/Interpolate.cs b/src/Numerics/Interpolate.cs index f3c4250b..f074b6cb 100644 --- a/src/Numerics/Interpolate.cs +++ b/src/Numerics/Interpolate.cs @@ -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 /// /// Creates an interpolation based on arbitrary points. /// - /// The sample points t. Supports both lists and arrays. - /// The sample point values x(t). Supports both lists and arrays. + /// The sample points t. + /// The sample point values x(t). /// /// An interpolation scheme optimized for the given sample points and values, /// which can then be used to compute interpolations and extrapolations /// on arbitrary points. /// + /// + /// if your data is already sorted in arrays, consider to use + /// MathNet.Numerics.Interpolation.Barycentric.InterpolateRationalFloaterHormannSorted + /// instead, which is more efficient. + /// public static IInterpolation Common(IEnumerable points, IEnumerable values) { return Barycentric.InterpolateRationalFloaterHormann(points, values); @@ -56,46 +61,57 @@ namespace MathNet.Numerics /// /// Create a floater hormann rational pole-free interpolation based on arbitrary points. /// - /// The sample points t. Supports both lists and arrays. - /// The sample point values x(t). Supports both lists and arrays. + /// The sample points t. + /// The sample point values x(t). /// /// An interpolation scheme optimized for the given sample points and values, /// which can then be used to compute interpolations and extrapolations /// on arbitrary points. /// + /// + /// if your data is already sorted in arrays, consider to use + /// MathNet.Numerics.Interpolation.Barycentric.InterpolateRationalFloaterHormannSorted + /// instead, which is more efficient. + /// public static IInterpolation RationalWithoutPoles(IEnumerable points, IEnumerable values) { return Barycentric.InterpolateRationalFloaterHormann(points, values); } /// - /// Create a burlisch stoer rational interpolation based on arbitrary points. + /// Create a Bulirsch Stoer rational interpolation based on arbitrary points. /// - /// The sample points t. Optimized for arrays.. - /// The sample point values x(t). Optimized for arrays. + /// The sample points t. + /// The sample point values x(t). /// /// An interpolation scheme optimized for the given sample points and values, /// which can then be used to compute interpolations and extrapolations /// on arbitrary points. /// + /// + /// if your data is already sorted in arrays, consider to use + /// MathNet.Numerics.Interpolation.BulirschStoerRationalInterpolation.InterpolateSorted + /// instead, which is more efficient. + /// public static IInterpolation RationalWithPoles(IEnumerable points, IEnumerable values) { - return new BulirschStoerRationalInterpolation(points, values); + return BulirschStoerRationalInterpolation.Interpolate(points, values); } /// /// Create a barycentric polynomial interpolation where the given sample points are equidistant. /// - /// The sample points t, must be equidistant. Optimized for arrays. - /// The sample point values x(t). Optimized for arrays. + /// The sample points t, must be equidistant. + /// The sample point values x(t). /// /// An interpolation scheme optimized for the given sample points and values, /// which can then be used to compute interpolations and extrapolations /// on arbitrary points. /// /// - /// 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. /// public static IInterpolation PolynomialEquidistant(IEnumerable points, IEnumerable values) { @@ -103,35 +119,41 @@ namespace MathNet.Numerics } /// - /// 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. /// - /// The sample points t. Optimized for arrays. - /// The sample point values x(t). Optimized for arrays. + /// The sample points t. + /// The sample point values x(t). /// /// An interpolation scheme optimized for the given sample points and values, /// which can then be used to compute interpolations and extrapolations /// on arbitrary points. /// + /// + /// if your data is already sorted in arrays, consider to use + /// MathNet.Numerics.Interpolation.NevillePolynomialInterpolation.InterpolateSorted + /// instead, which is more efficient. + /// public static IInterpolation Polynomial(IEnumerable points, IEnumerable values) { - return new NevillePolynomialInterpolation(points, values); + return NevillePolynomialInterpolation.Interpolate(points, values); } /// /// Create a piecewise linear spline interpolation based on arbitrary points. /// - /// The sample points t. Optimized for arrays. - /// The sample point values x(t). Optimized for arrays. + /// The sample points t. + /// The sample point values x(t). /// /// An interpolation scheme optimized for the given sample points and values, /// which can then be used to compute interpolations and extrapolations /// on arbitrary points. /// /// - /// 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. /// public static IInterpolation LinearSpline(IEnumerable points, IEnumerable values) { @@ -139,18 +161,20 @@ namespace MathNet.Numerics } /// - /// 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. /// - /// The sample points t. Optimized for arrays. - /// The sample point values x(t). Optimized for arrays. + /// The sample points t. + /// The sample point values x(t). /// /// An interpolation scheme optimized for the given sample points and values, /// which can then be used to compute interpolations and extrapolations /// on arbitrary points. /// /// - /// 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. /// public static IInterpolation CubicSpline(IEnumerable points, IEnumerable values) { @@ -158,18 +182,20 @@ namespace MathNet.Numerics } /// - /// 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. /// - /// The sample points t. Optimized for arrays. - /// The sample point values x(t). Optimized for arrays. + /// The sample points t. + /// The sample point values x(t). /// /// An interpolation scheme optimized for the given sample points and values, /// which can then be used to compute interpolations and extrapolations /// on arbitrary points. /// /// - /// 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. /// public static IInterpolation CubicSplineRobust(IEnumerable points, IEnumerable values) { @@ -177,10 +203,11 @@ namespace MathNet.Numerics } /// - /// 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. /// - /// The sample points t. Optimized for arrays. - /// The sample point values x(t). Optimized for arrays. + /// The sample points t. + /// The sample point values x(t). /// The slope at the sample points. Optimized for arrays. /// /// An interpolation scheme optimized for the given sample points and values, @@ -188,8 +215,9 @@ namespace MathNet.Numerics /// on arbitrary points. /// /// - /// 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. /// public static IInterpolation CubicSplineWithDerivatives(IEnumerable points, IEnumerable values, IEnumerable firstDerivatives) { diff --git a/src/Numerics/Interpolation/Barycentric.cs b/src/Numerics/Interpolation/Barycentric.cs index b5145b59..29369c5b 100644 --- a/src/Numerics/Interpolation/Barycentric.cs +++ b/src/Numerics/Interpolation/Barycentric.cs @@ -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; - /// Sample points (N), no sorting assumed. - /// Sample values (N). - /// Barycentric weights (N). + /// Sample points (N), sorted ascendingly. + /// Sample values (N), sorted ascendingly by x. + /// Barycentric weights (N), sorted ascendingly by x. 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 } /// - /// 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. /// - /// - /// 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. - /// - public static Barycentric InterpolatePolynomialEquidistant(IEnumerable x, IEnumerable 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); + } + + /// + /// 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. + /// + 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); + } + + /// + /// Create a barycentric polynomial interpolation from an unsorted set of (x,y) value pairs with equidistant x. + /// + public static Barycentric InterpolatePolynomialEquidistant(IEnumerable x, IEnumerable y) + { + // note: we must make a copy, even if the input was arrays already + return InterpolatePolynomialEquidistantInplace(x.ToArray(), y.ToArray()); } /// /// Create a barycentric polynomial interpolation from a set of values related to linearly/equidistant spaced points within an interval. /// - /// - /// 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. - /// public static Barycentric InterpolatePolynomialEquidistant(double leftBound, double rightBound, IEnumerable y) { var yy = (y as double[]) ?? y.ToArray(); var xx = Generate.LinearSpaced(yy.Length, leftBound, rightBound); - return InterpolatePolynomialEquidistant(xx, yy); + return InterpolatePolynomialEquidistantSorted(xx, yy); } /// /// 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. /// - /// Sample points (N), no sorting assumed. Optimized for arrays. - /// Sample values (N). Optimized for arrays. + /// Sample points (N), sorted ascendingly. + /// Sample values (N), sorted ascendingly by x. /// /// Order of the interpolation scheme, 0 <= order <= N. /// In most cases a value between 3 and 8 gives good results. /// - /// - /// 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. - /// - public static Barycentric InterpolateRationalFloaterHormann(IEnumerable x, IEnumerable 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 /// /// 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. /// - /// Sample points (N), no sorting assumed. Optimized for arrays. - /// Sample values (N). Optimized for arrays. - /// - /// 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. - /// + /// Sample points (N), no sorting assumed. + /// Sample values (N). + /// + /// Order of the interpolation scheme, 0 <= order <= N. + /// In most cases a value between 3 and 8 gives good results. + /// + 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); + } + + /// + /// Create a barycentric rational interpolation without poles, using Mike Floater and Kai Hormann's Algorithm. + /// + /// Sample points (N), no sorting assumed. + /// Sample values (N). + /// + /// Order of the interpolation scheme, 0 <= order <= N. + /// In most cases a value between 3 and 8 gives good results. + /// + public static Barycentric InterpolateRationalFloaterHormann(IEnumerable x, IEnumerable y, int order) + { + // note: we must make a copy, even if the input was arrays already + return InterpolateRationalFloaterHormannInplace(x.ToArray(), y.ToArray(), order); + } + + /// + /// 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. + /// + /// Sample points (N), sorted ascendingly. + /// Sample values (N), sorted ascendingly by x. + public static Barycentric InterpolateRationalFloaterHormannSorted(double[] x, double[] y) + { + return InterpolateRationalFloaterHormannSorted(x, y, Math.Min(3, x.Length - 1)); + } + + /// + /// 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. + /// + /// Sample points (N), no sorting assumed. + /// Sample values (N). + public static Barycentric InterpolateRationalFloaterHormannInplace(double[] x, double[] y) + { + return InterpolateRationalFloaterHormannInplace(x, y, Math.Min(3, x.Length - 1)); + } + + /// + /// Create a barycentric rational interpolation without poles, using Mike Floater and Kai Hormann's Algorithm. + /// + /// Sample points (N), no sorting assumed. + /// Sample values (N). public static Barycentric InterpolateRationalFloaterHormann(IEnumerable x, IEnumerable 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); } /// diff --git a/src/Numerics/Interpolation/BulirschStoerRationalInterpolation.cs b/src/Numerics/Interpolation/BulirschStoerRationalInterpolation.cs index a31d910b..bc4495d6 100644 --- a/src/Numerics/Interpolation/BulirschStoerRationalInterpolation.cs +++ b/src/Numerics/Interpolation/BulirschStoerRationalInterpolation.cs @@ -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; - /// - /// Initializes a new instance of the BulirschStoerRationalInterpolation class. - /// - /// Sample Points t - /// Sample Values x(t) - public BulirschStoerRationalInterpolation(IEnumerable x, IEnumerable y) + /// Sample Points t, sorted ascendingly. + /// Sample Values x(t), sorted ascendingly by x. + 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; + /// + /// Create a Bulirsch-Stoer rational interpolation from a set of (x,y) value pairs, sorted ascendingly by x. + /// + public static BulirschStoerRationalInterpolation InterpolateSorted(double[] x, double[] y) + { + return new BulirschStoerRationalInterpolation(x, y); + } + + /// + /// 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. + /// + 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); + } + + /// + /// Create a Bulirsch-Stoer rational interpolation from an unsorted set of (x,y) value pairs. + /// + public static BulirschStoerRationalInterpolation Interpolate(IEnumerable x, IEnumerable y) + { + // note: we must make a copy, even if the input was arrays already + return InterpolateInplace(x.ToArray(), y.ToArray()); } /// diff --git a/src/Numerics/Interpolation/CubicSpline.cs b/src/Numerics/Interpolation/CubicSpline.cs index 8c8ea2bf..5152cc43 100644 --- a/src/Numerics/Interpolation/CubicSpline.cs +++ b/src/Numerics/Interpolation/CubicSpline.cs @@ -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 } /// - /// 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. /// - /// - /// 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. - /// - public static CubicSpline InterpolateHermite(IEnumerable x, IEnumerable y, IEnumerable 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); } /// - /// 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. /// - /// - /// 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. - /// - public static CubicSpline InterpolateAkima(IEnumerable x, IEnumerable 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); + } + + /// + /// Create a hermite cubic spline interpolation from an unsorted set of (x,y) value pairs and their slope (first derivative). + /// + public static CubicSpline InterpolateHermite(IEnumerable x, IEnumerable y, IEnumerable firstDerivatives) + { + // note: we must make a copy, even if the input was arrays already + return InterpolateHermiteInplace(x.ToArray(), y.ToArray(), firstDerivatives.ToArray()); + } + + /// + /// Create an Akima cubic spline interpolation from a set of (x,y) value pairs, sorted ascendingly by x. + /// Akima splines are robust to outliers. + /// + 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); } /// - /// 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. /// - /// - /// 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. - /// - public static CubicSpline InterpolateNatural(IEnumerable x, IEnumerable 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); } /// - /// 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. /// - /// - /// 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. - /// - public static CubicSpline InterpolateBoundaries(IEnumerable x, IEnumerable y, + public static CubicSpline InterpolateAkima(IEnumerable x, IEnumerable y) + { + // note: we must make a copy, even if the input was arrays already + return InterpolateAkimaInplace(x.ToArray(), y.ToArray()); + } + + /// + /// Create a cubic spline interpolation from a set of (x,y) value pairs, sorted ascendingly by x, + /// and custom boundary/termination conditions. + /// + 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); + } + + /// + /// 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. + /// + 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); + } + + /// + /// Create a cubic spline interpolation from an unsorted set of (x,y) value pairs and custom boundary/termination conditions. + /// + public static CubicSpline InterpolateBoundaries(IEnumerable x, IEnumerable 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); + } + + /// + /// 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. + /// + public static CubicSpline InterpolateNaturalSorted(double[] x, double[] y) + { + return InterpolateBoundariesSorted(x, y, SplineBoundaryCondition.SecondDerivative, 0.0, SplineBoundaryCondition.SecondDerivative, 0.0); + } + + /// + /// 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. + /// + public static CubicSpline InterpolateNaturalInplace(double[] x, double[] y) + { + return InterpolateBoundariesInplace(x, y, SplineBoundaryCondition.SecondDerivative, 0.0, SplineBoundaryCondition.SecondDerivative, 0.0); + } + + /// + /// Create a natural cubic spline interpolation from an unsorted set of (x,y) value pairs + /// and zero second derivatives at the two boundaries. + /// + public static CubicSpline InterpolateNatural(IEnumerable x, IEnumerable y) + { + return InterpolateBoundaries(x, y, SplineBoundaryCondition.SecondDerivative, 0.0, SplineBoundaryCondition.SecondDerivative, 0.0); } /// diff --git a/src/Numerics/Interpolation/LinearSpline.cs b/src/Numerics/Interpolation/LinearSpline.cs index eaaaf1f5..fb9277f7 100644 --- a/src/Numerics/Interpolation/LinearSpline.cs +++ b/src/Numerics/Interpolation/LinearSpline.cs @@ -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 } /// - /// 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. /// - /// - /// 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. - /// - public static LinearSpline Interpolate(IEnumerable x, IEnumerable 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); + } + + /// + /// 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. + /// + 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); + } + + /// + /// Create a linear spline interpolation from an unsorted set of (x,y) value pairs. + /// + public static LinearSpline Interpolate(IEnumerable x, IEnumerable y) + { + // note: we must make a copy, even if the input was arrays already + return InterpolateInplace(x.ToArray(), y.ToArray()); } /// diff --git a/src/Numerics/Interpolation/NevillePolynomialInterpolation.cs b/src/Numerics/Interpolation/NevillePolynomialInterpolation.cs index baf69009..9f2ad234 100644 --- a/src/Numerics/Interpolation/NevillePolynomialInterpolation.cs +++ b/src/Numerics/Interpolation/NevillePolynomialInterpolation.cs @@ -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; - /// - /// Initializes a new instance of the NevillePolynomialInterpolation class. - /// - /// Sample Points t. Optimized for arrays. - /// Sample Values x(t). Optimized for arrays. - public NevillePolynomialInterpolation(IEnumerable x, IEnumerable y) + /// Sample Points t, sorted ascendingly. + /// Sample Values x(t), sorted ascendingly by x. + 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; + } + + /// + /// Create a Neville polynomial interpolation from a set of (x,y) value pairs, sorted ascendingly by x. + /// + public static NevillePolynomialInterpolation InterpolateSorted(double[] x, double[] y) + { + return new NevillePolynomialInterpolation(x, y); + } + + /// + /// 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. + /// + 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); + } + + /// + /// Create a Neville polynomial interpolation from an unsorted set of (x,y) value pairs. + /// + public static NevillePolynomialInterpolation Interpolate(IEnumerable x, IEnumerable y) + { + // note: we must make a copy, even if the input was arrays already + return InterpolateInplace(x.ToArray(), y.ToArray()); } /// diff --git a/src/UnitTests/InterpolationTests/BulirschStoerRationalTest.cs b/src/UnitTests/InterpolationTests/BulirschStoerRationalTest.cs index 3db86cd6..2370a46e 100644 --- a/src/UnitTests/InterpolationTests/BulirschStoerRationalTest.cs +++ b/src/UnitTests/InterpolationTests/BulirschStoerRationalTest.cs @@ -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); } diff --git a/src/UnitTests/InterpolationTests/CubicSplineTest.cs b/src/UnitTests/InterpolationTests/CubicSplineTest.cs index eba63f64..b12dbd2a 100644 --- a/src/UnitTests/InterpolationTests/CubicSplineTest.cs +++ b/src/UnitTests/InterpolationTests/CubicSplineTest.cs @@ -28,6 +28,9 @@ // OTHER DEALINGS IN THE SOFTWARE. // +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); + } } }