private void proc(UnivariateFunction fn) {
+ Calc calc = calcs.pop();// pop a calculation to be done.
+ double a = calc.a;
+ double b = calc.b;
+
+ // back off slightly from the edges (function evaluations typically go
haywire)
+ // scale for the edge approximation is set by the resolution.
+ if (a == edgeA) {
+ a += (b - a) * EPS;
+ }
+ if (b == edgeB) {
+ b -= (b - a) * EPS;
+ }
+
+ double delta = b - a;
+ double m = (a + b) / 2d;
+ double ma = (a + m) / 2d;
+ double mb = (b + m) / 2d;
+ double va = fn.value(a);
+ double vb = fn.value(b);
+ double vm = fn.value(m);
+ double Q1 = delta / 6 * (va + 4 * vm + vb);
+ double Q2 = delta / 12 * (va + 4 * fn.value(ma) + 2 * vm + 4 *
fn.value(mb) + vb);
+ if (Math.abs(Q2 - Q1) <= EPS) {
+ total += Q2 + (Q2 - Q1) / 15;
+ } else {
+ calcs.add(new Calc(a, m));
+ calcs.add(new Calc(m, b));
+ }
+ }
1) You form a new object Calc(a,m) and Calc(m,b), the m point is the same, yet
you evaluate it again when you compute vb and va, respectively, one iteration
down. I think it also repeats for a bunch of other points middle points.
2) I am of the opinion you shouldn't include special behavior unless you have
problem. If you want you can capture if eval returned a NaN and then move the
evaluation point.
3) I am not claiming you shouldn't replace what I wrote with a stack. For my
personal use it didn't matter.
On Jul 1, 2013, at 6:43 PM, Ajo Fod <ajo....@gmail.com> wrote:
> 1> Could you provide an example of redundant function evaluations in my
> code?
>
> 2> I made the adjustment in AQ because the size of the segment that was
> excluded from the integration could be determined adaptively to keep the
> error arbitrarily low. Note the size of the segment that is excluded (at
> the edges of the integration domain) can be defined as a function of the
> error tolerance. It is unlikely that that exclusion will cause problems for
> proper integrals. But it does save the day for when infinities are present
> ... not a bad trade off if you think about it.
>
> 3> BTW, I think your code is a lot like mine ... except that it uses the
> program stack ... so it will likely fail at the depth required for
> integration of the PDF with sigma=1000.
>
> Cheers,
> -Ajo
>
>
> On Mon, Jul 1, 2013 at 3:17 PM, Konstantin Berlin <kber...@gmail.com> wrote:
>
>> Just a few thoughts.
>>
>> It seems to me that this code suffers from redundant function evaluations.
>> I am not sure what to think about the movement from edges since it violates
>> proper behavior for well behaved functions in order to work for some
>> special cases. In case of infinite integrals it might be better to create
>> the initial offset there, rather than in the adaptive quadrature. In any
>> case the check should be moved out of the inside iteration.
>>
>> A while back I wrote an adaptive quadrature to evaluation vector
>> functions, some code can easily be modified into this class. Hope it helps.
>>
>> /**
>> * The Class SimpsonVectorQuad.
>> */
>> public final class SimpsonVectorQuad implements UnivariateVectorIntegrator
>> {
>>
>> /** The abs point tolerance. */
>> final double absPointTolerance;
>>
>> /** The real tolerance. */
>> final double realTolerance;
>>
>> /**
>> * Instantiates a new simpson vector quad.
>> *
>> * @param realTolerance
>> * the real tolerance
>> */
>> public SimpsonVectorQuad(double realTolerance)
>> {
>> this(realTolerance,1.0e-8);
>> }
>>
>> /**
>> * Instantiates a new simpson vector quad.
>> *
>> * @param realTolerance
>> * the real tolerance
>> * @param absPointTolerance
>> * the abs point tolerance
>> */
>> public SimpsonVectorQuad(double realTolerance, double
>> absPointTolerance)
>> {
>> this.realTolerance = realTolerance;
>> this.absPointTolerance = absPointTolerance;
>> }
>>
>> /* (non-Javadoc)
>> * @see
>> edu.umd.umiacs.armor.math.integration.UnivariateVectorIntegrator#integrate(edu.umd.umiacs.armor.math.func.UnivariateVectorFunction,
>> double, double, int)
>> */
>> @Override
>> public double[] integrate(UnivariateVectorFunction f, double a,
>> double b, int maxEvals)
>> {
>> IntegratorProperties property = new IntegratorProperties();
>> property.evalCount = 0;
>>
>> double h = 0.13579*(b-a);
>> double[] x = {a, a+h, a+2.0*h, (a+b)*0.5, b-2.0*h, b-h, b};
>>
>> double[][] y = new double[7][];
>> for (int iter=0; iter<x.length; iter++)
>> {
>> y[iter] = f.value(x[iter]);
>>
>> for (int nanLoop=0; nanLoop<y[iter].length; nanLoop++)
>> if (Double.isNaN(y[iter][nanLoop]))
>> throw new MathRuntimeException("Function
>> evaluation return NaN.");
>> }
>> property.evalCount += 7;
>>
>> double[] Q1 = new double[y[0].length];
>> double[] Q2 = new double[y[0].length];
>> double[] Q3 = new double[y[0].length];
>>
>> int depth = 0;
>>
>> quadstep(f,Q1,x[0],x[2],y[0],y[1],y[2],this.realTolerance/3.0,property,maxEvals,depth+1);
>>
>> quadstep(f,Q2,x[2],x[4],y[2],y[3],y[4],this.realTolerance/3.0,property,maxEvals,depth+1);
>>
>> quadstep(f,Q3,x[4],x[6],y[4],y[5],y[6],this.realTolerance/3.0,property,maxEvals,depth+1);
>>
>> //sum up the results
>> for (int iter=0; iter<Q1.length; iter++)
>> Q1[iter] = Q1[iter]+Q2[iter]+Q3[iter];
>>
>> return Q1;
>> }
>>
>> /**
>> * Quadstep.
>> *
>> * @param f
>> * the f
>> * @param Q
>> * the q
>> * @param a
>> * the a
>> * @param b
>> * the b
>> * @param fa
>> * the fa
>> * @param fc
>> * the fc
>> * @param fb
>> * the fb
>> * @param tol
>> * the tol
>> * @param property
>> * the property
>> * @param maxEvals
>> * the max evals
>> * @param depth
>> * the depth
>> */
>> private void quadstep(UnivariateVectorFunction f, double[] Q,
>> double a, double b,
>>
>> double[] fa, double[] fc, double[] fb,
>>
>> double tol, IntegratorProperties property, double
>> maxEvals, int depth)
>> {
>> // Evaluate integrant twice in interior of subinterval
>> [a,b].
>> double h = b - a;
>> double c = (a + b)*0.5;
>> double d = (a + c)*0.5;
>> double e = (c + b)*0.5;
>>
>> double[] fd = f.value(d);
>> for (int nanLoop=0; nanLoop<fd.length; nanLoop++)
>> if (Double.isNaN(fd[nanLoop]))
>> throw new ArithmeticException("Function evaluation
>> returned NaN.");
>>
>>
>> double[] fe = f.value(e);
>> for (int nanLoop=0; nanLoop<fe.length; nanLoop++)
>> if (Double.isNaN(fe[nanLoop]))
>> throw new ArithmeticException("Function evaluation
>> returned NaN.");
>>
>> property.evalCount += 2;
>>
>> double[] Q1 = new double[fd.length];
>> double[] Q2 = new double[fd.length];
>>
>> for (int iter=0; iter<Q.length; iter++)
>> {
>> // Three point Simpson's rule.
>> Q1[iter] = (h/6.0)*(fa[iter] + 4.0*fc[iter] +
>> fb[iter]);
>>
>> // Five point double Simpson's rule.
>> Q2[iter] = (h/12.0)*(fa[iter] + 4.0*fd[iter] +
>> 2.0*fc[iter] + 4.0*fe[iter] + fb[iter]);
>>
>> // One step of Romberg extrapolation.
>> Q[iter] = Q2[iter] + (Q2[iter] - Q1[iter])/15.0;
>> }
>>
>> // Maximum function count exceeded; singularity likely.
>> if (property.evalCount > maxEvals)
>> //return;
>> throw new MathRuntimeException("Integrator
>> exceeded maximum number of function evaluations. Singularity likely.");
>>
>> // Maximum depth count exceeded; singularity likely.
>> if (h<this.absPointTolerance)
>> return;
>>
>> //found the right value
>> boolean allBelow = true;
>> for (int iter=0; iter<Q.length; iter++)
>> if (Math.abs(Q2[iter] - Q[iter]) > tol)
>> allBelow = false;
>>
>> if (allBelow)
>> return;
>>
>>
>> // Subdivide into two subintervals.
>>
>> quadstep(f,Q1,a,c,fa,fd,fc,tol*0.5,property,maxEvals,depth+1);
>>
>> quadstep(f,Q2,c,b,fc,fe,fb,tol*0.5,property,maxEvals,depth+1);
>>
>> for (int iter=0; iter<Q.length; iter++)
>> Q[iter] = Q1[iter] + Q2[iter];
>>
>> return;
>> }
>> }
>>
>>
>> On Jul 1, 2013, at 1:50 PM, Ajo Fod <ajo....@gmail.com> wrote:
>>
>>> If you wanted to use the Math 3 codebase in AdaptiveQuadrature, you'd
>>> compute the calculations of Q1 and Q2 with something else. I'm not
>> entirely
>>> familiar with the apache Math codebase so my guess would be that you can
>>> replace the following line in AdaptiveQuadrature.proc():
>>>
>>> double Q1 = delta / 6 * (va + 4 * vm + vb);
>>>
>>> with a modification adapted from
>> IterativeLegendreGaussIntegrator.stage(n)
>>> line for integration as follows:
>>>
>>> Q1 = delta * FACTORY.legendreHighPrecision(numberOfPoints, a,
>>> b).integrate(fn);
>>> Q2 = delta * FACTORY.legendreHighPrecision(numberOfPoints+1, a,
>>> b).integrate(fn);
>>>
>>> ... or some slight modification of the above.
>>>
>>> Each of the tests in the patch is integrating a UnivariateFunction in
>>> [-1,1]. Infinity.wrap(fn) just provides that UnivariateFunction. [In the
>>> patches for MATH-995 the InfiniteIntegral was replaced by Infinity.wrap()
>>> ]. So, if you are saying that the intent of
>>> IterativeLegendreGaussIntegrator (refered to as LGQ) was not to integrate
>>> this kind of UnivariateFunction in [-1,1], ... what kind of univariate
>>> function would that be? If it is indeed supposed to do the integration,
>>> then AQ clearly does a better job.
>>>
>>> So, why does LGQ fail here? It is probably that the Adaptive division of
>>> the integration domain (as opposed to the uniform division with LGQ)
>> gives
>>> AQ the critical edge. The test you have for LGQ so far are pretty well
>>> behaved.
>>>
>>> Summary: I'm demonstrating a clear bug/inefficiency with LGQ and
>> providing
>>> you with an alternative that is more accurate.
>>>
>>> Cheers,
>>> Ajo.
>>>
>>>
>>> On Mon, Jul 1, 2013 at 8:22 AM, Gilles <gil...@harfang.homelinux.org>
>> wrote:
>>>
>>>> Hi.
>>>>
>>>>
>>>>
>>>>> I just noticed your request to write the algorithm along the lines of
>> the
>>>>> wikipedia article.
>>>>>
>>>>> The only major difference between my code and the article on Wikipedia
>> is
>>>>> that I found it necessary to move the recursive stack in into a data
>>>>> structure to avoid a StackOverflowException when the non polynomial
>>>>> curvature is concentrated in a corner of the domain of integration.
>> Notice
>>>>> that the Stack objects stores a Stack of limits of integration.
>>>>>
>>>>
>>>> There is a misunderstanding: I'm referring to the "high-level"
>>>> description of the algorithm that is the separation of concerns
>>>> between the quadrature method and the adaptive process. Your code
>>>> mixes the two. Moreover, it does not reuse any of the quadrature
>>>> schemes already implemented in CM, but implements a (new?) one
>>>> without any reference or comments.
>>>> [And this is even without delving into remarks concerning the
>>>> code structure itself.]
>>>>
>>>> Additionally, your patch also mixes two concepts: Adaptive
>>>> quadrature vs improper integral (which is also MATH-994); it is
>>>> hard to follow what problem this issue is supposed to point to,
>>>> and how the patch solves it. Indeed your unit tests shows a
>>>> problem with improper integrals which the class
>>>> "**IterativeGaussLegendreIntegrat**or" is _not_ meant to handle.[1]
>>>>
>>>> To be clear, hopefully, you are demonstrating a problem that
>>>> occurs when combining Commons Math code with code which you
>>>> created.
>>>> The first step is to create a unit test demonstrating whether
>>>> an issue exists with "**IterativeGaussLegendreIntegrat**or" code
>>>> only (i.e. without relying on your "InfiniteIntegral" class).[1]
>>>> If no independent issue exist, then MATH-995 should be replaced
>>>> by an appropriate feature request.
>>>> Also, it would certainly be helpful to pinpoint the reason why
>>>> the combination of "**IterativeGaussLegendreIntegrat**or" and
>>>> "InfiniteIntegral" is not legitimate (if that's the case).
>>>>
>>>>
>>>> Regards,
>>>> Gilles
>>>>
>>>> [1] Cf. also my latest comment on the MATH-995 page.
>>>>
>>>>
>>>>
>>>>> Cheers,
>>>>> Ajo.
>>>>>
>>>>>
>>>>> On Fri, Jun 28, 2013 at 11:07 AM, Ajo Fod <ajo....@gmail.com> wrote:
>>>>>
>>>>> BTW, it is possible that I'm not using LGQ correctly. If so, please
>> show
>>>>>> how to pass the tests I've added. I'd much rather use something that
>> is
>>>>>> better tested than my personal code.
>>>>>>
>>>>>> -Ajo.
>>>>>>
>>>>>>
>>>>>> On Fri, Jun 28, 2013 at 11:04 AM, Ajo Fod <ajo....@gmail.com> wrote:
>>>>>>
>>>>>> I just posted a patch on this issue. Feel free to edit as necessary to
>>>>>>> match your standards. There is a clear issue with LGQ.
>>>>>>>
>>>>>>> Cheers,
>>>>>>> Ajo.
>>>>>>>
>>>>>>>
>>>>>>> On Fri, Jun 28, 2013 at 10:54 AM, Gilles <
>> gil...@harfang.homelinux.org>
>>>>>>> **wrote:
>>>>>>>
>>>>>>> Ted,
>>>>>>>>
>>>>>>>>
>>>>>>>>
>>>>>>>> Did you read my other (rather more lengthy) post? Is that
>> "jumping"?
>>>>>>>>>
>>>>>>>>>>
>>>>>>>>>>
>>>>>>>>>> Yes. You jumped on him rather than helped him be productive. The
>>>>>>>>> general
>>>>>>>>> message is "we have something in the works, don't bother us with
>> your
>>>>>>>>> ideas".
>>>>>>>>>
>>>>>>>>>
>>>>>>>> Then please read all the messages pertaining to those issues more
>>>>>>>> carefully:
>>>>>>>> I never wrote such a thing (neither now nor in the past).
>>>>>>>> I pointed to a potential problem in the usage of the CM code.
>>>>>>>> I pointed (several times and in details) to problems in candidate
>>>>>>>> contributions,
>>>>>>>> with arguments that go well beyond "bad formatting".
>>>>>>>> I pointed out how we could improve the functionality _together_
>> (i.e.
>>>>>>>> by
>>>>>>>> using
>>>>>>>> what we have, instead of throwing it out without even trying to
>> figure
>>>>>>>> out how
>>>>>>>> good or bad it is).
>>>>>>>>
>>>>>>>> IMHO, these were all valid suggestions to be productive in helping
>> CM
>>>>>>>> to
>>>>>>>> become
>>>>>>>> better, instead of merely larger. The former indeed requires more
>>>>>>>> effort
>>>>>>>> than
>>>>>>>> the latter.
>>>>>>>>
>>>>>>>>
>>>>>>>>
>>>>>>>> Gilles
>>>>>>>>
>>>>>>>
>>>>
>>>>
>> ------------------------------**------------------------------**---------
>>>> To unsubscribe, e-mail: dev-unsubscribe@commons.**apache.org<
>> dev-unsubscr...@commons.apache.org>
>>>> For additional commands, e-mail: dev-h...@commons.apache.org
>>>>
>>>>
>>
>>
>> ---------------------------------------------------------------------
>> To unsubscribe, e-mail: dev-unsubscr...@commons.apache.org
>> For additional commands, e-mail: dev-h...@commons.apache.org
>>
>>
---------------------------------------------------------------------
To unsubscribe, e-mail: dev-unsubscr...@commons.apache.org
For additional commands, e-mail: dev-h...@commons.apache.org
