Merge branch 'fix_MATH-1458'

Closes #85

src/main/java/org/apache/commons/math4/analysis/integration/SimpsonIntegrator.java

@@ -37,7 +37,7 @@
 public class SimpsonIntegrator extends BaseAbstractUnivariateIntegrator {
 
     /** Maximal number of iterations for Simpson. */
-    public static final int SIMPSON_MAX_ITERATIONS_COUNT = 64;
+    public static final int SIMPSON_MAX_ITERATIONS_COUNT = 63;
 
     /**
      * Build a Simpson integrator with given accuracies and iterations counts.
@@ -97,33 +97,30 @@ public SimpsonIntegrator() {
 
     /** {@inheritDoc} */
     @Override
-    protected double doIntegrate()
-        throws TooManyEvaluationsException, MaxCountExceededException {
-
-        TrapezoidIntegrator qtrap = new TrapezoidIntegrator();
-        if (getMinimalIterationCount() == 1) {
-            return (4 * qtrap.stage(this, 1) - qtrap.stage(this, 0)) / 3.0;
-        }
-
+    protected double doIntegrate() {
         // Simpson's rule requires at least two trapezoid stages.
-        double olds = 0;
-        double oldt = qtrap.stage(this, 0);
+        // So we set the first sum using two trapezoid stages.
+        final TrapezoidIntegrator qtrap = new TrapezoidIntegrator();
+
+        final double s0 = qtrap.stage(this, 0);
+        double oldt = qtrap.stage(this, 1);
+        double olds = (4 * oldt - s0) / 3.0;
         while (true) {
-            final double t = qtrap.stage(this, iterations.getCount());
+            // The first iteration is the first refinement of the sum.
             iterations.incrementCount();
+            final int i = getIterations();
+            final double t = qtrap.stage(this, i + 1); // 1-stage ahead of the iteration
             final double s = (4 * t - oldt) / 3.0;
-            if (iterations.getCount() >= getMinimalIterationCount()) {
+            if (i >= getMinimalIterationCount()) {
                 final double delta = FastMath.abs(s - olds);
-                final double rLimit =
-                    getRelativeAccuracy() * (FastMath.abs(olds) + FastMath.abs(s)) * 0.5;
-                if ((delta <= rLimit) || (delta <= getAbsoluteAccuracy())) {
+                final double rLimit = getRelativeAccuracy() * (FastMath.abs(olds) + FastMath.abs(s)) * 0.5;
+                if (delta <= rLimit ||
+                    delta <= getAbsoluteAccuracy()) {
                     return s;
                 }
             }
             olds = s;
             oldt = t;
         }
-
     }
-
 }

src/test/java/org/apache/commons/math4/analysis/integration/SimpsonIntegratorTest.java

@@ -18,6 +18,8 @@
 
 import org.apache.commons.math4.analysis.QuinticFunction;
 import org.apache.commons.math4.analysis.UnivariateFunction;
+import org.apache.commons.math4.analysis.function.Identity;
+import org.apache.commons.math4.analysis.function.Inverse;
 import org.apache.commons.math4.analysis.function.Sin;
 import org.apache.commons.math4.analysis.integration.SimpsonIntegrator;
 import org.apache.commons.math4.analysis.integration.UnivariateIntegrator;
@@ -114,10 +116,272 @@ public void testParameters() {
         }
         try {
             // bad iteration limits
-            new SimpsonIntegrator(10, 99);
+            new SimpsonIntegrator(10, SimpsonIntegrator.SIMPSON_MAX_ITERATIONS_COUNT + 1);
             Assert.fail("Expecting NumberIsTooLargeException - bad iteration limits");
         } catch (NumberIsTooLargeException ex) {
             // expected
         }
     }
+
+    // Tests for MATH-1458:
+    // The SimpsonIntegrator had the following bugs:
+    // - minimalIterationCount==1 results in no possible iteration
+    // - minimalIterationCount==1 computes incorrect Simpson sum (following no iteration)
+    // - minimalIterationCount>1 computes the first iteration sum as the Trapezoid sum
+    // - minimalIterationCount>1 computes the second iteration sum as the first Simpson sum
+
+    /**
+     * Test iteration is possible when minimalIterationCount==1.
+     * <br/>
+     * MATH-1458: No iterations were performed when minimalIterationCount==1.
+     */
+    @Test
+    public void testIterationIsPossibleWhenMinimalIterationCountIs1() {
+        UnivariateFunction f = new Sin();
+        UnivariateIntegrator integrator = new SimpsonIntegrator(1,
+                SimpsonIntegrator.SIMPSON_MAX_ITERATIONS_COUNT);
+        // The range or result is not relevant.
+        // This sum should not converge at 1 iteration.
+        // This tests iteration occurred.
+        integrator.integrate(1000, f, 0, 1);
+        // MATH-1458: No iterations were performed when minimalIterationCount==1
+        Assert.assertTrue("Iteration is not above 1",
+                integrator.getIterations() > 1);
+    }
+
+    /**
+     * Test convergence at iteration 1 when minimalIterationCount==1.
+     * <br/>
+     * MATH-1458: No iterations were performed when minimalIterationCount==1.
+     */
+    @Test
+    public void testConvergenceIsPossibleAtIteration1() {
+    	// A linear function y=x should converge immediately
+        UnivariateFunction f = new Identity();
+        UnivariateIntegrator integrator = new SimpsonIntegrator(1,
+                SimpsonIntegrator.SIMPSON_MAX_ITERATIONS_COUNT);
+
+        double min, max, expected, result, tolerance;
+
+        min = 0; max = 1; expected = 0.5;
+        tolerance = FastMath.abs(expected * integrator.getRelativeAccuracy());
+        result = integrator.integrate(1000, f, min, max);
+        // MATH-1458: No iterations were performed when minimalIterationCount==1
+        Assert.assertTrue("Iteration is not above 0",
+                integrator.getIterations()  > 0);
+        // This should converge immediately
+        Assert.assertEquals("Iteration", integrator.getIterations(), 1);
+        Assert.assertEquals("Result", expected, result, tolerance);
+    }
+
+    /**
+     * Compute the integral using the composite Simpson's rule.
+     *
+     * @param f the function
+     * @param a the lower limit
+     * @param b the upper limit
+     * @param n the number of intervals (must be even)
+     * @return the integral between a and b
+     * @see <a href="https://en.wikipedia.org/wiki/Simpson%27s_rule#Composite_Simpson's_rule">
+     *       Composite_Simpson's_rule</a>
+     */
+    private static double compositeSimpsonsRule(UnivariateFunction f, double a,
+            double b, int n)
+    {
+        // Sum interval [a,b] split into n subintervals, with n an even number:
+        // sum ~ h/3 * [ f(x0) + 4f(x1) + 2f(x2) + 4f(x3) + 2f(x4) ... + 4f(xn-1) + f(xn) ]
+        // h = (b-a)/n
+        // f(xi) = f(a + i*h)
+        assert n > 0 && n % 2 == 0 : "n must be strictly positive and even";
+        final double h = (b - a) / n;
+        double sum4 = 0;
+        double sum2 = 0;
+        for (int i = 1; i < n; i++) {
+            // Alternate sums that are multiplied by 4 and 2
+            final double fxi = f.value(a + i * h);
+            if (i % 2 == 0)
+                sum2 += fxi;
+            else
+                sum4 += fxi;
+        }
+        return (h / 3) * (f.value(a) + 4 * sum4 + 2 * sum2 + f.value(b));
+    }
+
+    /**
+     * Compute the iteration of Simpson's rule.
+     *
+     * @param f the function
+     * @param a the lower limit
+     * @param b the upper limit
+     * @param iteration the refinement iteration
+     * @return the integral between a and b
+     */
+    private static double computeSimpsonIteration(UnivariateFunction f, double a,
+            double b, int iteration)
+    {
+        // The first possible Simpson's sum uses n=2.
+        // The next uses n=4. This is the 1st refinement expected when the
+        // integrator has performed 1 iteration.
+        final int n = 2 << iteration;
+        return compositeSimpsonsRule(f, a, b, n);
+    }
+
+    /**
+     * Test the reference Simpson integration is doing what is expected
+     */
+    @Test
+    public void testReferenceSimpsonItegrationIsCorrect() {
+        UnivariateFunction f = new Sin();
+
+        double a, b, h, expected, result, tolerance;
+
+        a = 0.5;
+        b = 1;
+
+        double b_a = b - a;
+
+        // First Simpson sum. 1 midpoint evaluation:
+        h = b_a / 2;
+        double f00 = f.value(a);
+        double f01 = f.value(a + 1 * h);
+        double f0n = f.value(b);
+        expected = (b_a / 6) * (f00 + 4 * f01 + f0n);
+        tolerance = FastMath.abs(expected * SimpsonIntegrator.DEFAULT_RELATIVE_ACCURACY);
+        result = computeSimpsonIteration(f, a, b, 0);
+        Assert.assertEquals("Result", expected, result, tolerance);
+
+        // Second Simpson sum: 2 more evaluations:
+        h = b_a / 4;
+        double f11 = f.value(a + 1 * h);
+        double f13 = f.value(a + 3 * h);
+        expected = (h / 3) * (f00 + 4 * f11 + 2 * f01 + 4 * f13 + f0n);
+        tolerance = FastMath.abs(expected * SimpsonIntegrator.DEFAULT_RELATIVE_ACCURACY);
+        result = computeSimpsonIteration(f, a, b, 1);
+        Assert.assertEquals("Result", expected, result, tolerance);
+
+        // Third Simpson sum: 4 more evaluations:
+        h = b_a / 8;
+        double f21 = f.value(a + 1 * h);
+        double f23 = f.value(a + 3 * h);
+        double f25 = f.value(a + 5 * h);
+        double f27 = f.value(a + 7 * h);
+        expected = (h / 3) * (f00 + 4 * f21 + 2 * f11 + 4 * f23 + 2 * f01 + 4 * f25 +
+                2 * f13 + 4 * f27 + f0n);
+        tolerance = FastMath.abs(expected * SimpsonIntegrator.DEFAULT_RELATIVE_ACCURACY);
+        result = computeSimpsonIteration(f, a, b, 2);
+        Assert.assertEquals("Result", expected, result, tolerance);
+    }
+
+    /**
+     * Test iteration 1 returns the expected sum when minimalIterationCount==1.
+     * <br/>
+     * MATH-1458: minimalIterationCount==1 computes incorrect Simpson sum
+     * (following no iteration).
+     */
+    @Test
+    public void testIteration1ComputesTheExpectedSimpsonSum() {
+        UnivariateFunction f = new Sin();
+        // Set convergence criteria to force immediate convergence
+        UnivariateIntegrator integrator = new SimpsonIntegrator(
+                0, Double.POSITIVE_INFINITY,
+                1, SimpsonIntegrator.SIMPSON_MAX_ITERATIONS_COUNT);
+        double min, max, expected, result, tolerance;
+
+        // MATH-1458: minimalIterationCount==1 computes incorrect
+        // Simpson sum (following no iteration)
+        min = 0;
+        max = 1;
+        result = integrator.integrate(1000, f, min, max);
+        // Immediate convergence
+        Assert.assertEquals("Iteration", 1, integrator.getIterations());
+
+        // Check the sum is as expected
+        expected = computeSimpsonIteration(f, min, max, 1);
+        tolerance = FastMath.abs(expected * SimpsonIntegrator.DEFAULT_RELATIVE_ACCURACY);
+        Assert.assertEquals("Result", expected, result, tolerance);
+    }
+
+    /**
+     * Test iteration N returns the expected sum when minimalIterationCount==1.
+     * <br/>
+     * MATH-1458: minimalIterationCount>1 computes the second iteration sum as
+     * the first Simpson sum.
+     */
+    @Test
+    public void testIterationNComputesTheExpectedSimpsonSum() {
+        // Use 1/x as the function as the sum will asymptote in a monotonic
+        // series. The convergence can then be controlled.
+        UnivariateFunction f = new Inverse();
+
+        double min, max, expected, result, tolerance;
+        int minIteration, maxIteration;
+
+        // Range for integration
+        min = 1;
+        max = 2;
+
+        // This is the expected sum.
+        // Each iteration will monotonically converge to this.
+        expected = FastMath.log(max) - FastMath.log(min);
+
+        // Test convergence at the given iteration
+        minIteration = 2;
+        maxIteration = 4;
+
+        // Compute the sums expected for different iterations.
+        // Add an additional sum so that the test can compare to the next value.
+        double[] sums = new double[maxIteration + 2];
+        for (int i = 0; i < sums.length; i++) {
+            sums[i] = computeSimpsonIteration(f, min, max, i);
+            // Check monotonic
+            if (i > 0) {
+                Assert.assertTrue("Expected series not monotonic descending",
+                        sums[i] < sums[i - 1]);
+                // Check monotonic difference
+                if (i > 1) {
+                    Assert.assertTrue("Expected convergence not monotonic descending",
+                           sums[i - 1] - sums[i] < sums[i - 2] - sums[i - 1]);
+                }
+            }
+        }
+
+        // Check the test function is correct.
+        tolerance = FastMath.abs(expected * SimpsonIntegrator.DEFAULT_RELATIVE_ACCURACY);
+        Assert.assertEquals("Expected result", expected, sums[maxIteration], tolerance);
+
+        // Set-up to test convergence at a specific iteration.
+        // Allow enough function evaluations.
+        // Iteration 0 = 3 evaluations
+        // Iteration 1 = 5 evaluations
+        // Iteration n = 2^(n+1)+1 evaluations
+        int evaluations = 2 << (maxIteration + 1) + 1;
+
+        for (int i = minIteration; i <= maxIteration; i++) {
+            // Create convergence criteria.
+            // (sum - previous) is monotonic descending.
+            // So use a point half-way between them:
+            // ((sums[i-1] - sums[i]) + (sums[i-2] - sums[i-1])) / 2
+            final double absoluteAccuracy = (sums[i - 2] - sums[i]) / 2;
+
+            // Use minimalIterationCount>1
+            UnivariateIntegrator integrator = new SimpsonIntegrator(
+                    0, absoluteAccuracy,
+                    2, SimpsonIntegrator.SIMPSON_MAX_ITERATIONS_COUNT);
+
+            result = integrator.integrate(evaluations, f, min, max);
+
+            // Check the iteration is as expected
+            Assert.assertEquals("Test failed to control iteration", i, integrator.getIterations());
+
+            // MATH-1458: minimalIterationCount>1 computes incorrect Simpson sum
+            // for the iteration. Check it is the correct sum.
+            // It should be closer to this one than the previous or next.
+            final double dp = FastMath.abs(sums[i-1] - result);
+            final double d  = FastMath.abs(sums[i]   - result);
+            final double dn = FastMath.abs(sums[i+1] - result);
+
+            Assert.assertTrue("Result closer to sum expected from previous iteration: " + i, d < dp);
+            Assert.assertTrue("Result closer to sum expected from next iteration: " + i, d < dn);
+        }
+    }
 }