org.apache.commons.math3.analysis.polynomials.PolynomialFunction Java Examples

@Test
public void testInterpolateLinear() {
    double x[] = { 0.0, 0.5, 1.0 };
    double y[] = { 0.0, 0.5, 0.0 };
    UnivariateInterpolator i = new SplineInterpolator();
    UnivariateFunction f = i.interpolate(x, y);
    verifyInterpolation(f, x, y);
    verifyConsistency((PolynomialSplineFunction) f, x);

    // Verify coefficients using analytical values
    PolynomialFunction polynomials[] = ((PolynomialSplineFunction) f).getPolynomials();
    double target[] = {y[0], 1.5d, 0d, -2d};
    TestUtils.assertEquals(polynomials[0].getCoefficients(), target, coefficientTolerance);
    target = new double[]{y[1], 0d, -3d, 2d};
    TestUtils.assertEquals(polynomials[1].getCoefficients(), target, coefficientTolerance);
}
 

@Test
public void testSmallError() {
    Random randomizer = new Random(53882150042l);
    double maxError = 0;
    for (int degree = 0; degree < 10; ++degree) {
        PolynomialFunction p = buildRandomPolynomial(degree, randomizer);

        PolynomialFitter fitter = new PolynomialFitter(new LevenbergMarquardtOptimizer());
        for (double x = -1.0; x < 1.0; x += 0.01) {
            fitter.addObservedPoint(1.0, x,
                                    p.value(x) + 0.1 * randomizer.nextGaussian());
        }

        final double[] init = new double[degree + 1];
        PolynomialFunction fitted = new PolynomialFunction(fitter.fit(init));

        for (double x = -1.0; x < 1.0; x += 0.01) {
            double error = FastMath.abs(p.value(x) - fitted.value(x)) /
                          (1.0 + FastMath.abs(p.value(x)));
            maxError = FastMath.max(maxError, error);
            Assert.assertTrue(FastMath.abs(error) < 0.1);
        }
    }
    Assert.assertTrue(maxError > 0.01);
}
 

@Test
public void testInterpolateLinear() {
    double x[] = { 0.0, 0.5, 1.0 };
    double y[] = { 0.0, 0.5, 0.0 };
    UnivariateInterpolator i = new SplineInterpolator();
    UnivariateFunction f = i.interpolate(x, y);
    verifyInterpolation(f, x, y);
    verifyConsistency((PolynomialSplineFunction) f, x);

    // Verify coefficients using analytical values
    PolynomialFunction polynomials[] = ((PolynomialSplineFunction) f).getPolynomials();
    double target[] = {y[0], 1.5d, 0d, -2d};
    TestUtils.assertEquals(polynomials[0].getCoefficients(), target, coefficientTolerance);
    target = new double[]{y[1], 0d, -3d, 2d};
    TestUtils.assertEquals(polynomials[1].getCoefficients(), target, coefficientTolerance);
}
 

@Test
public void testInterpolateLinear() {
    double x[] = { 0.0, 0.5, 1.0 };
    double y[] = { 0.0, 0.5, 0.0 };
    UnivariateInterpolator i = new SplineInterpolator();
    UnivariateFunction f = i.interpolate(x, y);
    verifyInterpolation(f, x, y);
    verifyConsistency((PolynomialSplineFunction) f, x);

    // Verify coefficients using analytical values
    PolynomialFunction polynomials[] = ((PolynomialSplineFunction) f).getPolynomials();
    double target[] = {y[0], 1.5d, 0d, -2d};
    TestUtils.assertEquals(polynomials[0].getCoefficients(), target, coefficientTolerance);
    target = new double[]{y[1], 0d, -3d, 2d};
    TestUtils.assertEquals(polynomials[1].getCoefficients(), target, coefficientTolerance);
}
 

@Test
public void testMixedDerivatives() {
    HermiteInterpolator interpolator = new HermiteInterpolator();
    interpolator.addSamplePoint(0.0, new double[] { 1.0 }, new double[] { 2.0 });
    interpolator.addSamplePoint(1.0, new double[] { 4.0 });
    interpolator.addSamplePoint(2.0, new double[] { 5.0 }, new double[] { 2.0 });
    Assert.assertEquals(4, interpolator.getPolynomials()[0].degree());
    DerivativeStructure y0 = interpolator.value(new DerivativeStructure(1, 1, 0, 0.0))[0];
    Assert.assertEquals(1.0, y0.getValue(), 1.0e-15);
    Assert.assertEquals(2.0, y0.getPartialDerivative(1), 1.0e-15);
    Assert.assertEquals(4.0, interpolator.value(1.0)[0], 1.0e-15);
    DerivativeStructure y2 = interpolator.value(new DerivativeStructure(1, 1, 0, 2.0))[0];
    Assert.assertEquals(5.0, y2.getValue(), 1.0e-15);
    Assert.assertEquals(2.0, y2.getPartialDerivative(1), 1.0e-15);
    checkPolynomial(new PolynomialFunction(new double[] { 1.0, 2.0, 4.0, -4.0, 1.0 }),
                    interpolator.getPolynomials()[0]);
}
 

@Test
public void testWikipedia() {
    // this test corresponds to the example from Wikipedia page:
    // http://en.wikipedia.org/wiki/Hermite_interpolation
    HermiteInterpolator interpolator = new HermiteInterpolator();
    interpolator.addSamplePoint(-1, new double[] { 2 }, new double[] { -8 }, new double[] { 56 });
    interpolator.addSamplePoint( 0, new double[] { 1 }, new double[] {  0 }, new double[] {  0 });
    interpolator.addSamplePoint( 1, new double[] { 2 }, new double[] {  8 }, new double[] { 56 });
    for (double x = -1.0; x <= 1.0; x += 0.125) {
        DerivativeStructure y = interpolator.value(new DerivativeStructure(1, 1, 0, x))[0];
        double x2 = x * x;
        double x4 = x2 * x2;
        double x8 = x4 * x4;
        Assert.assertEquals(x8 + 1, y.getValue(), 1.0e-15);
        Assert.assertEquals(8 * x4 * x2 * x, y.getPartialDerivative(1), 1.0e-15);
    }
    checkPolynomial(new PolynomialFunction(new double[] { 1, 0, 0, 0, 0, 0, 0, 0, 1 }),
                    interpolator.getPolynomials()[0]);
}
 

@Test
public void testFit() {
    final RealDistribution rng = new UniformRealDistribution(-100, 100);
    rng.reseedRandomGenerator(64925784252L);

    final LevenbergMarquardtOptimizer optim = new LevenbergMarquardtOptimizer();
    final PolynomialFitter fitter = new PolynomialFitter(optim);
    final double[] coeff = { 12.9, -3.4, 2.1 }; // 12.9 - 3.4 x + 2.1 x^2
    final PolynomialFunction f = new PolynomialFunction(coeff);

    // Collect data from a known polynomial.
    for (int i = 0; i < 100; i++) {
        final double x = rng.sample();
        fitter.addObservedPoint(x, f.value(x));
    }

    // Start fit from initial guesses that are far from the optimal values.
    final double[] best = fitter.fit(new double[] { -1e-20, 3e15, -5e25 });

    TestUtils.assertEquals("best != coeff", coeff, best, 1e-12);
}
 

@Test
public void testMixedDerivatives() {
    HermiteInterpolator interpolator = new HermiteInterpolator();
    interpolator.addSamplePoint(0.0, new double[] { 1.0 }, new double[] { 2.0 });
    interpolator.addSamplePoint(1.0, new double[] { 4.0 });
    interpolator.addSamplePoint(2.0, new double[] { 5.0 }, new double[] { 2.0 });
    Assert.assertEquals(4, interpolator.getPolynomials()[0].degree());
    DerivativeStructure y0 = interpolator.value(new DerivativeStructure(1, 1, 0, 0.0))[0];
    Assert.assertEquals(1.0, y0.getValue(), 1.0e-15);
    Assert.assertEquals(2.0, y0.getPartialDerivative(1), 1.0e-15);
    Assert.assertEquals(4.0, interpolator.value(1.0)[0], 1.0e-15);
    DerivativeStructure y2 = interpolator.value(new DerivativeStructure(1, 1, 0, 2.0))[0];
    Assert.assertEquals(5.0, y2.getValue(), 1.0e-15);
    Assert.assertEquals(2.0, y2.getPartialDerivative(1), 1.0e-15);
    checkPolynomial(new PolynomialFunction(new double[] { 1.0, 2.0, 4.0, -4.0, 1.0 }),
                    interpolator.getPolynomials()[0]);
}
 

@Test
public void testLargeSample() {
    Random randomizer = new Random(0x5551480dca5b369bl);
    double maxError = 0;
    for (int degree = 0; degree < 10; ++degree) {
        PolynomialFunction p = buildRandomPolynomial(degree, randomizer);

        PolynomialFitter fitter = new PolynomialFitter(new LevenbergMarquardtOptimizer());
        for (int i = 0; i < 40000; ++i) {
            double x = -1.0 + i / 20000.0;
            fitter.addObservedPoint(1.0, x,
                                    p.value(x) + 0.1 * randomizer.nextGaussian());
        }

        final double[] init = new double[degree + 1];
        PolynomialFunction fitted = new PolynomialFunction(fitter.fit(init));

        for (double x = -1.0; x < 1.0; x += 0.01) {
            double error = FastMath.abs(p.value(x) - fitted.value(x)) /
                          (1.0 + FastMath.abs(p.value(x)));
            maxError = FastMath.max(maxError, error);
            Assert.assertTrue(FastMath.abs(error) < 0.01);
        }
    }
    Assert.assertTrue(maxError > 0.001);
}
 

/**
 * Test of solver for the quadratic function.
 */
@Test
public void testQuadraticFunction() {
    double min, max, expected, result, tolerance;

    // p(x) = 2x^2 + 5x - 3 = (x+3)(2x-1)
    double coefficients[] = { -3.0, 5.0, 2.0 };
    PolynomialFunction f = new PolynomialFunction(coefficients);
    LaguerreSolver solver = new LaguerreSolver();

    min = 0.0; max = 2.0; expected = 0.5;
    tolerance = FastMath.max(solver.getAbsoluteAccuracy(),
                FastMath.abs(expected * solver.getRelativeAccuracy()));
    result = solver.solve(100, f, min, max);
    Assert.assertEquals(expected, result, tolerance);

    min = -4.0; max = -1.0; expected = -3.0;
    tolerance = FastMath.max(solver.getAbsoluteAccuracy(),
                FastMath.abs(expected * solver.getRelativeAccuracy()));
    result = solver.solve(100, f, min, max);
    Assert.assertEquals(expected, result, tolerance);
}
 

@Test
public void testSmallError() {
    Random randomizer = new Random(53882150042l);
    double maxError = 0;
    for (int degree = 0; degree < 10; ++degree) {
        PolynomialFunction p = buildRandomPolynomial(degree, randomizer);

        PolynomialFitter fitter = new PolynomialFitter(new LevenbergMarquardtOptimizer());
        for (double x = -1.0; x < 1.0; x += 0.01) {
            fitter.addObservedPoint(1.0, x,
                                    p.value(x) + 0.1 * randomizer.nextGaussian());
        }

        final double[] init = new double[degree + 1];
        PolynomialFunction fitted = new PolynomialFunction(fitter.fit(init));

        for (double x = -1.0; x < 1.0; x += 0.01) {
            double error = FastMath.abs(p.value(x) - fitted.value(x)) /
                          (1.0 + FastMath.abs(p.value(x)));
            maxError = FastMath.max(maxError, error);
            Assert.assertTrue(FastMath.abs(error) < 0.1);
        }
    }
    Assert.assertTrue(maxError > 0.01);
}
 

@Test
public void testInterpolateLinearDegenerateThreeSegment()
    {
    double x[] = { 0.0, 0.5, 1.0, 1.5 };
    double y[] = { 0.0, 0.5, 1.0, 1.5 };
    UnivariateInterpolator i = new SplineInterpolator();
    UnivariateFunction f = i.interpolate(x, y);
    verifyInterpolation(f, x, y);

    // Verify coefficients using analytical values
    PolynomialFunction polynomials[] = ((PolynomialSplineFunction) f).getPolynomials();
    double target[] = {y[0], 1d};
    TestUtils.assertEquals(polynomials[0].getCoefficients(), target, coefficientTolerance);
    target = new double[]{y[1], 1d};
    TestUtils.assertEquals(polynomials[1].getCoefficients(), target, coefficientTolerance);
    target = new double[]{y[2], 1d};
    TestUtils.assertEquals(polynomials[2].getCoefficients(), target, coefficientTolerance);

    // Check interpolation
    Assert.assertEquals(0,f.value(0), interpolationTolerance);
    Assert.assertEquals(1.4,f.value(1.4), interpolationTolerance);
    Assert.assertEquals(1.5,f.value(1.5), interpolationTolerance);
}
 

@Test
public void testInterpolateLinearDegenerateThreeSegment()
    {
    double x[] = { 0.0, 0.5, 1.0, 1.5 };
    double y[] = { 0.0, 0.5, 1.0, 1.5 };
    UnivariateInterpolator i = new SplineInterpolator();
    UnivariateFunction f = i.interpolate(x, y);
    verifyInterpolation(f, x, y);

    // Verify coefficients using analytical values
    PolynomialFunction polynomials[] = ((PolynomialSplineFunction) f).getPolynomials();
    double target[] = {y[0], 1d};
    TestUtils.assertEquals(polynomials[0].getCoefficients(), target, coefficientTolerance);
    target = new double[]{y[1], 1d};
    TestUtils.assertEquals(polynomials[1].getCoefficients(), target, coefficientTolerance);
    target = new double[]{y[2], 1d};
    TestUtils.assertEquals(polynomials[2].getCoefficients(), target, coefficientTolerance);

    // Check interpolation
    Assert.assertEquals(0,f.value(0), interpolationTolerance);
    Assert.assertEquals(1.4,f.value(1.4), interpolationTolerance);
    Assert.assertEquals(1.5,f.value(1.5), interpolationTolerance);
}
 

private void checkUnsolvableProblem(MultivariateVectorOptimizer optimizer,
                                    boolean solvable) {
    Random randomizer = new Random(1248788532l);
    for (int degree = 0; degree < 10; ++degree) {
        PolynomialFunction p = buildRandomPolynomial(degree, randomizer);

        PolynomialFitter fitter = new PolynomialFitter(optimizer);

        // reusing the same point over and over again does not bring
        // information, the problem cannot be solved in this case for
        // degrees greater than 1 (but one point is sufficient for
        // degree 0)
        for (double x = -1.0; x < 1.0; x += 0.01) {
            fitter.addObservedPoint(1.0, 0.0, p.value(0.0));
        }

        try {
            final double[] init = new double[degree + 1];
            fitter.fit(init);
            Assert.assertTrue(solvable || (degree == 0));
        } catch(ConvergenceException e) {
            Assert.assertTrue((! solvable) && (degree > 0));
        }
    }
}
 

@Test
public void testSmallError() {
    Random randomizer = new Random(53882150042l);
    double maxError = 0;
    for (int degree = 0; degree < 10; ++degree) {
        PolynomialFunction p = buildRandomPolynomial(degree, randomizer);

        PolynomialFitter fitter = new PolynomialFitter(new LevenbergMarquardtOptimizer());
        for (double x = -1.0; x < 1.0; x += 0.01) {
            fitter.addObservedPoint(1.0, x,
                                    p.value(x) + 0.1 * randomizer.nextGaussian());
        }

        final double[] init = new double[degree + 1];
        PolynomialFunction fitted = new PolynomialFunction(fitter.fit(init));

        for (double x = -1.0; x < 1.0; x += 0.01) {
            double error = FastMath.abs(p.value(x) - fitted.value(x)) /
                          (1.0 + FastMath.abs(p.value(x)));
            maxError = FastMath.max(maxError, error);
            Assert.assertTrue(FastMath.abs(error) < 0.1);
        }
    }
    Assert.assertTrue(maxError > 0.01);
}
 

private void checkUnsolvableProblem(DifferentiableMultivariateVectorOptimizer optimizer,
                                    boolean solvable) {
    Random randomizer = new Random(1248788532l);
    for (int degree = 0; degree < 10; ++degree) {
        PolynomialFunction p = buildRandomPolynomial(degree, randomizer);

        PolynomialFitter fitter = new PolynomialFitter(optimizer);

        // reusing the same point over and over again does not bring
        // information, the problem cannot be solved in this case for
        // degrees greater than 1 (but one point is sufficient for
        // degree 0)
        for (double x = -1.0; x < 1.0; x += 0.01) {
            fitter.addObservedPoint(1.0, 0.0, p.value(0.0));
        }

        try {
            final double[] init = new double[degree + 1];
            fitter.fit(init);
            Assert.assertTrue(solvable || (degree == 0));
        } catch(ConvergenceException e) {
            Assert.assertTrue((! solvable) && (degree > 0));
        }
    }
}
 

@Test
public void testSmallError() {
    Random randomizer = new Random(53882150042l);
    double maxError = 0;
    for (int degree = 0; degree < 10; ++degree) {
        PolynomialFunction p = buildRandomPolynomial(degree, randomizer);

        PolynomialFitter fitter = new PolynomialFitter(new LevenbergMarquardtOptimizer());
        for (double x = -1.0; x < 1.0; x += 0.01) {
            fitter.addObservedPoint(1.0, x,
                                    p.value(x) + 0.1 * randomizer.nextGaussian());
        }

        final double[] init = new double[degree + 1];
        PolynomialFunction fitted = new PolynomialFunction(fitter.fit(init));

        for (double x = -1.0; x < 1.0; x += 0.01) {
            double error = FastMath.abs(p.value(x) - fitted.value(x)) /
                          (1.0 + FastMath.abs(p.value(x)));
            maxError = FastMath.max(maxError, error);
            Assert.assertTrue(FastMath.abs(error) < 0.1);
        }
    }
    Assert.assertTrue(maxError > 0.01);
}
 

@Test
public void testInterpolateLinear() throws Exception {
    double x[] = { 0.0, 0.5, 1.0 };
    double y[] = { 0.0, 0.5, 0.0 };
    UnivariateInterpolator i = new SplineInterpolator();
    UnivariateFunction f = i.interpolate(x, y);
    verifyInterpolation(f, x, y);
    verifyConsistency((PolynomialSplineFunction) f, x);

    // Verify coefficients using analytical values
    PolynomialFunction polynomials[] = ((PolynomialSplineFunction) f).getPolynomials();
    double target[] = {y[0], 1.5d, 0d, -2d};
    TestUtils.assertEquals(polynomials[0].getCoefficients(), target, coefficientTolerance);
    target = new double[]{y[1], 0d, -3d, 2d};
    TestUtils.assertEquals(polynomials[1].getCoefficients(), target, coefficientTolerance);
}
 

private void checkUnsolvableProblem(MultivariateVectorOptimizer optimizer,
                                    boolean solvable) {
    Random randomizer = new Random(1248788532l);
    for (int degree = 0; degree < 10; ++degree) {
        PolynomialFunction p = buildRandomPolynomial(degree, randomizer);

        PolynomialFitter fitter = new PolynomialFitter(optimizer);

        // reusing the same point over and over again does not bring
        // information, the problem cannot be solved in this case for
        // degrees greater than 1 (but one point is sufficient for
        // degree 0)
        for (double x = -1.0; x < 1.0; x += 0.01) {
            fitter.addObservedPoint(1.0, 0.0, p.value(0.0));
        }

        try {
            final double[] init = new double[degree + 1];
            fitter.fit(init);
            Assert.assertTrue(solvable || (degree == 0));
        } catch(ConvergenceException e) {
            Assert.assertTrue((! solvable) && (degree > 0));
        }
    }
}
 

@Test
public void testInterpolateLinearDegenerateTwoSegment()
    {
    double x[] = { 0.0, 0.5, 1.0 };
    double y[] = { 0.0, 0.5, 1.0 };
    UnivariateInterpolator i = new SplineInterpolator();
    UnivariateFunction f = i.interpolate(x, y);
    verifyInterpolation(f, x, y);
    verifyConsistency((PolynomialSplineFunction) f, x);

    // Verify coefficients using analytical values
    PolynomialFunction polynomials[] = ((PolynomialSplineFunction) f).getPolynomials();
    double target[] = {y[0], 1d};
    TestUtils.assertEquals(polynomials[0].getCoefficients(), target, coefficientTolerance);
    target = new double[]{y[1], 1d};
    TestUtils.assertEquals(polynomials[1].getCoefficients(), target, coefficientTolerance);

    // Check interpolation
    Assert.assertEquals(0.0,f.value(0.0), interpolationTolerance);
    Assert.assertEquals(0.4,f.value(0.4), interpolationTolerance);
    Assert.assertEquals(1.0,f.value(1.0), interpolationTolerance);
}
 

@Test
public void testInterpolateLinear() {
    double x[] = { 0.0, 0.5, 1.0 };
    double y[] = { 0.0, 0.5, 0.0 };
    UnivariateInterpolator i = new SplineInterpolator();
    UnivariateFunction f = i.interpolate(x, y);
    verifyInterpolation(f, x, y);
    verifyConsistency((PolynomialSplineFunction) f, x);

    // Verify coefficients using analytical values
    PolynomialFunction polynomials[] = ((PolynomialSplineFunction) f).getPolynomials();
    double target[] = {y[0], 1.5d, 0d, -2d};
    TestUtils.assertEquals(polynomials[0].getCoefficients(), target, coefficientTolerance);
    target = new double[]{y[1], 0d, -3d, 2d};
    TestUtils.assertEquals(polynomials[1].getCoefficients(), target, coefficientTolerance);
}
 

@Test
public void testInterpolateLinearDegenerateTwoSegment()
    {
    double x[] = { 0.0, 0.5, 1.0 };
    double y[] = { 0.0, 0.5, 1.0 };
    UnivariateInterpolator i = new LinearInterpolator();
    UnivariateFunction f = i.interpolate(x, y);
    verifyInterpolation(f, x, y);

    // Verify coefficients using analytical values
    PolynomialFunction polynomials[] = ((PolynomialSplineFunction) f).getPolynomials();
    double target[] = {y[0], 1d};
    TestUtils.assertEquals(polynomials[0].getCoefficients(), target, coefficientTolerance);
    target = new double[]{y[1], 1d};
    TestUtils.assertEquals(polynomials[1].getCoefficients(), target, coefficientTolerance);

    // Check interpolation
    Assert.assertEquals(0.0,f.value(0.0), interpolationTolerance);
    Assert.assertEquals(0.4,f.value(0.4), interpolationTolerance);
    Assert.assertEquals(1.0,f.value(1.0), interpolationTolerance);
}
 

@Test
public void testNoError() {
    Random randomizer = new Random(64925784252l);
    for (int degree = 1; degree < 10; ++degree) {
        PolynomialFunction p = buildRandomPolynomial(degree, randomizer);

        PolynomialFitter fitter = new PolynomialFitter(new LevenbergMarquardtOptimizer());
        for (int i = 0; i <= degree; ++i) {
            fitter.addObservedPoint(1.0, i, p.value(i));
        }

        final double[] init = new double[degree + 1];
        PolynomialFunction fitted = new PolynomialFunction(fitter.fit(init));

        for (double x = -1.0; x < 1.0; x += 0.01) {
            double error = FastMath.abs(p.value(x) - fitted.value(x)) /
                           (1.0 + FastMath.abs(p.value(x)));
            Assert.assertEquals(0.0, error, 1.0e-6);
        }
    }
}
 

/** Compute the interpolation polynomials.
 * @return interpolation polynomials array
 * @exception NoDataException if sample is empty
 */
public PolynomialFunction[] getPolynomials()
    throws NoDataException {

    // safety check
    checkInterpolation();

    // iteration initialization
    final PolynomialFunction zero = polynomial(0);
    PolynomialFunction[] polynomials = new PolynomialFunction[topDiagonal.get(0).length];
    for (int i = 0; i < polynomials.length; ++i) {
        polynomials[i] = zero;
    }
    PolynomialFunction coeff = polynomial(1);

    // build the polynomials by iterating on the top diagonal of the divided differences array
    for (int i = 0; i < topDiagonal.size(); ++i) {
        double[] tdi = topDiagonal.get(i);
        for (int k = 0; k < polynomials.length; ++k) {
            polynomials[k] = polynomials[k].add(coeff.multiply(polynomial(tdi[k])));
        }
        coeff = coeff.multiply(polynomial(-abscissae.get(i), 1.0));
    }

    return polynomials;

}
 

@Test
public void testCompose() {
    double[] epsilon = new double[] { 1.0e-20, 5.0e-14, 2.0e-13, 3.0e-13, 2.0e-13, 1.0e-20 };
    PolynomialFunction poly =
            new PolynomialFunction(new double[] { 1.0, 2.0, 3.0, 4.0, 5.0, 6.0 });
    for (int maxOrder = 0; maxOrder < 6; ++maxOrder) {
        PolynomialFunction[] p = new PolynomialFunction[maxOrder + 1];
        p[0] = poly;
        for (int i = 1; i <= maxOrder; ++i) {
            p[i] = p[i - 1].polynomialDerivative();
        }
        for (double x = 0.1; x < 1.2; x += 0.001) {
            DerivativeStructure dsX = new DerivativeStructure(1, maxOrder, 0, x);
            DerivativeStructure dsY1 = dsX.getField().getZero();
            for (int i = poly.degree(); i >= 0; --i) {
                dsY1 = dsY1.multiply(dsX).add(poly.getCoefficients()[i]);
            }
            double[] f = new double[maxOrder + 1];
            for (int i = 0; i < f.length; ++i) {
                f[i] = p[i].value(x);
            }
            DerivativeStructure dsY2 = dsX.compose(f);
            DerivativeStructure zero = dsY1.subtract(dsY2);
            for (int n = 0; n <= maxOrder; ++n) {
                Assert.assertEquals(0.0, zero.getPartialDerivative(n), epsilon[n]);
            }
        }
    }
}
 

private PolynomialFunction buildRandomPolynomial(int degree, Random randomizer) {
    final double[] coefficients = new double[degree + 1];
    for (int i = 0; i <= degree; ++i) {
        coefficients[i] = randomizer.nextGaussian();
    }
    return new PolynomialFunction(coefficients);
}
 

private double exactIntegration(PolynomialFunction p, double a, double b) {
    final double[] coeffs = p.getCoefficients();
    double yb = coeffs[coeffs.length - 1] / coeffs.length;
    double ya = yb;
    for (int i = coeffs.length - 2; i >= 0; --i) {
        yb = yb * b + coeffs[i] / (i + 1);
        ya = ya * a + coeffs[i] / (i + 1);
    }
    return yb * b - ya * a;
}
 

@Test
public void testExactIntegration() {
    Random random = new Random(86343623467878363l);
    for (int n = 2; n < 6; ++n) {
        LegendreGaussIntegrator integrator =
            new LegendreGaussIntegrator(n,
                                        BaseAbstractUnivariateIntegrator.DEFAULT_RELATIVE_ACCURACY,
                                        BaseAbstractUnivariateIntegrator.DEFAULT_ABSOLUTE_ACCURACY,
                                        BaseAbstractUnivariateIntegrator.DEFAULT_MIN_ITERATIONS_COUNT,
                                        64);

        // an n points Gauss-Legendre integrator integrates 2n-1 degree polynoms exactly
        for (int degree = 0; degree <= 2 * n - 1; ++degree) {
            for (int i = 0; i < 10; ++i) {
                double[] coeff = new double[degree + 1];
                for (int k = 0; k < coeff.length; ++k) {
                    coeff[k] = 2 * random.nextDouble() - 1;
                }
                PolynomialFunction p = new PolynomialFunction(coeff);
                double result    = integrator.integrate(10000, p, -5.0, 15.0);
                double reference = exactIntegration(p, -5.0, 15.0);
                Assert.assertEquals(n + " " + degree + " " + i, reference, result, 1.0e-12 * (1.0 + FastMath.abs(reference)));
            }
        }

    }
}
 

private PolynomialFunction buildRandomPolynomial(int degree, Random randomizer) {
    final double[] coefficients = new double[degree + 1];
    for (int i = 0; i <= degree; ++i) {
        coefficients[i] = randomizer.nextGaussian();
    }
    return new PolynomialFunction(coefficients);
}
 

/** Compute the interpolation polynomials.
 * @return interpolation polynomials array
 * @exception NoDataException if sample is empty
 */
public PolynomialFunction[] getPolynomials()
    throws NoDataException {

    // safety check
    checkInterpolation();

    // iteration initialization
    final PolynomialFunction zero = polynomial(0);
    PolynomialFunction[] polynomials = new PolynomialFunction[topDiagonal.get(0).length];
    for (int i = 0; i < polynomials.length; ++i) {
        polynomials[i] = zero;
    }
    PolynomialFunction coeff = polynomial(1);

    // build the polynomials by iterating on the top diagonal of the divided differences array
    for (int i = 0; i < topDiagonal.size(); ++i) {
        double[] tdi = topDiagonal.get(i);
        for (int k = 0; k < polynomials.length; ++k) {
            polynomials[k] = polynomials[k].add(coeff.multiply(polynomial(tdi[k])));
        }
        coeff = coeff.multiply(polynomial(-abscissae.get(i), 1.0));
    }

    return polynomials;

}