commons-dev mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From Konstantin Berlin <kber...@gmail.com>
Subject Re: [math] On MATH-995: Problems with LegendreGaussQuadrature class.
Date Mon, 01 Jul 2013 23:57:04 GMT
  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.fod@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 <kberlin@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.fod@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 <gilles@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.fod@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.fod@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 <
>> gilles@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-unsubscribe@commons.apache.org>
>>>> For additional commands, e-mail: dev-help@commons.apache.org
>>>> 
>>>> 
>> 
>> 
>> ---------------------------------------------------------------------
>> To unsubscribe, e-mail: dev-unsubscribe@commons.apache.org
>> For additional commands, e-mail: dev-help@commons.apache.org
>> 
>> 


---------------------------------------------------------------------
To unsubscribe, e-mail: dev-unsubscribe@commons.apache.org
For additional commands, e-mail: dev-help@commons.apache.org


Mime
View raw message