diff --git a/src/Numerics/Interpolate.cs b/src/Numerics/Interpolate.cs index 202e9f50..08043329 100644 --- a/src/Numerics/Interpolate.cs +++ b/src/Numerics/Interpolate.cs @@ -142,7 +142,7 @@ namespace MathNet.Numerics } /// - /// Create a piecewise linear spline interpolation based on arbitrary points. + /// Create a piecewise linear interpolation based on arbitrary points. /// /// The sample points t. /// The sample point values x(t). @@ -156,30 +156,35 @@ namespace MathNet.Numerics /// MathNet.Numerics.Interpolation.LinearSpline.InterpolateSorted /// instead, which is more efficient. /// - public static IInterpolation LinearSpline(IEnumerable points, IEnumerable values) + public static IInterpolation Linear(IEnumerable points, IEnumerable values) { return Interpolation.LinearSpline.Interpolate(points, values); } /// - /// Create log linear spline interpolation based on arbitrary points. + /// Create piecewise log-linear 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 the values are passed as an array, they will be modified inplace, even it is already sorted. + /// if your data is already sorted in arrays, consider to use + /// MathNet.Numerics.Interpolation.LogLinear.InterpolateSorted + /// instead, which is more efficient. /// - public static IInterpolation LogLinearSpline(IEnumerable points, IEnumerable values) + public static IInterpolation LogLinear(IEnumerable points, IEnumerable values) + { + return Interpolation.LogLinear.Interpolate(points, values); + } + + [Obsolete("Use Linear instead. Will be removed in the next major version.")] + public static IInterpolation LinearSpline(IEnumerable points, IEnumerable values) { - return new LogLinearSpline(points, values); + return Interpolation.LinearSpline.Interpolate(points, values); } /// @@ -247,24 +252,23 @@ namespace MathNet.Numerics } /// - /// Create log linear spline interpolation based on arbitrary points. + /// Create a step-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 the values are passed as an array, they will be modified inplace, even it is already sorted. + /// if your data is already sorted in arrays, consider to use + /// MathNet.Numerics.Interpolation.StepInterpolation.InterpolateSorted + /// instead, which is more efficient. /// public static IInterpolation Step(IEnumerable points, IEnumerable values) { - return new StepInterpolation(points, values); + return StepInterpolation.Interpolate(points, values); } } } diff --git a/src/Numerics/Interpolation/Barycentric.cs b/src/Numerics/Interpolation/Barycentric.cs index 29369c5b..495e4d50 100644 --- a/src/Numerics/Interpolation/Barycentric.cs +++ b/src/Numerics/Interpolation/Barycentric.cs @@ -278,6 +278,11 @@ namespace MathNet.Numerics.Interpolation get { return false; } } + /// + /// Interpolate at point t. + /// + /// Point t to interpolate at. + /// Interpolated value x(t). public double Interpolate(double t) { // trivial case: only one sample? diff --git a/src/Numerics/Interpolation/LinearSpline.cs b/src/Numerics/Interpolation/LinearSpline.cs index fb9277f7..882f8aa6 100644 --- a/src/Numerics/Interpolation/LinearSpline.cs +++ b/src/Numerics/Interpolation/LinearSpline.cs @@ -36,7 +36,7 @@ using MathNet.Numerics.Properties; namespace MathNet.Numerics.Interpolation { /// - /// Linear Spline Interpolation. + /// Piece-wise Linear Interpolation. /// /// Supports both differentiation and integration. public class LinearSpline : IInterpolation diff --git a/src/Numerics/Interpolation/LogLinearSpline.cs b/src/Numerics/Interpolation/LogLinear.cs similarity index 58% rename from src/Numerics/Interpolation/LogLinearSpline.cs rename to src/Numerics/Interpolation/LogLinear.cs index 6ba0ea8e..e0f57f87 100644 --- a/src/Numerics/Interpolation/LogLinearSpline.cs +++ b/src/Numerics/Interpolation/LogLinear.cs @@ -1,10 +1,10 @@ -// +// // 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 +// 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 @@ -32,39 +32,80 @@ using MathNet.Numerics.Properties; using System; using System.Collections.Generic; using System.Linq; +using MathNet.Numerics.Threading; namespace MathNet.Numerics.Interpolation { /// - /// Log Linear Spline Interpolation + /// Piece-wise Log-Linear Interpolation /// /// This algorithm supports differentiation, not integration. - public class LogLinearSpline : IInterpolation + public class LogLinear : IInterpolation { /// /// Internal Spline Interpolation /// - private readonly LinearSpline _spline; + readonly LinearSpline _spline; + + /// Sample points (N), sorted ascending + /// Natural logarithm of the sample values (N) at the corresponding points + public LogLinear(double[] x, double[] logy) + { + _spline = LinearSpline.InterpolateSorted(x, logy); + } /// - /// Creates a log linear interpolation based on input data + /// Create a piecewise log-linear interpolation from a set of (x,y) value pairs, sorted ascendingly by x. /// - /// - /// - public LogLinearSpline(IEnumerable x, IEnumerable y) + public static LogLinear InterpolateSorted(double[] x, double[] y) { - var xx = (x as double[]) ?? x.ToArray(); - var yy = (y as double[]) ?? y.ToArray(); + if (x.Length != y.Length) + { + throw new ArgumentException(Resources.ArgumentVectorsSameLength); + } - if (xx.Length != yy.Length) + var logy = new double[y.Length]; + CommonParallel.For(0, y.Length, 4096, (a, b) => + { + for (int i = a; i < b; i++) + { + logy[i] = Math.Log(y[i]); + } + }); + + return new LogLinear(x, logy); + } + + /// + /// Create a piecewise log-linear interpolation from an unsorted set of (x,y) value pairs. + /// WARNING: Works in-place and can thus causes the data array to be reordered and modified. + /// + public static LogLinear InterpolateInplace(double[] x, double[] y) + { + if (x.Length != y.Length) { throw new ArgumentException(Resources.ArgumentVectorsSameLength); } - for (int i = 0; i < yy.Length; i++) - yy[i] = Math.Log(yy[i]); + Sorting.Sort(x, y); + CommonParallel.For(0, y.Length, 4096, (a, b) => + { + for (int i = a; i < b; i++) + { + y[i] = Math.Log(y[i]); + } + }); - _spline = LinearSpline.Interpolate(xx, yy); + return new LogLinear(x, y); + } + + /// + /// Create a piecewise log-linear interpolation from an unsorted set of (x,y) value pairs. + /// + public static LogLinear Interpolate(IEnumerable x, IEnumerable y) + { + // note: we must make a copy, even if the input was arrays already + return InterpolateInplace(x.ToArray(), y.ToArray()); } /// @@ -100,7 +141,7 @@ namespace MathNet.Numerics.Interpolation /// Interpolated first derivative at point t. public double Differentiate(double t) { - return Interpolate(t) * _spline.Differentiate(t); + return Interpolate(t)*_spline.Differentiate(t); } /// @@ -113,8 +154,8 @@ namespace MathNet.Numerics.Interpolation var linearFirstDerivative = _spline.Differentiate(t); var linearSecondDerivative = _spline.Differentiate2(t); - var secondDerivative = Differentiate(t) * linearFirstDerivative + - Interpolate(t) * linearSecondDerivative; + var secondDerivative = Differentiate(t)*linearFirstDerivative + + Interpolate(t)*linearSecondDerivative; return secondDerivative; } @@ -123,9 +164,9 @@ namespace MathNet.Numerics.Interpolation /// Indefinite integral at point t. /// /// Point t to integrate at. - public double Integrate(double t) + double IInterpolation.Integrate(double t) { - throw new NotImplementedException(); + throw new NotSupportedException(); } /// @@ -133,9 +174,9 @@ namespace MathNet.Numerics.Interpolation /// /// Left bound of the integration interval [a,b]. /// Right bound of the integration interval [a,b]. - public double Integrate(double a, double b) + double IInterpolation.Integrate(double a, double b) { - throw new NotImplementedException(); + throw new NotSupportedException(); } } } diff --git a/src/Numerics/Interpolation/StepInterpolation.cs b/src/Numerics/Interpolation/StepInterpolation.cs index 18892568..7f90d7e4 100644 --- a/src/Numerics/Interpolation/StepInterpolation.cs +++ b/src/Numerics/Interpolation/StepInterpolation.cs @@ -36,52 +36,104 @@ using MathNet.Numerics.Properties; namespace MathNet.Numerics.Interpolation { /// - /// A step function where the start of each segment is included, and last is excluded. Segment i is [x_i, x_i+1). + /// A step function where the start of each segment is included, and last is excluded. + /// Segment i is [x_i, x_i+1). /// The domain of the function is all real numbers, such that y = 0 where x < x_0 or x gt; x_n /// + /// Supports both differentiation and integration. public class StepInterpolation : IInterpolation { readonly double[] _x; readonly double[] _y; readonly Lazy _indefiniteIntegral; - /// Sample points (N+1) sorted in ascending order - /// Functional value (N) of each segment. - public StepInterpolation(IEnumerable x, IEnumerable y) + /// Sample points (N+1), sorted ascending + /// Functional values (N) of each segment + public StepInterpolation(double[] x, double[] sy) { - var xx = (x as double[]) ?? x.ToArray(); - var yy = (y as double[]) ?? y.ToArray(); - - if (xx.Length != yy.Length + 1) + if (x.Length != sy.Length + 1) { throw new ArgumentException(Resources.ArgumentVectorsSameLength); } - _x = xx; - _y = yy; + _x = x; + _y = sy; _indefiniteIntegral = new Lazy(ComputeIndefiniteIntegral); } - double[] ComputeIndefiniteIntegral() + /// + /// Create a linear spline interpolation from a set of (x,y) value pairs, sorted ascendingly by x. + /// The y-value corresponding to the largest x-sample is ignored. + /// + public static StepInterpolation InterpolateSorted(double[] x, double[] y) { - var integral = new double[_x.Length]; - for (int i = 0; i < integral.Length - 1; i++) + if (x.Length != y.Length) { - integral[i + 1] = integral[i] + (_x[i + 1] - _x[i])*_y[i]; + throw new ArgumentException(Resources.ArgumentVectorsSameLength); } - return integral; + + // drop the last value which is not part of any segment. + var segmentValues = new double[x.Length - 1]; + Array.Copy(y, 0, segmentValues, 0, segmentValues.Length); + + return new StepInterpolation(x, segmentValues); } - public bool SupportsDifferentiation + /// + /// Create a linear spline interpolation from an unsorted set of (x,y) value pairs. + /// The y-value corresponding to the largest x-sample is ignored. + /// WARNING: Works in-place and can thus causes the data array to be reordered. + /// + public static StepInterpolation 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 linear spline interpolation from an unsorted set of (x,y) value pairs. + /// The y-value corresponding to the largest x-sample is ignored. + /// + public static StepInterpolation Interpolate(IEnumerable x, IEnumerable y) + { + // note: we must make a copy, even if the input was arrays already + return InterpolateInplace(x.ToArray(), y.ToArray()); + } + + bool IInterpolation.SupportsDifferentiation { get { return true; } } - public bool SupportsIntegration + bool IInterpolation.SupportsIntegration { get { return true; } } + /// + /// Interpolate at point t. + /// + /// Point t to interpolate at. + /// Interpolated value x(t). + public double Interpolate(double t) + { + if (t < _x[0] || t >= _x[_x.Length - 1]) + return 0.0; + + int k = LeftBracketIndex(t); + return _y[k]; + } + + /// + /// Differentiate at point t. + /// + /// Point t to interpolate at. + /// Interpolated first derivative at point t. public double Differentiate(double t) { int index = Array.BinarySearch(_x, t); @@ -90,6 +142,11 @@ namespace MathNet.Numerics.Interpolation return 0d; } + /// + /// Differentiate twice at point t. + /// + /// Point t to interpolate at. + /// Interpolated second derivative at point t. public double Differentiate2(double t) { return Differentiate(t); @@ -122,13 +179,14 @@ namespace MathNet.Numerics.Interpolation return Integrate(b) - Integrate(a); } - public double Interpolate(double t) + double[] ComputeIndefiniteIntegral() { - if (t < _x[0] || t >= _x[_x.Length - 1]) - return 0.0; - - int k = LeftBracketIndex(t); - return _y[k]; + var integral = new double[_x.Length]; + for (int i = 0; i < integral.Length - 1; i++) + { + integral[i + 1] = integral[i] + (_x[i + 1] - _x[i])*_y[i]; + } + return integral; } /// diff --git a/src/Numerics/Interpolation/TransformedInterpolation.cs b/src/Numerics/Interpolation/TransformedInterpolation.cs new file mode 100644 index 00000000..8f2b918a --- /dev/null +++ b/src/Numerics/Interpolation/TransformedInterpolation.cs @@ -0,0 +1,175 @@ +// +// 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-2014 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. +// + +using MathNet.Numerics.Properties; +using System; +using System.Collections.Generic; +using System.Linq; +using MathNet.Numerics.Threading; + +namespace MathNet.Numerics.Interpolation +{ + public class TransformedInterpolation : IInterpolation + { + readonly IInterpolation _interpolation; + readonly Func _transform; + + public TransformedInterpolation(IInterpolation interpolation, Func transform) + { + _interpolation = interpolation; + _transform = transform; + } + + /// + /// Create a linear spline interpolation from a set of (x,y) value pairs, sorted ascendingly by x. + /// + public static TransformedInterpolation InterpolateSorted( + Func transform, Func transformInverse, + double[] x, double[] y) + { + if (x.Length != y.Length) + { + throw new ArgumentException(Resources.ArgumentVectorsSameLength); + } + + var yhat = new double[y.Length]; + CommonParallel.For(0, y.Length, 4096, (a, b) => + { + for (int i = a; i < b; i++) + { + yhat[i] = transformInverse(y[i]); + } + }); + + return new TransformedInterpolation(LinearSpline.InterpolateSorted(x, yhat), transform); + } + + /// + /// 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 and modified. + /// + public static TransformedInterpolation InterpolateInplace( + Func transform, Func transformInverse, + double[] x, double[] y) + { + if (x.Length != y.Length) + { + throw new ArgumentException(Resources.ArgumentVectorsSameLength); + } + + Sorting.Sort(x, y); + CommonParallel.For(0, y.Length, 4096, (a, b) => + { + for (int i = a; i < b; i++) + { + y[i] = transformInverse(y[i]); + } + }); + + return new TransformedInterpolation(LinearSpline.InterpolateSorted(x, y), transform); + } + + /// + /// Create a linear spline interpolation from an unsorted set of (x,y) value pairs. + /// + public static TransformedInterpolation Interpolate( + Func transform, Func transformInverse, + IEnumerable x, IEnumerable y) + { + // note: we must make a copy, even if the input was arrays already + return InterpolateInplace(transform, transformInverse, x.ToArray(), y.ToArray()); + } + + /// + /// Gets a value indicating whether the algorithm supports differentiation (interpolated derivative). + /// + bool IInterpolation.SupportsDifferentiation + { + get { return false; } + } + + /// + /// Gets a value indicating whether the algorithm supports integration (interpolated quadrature). + /// + bool IInterpolation.SupportsIntegration + { + get { return false; } + } + + /// + /// Interpolate at point t. + /// + /// Point t to interpolate at. + /// Interpolated value x(t). + public double Interpolate(double t) + { + return _transform(_interpolation.Interpolate(t)); + } + + /// + /// Differentiate at point t. NOT SUPPORTED. + /// + /// Point t to interpolate at. + /// Interpolated first derivative at point t. + double IInterpolation.Differentiate(double t) + { + throw new NotSupportedException(); + } + + /// + /// Differentiate twice at point t. NOT SUPPORTED. + /// + /// Point t to interpolate at. + /// Interpolated second derivative at point t. + double IInterpolation.Differentiate2(double t) + { + throw new NotSupportedException(); + } + + /// + /// Indefinite integral at point t. NOT SUPPORTED. + /// + /// Point t to integrate at. + double IInterpolation.Integrate(double t) + { + throw new NotSupportedException(); + } + + /// + /// Definite integral between points a and b. NOT SUPPORTED. + /// + /// Left bound of the integration interval [a,b]. + /// Right bound of the integration interval [a,b]. + double IInterpolation.Integrate(double a, double b) + { + throw new NotSupportedException(); + } + } +} diff --git a/src/Numerics/Interpolation/TransformedSpline.cs b/src/Numerics/Interpolation/TransformedSpline.cs deleted file mode 100644 index 03147afa..00000000 --- a/src/Numerics/Interpolation/TransformedSpline.cs +++ /dev/null @@ -1,76 +0,0 @@ -using MathNet.Numerics.Properties; -using System; -using System.Collections.Generic; -using System.Linq; -using System.Text; - -namespace MathNet.Numerics.Interpolation -{ - public class TransformedInterpolation : IInterpolation - { - private IInterpolation _baseInterpolation; - private Func _transformer; - - public TransformedInterpolation(IInterpolation baseInterp, Func transformer) - { - _baseInterpolation = baseInterp; - _transformer = transformer; - } - - public static TransformedInterpolation Interpolate( - Func transformer, Func transformerInverse, - IEnumerable x, IEnumerable y) - { - var xx = (x as double[]) ?? x.ToArray(); - var yy = (y as double[]) ?? y.ToArray(); - - if (xx.Length != yy.Length) - { - throw new ArgumentException(Resources.ArgumentVectorsSameLength); - } - - for (int i = 0; i < yy.Length; i++) - yy[i] = transformerInverse(yy[i]); - - var baseInterp = LinearSpline.Interpolate(xx, yy); - var interp = new TransformedInterpolation(baseInterp, transformer); - - return interp; - } - - public bool SupportsDifferentiation - { - get { return false; } - } - - public bool SupportsIntegration - { - get { return false; } - } - - public double Differentiate(double t) - { - throw new NotImplementedException(); - } - - public double Differentiate2(double t) - { - throw new NotImplementedException(); - } - - public double Integrate(double t) - { - throw new NotImplementedException(); - } - - public double Integrate(double a, double b) - { - throw new NotImplementedException(); - } - - public double Interpolate(double t) - { - return _transformer(_baseInterpolation.Interpolate(t)); - } - } -} diff --git a/src/Numerics/Numerics.csproj b/src/Numerics/Numerics.csproj index 5d75eb5e..a766e8ce 100644 --- a/src/Numerics/Numerics.csproj +++ b/src/Numerics/Numerics.csproj @@ -91,10 +91,10 @@ - + - +