commons-dev mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From Ajo Fod <ajo....@gmail.com>
Subject Re: [math] On MATH-995: Problems with LegendreGaussQuadrature class.
Date Mon, 01 Jul 2013 22:43:13 GMT
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
>
>

Mime
  • Unnamed multipart/alternative (inline, None, 0 bytes)
View raw message