[ https://issues.apache.org/jira/browse/MATH1458?page=com.atlassian.jira.plugin.system.issuetabpanels:alltabpanel
]
Alex D Herbert updated MATH1458:

Description:
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();
// 1stage 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.
was:
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*.
* The SimpsonIntegrator uses getIterations() to mark the current iteration and calls incrementCount()
immediately after.
* 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();
// 1stage 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.
> Simpson Integrator computes incorrect value at minimum iterations=1 and wastes an iteration
> 
>
> Key: MATH1458
> URL: https://issues.apache.org/jira/browse/MATH1458
> 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_1628u162b120ubuntu0.16.04.2b12)
> OpenJDK 64Bit Server VM (build 25.162b12, 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();
> // 1stage 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)
