Repository: Hipparchus-Math/hipparchus (java)

Issue:
Add Field equivalent to PolynomialFunctionLagrangeForm
There is currently no `Field` implementation of Lagrange interpolation

Pull request:
add Field version of Lagrange Polynomials
Closes #260 

Patch:
diff --git a/hipparchus-core/src/changes/changes.xml b/hipparchus-core/src/changes/changes.xml
index 9d94f8383..f095cc617 100644
--- a/hipparchus-core/src/changes/changes.xml
+++ b/hipparchus-core/src/changes/changes.xml
@@ -50,6 +50,9 @@ If the output is not quite correct, check for invisible trailing spaces!
   </properties>
   <body>
     <release version="4.0" date="TBD" description="TBD">
+      <action dev="serrof" type="add" issue="issues/260">
+        Added FieldPolynomialFunctionLagrangeForm.
+      </action>
       <action dev="serrof" type="add" issue="issues/370">
         Added native way to stack an independent variable in (Field)Gradient.
       </action>
diff --git a/hipparchus-core/src/main/java/org/hipparchus/analysis/polynomials/FieldPolynomialFunctionLagrangeForm.java b/hipparchus-core/src/main/java/org/hipparchus/analysis/polynomials/FieldPolynomialFunctionLagrangeForm.java
new file mode 100644
index 000000000..76b7d4394
--- /dev/null
+++ b/hipparchus-core/src/main/java/org/hipparchus/analysis/polynomials/FieldPolynomialFunctionLagrangeForm.java
@@ -0,0 +1,232 @@
+/*
+ * Licensed to the Hipparchus project under one or more
+ * contributor license agreements.  See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * The Hipparchus project licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License.  You may obtain a copy of the License at
+ *
+ *      https://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+package org.hipparchus.analysis.polynomials;
+
+import org.hipparchus.CalculusFieldElement;
+import org.hipparchus.Field;
+import org.hipparchus.analysis.CalculusFieldUnivariateFunction;
+import org.hipparchus.exception.LocalizedCoreFormats;
+import org.hipparchus.exception.MathIllegalArgumentException;
+import org.hipparchus.util.FastMath;
+import org.hipparchus.util.MathArrays;
+
+/**
+ * Implements the representation of a real polynomial function in
+ * <a href="http://mathworld.wolfram.com/LagrangeInterpolatingPolynomial.html">
+ * Lagrange Form</a>. For reference, see <b>Introduction to Numerical
+ * Analysis</b>, ISBN 038795452X, chapter 2.
+ * <p>
+ * The approximated function should be smooth enough for Lagrange polynomial
+ * to work well. Otherwise, consider using splines instead.</p>
+ * @see PolynomialFunctionLagrangeForm
+ * @since 4.0
+ *
+ */
+public class FieldPolynomialFunctionLagrangeForm<T extends CalculusFieldElement<T>>
+        implements CalculusFieldUnivariateFunction<T> {
+    /**
+     * The coefficients of the polynomial, ordered by degree -- i.e.
+     * coefficients[0] is the constant term and coefficients[n] is the
+     * coefficient of x^n where n is the degree of the polynomial.
+     */
+    private T[] coefficients;
+    /**
+     * Interpolating points (abscissas).
+     */
+    private final T[] x;
+    /**
+     * Function values at interpolating points.
+     */
+    private final T[] y;
+    /**
+     * Whether the polynomial coefficients are available.
+     */
+    private boolean coefficientsComputed;
+
+    /**
+     * Construct a Lagrange polynomial with the given abscissas and function
+     * values. The order of interpolating points is important.
+     * <p>
+     * The constructor makes copy of the input arrays and assigns them.</p>
+     *
+     * @param x interpolating points
+     * @param y function values at interpolating points
+     * @throws MathIllegalArgumentException if the array lengths are different.
+     * @throws MathIllegalArgumentException if the number of points is less than 2.
+     * @throws MathIllegalArgumentException if two abscissae have the same value.
+     * @throws MathIllegalArgumentException if the abscissae are not sorted.
+     */
+    public FieldPolynomialFunctionLagrangeForm(final T[] x, final T[] y)
+        throws MathIllegalArgumentException {
+        this.x = x.clone();
+        this.y = y.clone();
+        coefficientsComputed = false;
+
+        MathArrays.checkEqualLength(x, y);
+        if (x.length < 2) {
+            throw new MathIllegalArgumentException(LocalizedCoreFormats.WRONG_NUMBER_OF_POINTS, 2, x.length, true);
+        }
+        MathArrays.checkOrder(x, MathArrays.OrderDirection.INCREASING, true, true);
+    }
+
+    /**
+     * Calculate the function value at the given point.
+     *
+     * @param z Point at which the function value is to be computed.
+     * @return the function value.
+     * @throws MathIllegalArgumentException if {@code x} and {@code y} have
+     * different lengths.
+     * @throws MathIllegalArgumentException
+     * if {@code x} is not sorted in strictly increasing order.
+     * @throws MathIllegalArgumentException if the size of {@code x} is less
+     * than 2.
+     */
+    @Override
+    public T value(final T z) {
+        int nearest = 0;
+        final int n = x.length;
+        final T[] c = y.clone();
+        final T[] d = c.clone();
+        double minDist = Double.POSITIVE_INFINITY;
+        for (int i = 0; i < n; i++) {
+            // find out the abscissa closest to z
+            final double dist = FastMath.abs(z.subtract(x[i])).getReal();
+            if (dist < minDist) {
+                nearest = i;
+                minDist = dist;
+            }
+        }
+
+        // initial approximation to the function value at z
+        T value = y[nearest];
+
+        for (int i = 1; i < n; i++) {
+            for (int j = 0; j < n-i; j++) {
+                final T tc = x[j].subtract(z);
+                final T td = x[i+j].subtract(z);
+                final T divider = x[j].subtract(x[i+j]);
+                // update the difference arrays
+                final T w = (c[j+1].subtract(d[j])).divide(divider);
+                c[j] = tc.multiply(w);
+                d[j] = td.multiply(w);
+            }
+            // sum up the difference terms to get the final value
+            if (nearest < 0.5*(n-i+1)) {
+                value = value.add(c[nearest]);    // fork down
+            } else {
+                nearest--;
+                value = value.add(d[nearest]);    // fork up
+            }
+        }
+
+        return value;
+    }
+
+    /**
+     * Returns the degree of the polynomial.
+     *
+     * @return the degree of the polynomial
+     */
+    public int degree() {
+        return x.length - 1;
+    }
+
+    /**
+     * Returns a copy of the interpolating points array.
+     * <p>
+     * Changes made to the returned copy will not affect the polynomial.</p>
+     *
+     * @return a fresh copy of the interpolating points array
+     */
+    public T[] getInterpolatingPoints() {
+        return x.clone();
+    }
+
+    /**
+     * Returns a copy of the interpolating values array.
+     * <p>
+     * Changes made to the returned copy will not affect the polynomial.</p>
+     *
+     * @return a fresh copy of the interpolating values array
+     */
+    public T[] getInterpolatingValues() {
+        return y.clone();
+    }
+
+    /**
+     * Returns a copy of the coefficients array.
+     * <p>
+     * Changes made to the returned copy will not affect the polynomial.</p>
+     * <p>
+     * Note that coefficients computation can be ill-conditioned. Use with caution
+     * and only when it is necessary.</p>
+     *
+     * @return a fresh copy of the coefficients array
+     */
+    public T[] getCoefficients() {
+        if (!coefficientsComputed) {
+            computeCoefficients();
+        }
+        return coefficients.clone();
+    }
+
+    /**
+     * Calculate the coefficients of Lagrange polynomial from the
+     * interpolation data. It takes O(n^2) time.
+     * Note that this computation can be ill-conditioned: Use with caution
+     * and only when it is necessary.
+     */
+    protected void computeCoefficients() {
+        final int n = degree() + 1;
+        final Field<T> field = x[0].getField();
+        coefficients = MathArrays.buildArray(field, n);
+
+        // c[] are the coefficients of P(x) = (x-x[0])(x-x[1])...(x-x[n-1])
+        final T[] c = MathArrays.buildArray(field, n + 1);
+        c[0] = field.getOne();
+        for (int i = 0; i < n; i++) {
+            for (int j = i; j > 0; j--) {
+                c[j] = c[j-1].subtract(c[j].multiply(x[i]));
+            }
+            c[0] = c[0].multiply(x[i].negate());
+            c[i+1] = field.getOne();
+        }
+
+        final T[] tc = MathArrays.buildArray(field, n);
+        for (int i = 0; i < n; i++) {
+            // d = (x[i]-x[0])...(x[i]-x[i-1])(x[i]-x[i+1])...(x[i]-x[n-1])
+            T d = field.getOne();
+            for (int j = 0; j < n; j++) {
+                if (i != j) {
+                    d = d.multiply(x[i].subtract(x[j]));
+                }
+            }
+            final T t = y[i].divide(d);
+            // Lagrange polynomial is the sum of n terms, each of which is a
+            // polynomial of degree n-1. tc[] are the coefficients of the i-th
+            // numerator Pi(x) = (x-x[0])...(x-x[i-1])(x-x[i+1])...(x-x[n-1]).
+            tc[n-1] = c[n];     // actually c[n] = 1
+            coefficients[n-1] = coefficients[n-1].add(t.multiply(tc[n-1]));
+            for (int j = n-2; j >= 0; j--) {
+                tc[j] = c[j+1].add(tc[j+1].multiply(x[i]));
+                coefficients[j] = coefficients[j].add(t.multiply(tc[j]));
+            }
+        }
+
+        coefficientsComputed = true;
+    }
+}
diff --git a/hipparchus-core/src/main/java/org/hipparchus/analysis/polynomials/PolynomialFunctionLagrangeForm.java b/hipparchus-core/src/main/java/org/hipparchus/analysis/polynomials/PolynomialFunctionLagrangeForm.java
index e30b672bf..d1f962126 100644
--- a/hipparchus-core/src/main/java/org/hipparchus/analysis/polynomials/PolynomialFunctionLagrangeForm.java
+++ b/hipparchus-core/src/main/java/org/hipparchus/analysis/polynomials/PolynomialFunctionLagrangeForm.java
@@ -211,16 +211,16 @@ public class PolynomialFunctionLagrangeForm implements UnivariateFunction {
         final int n = x.length;
         final double[] c = new double[n];
         final double[] d = new double[n];
-        double min_dist = Double.POSITIVE_INFINITY;
+        double minDist = Double.POSITIVE_INFINITY;
         for (int i = 0; i < n; i++) {
             // initialize the difference arrays
             c[i] = y[i];
             d[i] = y[i];
             // find out the abscissa closest to z
             final double dist = FastMath.abs(z - x[i]);
-            if (dist < min_dist) {
+            if (dist < minDist) {
                 nearest = i;
-                min_dist = dist;
+                minDist = dist;
             }
         }
 
diff --git a/src/changes/changes.xml b/src/changes/changes.xml
index 041fe78e1..417b01b61 100644
--- a/src/changes/changes.xml
+++ b/src/changes/changes.xml
@@ -50,6 +50,9 @@ If the output is not quite correct, check for invisible trailing spaces!
   </properties>
   <body>
     <release version="4.0" date="TBD" description="TBD">
+      <action dev="serrof" type="add" issue="issues/260">
+        Added FieldPolynomialFunctionLagrangeForm.
+      </action>
       <action dev="serrof" type="update" issue="issues/372">
         Rename RungeKutta(Field)Integrator as FixedStepRungeKutta(Field)Integrator.
       </action>


Test patch:
diff --git a/hipparchus-core/src/test/java/org/hipparchus/analysis/polynomials/FieldPolynomialFunctionLagrangeFormTest.java b/hipparchus-core/src/test/java/org/hipparchus/analysis/polynomials/FieldPolynomialFunctionLagrangeFormTest.java
new file mode 100644
index 000000000..76aba7e31
--- /dev/null
+++ b/hipparchus-core/src/test/java/org/hipparchus/analysis/polynomials/FieldPolynomialFunctionLagrangeFormTest.java
@@ -0,0 +1,131 @@
+/*
+ * Licensed to the Hipparchus project under one or more
+ * contributor license agreements.  See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * The Hipparchus project licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License.  You may obtain a copy of the License at
+ *
+ *      https://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+package org.hipparchus.analysis.polynomials;
+
+import org.hipparchus.exception.MathIllegalArgumentException;
+import org.hipparchus.util.Binary64;
+import org.hipparchus.util.Binary64Field;
+import org.hipparchus.util.MathArrays;
+import org.junit.jupiter.api.Test;
+import org.junit.jupiter.params.ParameterizedTest;
+import org.junit.jupiter.params.provider.ValueSource;
+
+import static org.junit.jupiter.api.Assertions.*;
+
+class FieldPolynomialFunctionLagrangeFormTest {
+
+    @Test
+    void testExceptionSize() {
+        // GIVEN
+        final Binary64Field field = Binary64Field.getInstance();
+        final Binary64[] x = MathArrays.buildArray(field, 1);
+        final Binary64[] y = x.clone();
+        // WHEN & THEN
+        assertThrows(MathIllegalArgumentException.class, () -> new FieldPolynomialFunctionLagrangeForm<>(x, y));
+    }
+
+    @Test
+    void testExceptionOrder() {
+        // GIVEN
+        final Binary64Field field = Binary64Field.getInstance();
+        final Binary64[] x = MathArrays.buildArray(field, 2);
+        x[0] = field.getOne();
+        x[1] = field.getZero();
+        final Binary64[] y = x.clone();
+        // WHEN & THEN
+        assertThrows(MathIllegalArgumentException.class, () -> new FieldPolynomialFunctionLagrangeForm<>(x, y));
+    }
+
+    @Test
+    void testGetter() {
+        // GIVEN
+        final Binary64Field field = Binary64Field.getInstance();
+        Binary64[] x = { field.getOne(), field.getOne().multiply(2) };
+        Binary64[] y = { field.getOne().negate(), field.getOne().multiply(4.) };
+        final FieldPolynomialFunctionLagrangeForm<Binary64> lagrangeForm = new FieldPolynomialFunctionLagrangeForm<>(x, y);
+        // WHEN & THEN
+        assertArrayEquals(x, lagrangeForm.getInterpolatingPoints());
+        assertArrayEquals(y, lagrangeForm.getInterpolatingValues());
+        final Binary64[] expected = MathArrays.buildArray(field, 2);
+        expected[0] = field.getOne().multiply(-6.);
+        expected[1] = field.getOne().multiply(5.);
+        assertArrayEquals(expected, lagrangeForm.getCoefficients());
+    }
+
+    @Test
+    void testGetCoefficients() {
+        // GIVEN
+        final Binary64Field field = Binary64Field.getInstance();
+        Binary64[] x = { field.getOne(), field.getOne().multiply(2) };
+        Binary64[] y = { field.getOne().negate(), field.getOne().multiply(4.) };
+        final FieldPolynomialFunctionLagrangeForm<Binary64> lagrangeForm = new FieldPolynomialFunctionLagrangeForm<>(x, y);
+        // WHEN
+        final Binary64[] coefficients = lagrangeForm.getCoefficients();
+        // WHEN & THEN
+        assertArrayEquals(lagrangeForm.getCoefficients(), coefficients);
+    }
+
+    @ParameterizedTest
+    @ValueSource(ints = {2, 3, 4, 5, 6, 7, 8, 9, 10})
+    void testValueDegree(final int n) {
+        // GIVEN
+        final Binary64Field field = Binary64Field.getInstance();
+        final Binary64[] x = MathArrays.buildArray(field, n);
+        final Binary64[] y = MathArrays.buildArray(field, n);
+        for (int i = 0; i < n; i++) {
+            x[i] = field.getOne().multiply(i);
+            y[i] = field.getOne().multiply(2 * i + 1);
+        }
+        final FieldPolynomialFunctionLagrangeForm<Binary64> lagrangeForm = new FieldPolynomialFunctionLagrangeForm<>(x, y);
+        final Binary64 five = field.getOne().multiply(5);
+        // WHEN
+        final Binary64 fieldValue = lagrangeForm.value(five);
+        // WHEN & THEN
+        final double[] xDouble = new double[n];
+        final double[] yDouble = new double[n];
+        for (int i = 0; i < n; i++) {
+            xDouble[i] = x[i].getReal();
+            yDouble[i] = y[i].getReal();
+        }
+        assertEquals(PolynomialFunctionLagrangeForm.evaluate(xDouble, yDouble, five.getReal()),
+                fieldValue.getReal());
+    }
+
+    @ParameterizedTest
+    @ValueSource(doubles = {2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.6, 11})
+    void testValuePoint(final double z) {
+        // GIVEN
+        final Binary64Field field = Binary64Field.getInstance();
+        final Binary64[] x = MathArrays.buildArray(field, 10);
+        final Binary64[] y = MathArrays.buildArray(field, 10);
+        for (int i = 0; i < x.length; i++) {
+            x[i] = field.getOne().multiply(i);
+            y[i] = field.getOne().multiply(2 * i + 1);
+        }
+        final FieldPolynomialFunctionLagrangeForm<Binary64> lagrangeForm = new FieldPolynomialFunctionLagrangeForm<>(x, y);
+        // WHEN
+        final Binary64 fieldValue = lagrangeForm.value(new Binary64(z));
+        // WHEN & THEN
+        final double[] xDouble = new double[x.length];
+        final double[] yDouble = xDouble.clone();
+        for (int i = 0; i < xDouble.length; i++) {
+            xDouble[i] = x[i].getReal();
+            yDouble[i] = y[i].getReal();
+        }
+        assertEquals(PolynomialFunctionLagrangeForm.evaluate(xDouble, yDouble, z), fieldValue.getReal());
+    }
+}
diff --git a/hipparchus-core/src/test/java/org/hipparchus/analysis/polynomials/PolynomialFunctionLagrangeFormTest.java b/hipparchus-core/src/test/java/org/hipparchus/analysis/polynomials/PolynomialFunctionLagrangeFormTest.java
index 8f0a8e809..7c09a0754 100644
--- a/hipparchus-core/src/test/java/org/hipparchus/analysis/polynomials/PolynomialFunctionLagrangeFormTest.java
+++ b/hipparchus-core/src/test/java/org/hipparchus/analysis/polynomials/PolynomialFunctionLagrangeFormTest.java
@@ -24,6 +24,7 @@ package org.hipparchus.analysis.polynomials;
 import org.hipparchus.exception.MathIllegalArgumentException;
 import org.junit.jupiter.api.Test;
 
+import static org.junit.jupiter.api.Assertions.assertArrayEquals;
 import static org.junit.jupiter.api.Assertions.assertEquals;
 import static org.junit.jupiter.api.Assertions.fail;
 
@@ -37,6 +38,43 @@ import static org.junit.jupiter.api.Assertions.fail;
  */
 final class PolynomialFunctionLagrangeFormTest {
 
+    @Test
+    void testGetter() {
+        // GIVEN
+        double[] x = { 1., 2. };
+        double[] y = { -1., 4. };
+        final PolynomialFunctionLagrangeForm lagrangeForm = new PolynomialFunctionLagrangeForm(x, y);
+        // WHEN & THEN
+        assertArrayEquals(x, lagrangeForm.getInterpolatingPoints());
+        assertArrayEquals(y, lagrangeForm.getInterpolatingValues());
+        assertArrayEquals(new double[] { -6., 5. }, lagrangeForm.getCoefficients());
+    }
+
+    @Test
+    void testGetCoefficients() {
+        // GIVEN
+        double[] x = { 1., 2. };
+        double[] y = { -1., 4. };
+        final PolynomialFunctionLagrangeForm lagrangeForm = new PolynomialFunctionLagrangeForm(x, y);
+        // WHEN
+        final double[] coefficients = lagrangeForm.getCoefficients();
+        // WHEN & THEN
+        assertArrayEquals(lagrangeForm.getCoefficients(), coefficients);
+    }
+
+    @Test
+    void testEvaluate() {
+        // GIVEN
+        double[] x = { 1., 2. };
+        double[] y = { -1., 4. };
+        final double point = 3.;
+        // WHEN
+        final double value = PolynomialFunctionLagrangeForm.evaluate(x, y, point);
+        // THEN
+        final PolynomialFunctionLagrangeForm lagrangeForm = new PolynomialFunctionLagrangeForm(x, y);
+        assertEquals(lagrangeForm.value(point), value);
+    }
+
     /**
      * Test of polynomial for the linear function.
      */