From 54a63dc0d536884b335ad56c08592d93359b1a9e Mon Sep 17 00:00:00 2001 From: Christoph Ruegg Date: Sat, 21 Dec 2013 12:25:14 +0100 Subject: [PATCH] Interpolation: migrate rational floater hormann to Barycentric class --- src/Numerics/Interpolate.cs | 8 +- src/Numerics/Interpolation/Barycentric.cs | 78 ++++++++++++++++++- .../FloaterHormannRationalTest.cs | 33 +++----- 3 files changed, 91 insertions(+), 28 deletions(-) diff --git a/src/Numerics/Interpolate.cs b/src/Numerics/Interpolate.cs index b799d9c9..cc2ec45e 100644 --- a/src/Numerics/Interpolate.cs +++ b/src/Numerics/Interpolate.cs @@ -48,9 +48,9 @@ namespace MathNet.Numerics /// which can then be used to compute interpolations and extrapolations /// on arbitrary points. /// - public static IInterpolation Common(IList points, IList values) + public static IInterpolation Common(IEnumerable points, IEnumerable values) { - return new FloaterHormannRationalInterpolation(points, values); + return Interpolation.Barycentric.InterpolateRationalFloaterHormann(points, values); } /// @@ -63,9 +63,9 @@ namespace MathNet.Numerics /// which can then be used to compute interpolations and extrapolations /// on arbitrary points. /// - public static IInterpolation RationalWithoutPoles(IList points, IList values) + public static IInterpolation RationalWithoutPoles(IEnumerable points, IEnumerable values) { - return new FloaterHormannRationalInterpolation(points, values); + return Interpolation.Barycentric.InterpolateRationalFloaterHormann(points, values); } /// diff --git a/src/Numerics/Interpolation/Barycentric.cs b/src/Numerics/Interpolation/Barycentric.cs index c3e0b70d..b5145b59 100644 --- a/src/Numerics/Interpolation/Barycentric.cs +++ b/src/Numerics/Interpolation/Barycentric.cs @@ -113,6 +113,82 @@ namespace MathNet.Numerics.Interpolation return InterpolatePolynomialEquidistant(xx, yy); } + /// + /// Create a barycentric rational interpolation without poles, using Mike Floater and Kai Hormann's Algorithm. + /// + /// Sample points (N), no sorting assumed. Optimized for arrays. + /// Sample values (N). Optimized for arrays. + /// + /// 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) + { + 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 (0 > order || xx.Length <= order) + { + throw new ArgumentOutOfRangeException("order"); + } + + Sorting.Sort(xx, yy); + + var weights = new double[xx.Length]; + + // order: odd -> negative, even -> positive + double sign = ((order & 0x1) == 0x1) ? -1.0 : 1.0; + + // compute barycentric weights + for (int k = 0; k < xx.Length; k++) + { + double s = 0; + for (int i = Math.Max(k - order, 0); i <= Math.Min(k, weights.Length - 1 - order); i++) + { + double v = 1; + for (int j = i; j <= i + order; j++) + { + if (j != k) + { + v = v/Math.Abs(xx[k] - xx[j]); + } + } + + s = s + v; + } + + weights[k] = sign*s; + sign = -sign; + } + + return new Barycentric(xx, yy, weights); + } + + /// + /// Create a barycentric rational interpolation without poles, using Mike Floater and Kai Hormann's Algorithm. + /// + /// 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. + /// + public static Barycentric InterpolateRationalFloaterHormann(IEnumerable x, IEnumerable y) + { + var xx = (x as double[]) ?? x.ToArray(); + var order = Math.Min(3, xx.Length - 1); + return InterpolateRationalFloaterHormann(xx, y, order); + } + /// /// Gets a value indicating whether the algorithm supports differentiation (interpolated derivative). /// @@ -152,7 +228,7 @@ namespace MathNet.Numerics.Interpolation // trivial case: on a known sample point? if (offset == 0.0) { - // NOTE (cdrnet, 200908) not offset.AlmostZero() by design + // NOTE (cdrnet, 2009-08) not offset.AlmostZero() by design return _y[closestPoint]; } diff --git a/src/UnitTests/InterpolationTests/FloaterHormannRationalTest.cs b/src/UnitTests/InterpolationTests/FloaterHormannRationalTest.cs index b65da7ba..1aa9380c 100644 --- a/src/UnitTests/InterpolationTests/FloaterHormannRationalTest.cs +++ b/src/UnitTests/InterpolationTests/FloaterHormannRationalTest.cs @@ -33,21 +33,11 @@ using NUnit.Framework; namespace MathNet.Numerics.UnitTests.InterpolationTests { - /// - /// FloaterHormannRational test case. - /// [TestFixture, Category("Interpolation")] public class FloaterHormannRationalTest { - /// - /// Sample points. - /// readonly double[] _t = { -2.0, -1.0, 0.0, 1.0, 2.0 }; - - /// - /// Sample values. - /// - readonly double[] _x = { 1.0, 2.0, -1.0, 0.0, 1.0 }; + readonly double[] _y = { 1.0, 2.0, -1.0, 0.0, 1.0 }; /// /// Verifies that the interpolation matches the given value at all the provided polynomial sample points. @@ -55,11 +45,10 @@ namespace MathNet.Numerics.UnitTests.InterpolationTests [Test] public void PolyomnialFitsAtSamplePoints() { - IInterpolation interpolation = new FloaterHormannRationalInterpolation(_t, _x); - - for (int i = 0; i < _x.Length; i++) + IInterpolation it = Barycentric.InterpolateRationalFloaterHormann(_t, _y); + for (int i = 0; i < _y.Length; i++) { - Assert.AreEqual(_x[i], interpolation.Interpolate(_t[i]), "A Exact Point " + i); + Assert.AreEqual(_y[i], it.Interpolate(_t[i]), "A Exact Point " + i); } } @@ -85,9 +74,8 @@ namespace MathNet.Numerics.UnitTests.InterpolationTests [TestCase(-10.0, -5071.0, 1e-8)] public void PolynomialFitsAtArbitraryPointsWithMaple(double t, double x, double maxAbsoluteError) { - IInterpolation interpolation = new FloaterHormannRationalInterpolation(_t, _x); - - Assert.AreEqual(x, interpolation.Interpolate(t), maxAbsoluteError, "Interpolation at {0}", t); + IInterpolation it = Barycentric.InterpolateRationalFloaterHormann(_t, _y); + Assert.AreEqual(x, it.Interpolate(t), maxAbsoluteError, "Interpolation at {0}", t); } /// @@ -107,11 +95,10 @@ namespace MathNet.Numerics.UnitTests.InterpolationTests x[i] = 1.0/(1.0 + (tt*tt)); } - IInterpolation interpolation = new FloaterHormannRationalInterpolation(t, x); - + IInterpolation it = Barycentric.InterpolateRationalFloaterHormann(t, x); for (int i = 0; i < x.Length; i++) { - Assert.AreEqual(x[i], interpolation.Interpolate(t[i]), "A Exact Point " + i); + Assert.AreEqual(x[i], it.Interpolate(t[i]), "A Exact Point " + i); } } @@ -126,10 +113,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 FloaterHormannRationalInterpolation(x, y); + IInterpolation it = Barycentric.InterpolateRationalFloaterHormann(x, y); for (int i = 0; i < xtest.Length; i++) { - Assert.AreEqual(ytest[i], interpolation.Interpolate(xtest[i]), 1e-14, "Linear with {0} samples, sample {1}", samples, i); + Assert.AreEqual(ytest[i], it.Interpolate(xtest[i]), 1e-14, "Linear with {0} samples, sample {1}", samples, i); } } }