Adds an alternative method to compute eta_function in CustomSD to deal with large memory-lengths

[emile-cochin]

Summary

This adds an alt_integrator parameter to PowerLawSD/CustomSD and to eta_function/correlation_2d_integral that dictates whether or not to use an alternative method to numerically integrate eta_function (splitting the integrals into oscillating terms and using cos/sin weighted integrals instead as described more precisely in #150). Also adds corresponding coverage tests. It adds a new physics test bath_correlations_test.py which compares exact values of these 2d integrals (for an Ohmic SD with exponential cutoff at 0 temperature) to both versions of the integrators, and for short and long memory-lengths. More specifically, this test expects either a sufficiently good agreement or at least an IntegrationWarning being raised if the results don't match. For now it passes in about 2s with the default integrator issuing a warning for long memory-lengths.

This is a draft since there are a few choices to discuss before this is potentially publishable. I'm notably thinking about:

• Whether giving the choice through this additional parameter is the way to go (or if its name is sufficiently clear).
• What choice of $\tilde{\omega}$ (omega_tilde) is best suited for each $\tau$ (tau). For now the function implemented is omega_tilde = omega_cutoff/int(2+omega_cutoff*tau/10) which seemed to give convincing results for a large number of spectral densities and parameter choices.
  It is heuristically motivated since the idea is to have a smaller first interval for larger $\tau$ to limit the number of oscillations that the non-specialized integrator has to deal with. Still, it looks a bit arbitrary and it might be nice to let the user be able to play around with this for their particular use case (in which case one could pass a function to alt_integrator or alternatively rewrite all of this to pass an integrator_options dictionary instead).
• If there should be a way to visualize the results of correlation_2d_integral in the helpers submodule and, if so, if it should let the user plot any correlation or specifically the ones as used in the influence_matrix function.
• If I should add an example of how to use this feature or if it is small enough not to be necessary.

Pull Request Check List

This checklist serves as a guide for contributors and maintainers to assess
whether a contribution can be merged into the main branch and thus serves
to guarantee the quality of the package. Please read the
contribution guideline.

Before you submit a pull request, please go through this checklist once
more, as it will decrease the number of unnecessary review cycles.
However, this list is there to help all contributors to produce high
quality code, not to deter you from contributing. If you can't tick some of the
boxes, feel free to submit the pull request anyway and report about it in the
pull request message. Also, some of the boxes might not apply to the specific
type of your contribution.

• The contribution has been discussed and agreed on in the Issue section.
• Code contributions do its best to follow the zen of python.
• The automated test are all positive:
  • tox -e py36 (to run pytest) the code tests.
  • tox -e style (to run pylint) the code style tests.
• Added test for changed/added code.
• The documentation has been updated:
  • docstring for all new functions/methods/classes/modules.
  • consistent style of all docstrings.
  • for new modules: /docs/pages/modules.rst has been updated.
  • for new functionality: /docs/pages/api.rst has been updated.
  • for new functionality: tutorials and examples have been updated.

[piperfw]

Some comments on docstrings relevant to ongoing Issue discussion

[piperfw]

Suggest small revision to mention scipy,

"Whether to use an alternative integration scheme leveraging sine/cosine weighted versions of scipy.integrate.quad to handle rapid oscillations. This may improve numerical accuracy and stability, especially for long memory times.
[Blank line]
The integral...
"

[piperfw]

Do we want to assert that, if num_oscillations is provided, alt_integrator should be true? Or at least provide a warning that it will be discarded when alt_integrator is False?

[piperfw]

Make consistent with above

[piperfw]

I think in docstrings other than CustomSD (or whatever you deem the 'parent' one to be) we can either a) copy the entire docstring or b) (probably neater) reference CustomSD, e.g.

Whether to use an alternative integration scheme leveraging sine/cosine weighted versions of scipy.integrate.quad to handle rapid oscillations. See...

[gefux]

Hi @emile-cochin, I am sorry for this massively delayed reply - I simply didn't get round to having a close enough look at your code until now. Thank you very much for this contribution, including the tests that compare the result with exact analytical results. This is very helpful.
Concerning the structure of the code, and the parameters I think this can be streamlined without all too much effort:

Currently the logic in the CustomSD.eta_function() is quite a bit tangled through different combinations of matsubara=True/False, alt_integrator=True/False, and T=0/T>0. Although the code gives the right results it is tough to understand for anyone joining the project and thus prone to produce errors/confusion in the future. Here are a few suggestions to make things simpler:

1. As I understand it, there is no known serious disadvantage in using your alternative method over the previous one. I therefore suggest we simply replace the previous method with yours.
2. Either the imaginary time computation can be brought into the exact same form as the real time computation (as it has been before -- i.e. without any if-blocks that are longer then a single line), or we should put them into different methods. Because the existence of an eta function is anyways special to the CustomSD (i.e. doesn't exist in BaseCorrelation), I don't see any problem with putting the two versions in separate specialized methods.
3. As the integration parameters now start to pile up, I agree that we should give the user the option to specify the integration parameters through a dictionary, that can be set in the constructor or the call of the CustomSD.eta_function(). This dictionary should have the fields epsilon, subdiv_limit, plus whatever you think are good parameters for your integration method. In the CustomSD.__init__() the integration_parameters should default to a dictionary set in config.py. [note that this needs to be done through using integration_params: Dict = None because dictionaries, like lists, are mutable (https://dnmtechs.com/best-practice-for-setting-default-value-of-list-parameter-in-python-3/)]. In the CustomSD.eta_function() it should then default to the values set in CustomSD.__init__(). Of course the description of all parameters would need to go into the single integration_parameters parameter docstring.

@emile-cochin @piperfw : Let me know what you think of this! --- Best, Gerald

[piperfw]

Great input on streamlining @gefux, I agree simplifying the logic we have here would be beneficial. Just on your points

1. There is one known issue that can affect us, that the weighted integrals do not handle small wvars incorrectly. Indeed Emile mentions problems when $\tau$ is small in Adding an alternative integrator to deal with fast oscillations in long memory-length numerical integrations for CustomSD #150. In fact, I think we need to check for and protect against this case whether we keep or replace the previous method. The interface simplification from having a single method is appealing.

A related question, @emile-cochin have you done any performance testing? Not that we prioritise fast code, but it would be interesting to know whether the gain from using the weighted version outweighs the overhead from replacing two calls to scipy.quad with (I think) six calls, for every eta integral.

2. Two separate versions sounds good unless Emile particularly wants to work on unifying them.

3. A dictionary would work well, subject to 1.: if we keep both methods, then the flexibility Emile included by allowing method parameters to override the value of alt_integrator passed to the customSD constructor may be useful e.g. to handle small tau?

[emile-cochin]

Hi @gefux, thank you very much for your inputs and very sorry for the delay, I'm back to working on OQuPy for the beginning of my PhD which started today. Here are my thoughts on these issues:

Just to answer @piperfw about performance: the cases where the memory length is low are basically instant in all cases and for long memory lengths, the new code performs a lot faster (sometimes up to x50 speedup) since it doesn't struggle with the oscillations like the previous one. If you want a proper benchmark I can also try to measure it properly, but since the difference was in favor of the new method, I didn't go beyond this.

1. For this point there is indeed the issue of very small $\tau$ for which the oscillatory integrals won't work. The new code deals with it by automatically switching to the old method when $\tau$ is too small (namely when $\tau < 1/\omega_c$ where there are no issues with oscillations anyway). For now, the old method was kept for this reason, but I agree that it can be pretty ugly and hard to understand for anyone joining in. One way to simplify it would be to completely replace the previous method by this one by using this 6-way splitting of the integral for everything and only deciding on which quad integrator (cos/sin/default) to use for the parts that can use it/when the conditions on $\tau$ are met.
2. This I don't know if I've understood well but if you're talking about the matsubara case, I don't believe the new code does it any differently than the previous one (apart from checking if alt_integrator and matsubara aren't both true at the same time). However, if we decide to remove the previous code, it might indeed be easier to have a separate method (I'll try and see if it is possible to fit it in no more than one if-block).
3. A dictionary seems like a good idea to me too, I'll implement this!