[ https://issues.apache.org/jira/browse/MATH-1458?page=com.atlassian.jira.plugin.system.issuetabpanels:comment-tabpanel&focusedCommentId=16460872#comment-16460872 ]
Gilles commented on MATH-1458: ------------------------------ Thanks a lot for the report and thorough analysis. The next step would be to write a unit test that displays the bug. And then a patch/PR to fix it. I don't follow the argument about {{Incrementor}} but please note that it is _deprecated_ in the development version of the library: updates/fixes should be performed against the "master" branch in the code repository. > Simpson Integrator computes incorrect value at minimum iterations=1 and > wastes an iteration > ------------------------------------------------------------------------------------------- > > Key: MATH-1458 > URL: https://issues.apache.org/jira/browse/MATH-1458 > Project: Commons Math > Issue Type: Bug > Affects Versions: 3.6.1 > Environment: openjdk version "1.8.0_162" > > OpenJDK Runtime Environment (build 1.8.0_162-8u162-b12-0ubuntu0.16.04.2-b12) > > OpenJDK 64-Bit Server VM (build 25.162-b12, mixed mode) > Reporter: Alex D Herbert > Priority: Minor > Labels: documentation, easyfix, newbie, patch > > org.apache.commons.math3.analysis.integration.SimpsonIntergrator > When used with minimalIterationCount == 1 the integrator computes the wrong > value due to the inlining of computation of stage 1 and stage 0 of the > TrapezoidIntegrator. Each stage is successive since it relies on the result > of the previous stage. So stage 0 must be computed first. This inlining > causes stage 1 to be computed before stage 0: > {code:java} > return (4 * qtrap.stage(this, 1) - qtrap.stage(this, 0)) / 3.0; > {code} > This can be fixed using: > {code:java} > final double s0 = qtrap.stage(this, 0); > return (4 * qtrap.stage(this, 1) - s0) / 3.0;{code} > What is not clear is why setting minimum iterations to 1 results in no > iteration. This is not documented. I would expect setting it to 1 would > compute the first Simpson sum and then perform 1 refinement. This would make > it functionality equivalent to the other Integrator classes which compute two > sums for the first iteration and allow them to be compared if minimum > iterations = 1. If convergence fails then each additional iteration computes > an additional sum. > Note when used with minimalIterationCount > 1 the SimpsonIntegrator wastes a > stage since it computes the following stages: 0, 0, 1, 2, 3. i.e. stage 0 is > computed twice. This is because the iteration is incremented after the stage > is computed: > {code:java} > final double t = qtrap.stage(this, getIterations()); > incrementCount(); > {code} > This should be: > {code:java} > incrementCount(); > final double t = qtrap.stage(this, getIterations()); > {code} > On the first iteration it thus computes the trapezoid sum and compares it to > zero. This would result in a bad computation if integrating a function whose > trapezoid sum is zero (e.g. y=x^3 in the range -1 to 1). However since > iteration only occurs for minimalIterationCount>1 no termination comparison > is made on the first loop. The first termination comparison can be made at > iteration=2 where the comparison will be between the Trapezoid sum and the > first Simpson sum. This is a bug. > However I do not want to submit a formal patch as there is a lack of > consistency across all the integrators in their doIntegrate() method with the > use of incrementCount() and what the iteration number should be at the start > of the while (true) loop: > * IterativeLegendreGauss integrator uses getIterations()+1 to mark the > current iteration inside the loop and calls incrementCount() at the end. > * TrapezoidIntegrator calls incrementCount() outside the loop, uses > getIterations() to mark the current iteration and calls incrementCount() at > the end. > * The MidpointIntegrator calls incrementCount() at the start of the loop and > uses getIterations() to mark the current iteration. > * The RombergIntegrator calls incrementCount() outside the loop, uses > getIterations() to mark the current iteration and calls incrementCount() in > the middle of the loop before termination conditions have been checked. This > allows it to fail when the iterations are equal to the maximum iterations > even if convergence has been achieved (see Note*). > * The SimpsonIntegrator uses getIterations() to mark the current iteration > and calls incrementCount() immediately after. > Note*: This may not be discovered in a unit test since the incrementCount() > uses a backing Incrementor where the Incrementor.increment() method calls > Incrementor.increment(1) which ends up calling canIncrement(0) \{ instead of > canIncrement(nTimes) } to check if the maxCountCallback should be triggered. > I expect that all uses of the Incrementor actually trigger the > maxCountCallback when the count has actually exceeded the maximalCount. I > don't want to get into debugging that class since it also has an iterator > using hasNext() with a call to canIncrement(0) and I do not know the contract > that the iterator is working under. > A consistent approach would be: > * Compute the first sum before the loop > * Enter the loop and increment the iteration (so the first loop execution > would be iteration 1) > * Compute the next sum > * Check termination conditions > An example for the SimpsonIntegrator is below: > {code:java} > protected double doIntegrate() throws > TooManyEvaluationsException, MaxCountExceededException > { > // This is a modification from the default SimpsonIntegrator. > // That only computed a single iteration if > // getMinimalIterationCount() == 1. > // Simpson's rule requires at least two trapezoid stages. > // So we set the first sum using two trapezoid stages. > 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) > { > // The first iteration is now the first refinement of the sum. > // This matches how the MidPointIntegrator works. > incrementCount(); > final int i = getIterations(); > // 1-stage ahead of the iteration > final double t = qtrap.stage(this, i + 1); > final double s = (4 * t - oldt) / 3.0; > 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())) > { > return s; > } > } > olds = s; > oldt = t; > } > } > {code} > Note: If this method is accepted then the SIMPSON_MAX_ITERATIONS_COUNT must > be reduced by 1 to 63, since the stage method of the TrapezoidIntegrator has > a maximum valid input argument of 64. > -- This message was sent by Atlassian JIRA (v7.6.3#76005)