Vegas+

[scarrazza]

This PR implements Vegas+ from paper.

[scarlehoff]

Nice! Are there benchmarks for this?

[scarrazza]

Very good! Thanks for this great implementation.

Todo / summary:

• cross-check d_h definition.
• in order to include this module please edit the src/vegasflow/__init__.py.
• extend tests and examples.
• profile (performance) and compare to vegas importance sampling
• tests on multiple hardware (CPU, GPU).

[andrea-pasquale]

Nice! Are there benchmarks for this?

Benchmark using StratifiedFlow on my local machine (4 threads)

Gauss 4-d, 1e5 samples per iteration, rtol = 1e-3
Integrator VegasFlow converged in 4 iterations and took time 6.052872896194458 s
RESULT: 1.0003468042368207 +/- 0.0009635567255055014
Integrator PlainFlow cannot converge to the required percent uncertainty
Integrator StratifiedFlow cannot converge to the required percent uncertainty

Gauss 4-d, 1e5 samples per iteration, rtol = 1e-2
Integrator VegasFlow converged in 2 iterations and took time 5.815392017364502 s
RESULT: 1.004303283394375 +/- 0.00301380523785333
Integrator PlainFlow converged in 26 iterations and took time 0.5862269401550293 s
RESULT: 0.9950565948150915 +/- 0.009786743905761055
Integrator StratifiedFlow converged in 3 iterations and took time 1.4338939189910889 s
RESULT: 0.9842667020935826 +/- 0.009709536189911981

Gauss 4-d, 1e6 samples per iteration, rtol = 1e-3
Integrator VegasFlow converged in 2 iterations and took time 10.169922113418579 s
RESULT: 0.9999929734060133 +/- 0.0009445809568101617
Integrator PlainFlow cannot converge to the required percent uncertainty
Integrator StratifiedFlow converged in 54 iterations and took time 7.24327278137207 s
RESULT: 1.0006502778503834 +/- 0.0009941754908833862

Gauss 4-d, 1e6 samples per iteration, rtol = 1e-2
Integrator VegasFlow converged in 2 iterations and took time 10.056797981262207 s
RESULT: 1.0015330194146999 +/- 0.0009388122928216686
Integrator PlainFlow converged in 3 iterations and took time 0.5274820327758789 s
RESULT: 0.996918300290114 +/- 0.009100752007738737
Integrator StratifiedFlow converged in 1 iterations and took time 0.9030961990356445 s
RESULT: 0.9979067477189824 +/- 0.0020800953035176426

Gauss 8-d, 1e5 samples per iteration, rtol = 1e-2
Integrator VegasFlow converged in 4 iterations and took time 9.744021892547607 s
RESULT: 1.0010831277779253 +/- 0.004278006359363039
Integrator PlainFlow cannot converge to the required percent uncertainty
Integrator StratifiedFlow cannot converge to the required percent uncertainty

Gauss 8-d, 1e5 samples per iteration, rtol = 1e-3
Integrator VegasFlow converged in 9 iterations and took time 7.46278190612793 s
RESULT: 0.9999209309131106 +/- 0.0009953742839453855
Integrator PlainFlow cannot converge to the required percent uncertainty
Integrator StratifiedFlow cannot converge to the required percent uncertainty

Gauss 8-d, 1e6 samples per iteration, rtol = 1e-2
Integrator VegasFlow converged in 2 iterations and took time 13.53518295288086 s
RESULT: 1.003553386021623 +/- 0.0036343758510322724
Integrator PlainFlow cannot converge to the required percent uncertainty
Integrator StratifiedFlow converged in 35 iterations and took time 10.357421159744263 s
RESULT: 0.8813677900006324 +/- 0.00870484845558766

Gauss 8-d, 1e6 samples per iteration, rtol = 1e-3
Integrator VegasFlow converged in 3 iterations and took time 13.556339979171753 s
RESULT: 1.0007357025364503 +/- 0.0009060542874034692
Integrator PlainFlow cannot converge to the required percent uncertainty
Integrator StratifiedFlow cannot converge to the required percent uncertainty

[andrea-pasquale]

I tried to fix the retracing problem moving the generation of the random numbers between 0 and 1 from _run_event to generate_samples_in_hypercubes and adding an input signature to the function generate_samples_in_hypercubes but it didn't work.
I also added an example to show the retracing problem.
The function that triggers the retracing is <function MonteCarloFlow.compile.<locals>.run_event
What is the best way to solve this issue?

[scarrazza]

@scarlehoff could you please help here?

[scarlehoff]

My two cents would be to not compile at all _run_event. Using a variable number of events is tricky by nature and the content of _run_event is mostly summations so I don't think it would be that terrible.

Also, I think the retracing is "good" here because you are changing the self.n_ev so you need to either reassign it as it is done with the self.divisions or pass it through. I've chosen the pass it through because wasn't sure whether the number of elements would be changing.

I'm going to open a branch with suggested changes, you can either merge it if you like it or select only the changes you agree with.

btw I would take the generation of the random numbers outside of generate_samples_in_hypercubes (like in the commented code).

Edit: got confused about my own edit.

[andrea-pasquale]

Sorry I misclick

[scarlehoff]

I will be merging this during this week, there are quite a few conflicts so it might be non trivial though... If anyone has commits or changes they haven't pushed yet please let me know.

[andrea-pasquale]

I will be merging this during this week, there are quite a few conflicts so it might be non trivial though... If anyone has commits or changes they haven't pushed yet please let me know.

I have no commits to push. I won't be able to review the code until saturday since I am currently away from home.

[scarlehoff]

No worries, there's no urgency at all ^__^ thanks for confirming

[scarlehoff]

Ok, done. I've done some benchmarks and everything seems ok in terms of speed so I haven't broken anything I believe. I'll play a bit more of course, maybe even add the option to use Vegas+ to the script in https://github.com/N3PDF/madflow

TODOs:

Some documentation (just saying the algorithm is implemented and maybe linking to Lepage's paper) would be nice.

A benchmark against the implementation in Torchquad (https://github.com/esa/torchquad/) would be great to have as well.

However, these things might require time we might not have right now (from tomorrow and in the next few weeks I'll be busy) so I'll open separate issues for that (it's "out of scope" for this PR)

[scarlehoff]

Reading the changelogs for Tensorflow 2.6, it is possible that the problem we were seeing with tf.repeat are fixed now (#64 (comment))

Anyway, I still think the limit of 13 dimensions is a reasonable one.

[andrea-pasquale]

Reading the changelogs for Tensorflow 2.6, it is possible that the problem we were seeing with tf.repeat are fixed now (#64 (comment))
Anyway, I still think the limit of 13 dimensions is a reasonable one.

Nice.
I've just done some benchmarks and I can confirm that everything works fine.

[andrea-pasquale]

Right now, vegas+ cannot break the integration, correct?
Maybe it is worth considering for later developments of vegas+.

[scarlehoff]

Yes, and actually we should add a warning for now instead of doing it silently.
It should be indeed a priority.

[scarlehoff]

Reading the changelogs for Tensorflow 2.6, it is possible that the problem we were seeing with tf.repeat are fixed now (#64 (comment))
Anyway, I still think the limit of 13 dimensions is a reasonable one.

Nice.
I've just done some benchmarks and I can confirm that everything works fine.

Should I merge this then?

[andrea-pasquale]

Reading the changelogs for Tensorflow 2.6, it is possible that the problem we were seeing with tf.repeat are fixed now (#64 (comment))
Anyway, I still think the limit of 13 dimensions is a reasonable one.

Nice.
I've just done some benchmarks and I can confirm that everything works fine.

Should I merge this then?

Yes, I have no other changes to add.