diff --git a/src/Numerics/Interpolate.cs b/src/Numerics/Interpolate.cs index cffc0446..b799d9c9 100644 --- a/src/Numerics/Interpolate.cs +++ b/src/Numerics/Interpolate.cs @@ -86,16 +86,20 @@ namespace MathNet.Numerics /// /// Create a barycentric polynomial interpolation where the given sample points are equidistant. /// - /// The sample points t, must be equidistant. Supports both lists and arrays. - /// The sample point values x(t). Supports both lists and arrays. + /// The sample points t, must be equidistant. Optimized for arrays. + /// The sample point values x(t). Optimized for arrays. /// /// An interpolation scheme optimized for the given sample points and values, /// which can then be used to compute interpolations and extrapolations /// on arbitrary points. /// - public static IInterpolation PolynomialEquidistant(IList points, IList values) + /// + /// 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 IInterpolation PolynomialEquidistant(IEnumerable points, IEnumerable values) { - return new EquidistantPolynomialInterpolation(points, values); + return Interpolation.Barycentric.InterpolatePolynomialEquidistant(points, values); } /// diff --git a/src/Numerics/Interpolation/Barycentric.cs b/src/Numerics/Interpolation/Barycentric.cs new file mode 100644 index 00000000..c3e0b70d --- /dev/null +++ b/src/Numerics/Interpolation/Barycentric.cs @@ -0,0 +1,226 @@ +// +// Math.NET Numerics, part of the Math.NET Project +// http://numerics.mathdotnet.com +// http://github.com/mathnet/mathnet-numerics +// http://mathnetnumerics.codeplex.com +// +// Copyright (c) 2009-2013 Math.NET +// +// Permission is hereby granted, free of charge, to any person +// obtaining a copy of this software and associated documentation +// files (the "Software"), to deal in the Software without +// restriction, including without limitation the rights to use, +// copy, modify, merge, publish, distribute, sublicense, and/or sell +// copies of the Software, and to permit persons to whom the +// Software is furnished to do so, subject to the following +// conditions: +// +// The above copyright notice and this permission notice shall be +// included in all copies or substantial portions of the Software. +// +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, +// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES +// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND +// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT +// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, +// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING +// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR +// OTHER DEALINGS IN THE SOFTWARE. +// + +using System; +using System.Collections.Generic; +using System.Linq; +using MathNet.Numerics.Properties; + +namespace MathNet.Numerics.Interpolation +{ + /// + /// Barycentric Interpolation Algorithm. + /// + /// Supports neither differentiation nor integration. + public class Barycentric : IInterpolation + { + readonly double[] _x; + readonly double[] _y; + readonly double[] _w; + + /// Sample points (N), no sorting assumed. + /// Sample values (N). + /// Barycentric weights (N). + public Barycentric(double[] x, double[] y, double[] w) + { + if (x.Length != y.Length || x.Length != w.Length) + { + throw new ArgumentException(Resources.ArgumentVectorsSameLength); + } + + if (x.Length < 1) + { + throw new ArgumentOutOfRangeException("x"); + } + + _x = x; + _y = y; + _w = w; + } + + /// + /// Create a barycentric polynomial interpolation from a set of (x,y) value pairs with equidistant x. No sorting is assumed. + /// + /// + /// 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) + { + var xx = (x as double[]) ?? x.ToArray(); + var yy = (y as double[]) ?? y.ToArray(); + + if (xx.Length != yy.Length) + { + throw new ArgumentException(Resources.ArgumentVectorsSameLength); + } + + if (xx.Length < 1) + { + throw new ArgumentOutOfRangeException("x"); + } + + Sorting.Sort(xx, yy); + + var weights = new double[xx.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); + } + + /// + /// 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); + } + + /// + /// 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; } + } + + public double Interpolate(double t) + { + // trivial case: only one sample? + if (_x.Length == 1) + { + return _y[0]; + } + + // evaluate closest point and offset from that point (no sorting assumed) + int closestPoint = 0; + double offset = t - _x[0]; + for (int i = 1; i < _x.Length; i++) + { + if (Math.Abs(t - _x[i]) < Math.Abs(offset)) + { + offset = t - _x[i]; + closestPoint = i; + } + } + + // trivial case: on a known sample point? + if (offset == 0.0) + { + // NOTE (cdrnet, 200908) not offset.AlmostZero() by design + return _y[closestPoint]; + } + + if (Math.Abs(offset) > 1e-150) + { + // no need to guard against overflow, so use fast formula + closestPoint = -1; + offset = 1.0; + } + + double s1 = 0.0; + double s2 = 0.0; + for (int i = 0; i < _x.Length; i++) + { + if (i != closestPoint) + { + double v = offset*_w[i]/(t - _x[i]); + s1 = s1 + (v*_y[i]); + s2 = s2 + v; + } + else + { + double v = _w[i]; + s1 = s1 + (v*_y[i]); + s2 = s2 + v; + } + } + + return s1/s2; + } + + /// + /// 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/Numerics.csproj b/src/Numerics/Numerics.csproj index 0057f245..89d46a5a 100644 --- a/src/Numerics/Numerics.csproj +++ b/src/Numerics/Numerics.csproj @@ -90,6 +90,7 @@ + diff --git a/src/UnitTests/InterpolationTests/EquidistantPolynomialTest.cs b/src/UnitTests/InterpolationTests/EquidistantPolynomialTest.cs index e7c4d618..c1ef7f0e 100644 --- a/src/UnitTests/InterpolationTests/EquidistantPolynomialTest.cs +++ b/src/UnitTests/InterpolationTests/EquidistantPolynomialTest.cs @@ -33,26 +33,12 @@ using NUnit.Framework; namespace MathNet.Numerics.UnitTests.InterpolationTests { - /// - /// EquidistantPolynomial Test case. - /// [TestFixture, Category("Interpolation")] public class EquidistantPolynomialTest { - /// - /// Left bound; - /// const double Tmin = 0.0; - - /// - /// Right bound. - /// const double Tmax = 4.0; - - /// - /// Sample values. - /// - readonly double[] _x = { 0.0, 3.0, 2.5, 1.0, 3.0 }; + readonly double[] _y = { 0.0, 3.0, 2.5, 1.0, 3.0 }; /// /// Verifies that the interpolation matches the given value at all the provided sample points. @@ -60,11 +46,10 @@ namespace MathNet.Numerics.UnitTests.InterpolationTests [Test] public void FitsAtSamplePoints() { - IInterpolation interpolation = new EquidistantPolynomialInterpolation(Tmin, Tmax, _x); - - for (int i = 0; i < _x.Length; i++) + IInterpolation it = Barycentric.InterpolatePolynomialEquidistant(Tmin, Tmax, _y); + for (int i = 0; i < _y.Length; i++) { - Assert.AreEqual(_x[i], interpolation.Interpolate(i), "A Exact Point " + i); + Assert.AreEqual(_y[i], it.Interpolate(i), "A Exact Point " + i); } } @@ -86,11 +71,10 @@ namespace MathNet.Numerics.UnitTests.InterpolationTests [TestCase(4.5, 7.265625, 1e-14)] [TestCase(10.0, 592.5, 1e-10)] [TestCase(-10.0, 657.5, 1e-9)] - public void FitsAtArbitraryPointsWithMaple(double t, double x, double maxAbsoluteError) + public void FitsAtArbitraryPoints(double t, double x, double maxAbsoluteError) { - IInterpolation interpolation = new EquidistantPolynomialInterpolation(Tmin, Tmax, _x); - - Assert.AreEqual(x, interpolation.Interpolate(t), maxAbsoluteError, "Interpolation at {0}", t); + IInterpolation it = Barycentric.InterpolatePolynomialEquidistant(Tmin, Tmax, _y); + Assert.AreEqual(x, it.Interpolate(t), maxAbsoluteError, "Interpolation at {0}", t); } /// @@ -104,10 +88,10 @@ namespace MathNet.Numerics.UnitTests.InterpolationTests { double[] x, y, xtest, ytest; LinearInterpolationCase.Build(out x, out y, out xtest, out ytest, samples); - IInterpolation interpolation = new EquidistantPolynomialInterpolation(x, y); + IInterpolation it = Barycentric.InterpolatePolynomialEquidistant(x, y); for (int i = 0; i < xtest.Length; i++) { - Assert.AreEqual(ytest[i], interpolation.Interpolate(xtest[i]), 1e-12, "Linear with {0} samples, sample {1}", samples, i); + Assert.AreEqual(ytest[i], it.Interpolate(xtest[i]), 1e-12, "Linear with {0} samples, sample {1}", samples, i); } } } diff --git a/src/UnitTests/InterpolationTests/FloaterHormannRationalTest.cs b/src/UnitTests/InterpolationTests/FloaterHormannRationalTest.cs index 37f7be20..b65da7ba 100644 --- a/src/UnitTests/InterpolationTests/FloaterHormannRationalTest.cs +++ b/src/UnitTests/InterpolationTests/FloaterHormannRationalTest.cs @@ -81,11 +81,11 @@ namespace MathNet.Numerics.UnitTests.InterpolationTests [TestCase(0.1, -1.10805, 1e-15)] [TestCase(0.4, -1.1248, 1e-15)] [TestCase(1.2, 0.5392, 1e-15)] - [TestCase(10.0, -4431.0, 1e-9)] - [TestCase(-10.0, -5071.0, 1e-9)] + [TestCase(10.0, -4431.0, 1e-8)] + [TestCase(-10.0, -5071.0, 1e-8)] public void PolynomialFitsAtArbitraryPointsWithMaple(double t, double x, double maxAbsoluteError) { - IInterpolation interpolation = new EquidistantPolynomialInterpolation(_t, _x); + IInterpolation interpolation = new FloaterHormannRationalInterpolation(_t, _x); Assert.AreEqual(x, interpolation.Interpolate(t), maxAbsoluteError, "Interpolation at {0}", t); }