Mailing-List: contact dev-help@commons.apache.org; run by ezmlm
Precedence: bulk
Reply-To: "Commons Developers List" <dev@commons.apache.org>
Received-SPF: pass (nike.apache.org: domain of ajo.fod@gmail.com designates
 209.85.216.52 as permitted sender)
MIME-Version: 1.0
In-Reply-To: 
 <CAMFBX4hS5Oa0DMn0Z=wWiqRQqBzRv70zPCp61g3QH0f=cEKtvA@mail.gmail.com>
References: 
 <CAOOSspget8kcriyvBaGKqzedKtokeRSBWDXnj5ZO0T5HohSAYg@mail.gmail.com>
	<d94258cbd34131e5bee7706f16b27a52@scarlet.be>
	<CAJwFCa2uF7k_A9Z1O3_dDgZY+rT4bS5EmVhmG+EVHS1D6GQHUg@mail.gmail.com>
	<3ccc25c8612e4d60e47d038d4e1f2b8e@scarlet.be>
	<CAJwFCa0og_Zz8ayQSCCTL2nD2Nw4g--8bgskdEb39vSLE=2ngw@mail.gmail.com>
	<50ccb4d8edeb3d9f3075e346fb4f5eec@scarlet.be>
	<CAOOSspgahjqBLVcPziDeP6Te=h+dkqOWVHNxzmddU3gRY00gcA@mail.gmail.com>
	<CAOOSspjzwNe5YZTkrmQXgm_-9LZbd6s6Ryzts-U0RS=1g7cHgQ@mail.gmail.com>
	<CAOOSspieWEWSDFd6VRCggvnW4O1jjKZCAV0wq+gzwENRQ0N5yw@mail.gmail.com>
	<6a0d540f36cf14bc00ee71b52ca8a368@scarlet.be>
	<CAOOSspis2ZsbMYSMZTFQFVy24Oxf3fKwUM_BrDAsQv2F44ytTA@mail.gmail.com>
	<751027a8df1a5144ca7684c617c9db48@scarlet.be>
	<CAOOSspjwgOz99zHhfwJTqgqvY51XPe7BF0j7ha+XLcSa9ObcAw@mail.gmail.com>
	<-3037739403477643343@unknownmsgid>
	<CAOOSsphsSe_HjfG4cFZWLzv818sP5nvMhz4vOJgeju5k0a2B0Q@mail.gmail.com>
	<CAMFBX4hS5Oa0DMn0Z=wWiqRQqBzRv70zPCp61g3QH0f=cEKtvA@mail.gmail.com>
Date: Tue, 2 Jul 2013 09:51:28 -0700
Message-ID: 
 <CAOOSsphhqSE-XqubJ5w_254FDJrSs4BYznaFMBVN4x1+Bk8U_g@mail.gmail.com>
Subject: Re: [math] On MATH-995: Problems with LegendreGaussQuadrature class.
From: Ajo Fod <ajo.fod@gmail.com>
To: Commons Developers List <dev@commons.apache.org>
Content-Type: multipart/alternative; boundary=047d7bdca0ea0efd2704e08a286e

--047d7bdca0ea0efd2704e08a286e
Content-Type: text/plain; charset=ISO-8859-1

Konstantin,

In essence, we're both rooting for Adaptive Quadratures. IMHO
AdaptiveQuadrature (or some form of it that passes all Apache standards)
must replace IterateiveLegendreGaussIntegrator as soon as possible before
somebody else gets hit by a report of convergence to the wrong answer.
**

Cheers,
-Ajo


On Tue, Jul 2, 2013 at 9:33 AM, Konstantin Berlin <kberlin@gmail.com> wrote:

> At first it seems you are still compute redundant points. See my example
> that I posted, where I propagate 3 functional values not two.
>
> In regards to improvement. I am not an expert of different integration
> strategies but:
> The concept of adaptive quadrature is separate from how you integrate the
> subinterval. In your implementation (and mine) we use Simpson's rule. This
> is why putting comments on where your actual formulas come from is
> important for public code.
>         double Q1 = delta / 6 * (va + 4 * vm + vb); //this is simpson's
> rule of integration
> It does not make adaptive integration better or worse than Legendre-Gauss
> quadrature, rather they solve different problems.
> 2) The error tolerance should be an input to the quadrature, since
> tolerance is a problem dependent value. This values should be subdivided
> for each interval, such that the total sum would equal to the overall
> desired tolerance. Look at the code I provided. In any case, like Gilles
> said, the setup should be in the same spirit as other quadrature methods in
> commons.
> 3) Looking at Commons numerical integrations it is a little hard to read. I
> would have though that the basic quadratures, like trapezoidal or simpson's
> rule, are adaptive already. I think its rare to use these basic formulas
> without an adaptive quadrature. I would also think that Guass Kronrod would
> be a better way of going about the integration than what is there now for
> LGQ.
>  http://en.wikipedia.org/wiki/Gauss%E2%80%93Kronrod_quadrature_formula
>
>
> On Mon, Jul 1, 2013 at 11:37 PM, Ajo Fod <ajo.fod@gmail.com> wrote:
>
> > Hi Konstantin,
> >
> > Thanks for pointing out the inefficiency in AQ. I just improved the
> > efficiency of AQ to 1.41x that of LGQ (up from 1.05x) - measured in
> digits
> > of accuracy per evaluation for integral of normal with sigma 1000 in
> range
> > [-5000, 5000]
> >
> > Please let me know if this doesn't answer your question about the
> > discussion:
> >
> > In essence Gilles thinks there is no problem with LGQ because it
> integrates
> > the low frequency functions in his unit tests accurately. He thinks that
> > the problem is with the function I provided because it has "numerical
> > instabilities" while I think it is the high frequency nature of the
> > function that LGQ can't handle because it divides intervals
> > indiscriminately. So, it is not clear to me how Gilles explains why AQ
> > converges to the right answer in the presence of these "numerical
> > instabilities" ... after all, LGQ and AQ are being passed the same
> function
> > and identical limits.
> >
> > This problem should appear with any function that has sufficiently high
> > frequency components. It would be better if LGQ threw an exception when
> it
> > encounters a high frequency function. Instead, it "converges" confidently
> > to the wrong answer. I personally think that AQ will fail at some point
> at
> > a high enough frequency ... but that will be well beyond the point at
> which
> > LGQ fails.
> >
> > I've avoided the complex method of fetching weights used in the current
> > schemes because the  improvement in efficiency arises from the adaptive
> > nature of the AQ method. You may notice that I'm using the weights that
> you
> > use in your code. I think Gilles requires that I use the weight
> generation
> > scheme he has worked with in the codebase in order to consider the code
> > usable in Apache MATH.
> >
> > In summary, I feel the accuracy and versatility of AQ are being ignored
> in
> > favor of the familiarity of LGQ in the apache codebase. If there are
> tests
> > that AQ fails, I'll update my opinion.
> >
> > Cheers,
> > Ajo
> >
> >
> > On Mon, Jul 1, 2013 at 6:49 PM, Konstantin Berlin <kberlin@gmail.com>
> > wrote:
> >
> > > I am not understanding the discussion here. Adaptive integration is
> > > designed for functions that have different behavior in different
> > > regions. Some regions are smoother, some have higher frequeniesy. How
> > > you integrate a divided region, Simpson rule or whatever is a separate
> > > question.
> > >
> > > Adaptive integration saves of function evaluations by avoiding large
> > > series approximation in smooth regions. Nothing to do with how you
> > > compute the subdivided regions.
> > >
> > > On Jul 1, 2013, at 8:45 PM,  Fod <ajo.fod@gmail.com> wrote:
> > >
> > > > Hi Gilles,
> > > >
> > > > Your accuracy concern made me wonder. So, I dropped the
> > > > AdaptiveQuadrature.EPS to 1e-2 from 1e-9 in the code and ran the test
> > in
> > > > the patch.
> > > >
> > > > I computed the log of the error per evaluation ...i.e a measure of
> the
> > > > efficiency of the algorithm.
> > > > And wait for it ... AQ beats LGQ by about 5% for the particular
> > > formulation
> > > > of the problem.
> > > >
> > > > Your request to use the Math classes falls under "coding style" IMHO.
> > If
> > > it
> > > > doesn't satisfy your standards, feel free to modify. I'm happy with
> it.
> > > > Although as far as accuracy and convergence goes, I'd use AQ always.
> > > >
> > > >
> > > > Got to compare apples to apples Gilles !
> > > >
> > > > Cheers,
> > > > -Ajo
> > > >
> > > >
> > > >
> > > >
> > > > On Mon, Jul 1, 2013 at 4:16 PM, Gilles <gilles@harfang.homelinux.org
> >
> > > wrote:
> > > >
> > > >> Hi.
> > > >>
> > > >>
> > > >> On Mon, 1 Jul 2013 10:50:19 -0700, Ajo Fod 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 [...]
> > > >>
> > > >> You could file a "wish" request as a Commons Math's user. Then, if
> > > >> and when some regular contributor finds some time, he will try to
> > > >> implement the functionality.
> > > >>
> > > >> However, when you provide a patch for inclusion in the codebase, it
> > > >> is necessary to be more informed about similar functionality that
> > > >> would already exist in Commons Math, so that the contribution can be
> > > >> merged gracefully (i.e. with "minor" changes which committers will
> > > >> happily perform for you).
> > > >> You are welcome to ask questions in order to be able to contribute.
> > > >> As I tried to explain in more than one way, modifying your code is
> > > >> far from being trivial. If the committer has to figure out how to
> > > >> change/adapt/comment a significant part of the contribution, it
> > > >> ends up being easier to implement the feature from scratch!
> > > >>
> > > >>
> > > >>
> > > >>> 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
> > > >>> IterativeLegendreGaussIntegrat**or (refered to as LGQ) was not to
> > > >>> integrate
> > > >>> this kind of UnivariateFunction in [-1,1], ... what kind of
> > univariate
> > > >>> function would that be?
> > > >>
> > > >> Again, it is not just any UnivariateFunction, it is a function that
> > > >> maps the [-inf, +inf] interval into [-1, 1].
> > > >> It seems that the Gauss-Legendre quadrature is not appropriate for
> > > >> this. This is probably because the sample integration points do not
> > > >> cover the _whole_ interval: for the 10-point rule, the first point
> > > >> is at
> > > >> -0.9739065285171717
> > > >> and the last point is at
> > > >>  0.9994069665572084
> > > >> The interval in the original variable is thus [-18.908, 842.872].
> This
> > > >> is far from adequate for integrating a Gaussian function with
> > > sigma=1000.
> > > >> [And, as Phil pointed out from the outset, I suspect that the change
> > of
> > > >> variable also introduces numerical errors since the result becomes
> > worse
> > > >> when increasing the number of sample points. Increasing the
> requested
> > > >> precision leads to a prohibitive increase of the number of
> > evaluations,
> > > >> without improvement of the accuracy. In itself it is not sufficient
> to
> > > >> indicate a bug of "IterativeGaussLegendre"; it could simply be a
> > > >> limitation inherent to the algorithm.]
> > > >>
> > > >>
> > > >> If it is indeed supposed to do the integration,
> > > >>> then AQ clearly does a better job.
> > > >>
> > > >> Adaptive methods are certainly useful, but we need examples where
> > > >> its usage is appropriate. It is _not_ indicated for the improper
> > > >> integral of a Gaussian (even though it indeed performs better than
> > > >> Gauss-Legendre).
> > > >>
> > > >>
> > > >>
> > > >>> 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.
> > > >>
> > > >> Cf. above.
> > > >>
> > > >> AFAIU, the problem reported by MATH-995 is not a bug in
> > > >> "IterativeLegendreIntegrator": it correctly integrates a Gaussian
> > > >> with a large sigma _if_ the integration interval is "large enough"
> > > >> (cf. unit test referred to in my comment to MATH-995).
> > > >>
> > > >> Unless someone can point to something I'm missing in MATH-995, I'll
> > > >> close that issue.
> > > >>
> > > >>
> > > >>
> > > >>> Summary: I'm demonstrating a clear bug/inefficiency with LGQ
> > > >>
> > > >> I don't agree with that statement.
> > > >> "IterativeGaussLegendre" produces the correct answer (at 1e-6
> > > >> accuracy) in less than 60 function evaluations. To achieve the same,
> > > >> your code needs 995 evaluations.
> > > >>
> > > >>
> > > >> and providing
> > > >>> you with an alternative that is more accurate.
> > > >>
> > > >> Cf. above (and my previous post), about how to contribute to
> > > >> Commons Math.
> > > >> Please open a new feature request.
> > > >>
> > > >>
> > > >> Regards,
> > > >> Gilles
> > > >>
> > > >>
> > > >>> 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
> > >
> > >
> >
>

--047d7bdca0ea0efd2704e08a286e--