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 <[email protected]> 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 <[email protected]> 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 <[email protected]> > > 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 <[email protected]> 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 <[email protected] > > > > > 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 < > [email protected] > > > > > > >>> 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 <[email protected]> > > 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 <[email protected]> > > > 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 < > > > >>>>>>> [email protected]> > > > >>>>>>> **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< > > > [email protected]> > > > >> For additional commands, e-mail: [email protected] > > > >> > > > >> > > > > > > --------------------------------------------------------------------- > > > To unsubscribe, e-mail: [email protected] > > > For additional commands, e-mail: [email protected] > > > > > > > > >
