From issues-return-67546-archive-asf-public=cust-asf.ponee.io@commons.apache.org Tue May 1 15:59:05 2018 Return-Path: X-Original-To: archive-asf-public@cust-asf.ponee.io Delivered-To: archive-asf-public@cust-asf.ponee.io Received: from mail.apache.org (hermes.apache.org [140.211.11.3]) by mx-eu-01.ponee.io (Postfix) with SMTP id D0D17180645 for ; Tue, 1 May 2018 15:59:04 +0200 (CEST) Received: (qmail 41514 invoked by uid 500); 1 May 2018 13:59:03 -0000 Mailing-List: contact issues-help@commons.apache.org; run by ezmlm Precedence: bulk List-Help: List-Unsubscribe: List-Post: List-Id: Reply-To: issues@commons.apache.org Delivered-To: mailing list issues@commons.apache.org Received: (qmail 41502 invoked by uid 99); 1 May 2018 13:59:03 -0000 Received: from pnap-us-west-generic-nat.apache.org (HELO spamd1-us-west.apache.org) (209.188.14.142) by apache.org (qpsmtpd/0.29) with ESMTP; Tue, 01 May 2018 13:59:03 +0000 Received: from localhost (localhost [127.0.0.1]) by spamd1-us-west.apache.org (ASF Mail Server at spamd1-us-west.apache.org) with ESMTP id B6910C6E74 for ; Tue, 1 May 2018 13:59:02 +0000 (UTC) X-Virus-Scanned: Debian amavisd-new at spamd1-us-west.apache.org X-Spam-Flag: NO X-Spam-Score: -109.501 X-Spam-Level: X-Spam-Status: No, score=-109.501 tagged_above=-999 required=6.31 tests=[ENV_AND_HDR_SPF_MATCH=-0.5, KAM_ASCII_DIVIDERS=0.8, RCVD_IN_DNSWL_MED=-2.3, SPF_PASS=-0.001, USER_IN_DEF_SPF_WL=-7.5, USER_IN_WHITELIST=-100] autolearn=disabled Received: from mx1-lw-us.apache.org ([10.40.0.8]) by localhost (spamd1-us-west.apache.org [10.40.0.7]) (amavisd-new, port 10024) with ESMTP id 0YUXkE02TBXG for ; Tue, 1 May 2018 13:59:01 +0000 (UTC) Received: from mailrelay1-us-west.apache.org (mailrelay1-us-west.apache.org [209.188.14.139]) by mx1-lw-us.apache.org (ASF Mail Server at mx1-lw-us.apache.org) with ESMTP id D36EF5F1F0 for ; Tue, 1 May 2018 13:59:00 +0000 (UTC) Received: from jira-lw-us.apache.org (unknown [207.244.88.139]) by mailrelay1-us-west.apache.org (ASF Mail Server at mailrelay1-us-west.apache.org) with ESMTP id 5BEFEE01BF for ; Tue, 1 May 2018 13:59:00 +0000 (UTC) Received: from jira-lw-us.apache.org (localhost [127.0.0.1]) by jira-lw-us.apache.org (ASF Mail Server at jira-lw-us.apache.org) with ESMTP id 0B11A21297 for ; Tue, 1 May 2018 13:59:00 +0000 (UTC) Date: Tue, 1 May 2018 13:59:00 +0000 (UTC) From: "Alex D Herbert (JIRA)" To: issues@commons.apache.org Message-ID: In-Reply-To: References: Subject: [jira] [Updated] (MATH-1458) Simpson Integrator computes incorrect value at minimum iterations=1 and wastes an iteration MIME-Version: 1.0 Content-Type: text/plain; charset=utf-8 Content-Transfer-Encoding: quoted-printable X-JIRA-FingerPrint: 30527f35849b9dde25b450d4833f0394 [ https://issues.apache.org/jira/browse/MATH-1458?page=3Dcom.atlassian= .jira.plugin.system.issuetabpanels:all-tabpanel ] Alex D Herbert updated MATH-1458: --------------------------------- Description:=20 org.apache.commons.math3.analysis.integration.SimpsonIntergrator When used with minimalIterationCount =3D=3D 1 the integrator computes the w= rong value due to the inlining of computation of stage 1 and stage 0 of the= =C2=A0TrapezoidIntegrator. Each stage is successive since it relies on the = result of the previous stage. So stage 0 must be computed first. This inlin= ing 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 =3D 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 iter= ation. 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 fu= nctionality equivalent to the other Integrator classes which compute two su= ms for the first iteration and allow them to be compared if minimum iterati= ons =3D 1. If=C2=A0convergence fails then each additional iteration compute= s an additional sum. Note when used with=C2=A0minimalIterationCount > 1 the SimpsonIntegrator wa= stes a stage since it computes the following stages: 0, 0, 1, 2, 3. i.e. st= age 0 is computed twice. This is because the iteration is incremented after= the stage is computed: {code:java} final double t =3D qtrap.stage(this, getIterations()); incrementCount(); {code} This should be: {code:java} incrementCount(); final double t =3D qtrap.stage(this, getIterations()); {code} On the first iteration it thus computes the trapezoid sum and compares it t= o zero. This would=C2=A0 result in a bad computation if integrating a funct= ion whose trapezoid sum is zero (e.g. y=3Dx^3 in the range -1 to 1). Howeve= r since iteration only occurs for minimalIterationCount>1 no termination co= mparison is made on the first loop. The first termination comparison can be= made at iteration=3D2 where the comparison will be between the Trapezoid s= um 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 consis= tency across all the integrators in their doIntegrate() method with the use= of incrementCount() and what=C2=A0the iteration number should be at the st= art of the while (true) loop: * IterativeLegendreGauss integrator uses getIterations()+1 to mark the cur= rent iteration inside the loop and calls incrementCount() at the end.=C2=A0 * TrapezoidIntegrator calls incrementCount() outside the loop, uses getIte= rations() to mark the current iteration and calls=C2=A0incrementCount() at = the end. * The MidpointIntegrator=C2=A0calls incrementCount() at the start of the l= oop and uses getIterations() to mark the current iteration. * The=C2=A0RombergIntegrator=C2=A0calls incrementCount() outside the loop,= uses getIterations() to mark the current iteration and calls=C2=A0incremen= tCount() in the middle of the loop before termination conditions have been = checked. This allows it to fail when the iterations are equal to the maximu= m iterations even if convergence has been achieved (see Note*). * The SimpsonIntegrator=C2=A0uses getIterations() to mark the current iter= ation and calls=C2=A0incrementCount() 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= =C2=A0Incrementor.increment(1) which ends up calling canIncrement(0) \{ ins= tead of canIncrement(nTimes) } to check if the maxCountCallback should be t= riggered. I expect that all uses of the Incrementor actually trigger the ma= xCountCallback when=C2=A0the 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=C2=A0know 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=20 =C2=A0 TooManyEvaluationsException, MaxCountExceededException { // This is a modification from the default SimpsonIntegrator. // That only computed a single iteration if=20 // getMinimalIterationCount() =3D=3D 1. // Simpson's rule requires at least two trapezoid stages. // So we set the first sum using two trapezoid stages. TrapezoidIntegrator qtrap =3D new TrapezoidIntegrator(); =C2=A0 final double s0 =3D qtrap.stage(this, 0); double oldt =3D qtrap.stage(this, 1); =C2=A0 double olds =3D (4 * oldt - s0) / 3.0; =C2=A0 while (true) =C2=A0 { =C2=A0 // The first iteration is now the first refinement of the sum. =C2=A0 // This matches how the MidPointIntegrator works. =C2=A0 incrementCount(); =C2=A0 final int i =3D getIterations(); =C2=A0 // 1-stage ahead of the iteration =C2=A0 final double t =3D qtrap.stage(this, i + 1); =C2=A0 final double s =3D (4 * t - oldt) / 3.0; =C2=A0 if (i >=3D getMinimalIterationCount()) =C2=A0 { =C2=A0 final double delta =3D FastMath.abs(s - olds); =C2=A0 final double rLimit =3D getRelativeAccuracy() *=20 =C2=A0 =C2=A0 =C2=A0 =C2=A0 (FastMath.abs(olds) + FastMath.abs(s)) * 0.5; =C2=A0 if ((delta <=3D rLimit) || (delta <=3D getAbsoluteAccuracy())) =C2=A0 { =C2=A0 return s; =C2=A0 } =C2=A0 } =C2=A0 olds =3D s; =C2=A0 oldt =3D t; =C2=A0 } } {code} Note: If this method is accepted then the=C2=A0SIMPSON_MAX_ITERATIONS_COUNT= must be reduced by 1 to 63, since the stage method of the=C2=A0TrapezoidIn= tegrator has a maximum valid input argument of 64. =C2=A0 was: org.apache.commons.math3.analysis.integration.SimpsonIntergrator When used with minimalIterationCount =3D=3D 1 the integrator computes the w= rong value due to the inlining of computation of stage 1 and stage 0 of the= =C2=A0TrapezoidIntegrator. Each stage is successive since it relies on the = result of the previous stage. So stage 0 must be computed first. This inlin= ing 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 =3D 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 iter= ation. 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 fu= nctionality equivalent to the other Integrator classes which compute two su= ms for the first iteration and allow them to be compared if minimum iterati= ons =3D 1. If=C2=A0convergence fails then each additional iteration compute= s an additional sum. Note when used with=C2=A0minimalIterationCount > 1 the SimpsonIntegrator wa= stes a stage since it computes the following stages: 0, 0, 1, 2, 3. i.e. st= age 0 is computed twice. This is because the iteration is incremented after= the stage is computed: {code:java} final double t =3D qtrap.stage(this, getIterations()); incrementCount(); {code} This should be: {code:java} incrementCount(); final double t =3D qtrap.stage(this, getIterations()); {code} On the first iteration it thus computes the trapezoid sum and compares it t= o zero. This would=C2=A0 result in a bad computation if integrating a funct= ion whose trapezoid sum is zero (e.g. y=3Dx^3 in the range -1 to 1). Howeve= r since iteration only occurs for minimalIterationCount>1 no termination co= mparison is made on the first loop. The first termination comparison can be= made at iteration=3D2 where the comparison will be between the Trapezoid s= um 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 consis= tency across all the integrators in their doIntegrate() method with the use= of incrementCount() and what=C2=A0the iteration number should be at the st= art of the while (true) loop: * IterativeLegendreGauss integrator uses getIterations()+1 to mark the cur= rent iteration inside the loop and calls incrementCount() at the end.=C2=A0 * TrapezoidIntegrator calls incrementCount() outside the loop, uses getIte= rations() to mark the current iteration and calls=C2=A0incrementCount() at = the end. * The MidpointIntegrator=C2=A0calls incrementCount() at the start of the l= oop and uses getIterations() to mark the current iteration. * The=C2=A0RombergIntegrator=C2=A0calls incrementCount() outside the loop,= uses getIterations() to mark the current iteration and calls=C2=A0incremen= tCount() in the middle of the loop before termination conditions have been = checked. This allows it to fail when the iterations are equal to the maximu= m iterations even if convergence has been achieved*. * The SimpsonIntegrator=C2=A0uses getIterations() to mark the current iter= ation and calls=C2=A0incrementCount() immediately after. * This may not be discovered in a unit test since the incrementCount() uses= a backing Incrementor where the Incrementor.increment() method calls=C2=A0= Incrementor.increment(1) which ends up calling canIncrement(0) \{ instead o= f canIncrement(nTimes) } to check if the maxCountCallback should be trigger= ed. I expect that all uses of the Incrementor actually trigger the maxCount= Callback when=C2=A0the 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=C2=A0know the contra= ct 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=20 =C2=A0 TooManyEvaluationsException, MaxCountExceededException { // This is a modification from the default SimpsonIntegrator. // That only computed a single iteration if=20 // getMinimalIterationCount() =3D=3D 1. // Simpson's rule requires at least two trapezoid stages. // So we set the first sum using two trapezoid stages. TrapezoidIntegrator qtrap =3D new TrapezoidIntegrator(); =C2=A0 final double s0 =3D qtrap.stage(this, 0); double oldt =3D qtrap.stage(this, 1); =C2=A0 double olds =3D (4 * oldt - s0) / 3.0; =C2=A0 while (true) =C2=A0 { =C2=A0 // The first iteration is now the first refinement of the sum. =C2=A0 // This matches how the MidPointIntegrator works. =C2=A0 incrementCount(); =C2=A0 final int i =3D getIterations(); =C2=A0 // 1-stage ahead of the iteration =C2=A0 final double t =3D qtrap.stage(this, i + 1); =C2=A0 final double s =3D (4 * t - oldt) / 3.0; =C2=A0 if (i >=3D getMinimalIterationCount()) =C2=A0 { =C2=A0 final double delta =3D FastMath.abs(s - olds); =C2=A0 final double rLimit =3D getRelativeAccuracy() *=20 =C2=A0 =C2=A0 =C2=A0 =C2=A0 (FastMath.abs(olds) + FastMath.abs(s)) * 0.5; =C2=A0 if ((delta <=3D rLimit) || (delta <=3D getAbsoluteAccuracy())) =C2=A0 { =C2=A0 return s; =C2=A0 } =C2=A0 } =C2=A0 olds =3D s; =C2=A0 oldt =3D t; =C2=A0 } } {code} Note: If this method is accepted then the=C2=A0SIMPSON_MAX_ITERATIONS_COUNT= must be reduced by 1 to 63, since the stage method of the=C2=A0TrapezoidIn= tegrator has a maximum valid input argument of 64. =C2=A0 > Simpson Integrator computes incorrect value at minimum iterations=3D1 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" =C2=A0=C2=A0=C2=A0=C2=A0= =C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2= =A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0= =C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2= =A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0= =C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2= =A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0= =C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2= =A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0= =C2=A0=C2=A0 > OpenJDK Runtime Environment (build 1.8.0_162-8u162-b12-0ubuntu0.16.04.2-b= 12) =C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2= =A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0= =C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2= =A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0= =C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0 > OpenJDK 64-Bit Server VM (build 25.162-b12, mixed mode) =C2=A0=C2=A0=C2= =A0 > Reporter: Alex D Herbert > Priority: Minor > Labels: documentation, easyfix, newbie, patch > > org.apache.commons.math3.analysis.integration.SimpsonIntergrator > When used with minimalIterationCount =3D=3D 1 the integrator computes the= wrong value due to the inlining of computation of stage 1 and stage 0 of t= he=C2=A0TrapezoidIntegrator. Each stage is successive since it relies on th= e result of the previous stage. So stage 0 must be computed first. This inl= ining 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 =3D 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 it= eration. This is not documented. I would expect setting it to 1 would compu= te 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 itera= tions =3D 1. If=C2=A0convergence fails then each additional iteration compu= tes an additional sum. > Note when used with=C2=A0minimalIterationCount > 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 aft= er the stage is computed: > {code:java} > final double t =3D qtrap.stage(this, getIterations()); > incrementCount(); > {code} > This should be: > {code:java} > incrementCount(); > final double t =3D qtrap.stage(this, getIterations()); > {code} > On the first iteration it thus computes the trapezoid sum and compares it= to zero. This would=C2=A0 result in a bad computation if integrating a fun= ction whose trapezoid sum is zero (e.g. y=3Dx^3 in the range -1 to 1). Howe= ver 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=3D2 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 cons= istency across all the integrators in their doIntegrate() method with the u= se of incrementCount() and what=C2=A0the iteration number should be at the = start of the while (true) loop: > * IterativeLegendreGauss integrator uses getIterations()+1 to mark the c= urrent iteration inside the loop and calls incrementCount() at the end.=C2= =A0 > * TrapezoidIntegrator calls incrementCount() outside the loop, uses getI= terations() to mark the current iteration and calls=C2=A0incrementCount() a= t the end. > * The MidpointIntegrator=C2=A0calls incrementCount() at the start of the= loop and uses getIterations() to mark the current iteration. > * The=C2=A0RombergIntegrator=C2=A0calls incrementCount() outside the loo= p, uses getIterations() to mark the current iteration and calls=C2=A0increm= entCount() in the middle of the loop before termination conditions have bee= n checked. This allows it to fail when the iterations are equal to the maxi= mum iterations even if convergence has been achieved (see Note*). > * The SimpsonIntegrator=C2=A0uses getIterations() to mark the current it= eration and calls=C2=A0incrementCount() immediately after. > Note*: This may not be discovered in a unit test since the incrementCount= () uses a backing Incrementor where the Incrementor.increment() method call= s=C2=A0Incrementor.increment(1) which ends up calling canIncrement(0) \{ in= stead of canIncrement(nTimes) } to check if the maxCountCallback should be = triggered. I expect that all uses of the Incrementor actually trigger the m= axCountCallback when=C2=A0the count has actually exceeded the maximalCount.= I don't want to get into debugging that class since it also has an iterato= r using hasNext() with a call to canIncrement(0) and I do not=C2=A0know 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 executio= n would be iteration 1) > * Compute the next sum > * Check termination conditions > An example for the SimpsonIntegrator is below: > {code:java} > protected double doIntegrate() throws=20 > =C2=A0 TooManyEvaluationsException, MaxCountExceededException > { > // This is a modification from the default SimpsonIntegrator. > // That only computed a single iteration if=20 > // getMinimalIterationCount() =3D=3D 1. > // Simpson's rule requires at least two trapezoid stages. > // So we set the first sum using two trapezoid stages. > TrapezoidIntegrator qtrap =3D new TrapezoidIntegrator(); > =C2=A0 final double s0 =3D qtrap.stage(this, 0); > double oldt =3D qtrap.stage(this, 1); > =C2=A0 double olds =3D (4 * oldt - s0) / 3.0; > =C2=A0 while (true) > =C2=A0 { > =C2=A0 // The first iteration is now the first refinement of the sum. > =C2=A0 // This matches how the MidPointIntegrator works. > =C2=A0 incrementCount(); > =C2=A0 final int i =3D getIterations(); > =C2=A0 // 1-stage ahead of the iteration > =C2=A0 final double t =3D qtrap.stage(this, i + 1); > =C2=A0 final double s =3D (4 * t - oldt) / 3.0; > =C2=A0 if (i >=3D getMinimalIterationCount()) > =C2=A0 { > =C2=A0 final double delta =3D FastMath.abs(s - olds); > =C2=A0 final double rLimit =3D getRelativeAccuracy() *=20 > =C2=A0 =C2=A0 =C2=A0 =C2=A0 (FastMath.abs(olds) + FastMath.abs(s)) * 0.5; > =C2=A0 if ((delta <=3D rLimit) || (delta <=3D getAbsoluteAccuracy())) > =C2=A0 { > =C2=A0 return s; > =C2=A0 } > =C2=A0 } > =C2=A0 olds =3D s; > =C2=A0 oldt =3D t; > =C2=A0 } > } > {code} > Note: If this method is accepted then the=C2=A0SIMPSON_MAX_ITERATIONS_COU= NT must be reduced by 1 to 63, since the stage method of the=C2=A0Trapezoid= Integrator has a maximum valid input argument of 64. > =C2=A0 -- This message was sent by Atlassian JIRA (v7.6.3#76005)