DEP: integrate: deprecate romberg and quadrature

[mdhaber]

Reference issue

Toward gh-18574

What does this implement/fix?

scipy.integrate.quadrature and scipy.integrate.romberg do not reliably meet the accuracy target and are often slower the scipy.integrate.quad with the default vectorized=False; see #18574 (comment). This PR deprecates them in favor of scipy.integrate.quad.

Additional information

To do:

• Seek approval from mailing list
• Update tutorial
• Add to release notes
• Add to META: deprecation tracker #15765

[andyfaff]

I don't think there's a like to like equivalence for the replacement function. For example, I have a use case that requires adaptive quadrature. However, I want to start at an order of 20 and work upwards. I don't think I could do that with quad

[steppi]

Even without a use case, I think these still would still have value as reference implementations for classic algorithms.

[andyfaff]

@ev-br indicated +1 for deprecation.

Could we deprecate, but not remove?

[steppi]

Could we deprecate, but not remove?

I’d be ok with that. I think @mdhaber’s right that most users probably would be better served by using quad, and the deprecation warning could point them in the right direction. Someone who has a specific reason to use Romberg integration or Gaussian quadrature would still have the option at least and could just silence the warning.

[j-bowhay]

Could we deprecate, but not remove?

I would be -1 on deprecating but not removing. It creates confusion and the warnings are just annoying. Instead can we not just add a warning to the docs saying in most cases you would be better served by xyz.

Even without a use case, I think these still would still have value as reference implementations for classic algorithms.

I have to admit personally I do not think this is hugely valuable, most guass quadrature methods are trivial to implement given the collection of orthogonal polynomials in special. Plus there are other classics we are missing eg. other than the interface to quadpack we don't have a Clenshaw–Curtis implementation.

[mdhaber]

Could we deprecate, but not remove?

Another option is to declare them "legacy". It has been suggested that "legacy" code would not be maintained and would be released once as a separate package alongside SciPy 2.0, from which it would be removed. Documentation would steer new users away from legacy functions, but I don't think legacy code would emit a warning until a SciPy 2.0 is planned.

@h-vetinari and I had discussed reclaiming the name quadrature before that, though. The interface of quad has become a mess (gh-18574) and is difficult to fix in a backward compatible way, so we were looking for a new name.

I could live with declaring these as legacy if we can call a new function quadrate : )

That said, I would be surprised if there are many quadrature and romberg users out there. Otherwise, wouldn't we be flooded with bug reports about them and fewer about quad?

[andyfaff]

Is it possible to add new keywords to quad? e.g. Can we specify the minimum and maximum orders it'll try? For my personal projects I know that using a minimum of e.g. 21 is required, so there's no point in searching lower.

[steppi]

That said, I would be surprised if there are many quadrature and romberg users out there. Otherwise, wouldn't we be flooded with bug reports about them and fewer about quad?

quad must get more use, but Romberg integration and fixed tolerance Gaussian quadrature are pretty straightforward algorithms, and they're written in clean Python, so I don't think they have the kind of error surface QUADPACK has.

Also, I saw in your benchmarks above that romberg isn't returning error estimates. I'd consider that a bug. The absolute value of the difference between the values at the current and previous iterations gives a good error estimate.

I could live with declaring these as legacy if we can call a new function quadrate : )

Why not rename quadrature, gaussian_quadrature_fixed_tol or something, and deprecate and then reclaim quadrature?

[mdhaber]

so I don't think they have the kind of error surface QUADPACK has.

Right, they just give incorrect output and the user doesn't know, so they don't report an error : )

Also, I saw in your benchmarks above that romberg isn't returning error estimates. I'd consider that a bug.

The output is currently a NumPy float. We could add the error as an attribute (something I would actually like), but I get the impression that this would not be popular. So adding an error estimate would probably mean adding a keyword argument that changes the return type, and that's something we're trying to get away from.

Why not rename quadrature, gaussian_quadrature_fixed_tol...

That's always a possibility. I just prefer to trim it because I think there are better alternatives available.

Is it possible to add new keywords to quad?

Perhaps it would be possible when QUADPACK gets translated. I think @ilayn was interested in that eventually?

[steppi]

The output is currently a NumPy float. We could add the error as an attribute (something I would actually like), but I get the impression that this would not be popular. So adding an error estimate would probably mean adding a keyword argument that changes the return type, and that's something we're trying to get away from.

I see. Having romberg integration not return an error estimate is kind of odd, since one of its advantages is the quality of its error estimates. I don't think testing romberg over a large set of integrals of varying complexity is a good measure of whether it's worth keeping around. What it's suited for is quickly integrating over a finite interval a function which is smooth and well behaved. It does not do well with spikes and messy behavior, and I don't think it can handle singularities at all. There are particular cases where it can be a good choice though.

I think my preference would be to keep it around, but with the error estimates, and with documentation clearly stating the limitations of the method, and the situations where it may actually be useful.

[ilayn]

I think @ilayn was interested in that eventually?

Yes, intention is no unmaintained F77 code eventually. But I'd be really happy to not to translate if things are already not up to production level. Already feeling dizzy from cross referencing cryptic hieroglyphs. There is no point providing low-quality options to the users. If needed, anyone can download from Netlib and use it. Or alternatively we can leave it in a separate package as is.

Hence my vote is going for writing a bit of meson magic on top of things we deprecate and separate, so the legacy stuff is still pip installable and but leaving everything in the past like we did with scipy-weave.

This is all contingent upon having a better alternative of course.

[mdhaber]

Oops @ilayn Fortran quad is still the best general-purpose option, so it would be good to maintain that. The idea here is to deprecate two less general-purpose functions.

[j-bowhay]

It seems like the only objection is to deprecating quadrature not romberg. Could we not proceed with the deprecation of romberg and continue to discuss quadrature?

EDIT:
Searching reveals there is very little use of romberg: https://github.com/search?q=romberg+language%3APython+NOT+path%3A%2F%5Escipy%5C%2Fintegrate%5C%2F%2F+NOT+path%3A%2F%5Econtrib%5C%2Fpython%5C%2Fscipy%5C%2F%2F+NOT+path%3A%2F%5Evenv%5C%2FLib%5C%2Fsite-packages%5C%2Fscipy%5C%2Fintegrate%5C%2F%2F&type=code&l=Python

[ilayn]

The idea here is to deprecate two less general-purpose functions.

Then I don't know what I'm talking about 😃 so yes, at some point we will translate and for the rest I shut up.

[andyfaff]

With the help of Matt I've just been investigating the use of _tanhsinh for some of my personal projects. It's clear that that approach and indeed quad have much better than the existing quadrature. I'd like to change my -1 to +1 for deprecation, then reclaiming quadrature for a better performing algorithm.

[steppi]

I have no real objections to deprecating both since that seems to the consensus. I just want to make sure everyone’s aware that the benchmark set from #18574 (comment) contains examples which violate the assumptions of Romberg integration, e.g. singularities, points that evaluate to nan, rapidly oscillating or spiky functions. For smooth, well behaved functions, it converges rapidly with good error estimates. I think there are situations where a good Romberg implementation could be useful, but if someone was using it for these targeted cases, they would probably have asked for the error estimates to be returned.

[mdhaber]

Thanks @andyfaff!

It seems like the only objection is to deprecating quadrature not romberg.

@j-bowhay I think @steppi had reservations about removing the Romberg algorithm, although IIUC @steppi you would be OK deprecating if the algorithm were made available otherwise?

I thought about a possible compromise.

• deprecate romberg
• add a method argument to quadrature, and make 'Romberg', 'Gaussian', and 'tanh-sinh' available as methods. (This would allow the Romberg algorithm to return an error estimate.)
• warn that we are changing the following things about quadrature:
  • all arguments after b would become keyword only
  • name of tol would become atol
  • we might want to change the names of maxiter/miniter (or make them part of an options argument)
  • default method would become 'tanh-sinh'
• use make_tuple_bunch to turn the output tuple into a rich object in a backward-compatible way. New attributes would include iteration count, function evaluation count, and termination status.
• add the additional functionality of _tanhsinh to the other algorithms
  • vectorized support for elementwise callables (i.e. no need for for loops; just pass array limits of integration and args)
  • infinite limits of integration
  • iteration and function evaluation counts
  • termination status
  • "log-integration" (evaluate log of integral given callable that evaluates log of integrand)
  • callback interface
  • respect for dtype
  • (eventually) Array API compatibility so different backends can be used
• (eventually) add adaptive subdivision of the integral (basics would not be very hard)

This might work, but I see a few issues:

• It will be painful for us. The extra work tiptoeing around the existing interface may greatly exceed the reduction in impact to downstream users.
• If we don't change the output to a regular object now, we will probably end up supporting tuple unpacking of the integral and error estimate forever.
• (maybe most important) We will probably end up with bug reports from users expecting methods 'Gaussian' and 'Romberg' quadrature to work exactly as they have in the past. Try as we might, I don't think we will be able to move forward to improve these without little tweaks, accidental or desired.

Personally, I'd prefer a clean break, and once quadrature is available, we do all the good things with none of the pain points.

[steppi]

Personally, I'd prefer a clean break, and once quadrature is available, we do all the good things with none of the downsides.

I think I’m fine with the clean break then. There’d be opportunity cost with the extra work required. The effort would probably be better spent elsewhere.

[ev-br]

As a aside: is there a write-up somewhere about interface differences/ target use cases of tanhsinh vs quad_vec?

[mdhaber]

Not yet. Depending on how you write the callable, you can do integrals of vector-valued functions with _tanhsinh that you can do with quad_vec, but currently it is geared toward "elementwise" functions. The difference is that in an "elementwise" integrand (as I'm using the term), the expression used to evaluate the integrand (c*x**d) should be the same for different outputs - but the parameters (c, d) are allowed to be arrays. For a vector-valued integrand, the outputs could be evaluated by completely different functions, e.g. ([c*x**d, sin(x)]).

In any case, I plan to get rid of this limitations of _tanhsinh with an option. (Originally, it didn't have the limitations; I added it for performance in it's intended application.) _tanhsinh is also much, much faster, it is more reliable for general use (quad_vec is is not close to quad in this regard), and it has many other features mentioned in #19510 (comment).

[j-bowhay]

Think this is needed (as well as the fix in main) for the doc build

[mdhaber]

Thanks @j-bowhay. I merged main and added your suggestions:

[j-bowhay]

LGTM, happy to merge if that is the general consensus

[zachjweiner]

I just wanted to express (as a user) my strong opposition to removing algorithms, exactly per @steppi's original comment:

Even without a use case, I think these still would still have value as reference implementations for classic algorithms.

I'm sure at some point in the future I will use them, whether for sanity checks or teaching purposes.

I love the direction to unify the quadrature interfaces (and add better algorithms) - does that have to come at the cost of dropping other, already-implemented methods? Even if making Romberg and fixed-tol Gaussian options for the future interface is not feasible, why not just keep the old functions, sans deprecation (even if quadrature has to be renamed)? I would presume keeping them as-is incurs little maintenance burden. (And of course it makes sense to seclude them in the API reference and admonish their usage for general purposes.)

[mdhaber]

@j-bowhay @steppi @andyfaff I added gauss_quad as a replacement for quadrature, deprecated quadrature, and made gauss_quad and romberg legacy. Questions:

• Is the name gauss_quad OK? I considered gaussian to be like romberg, but went with gauss_quad along the lines of fixed_quad and qmc_quad.
• I'm not sure the admonition is supposed to be the first line. I think it looks a bit funny to replace the description in the TOC with an admonition box that is slightly below center. What do you think about swapping the admonition back to be the second line, since I've added a section to the TOC for "Legacy and Deprecated API"?

[j-bowhay]

@j-bowhay @steppi @andyfaff I added gauss_quad as a replacement for quadrature, deprecated quadrature, and made gauss_quad and romberg legacy. Questions:

I guess we should probably go back to the mailing list if we are going to add guass_quad to the public api

[mdhaber]

Sent

[mdhaber]

Sorry @zachjweiner, but I've had to revert back to deprecation due to feedback from the mailing list. These functions will always be available in the SciPy source for reference - git never erases anything from the history - and the plan is to replace them with better versions of themselves. You have them in the released SciPy for the next two versions, and they are short enough to copy-paste into your own code.

@j-bowhay I think this is back to the state of 5baf2f2 in which you approved it. I've seen support for this deprecation from six maintainers, which is more than we ever get, so I'd be happy with merging when tests pass.

[j-bowhay]

Sorry @zachjweiner, but I've had to revert back to deprecation due to feedback from the mailing list. These functions will always be available in the SciPy source for reference - git never erases anything from the history - and the plan is to replace them with better versions of themselves. You have them in the released SciPy for the next two versions, and they are short enough to copy-paste into your own code.

@j-bowhay I think this is back to the state of 5baf2f2 in which you approved it. I've seen support for this deprecation from six maintainers, which is more than we ever get, so I'd be happy with merging when tests pass.

Sounds good. @zachjweiner I think your concern would partially be addressed by the example for special.roots_legendre which shows how to write a simple code for guass quadrature

[zachjweiner]

and the plan is to replace them with better versions of themselves.

I got the opposite impression from the above discussion - that the plan was to drop the algorithms permanently rather than to reimplement them in a new interface - apologies if my perspective was misinformed! My opposition was only to reducing the number of implemented algorithms (I'm generally amenable to adapting personal code to API changes), and I thought keeping the old methods hanging around with whatever names would be the least burdensome way to avoid doing so. I appreciate the consideration either way.