Implement adaptive integration algorithms

[haasal]

I propose following and extending scipys quad function which is based on QUADPACKS adaptive QAGS algorithm. I also propose setting the default integrate function from the prelude to this one.

The reason behind this is as follows:

1. The current integrate is slow/inaccurate for big intervals -> QAGS is much more versatile
2. All other integrators beside GK (slow and inaccurate) failed for me for peaked (delta-distribution-like) functions (output=0)
3. There is a huge array of problems one can solve with QAGS (S stands for singularities)
4. Currently there is no way to easily do infinite integration so I propose adding QAGI too for that case.
5. Incredibly fast and relatively intuitive to tune (relative error and max iterations)
6. I haven't really found any bigger rust library supporting QAGS yet (besides the gsl bindings).

In short: there is a reason the scipy project decided on QAGS as their quadrature algo of choice.

I would take the implementation of the gsl project.

This probably wouldn't be too hard to implement as the original algorithm is based on GK which we already have.

I would of course help if this gets accepted!

[Axect]

Thank you for this excellent and well-thought-out proposal! You raise several compelling points, especially regarding the limitations of current integration methods. You're absolutely right about the challenges with Dirac delta-like distributions - most other integration techniques indeed fail in these cases, and the current implementation does struggle with both accuracy and performance over large intervals.

The benefits of QAGS sound quite promising based on your description. While I haven't had the chance to look into it in detail yet, the ability to handle singularities effectively is a significant advantage. Additionally, your suggestion to implement QAGI for infinite interval integration would fill an important gap in our current functionality, as we currently don't have a robust solution for this case.

However, I should be transparent about the timeline: due to current workload, I likely won't be able to start implementing this until December. That said, if you're interested in implementing this yourself, we would definitely welcome a PR! Your understanding of both the mathematical background and the practical benefits of QAGS/QAGI would be valuable in this implementation.

Thanks again for bringing this up and offering to help with the implementation. It's exactly these kinds of well-reasoned proposals that help improve the library's functionality and performance.

[haasal]

Very nicely written response ;)

I would use the gsl implementations as a reference. The python c-implementation is horrible to read. But pythons QAGS implemtentation is different from gsl's implementation. For example in gsl QAGS I can not integrate from -inf to inf (only using QAGI but this resulted for me in oscillations probably due to numerical errors). In scipy you can enter -inf to inf as integration boundries into quad(...).
This whole thing hast cost me countless hours over the last three days and I still haven't found a solution to my particular problem.
I'm trying to fit superconducting gaps which are notoriously difficult so maybe this is the issue.

My mathematical background in numerical integration isn't particularly versed but I will probably try to read up on these algorithms.

Reference implementation.
Note: In scipy they have default handlers for non-convergence or round-off errors and in these cases the integrator simply returns nan (I think). I propose simply returning a Result instead.

[haasal]

Another idea: I think for many users it will be intuitive to insert ±std::f64::INFINITY into the integral bounds. So how about only for the simple interface if someone does that we log a warning (or something similar) and then instead of using the default qags we use qagi automatically. We definitely will have to add examples for this case.

Reasoning: The default/simple interface will be used mainly by people who just want the job done so this should work out-of-the-box for most cases. Many people don't bother looking at the integrant first and infinite bounds is a pretty common case (at least for physics).

[Axect]

That's a great suggestion about automatically switching to qagi when ±std::f64::INFINITY is detected in the bounds! I agree that implementing this behavior in the prelude's integrate function as the default would provide a more intuitive user experience. Many users would naturally try to use INFINITY as bounds, so having it "just work" would be quite beneficial.

However, I have some concerns about including warnings in this case. While I understand the intent to inform users about the automatic switch to qagi, this could be problematic in certain use cases. For instance, when processing large numbers of parallel integrations (which is quite common in scientific computing), these warnings could create significant noise in the output and potentially impact performance.

Perhaps instead of warnings, we could document this behavior clearly in the API docs and provide comprehensive examples showing both finite and infinite integration cases? This would maintain the convenient automatic switching while avoiding potential issues with warning spam in production use cases.

What are your thoughts on handling this without warnings? We could explore other ways to make the behavior clear to users if you think some form of notification is important.