IterateiveLegendreGaussIntegrator should be replaced by adaptive guass-kronrod. The simpson and trapezoidal methods should be replaced by their adaptive versions, or at least adaptive wrapper provided for them.
On Tue, Jul 2, 2013 at 12:51 PM, Ajo Fod <ajo....@gmail.com> wrote: > 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 <kber...@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....@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 <kber...@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....@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 < > gil...@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 < > > 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 > > > > > > > > > > > > > >
