Interpolation: migrate rational floater hormann to Barycentric class

src/Numerics/Interpolate.cs

@@ -48,9 +48,9 @@ namespace MathNet.Numerics
         /// which can then be used to compute interpolations and extrapolations
         /// on arbitrary points.
         /// </returns>
-        public static IInterpolation Common(IList<double> points, IList<double> values)
+        public static IInterpolation Common(IEnumerable<double> points, IEnumerable<double> values)
         {
-            return new FloaterHormannRationalInterpolation(points, values);
+            return Interpolation.Barycentric.InterpolateRationalFloaterHormann(points, values);
         }
 
         /// <summary>
@@ -63,9 +63,9 @@ namespace MathNet.Numerics
         /// which can then be used to compute interpolations and extrapolations
         /// on arbitrary points.
         /// </returns>
-        public static IInterpolation RationalWithoutPoles(IList<double> points, IList<double> values)
+        public static IInterpolation RationalWithoutPoles(IEnumerable<double> points, IEnumerable<double> values)
         {
-            return new FloaterHormannRationalInterpolation(points, values);
+            return Interpolation.Barycentric.InterpolateRationalFloaterHormann(points, values);
         }
 
         /// <summary>

src/Numerics/Interpolation/Barycentric.cs

@@ -113,6 +113,82 @@ namespace MathNet.Numerics.Interpolation
             return InterpolatePolynomialEquidistant(xx, yy);
         }
 
+        /// <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. Optimized for arrays.</param>
+        /// <param name="y">Sample values (N). Optimized for arrays.</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)
+        {
+            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);
+        }
+
+        /// <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. 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>
+        public static Barycentric InterpolateRationalFloaterHormann(IEnumerable<double> x, IEnumerable<double> y)
+        {
+            var xx = (x as double[]) ?? x.ToArray();
+            var order = Math.Min(3, xx.Length - 1);
+            return InterpolateRationalFloaterHormann(xx, y, order);
+        }
+
         /// <summary>
         /// Gets a value indicating whether the algorithm supports differentiation (interpolated derivative).
         /// </summary>
@@ -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];
             }
 

src/UnitTests/InterpolationTests/FloaterHormannRationalTest.cs

@@ -33,21 +33,11 @@ using NUnit.Framework;
 
 namespace MathNet.Numerics.UnitTests.InterpolationTests
 {
-    /// <summary>
-    /// FloaterHormannRational test case.
-    /// </summary>
     [TestFixture, Category("Interpolation")]
     public class FloaterHormannRationalTest
     {
-        /// <summary>
-        /// Sample points.
-        /// </summary>
         readonly double[] _t = { -2.0, -1.0, 0.0, 1.0, 2.0 };
-
-        /// <summary>
-        /// Sample values.
-        /// </summary>
-        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 };
 
         /// <summary>
         /// 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);
         }
 
         /// <summary>
@@ -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);
             }
         }
     }