From 4a9933a56f131003d0c7f782a1a3739ead31b70d Mon Sep 17 00:00:00 2001 From: Andrea Pasquale Date: Tue, 9 Feb 2021 10:58:42 +0100 Subject: [PATCH 01/27] stratified NON adaptive sampling --- src/vegasflow/stratified.py | 154 ++++++++++++++++++++++++++++++++++++ 1 file changed, 154 insertions(+) create mode 100644 src/vegasflow/stratified.py diff --git a/src/vegasflow/stratified.py b/src/vegasflow/stratified.py new file mode 100644 index 0000000..74224f2 --- /dev/null +++ b/src/vegasflow/stratified.py @@ -0,0 +1,154 @@ +""" + Plain implementation of the plainest possible MonteCarlo +""" + +from vegasflow.configflow import DTYPE, fone, fzero, DTYPEINT,DTYPE +from vegasflow.monte_carlo import MonteCarloFlow, wrapper +import tensorflow as tf + +import numpy as np +from itertools import chain,repeat,product +N_STRAT_MIN = 4 + + +@tf.function +def generate_random_array(rnds,n_strat,n_ev,hypercubes): + """Receives an array of random numbers 0 and 1 and + distribute them in each hypercube according to the + number of sample in each hypercube specified by n_ev + + Parameters + ---------- + `rnds`: tensor of random number betweenn 0 and 1 + `n_strat`: tensor with number of stratifications in each dimension + `n_ev`: tensor containing number of sample per each hypercube + `hypercubes`: tensor containing all different hypercube + + Returns + ------- + `x` : random numbers collocated in hypercubes + `w` : weight of each event + `segm` : segmentantion for later computations + """ + + # stratification width + delta_y = tf.cast(1/n_strat,DTYPE) + + random = rnds*delta_y + points = tf.cast(tf.repeat(hypercubes,n_ev,axis=0),DTYPE) + n_dim = int(tf.shape(hypercubes)[1]) + + x = points*delta_y + random - tf.cast([delta_y] * n_dim,DTYPE) + w = tf.cast(tf.repeat(1/n_ev,n_ev),DTYPE) + segm = tf.cast(tf.repeat(tf.range(fzero,tf.shape(hypercubes)[0]),n_ev),DTYPEINT) + return x,w,segm + + +class StratifiedFlow(MonteCarloFlow): + """ + Simple Monte Carlo integrator with stratified sampling. + """ + + def __init__(self, n_dim, n_events, **kwargs): + super().__init__(n_dim, n_events, **kwargs) + + # Initialize stratifications + self.n_strat = tf.math.floor(tf.math.pow(n_events/N_STRAT_MIN, 1/n_dim)) + substratification_np = np.linspace(0,1,int(self.n_strat)+1) + stratifications_np = substratification_np.repeat(n_dim).reshape(-1, n_dim).T + self.stratifications = tf.Variable(stratifications_np, dtype=DTYPE) + + # Initialize hypercubes + hypercubes_one_dim = np.arange(1,int(self.n_strat)+1) + hypercubes = [list(p) for p in product(hypercubes_one_dim, repeat=int(n_dim))] + hypercubes_tf=tf.convert_to_tensor(hypercubes) + self.hypercubes = hypercubes_tf + + #Initialize n_ev + n_ev = tf.math.floor(n_events/tf.math.pow(self.n_strat,n_dim)) + self.n_ev=tf.fill([1 ,len(hypercubes)], int(n_ev)) + self.n_ev = tf.reshape(self.n_ev,[-1]) + + #correction of self.n_events due to samples per hypercube + self.n_events = int(n_ev*tf.math.pow(self.n_strat,n_dim)) + + + + + def _run_event(self, integrand, ncalls=None): + + if ncalls is None: + n_events = self.n_events + else: + n_events = ncalls + + # Generate all random numbers + rnds = tf.random.uniform((n_events, self.n_dim), minval=0, maxval=1, dtype=DTYPE) + + # Pass random numbers in hypercubes + x, w, segm = generate_random_array(rnds,self.n_strat,self.n_ev,self.hypercubes) + + # compute integrand + xjac = w + tmp = integrand(x, n_dim=self.n_dim, weight=xjac) * xjac + tmp2 = tf.square(tmp) + + # tensor containing resummed component for each hypercubes + ress = tf.math.segment_sum(tmp,segm) + ress2 = tf.math.segment_sum(tmp2,segm) + + return ress, ress2 + + + def _run_iteration(self): + + ress, raw_ress2 = self.run_event() + + # compute variance for each hypercube + ress2 = raw_ress2 * tf.cast(self.n_ev,DTYPE) + err_tmp2s = (ress2 - tf.square(ress))/(tf.cast(self.n_ev,DTYPE)-fone) + sigmas2 = tf.maximum(err_tmp2s,fzero) + + + res = tf.reduce_sum(ress) / tf.cast(tf.math.pow(self.n_strat,self.n_dim),DTYPE) + sigma2 = tf.reduce_sum(sigmas2)/tf.cast(tf.math.pow(self.n_strat,self.n_dim),DTYPE)/self.n_events + sigma = tf.sqrt(sigma2) + return res, sigma + + +def plain_wrapper(*args): + """ Wrapper around PlainFlow """ + return wrapper(StratifiedFlow, *args) + + + + +@tf.function +def symgauss(xarr, n_dim=None, **kwargs): + """symgauss test function""" + if n_dim is None: + n_dim = xarr.shape[-1] + a = tf.constant(0.1, dtype=DTYPE) + n100 = tf.cast(100 * n_dim, dtype=DTYPE) + pref = tf.pow(1.0 / a / np.sqrt(np.pi), n_dim) + coef = tf.reduce_sum(tf.range(n100 + 1)) + coef += tf.reduce_sum(tf.square((xarr - 1.0 / 2.0) / a), axis=1) + coef -= (n100 + 1) * n100 / 2.0 + return pref * tf.exp(-coef) + + +if __name__ == "__main__": + + import time + # MC integration setup + dim = 4 + ncalls = int(1e3) + n_iter = 5 + + start = time.time() + vegas_instance = StratifiedFlow(dim, ncalls,simplify_signature=True) + #print(int(tf.shape(vegas_instance.hypercubes)[1])) + vegas_instance.compile(symgauss) + result = vegas_instance.run_integration(n_iter) + end = time.time() + print(f"Vegas took: time (s): {end-start}") \ No newline at end of file From a7517fc69d71e1303e0cf33e647ece07bacc8248 Mon Sep 17 00:00:00 2001 From: Andrea Pasquale Date: Tue, 9 Feb 2021 12:02:11 +0100 Subject: [PATCH 02/27] added adaptive-stratified sampling --- src/vegasflow/stratified.py | 32 +++++++++++++++++++++++++------- 1 file changed, 25 insertions(+), 7 deletions(-) diff --git a/src/vegasflow/stratified.py b/src/vegasflow/stratified.py index 74224f2..adb51e5 100644 --- a/src/vegasflow/stratified.py +++ b/src/vegasflow/stratified.py @@ -9,6 +9,7 @@ import numpy as np from itertools import chain,repeat,product N_STRAT_MIN = 4 +BETA = 0.75 @tf.function @@ -46,18 +47,19 @@ def generate_random_array(rnds,n_strat,n_ev,hypercubes): class StratifiedFlow(MonteCarloFlow): """ - Simple Monte Carlo integrator with stratified sampling. + Simple Monte Carlo integrator with Stratified Sampling. """ - def __init__(self, n_dim, n_events, **kwargs): + def __init__(self, n_dim, n_events,adaptive=True, **kwargs): super().__init__(n_dim, n_events, **kwargs) + self.adaptive = adaptive + # Initialize stratifications self.n_strat = tf.math.floor(tf.math.pow(n_events/N_STRAT_MIN, 1/n_dim)) substratification_np = np.linspace(0,1,int(self.n_strat)+1) stratifications_np = substratification_np.repeat(n_dim).reshape(-1, n_dim).T self.stratifications = tf.Variable(stratifications_np, dtype=DTYPE) - # Initialize hypercubes hypercubes_one_dim = np.arange(1,int(self.n_strat)+1) hypercubes = [list(p) for p in product(hypercubes_one_dim, repeat=int(n_dim))] @@ -72,7 +74,13 @@ def __init__(self, n_dim, n_events, **kwargs): #correction of self.n_events due to samples per hypercube self.n_events = int(n_ev*tf.math.pow(self.n_strat,n_dim)) - + def redistribute_samples(self,arr_var): + """Receives an array with the variance of the integrand in each hypercube + and recalculate the samples per hypercube according to VEGAS+ algorithm""" + + damped_arr_var = tf.pow(arr_var,BETA) + new_n_ev = tf.maximum(2,damped_arr_var * self.n_events/tf.reduce_sum(damped_arr_var)) + self.n_ev = tf.math.floor(new_n_ev) def _run_event(self, integrand, ncalls=None): @@ -97,12 +105,18 @@ def _run_event(self, integrand, ncalls=None): ress = tf.math.segment_sum(tmp,segm) ress2 = tf.math.segment_sum(tmp2,segm) - return ress, ress2 + #if adaptive save variance of each hypercube + arr_var = None + if self.adaptive: + hypercube_volume = tf.cast(tf.math.pow(1/self.n_strat, self.n_dim),DTYPE) + arr_var = ress2 * tf.cast(self.n_ev,DTYPE)* tf.square(hypercube_volume) - tf.square(ress*hypercube_volume) + + return ress, ress2, arr_var def _run_iteration(self): - ress, raw_ress2 = self.run_event() + ress, raw_ress2, arr_var = self.run_event() # compute variance for each hypercube ress2 = raw_ress2 * tf.cast(self.n_ev,DTYPE) @@ -113,6 +127,11 @@ def _run_iteration(self): res = tf.reduce_sum(ress) / tf.cast(tf.math.pow(self.n_strat,self.n_dim),DTYPE) sigma2 = tf.reduce_sum(sigmas2)/tf.cast(tf.math.pow(self.n_strat,self.n_dim),DTYPE)/self.n_events sigma = tf.sqrt(sigma2) + + # If adaptive is active redistributes samples + if self.adaptive: + self.redistribute_samples(arr_var) + return res, sigma @@ -147,7 +166,6 @@ def symgauss(xarr, n_dim=None, **kwargs): start = time.time() vegas_instance = StratifiedFlow(dim, ncalls,simplify_signature=True) - #print(int(tf.shape(vegas_instance.hypercubes)[1])) vegas_instance.compile(symgauss) result = vegas_instance.run_integration(n_iter) end = time.time() From c9ec2fc922e2f6b11ed1ad8941ee6e1ee1064836 Mon Sep 17 00:00:00 2001 From: Andrea Pasquale Date: Sat, 13 Feb 2021 09:42:30 +0100 Subject: [PATCH 03/27] PEP8 style corrections --- src/vegasflow/stratified.py | 127 ++++++++++++++++++++---------------- 1 file changed, 69 insertions(+), 58 deletions(-) diff --git a/src/vegasflow/stratified.py b/src/vegasflow/stratified.py index adb51e5..6c85aca 100644 --- a/src/vegasflow/stratified.py +++ b/src/vegasflow/stratified.py @@ -1,21 +1,19 @@ """ - Plain implementation of the plainest possible MonteCarlo + Implementation of the stratified sampling """ - -from vegasflow.configflow import DTYPE, fone, fzero, DTYPEINT,DTYPE -from vegasflow.monte_carlo import MonteCarloFlow, wrapper +from itertools import chain, repeat, product import tensorflow as tf - import numpy as np -from itertools import chain,repeat,product +from vegasflow.configflow import DTYPE, fone, fzero, DTYPEINT +from vegasflow.monte_carlo import MonteCarloFlow, wrapper + N_STRAT_MIN = 4 BETA = 0.75 - @tf.function -def generate_random_array(rnds,n_strat,n_ev,hypercubes): - """Receives an array of random numbers 0 and 1 and - distribute them in each hypercube according to the +def generate_random_array(rnds, n_strat, n_ev, hypercubes): + """Receives an array of random numbers 0 and 1 and + distribute them in each hypercube according to the number of sample in each hypercube specified by n_ev Parameters @@ -29,20 +27,21 @@ def generate_random_array(rnds,n_strat,n_ev,hypercubes): ------- `x` : random numbers collocated in hypercubes `w` : weight of each event - `segm` : segmentantion for later computations + `segm` : segmentantion for later computations """ - # stratification width - delta_y = tf.cast(1/n_strat,DTYPE) + delta_y = tf.cast(1/n_strat, DTYPE) random = rnds*delta_y - points = tf.cast(tf.repeat(hypercubes,n_ev,axis=0),DTYPE) + points = tf.cast(tf.repeat(hypercubes, n_ev, axis=0), DTYPE) n_dim = int(tf.shape(hypercubes)[1]) - - x = points*delta_y + random - tf.cast([delta_y] * n_dim,DTYPE) - w = tf.cast(tf.repeat(1/n_ev,n_ev),DTYPE) - segm = tf.cast(tf.repeat(tf.range(fzero,tf.shape(hypercubes)[0]),n_ev),DTYPEINT) - return x,w,segm + x = points*delta_y + random - tf.cast([delta_y] * n_dim, DTYPE) + w = tf.cast(tf.repeat(1/n_ev, n_ev), DTYPE) + segm = tf.cast(tf.repeat(tf.range(fzero, + tf.shape(hypercubes)[0]), + n_ev), + DTYPEINT) + return x, w, segm class StratifiedFlow(MonteCarloFlow): @@ -50,39 +49,43 @@ class StratifiedFlow(MonteCarloFlow): Simple Monte Carlo integrator with Stratified Sampling. """ - def __init__(self, n_dim, n_events,adaptive=True, **kwargs): + def __init__(self, n_dim, n_events, adaptive=True, **kwargs): super().__init__(n_dim, n_events, **kwargs) self.adaptive = adaptive # Initialize stratifications - self.n_strat = tf.math.floor(tf.math.pow(n_events/N_STRAT_MIN, 1/n_dim)) - substratification_np = np.linspace(0,1,int(self.n_strat)+1) - stratifications_np = substratification_np.repeat(n_dim).reshape(-1, n_dim).T + self.n_strat = tf.math.floor(tf.math.pow(n_events/N_STRAT_MIN, + 1/n_dim)) + substratification_np = np.linspace(0, 1, int(self.n_strat)+1) + stratifications_np = substratification_np.repeat(n_dim).reshape(-1, + n_dim).T self.stratifications = tf.Variable(stratifications_np, dtype=DTYPE) # Initialize hypercubes - hypercubes_one_dim = np.arange(1,int(self.n_strat)+1) - hypercubes = [list(p) for p in product(hypercubes_one_dim, repeat=int(n_dim))] - hypercubes_tf=tf.convert_to_tensor(hypercubes) + hypercubes_one_dim = np.arange(1, int(self.n_strat)+1) + hypercubes = [list(p) for p in product(hypercubes_one_dim, + repeat=int(n_dim))] + hypercubes_tf = tf.convert_to_tensor(hypercubes) self.hypercubes = hypercubes_tf - #Initialize n_ev - n_ev = tf.math.floor(n_events/tf.math.pow(self.n_strat,n_dim)) - self.n_ev=tf.fill([1 ,len(hypercubes)], int(n_ev)) - self.n_ev = tf.reshape(self.n_ev,[-1]) + # Initialize n_ev + n_ev = tf.math.floor(n_events/tf.math.pow(self.n_strat, n_dim)) + self.n_ev = tf.fill([1, len(hypercubes)], int(n_ev)) + self.n_ev = tf.reshape(self.n_ev, [-1]) - #correction of self.n_events due to samples per hypercube - self.n_events = int(n_ev*tf.math.pow(self.n_strat,n_dim)) + # Correction of self.n_events due to samples per hypercube + self.n_events = int(n_ev*tf.math.pow(self.n_strat, n_dim)) - def redistribute_samples(self,arr_var): - """Receives an array with the variance of the integrand in each hypercube - and recalculate the samples per hypercube according to VEGAS+ algorithm""" + def redistribute_samples(self, arr_var): + """Receives an array with the variance of the integrand in each + hypercube and recalculate the samples per hypercube according + to VEGAS+ algorithm""" - damped_arr_var = tf.pow(arr_var,BETA) - new_n_ev = tf.maximum(2,damped_arr_var * self.n_events/tf.reduce_sum(damped_arr_var)) + damped_arr_var = tf.pow(arr_var, BETA) + new_n_ev = tf.maximum(2, + damped_arr_var * self.n_events/tf.reduce_sum(damped_arr_var)) self.n_ev = tf.math.floor(new_n_ev) - def _run_event(self, integrand, ncalls=None): if ncalls is None: @@ -90,11 +93,13 @@ def _run_event(self, integrand, ncalls=None): else: n_events = ncalls - # Generate all random numbers - rnds = tf.random.uniform((n_events, self.n_dim), minval=0, maxval=1, dtype=DTYPE) + # Generate all random numbers + rnds = tf.random.uniform((n_events, self.n_dim), minval=0, maxval=1, + dtype=DTYPE) # Pass random numbers in hypercubes - x, w, segm = generate_random_array(rnds,self.n_strat,self.n_ev,self.hypercubes) + x, w, segm = generate_random_array(rnds, self.n_strat, self.n_ev, + self.hypercubes) # compute integrand xjac = w @@ -102,30 +107,37 @@ def _run_event(self, integrand, ncalls=None): tmp2 = tf.square(tmp) # tensor containing resummed component for each hypercubes - ress = tf.math.segment_sum(tmp,segm) - ress2 = tf.math.segment_sum(tmp2,segm) + ress = tf.math.segment_sum(tmp, segm) + ress2 = tf.math.segment_sum(tmp2, segm) - #if adaptive save variance of each hypercube + # if adaptive save variance of each hypercube arr_var = None if self.adaptive: - hypercube_volume = tf.cast(tf.math.pow(1/self.n_strat, self.n_dim),DTYPE) - arr_var = ress2 * tf.cast(self.n_ev,DTYPE)* tf.square(hypercube_volume) - tf.square(ress*hypercube_volume) + hypercube_volume = tf.cast(tf.math.pow(1/self.n_strat, self.n_dim), + DTYPE) + arr_var = ((ress2 * + tf.cast(self.n_ev, DTYPE) * + tf.square(hypercube_volume)) - + tf.square(ress*hypercube_volume)) return ress, ress2, arr_var - def _run_iteration(self): ress, raw_ress2, arr_var = self.run_event() # compute variance for each hypercube - ress2 = raw_ress2 * tf.cast(self.n_ev,DTYPE) - err_tmp2s = (ress2 - tf.square(ress))/(tf.cast(self.n_ev,DTYPE)-fone) - sigmas2 = tf.maximum(err_tmp2s,fzero) - - - res = tf.reduce_sum(ress) / tf.cast(tf.math.pow(self.n_strat,self.n_dim),DTYPE) - sigma2 = tf.reduce_sum(sigmas2)/tf.cast(tf.math.pow(self.n_strat,self.n_dim),DTYPE)/self.n_events + ress2 = raw_ress2 * tf.cast(self.n_ev, DTYPE) + err_tmp2s = (ress2 - tf.square(ress))/(tf.cast(self.n_ev, DTYPE)-fone) + sigmas2 = tf.maximum(err_tmp2s, fzero) + + res = tf.reduce_sum(ress) / tf.cast(tf.math.pow(self.n_strat, + self.n_dim), + DTYPE) + + sigma2 = tf.reduce_sum(sigmas2)/tf.cast(tf.math.pow(self.n_strat, + self.n_dim), + DTYPE) / self.n_events sigma = tf.sqrt(sigma2) # If adaptive is active redistributes samples @@ -140,8 +152,6 @@ def plain_wrapper(*args): return wrapper(StratifiedFlow, *args) - - @tf.function def symgauss(xarr, n_dim=None, **kwargs): """symgauss test function""" @@ -165,8 +175,9 @@ def symgauss(xarr, n_dim=None, **kwargs): n_iter = 5 start = time.time() - vegas_instance = StratifiedFlow(dim, ncalls,simplify_signature=True) + vegas_instance = StratifiedFlow(dim, ncalls, simplify_signature=True) vegas_instance.compile(symgauss) result = vegas_instance.run_integration(n_iter) end = time.time() - print(f"Vegas took: time (s): {end-start}") \ No newline at end of file + print(f"Vegas took: time (s): {end-start}") + \ No newline at end of file From 883cd007d472be62e8c43fec321ccc37f703ec46 Mon Sep 17 00:00:00 2001 From: Andrea Pasquale Date: Sat, 13 Feb 2021 10:00:21 +0100 Subject: [PATCH 04/27] updated __init__.py and test_algs.py --- src/vegasflow/__init__.py | 1 + src/vegasflow/stratified.py | 3 +-- src/vegasflow/tests/test_algs.py | 6 ++++++ 3 files changed, 8 insertions(+), 2 deletions(-) diff --git a/src/vegasflow/__init__.py b/src/vegasflow/__init__.py index f4e0c11..fc8c171 100644 --- a/src/vegasflow/__init__.py +++ b/src/vegasflow/__init__.py @@ -4,5 +4,6 @@ # Expose the main interfaces from vegasflow.vflow import VegasFlow, vegas_wrapper from vegasflow.plain import PlainFlow, plain_wrapper +from vegasflow.stratified import StratifiedFlow, stratified_wrapper __version__ = "1.2.0" diff --git a/src/vegasflow/stratified.py b/src/vegasflow/stratified.py index 6c85aca..4042770 100644 --- a/src/vegasflow/stratified.py +++ b/src/vegasflow/stratified.py @@ -147,7 +147,7 @@ def _run_iteration(self): return res, sigma -def plain_wrapper(*args): +def stratified_wrapper(*args): """ Wrapper around PlainFlow """ return wrapper(StratifiedFlow, *args) @@ -180,4 +180,3 @@ def symgauss(xarr, n_dim=None, **kwargs): result = vegas_instance.run_integration(n_iter) end = time.time() print(f"Vegas took: time (s): {end-start}") - \ No newline at end of file diff --git a/src/vegasflow/tests/test_algs.py b/src/vegasflow/tests/test_algs.py index f2d0ad9..a2032d7 100644 --- a/src/vegasflow/tests/test_algs.py +++ b/src/vegasflow/tests/test_algs.py @@ -13,6 +13,7 @@ from vegasflow.configflow import DTYPE from vegasflow.vflow import VegasFlow from vegasflow.plain import PlainFlow +from vegasflow.stratified import StratifiedFlow # Test setup dim = 2 @@ -109,3 +110,8 @@ def test_PlainFlow(): plain_instance = instance_and_compile(PlainFlow) result = plain_instance.run_integration(n_iter) check_is_one(result) + +def test_StratifiedFlow(): + stratified_instance = instance_and_compile(StratifiedFlow) + result = stratified_instance.run_integration(n_iter) + check_is_one(result) From befb0322d39e5215dccb9bc38bda3f8e80c159b1 Mon Sep 17 00:00:00 2001 From: Andrea Pasquale Date: Sat, 13 Feb 2021 11:15:18 +0100 Subject: [PATCH 05/27] added example using StratifiedFlow --- examples/example_stratified.py | 42 ++++++++++++++++++++++++++++++++++ src/vegasflow/stratified.py | 31 +------------------------ 2 files changed, 43 insertions(+), 30 deletions(-) create mode 100644 examples/example_stratified.py diff --git a/examples/example_stratified.py b/examples/example_stratified.py new file mode 100644 index 0000000..894e05c --- /dev/null +++ b/examples/example_stratified.py @@ -0,0 +1,42 @@ +""" + Example: integration with stratified sampling + + Basic example using the stratified_wrapper helper +""" + +from vegasflow.configflow import DTYPE +import time +import numpy as np +import tensorflow as tf +from vegasflow.vflow import vegas_wrapper +from vegasflow.plain import plain_wrapper +from vegasflow.stratified import stratified_wrapper + +# MC integration setup +dim = 4 +ncalls = np.int32(1e5) +n_iter = 5 + + +@tf.function +def symgauss(xarr, n_dim=None, **kwargs): + """symgauss test function""" + if n_dim is None: + n_dim = xarr.shape[-1] + a = tf.constant(0.1, dtype=DTYPE) + n100 = tf.cast(100 * n_dim, dtype=DTYPE) + pref = tf.pow(1.0 / a / np.sqrt(np.pi), n_dim) + coef = tf.reduce_sum(tf.range(n100 + 1)) + coef += tf.reduce_sum(tf.square((xarr - 1.0 / 2.0) / a), axis=1) + coef -= (n100 + 1) * n100 / 2.0 + return pref * tf.exp(-coef) + + +if __name__ == "__main__": + """Testing several different integrations""" + print(f"VEGAS MC, ncalls={ncalls}:") + start = time.time() + ncalls = 10*ncalls + r = stratified_wrapper(symgauss, dim, n_iter, ncalls) + end = time.time() + print(f"Vegas took: time (s): {end-start}") \ No newline at end of file diff --git a/src/vegasflow/stratified.py b/src/vegasflow/stratified.py index 4042770..b640bc4 100644 --- a/src/vegasflow/stratified.py +++ b/src/vegasflow/stratified.py @@ -10,6 +10,7 @@ N_STRAT_MIN = 4 BETA = 0.75 + @tf.function def generate_random_array(rnds, n_strat, n_ev, hypercubes): """Receives an array of random numbers 0 and 1 and @@ -150,33 +151,3 @@ def _run_iteration(self): def stratified_wrapper(*args): """ Wrapper around PlainFlow """ return wrapper(StratifiedFlow, *args) - - -@tf.function -def symgauss(xarr, n_dim=None, **kwargs): - """symgauss test function""" - if n_dim is None: - n_dim = xarr.shape[-1] - a = tf.constant(0.1, dtype=DTYPE) - n100 = tf.cast(100 * n_dim, dtype=DTYPE) - pref = tf.pow(1.0 / a / np.sqrt(np.pi), n_dim) - coef = tf.reduce_sum(tf.range(n100 + 1)) - coef += tf.reduce_sum(tf.square((xarr - 1.0 / 2.0) / a), axis=1) - coef -= (n100 + 1) * n100 / 2.0 - return pref * tf.exp(-coef) - - -if __name__ == "__main__": - - import time - # MC integration setup - dim = 4 - ncalls = int(1e3) - n_iter = 5 - - start = time.time() - vegas_instance = StratifiedFlow(dim, ncalls, simplify_signature=True) - vegas_instance.compile(symgauss) - result = vegas_instance.run_integration(n_iter) - end = time.time() - print(f"Vegas took: time (s): {end-start}") From f0fd6098861363c5a71b4d5567fb2fdf80add711 Mon Sep 17 00:00:00 2001 From: Andrea Pasquale Date: Mon, 15 Feb 2021 11:53:30 +0100 Subject: [PATCH 06/27] fixed error in redistribute_samples --- src/vegasflow/stratified.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/vegasflow/stratified.py b/src/vegasflow/stratified.py index b640bc4..65f5837 100644 --- a/src/vegasflow/stratified.py +++ b/src/vegasflow/stratified.py @@ -82,9 +82,9 @@ def redistribute_samples(self, arr_var): hypercube and recalculate the samples per hypercube according to VEGAS+ algorithm""" - damped_arr_var = tf.pow(arr_var, BETA) + damped_arr_sdev = tf.pow(arr_var, BETA/2) new_n_ev = tf.maximum(2, - damped_arr_var * self.n_events/tf.reduce_sum(damped_arr_var)) + damped_arr_sdev * self.n_events/tf.reduce_sum(damped_arr_sdev)) self.n_ev = tf.math.floor(new_n_ev) def _run_event(self, integrand, ncalls=None): From bf138495222c0e37f08121fd50f158a426279681 Mon Sep 17 00:00:00 2001 From: Andrea Pasquale Date: Tue, 2 Mar 2021 11:30:21 +0100 Subject: [PATCH 07/27] first implementation of vegas+ in vegasflow --- src/vegasflow/__init__.py | 1 + src/vegasflow/vflowplus.py | 209 +++++++++++++++++++++++++++++++++++++ 2 files changed, 210 insertions(+) create mode 100644 src/vegasflow/vflowplus.py diff --git a/src/vegasflow/__init__.py b/src/vegasflow/__init__.py index fc8c171..7af16fb 100644 --- a/src/vegasflow/__init__.py +++ b/src/vegasflow/__init__.py @@ -5,5 +5,6 @@ from vegasflow.vflow import VegasFlow, vegas_wrapper from vegasflow.plain import PlainFlow, plain_wrapper from vegasflow.stratified import StratifiedFlow, stratified_wrapper +from vegasflow.vflowplus import VegasFlowPlus, vegasflowplus_wrapper __version__ = "1.2.0" diff --git a/src/vegasflow/vflowplus.py b/src/vegasflow/vflowplus.py new file mode 100644 index 0000000..6a540a2 --- /dev/null +++ b/src/vegasflow/vflowplus.py @@ -0,0 +1,209 @@ +#!/usr/env python +""" + Implementation of vegas+ algorithm: + adaptive importance sampling + adaptive stratified sampling + +""" +from itertools import chain, repeat, product +import tensorflow as tf +import numpy as np +from vegasflow.configflow import DTYPE, DTYPEINT, fone, fzero, float_me, ione, int_me, BINS_MAX +from vegasflow.monte_carlo import wrapper +from vegasflow.vflow import VegasFlow +from vegasflow.utils import consume_array_into_indices + + +N_STRAT_MIN = 4 +BETA = 0.75 +FBINS = float_me(BINS_MAX) + + +@tf.function +def generate_samples_in_hypercubes(randoms, n_strat, n_ev, hypercubes, divisions): + """Receives an array of random numbers 0 and 1 and + distribute them in each hypercube according to the + number of sample in each hypercube specified by n_ev + + Parameters + ---------- + `rnds`: tensor of random number betweenn 0 and 1 + `n_strat`: tensor with number of stratifications in each dimension + `n_ev`: tensor containing number of sample per each hypercube + `hypercubes`: tensor containing all different hypercube + `divisions`: vegas grid + + Returns + ------- + `x` : random numbers collocated in hypercubes + `w` : weight of each event + div_index: division index in which each (n_dim) set of random numbers fall + `segm` : segmentantion for later computations + """ + points = tf.repeat(hypercubes, n_ev, axis=0) + xn = tf.transpose((points+tf.transpose(randoms))*FBINS/float_me(n_strat)) + segm = tf.cast(tf.repeat(tf.range(fzero, + tf.shape(hypercubes)[0]), + n_ev), + DTYPEINT) + + # same as vegasflow generate_random_array + @tf.function + def digest(xn): + ind_i = tf.cast(xn, DTYPEINT) + # Get the value of the left and right sides of the bins + ind_f = ind_i + ione + x_ini = tf.gather(divisions, ind_i, batch_dims=1) + x_fin = tf.gather(divisions, ind_f, batch_dims=1) + # Compute the width of the bins + xdelta = x_fin - x_ini + return ind_i, x_ini, xdelta + + ind_i, x_ini, xdelta = digest(xn) + # Compute the random number between 0 and 1 + # This is the heavy part of the calc + + @tf.function + def compute_x(x_ini, xn, xdelta): + aux_rand = xn - tf.math.floor(xn) + return x_ini + xdelta * aux_rand + + x = compute_x(x_ini, xn, xdelta) + # Compute the random number between the limits + # x = reg_i + rand_x * (reg_f - reg_i) + # and the weight + weights = tf.reduce_prod(xdelta * FBINS, axis=0) + x_t = tf.transpose(x) + int_xn = tf.transpose(ind_i) + return x_t, int_xn, weights, segm + + +class VegasFlowPlus(VegasFlow): + """ + Implementation of the VEGAS+ algorithm + """ + + def __init__(self, n_dim, n_events, train=True, adaptive=True, **kwargs): + super().__init__(n_dim, n_events, train, **kwargs) + + self.adaptive = adaptive + + # Initialize stratifications + self.n_strat = tf.math.floor(tf.math.pow(self.n_events/N_STRAT_MIN, 1/n_dim)) + self.n_strat = int_me(self.n_strat) + + # Initialize hypercubes + hypercubes_one_dim = np.arange(0, int(self.n_strat)) + hypercubes = [list(p) for p in product(hypercubes_one_dim, + repeat=int(n_dim))] + self.hypercubes = tf.convert_to_tensor(hypercubes, dtype=DTYPE) + + if len(hypercubes) != int(tf.math.pow(self.n_strat, n_dim)): + raise ValueError("Hypercubes problem!") + + # Initialize n_ev + n_ev = int_me(tf.math.floordiv(self.n_events, tf.math.pow(self.n_strat, n_dim))) + n_ev = tf.math.maximum(n_ev, 2) + self.n_ev = tf.fill([1, len(hypercubes)], n_ev) + self.n_ev = tf.reshape(self.n_ev, [-1]) + + self.n_events = int(tf.reduce_sum(self.n_ev)) + self.xjac = 1 / self.n_events + + def redistribute_samples(self, arr_var): + """Receives an array with the variance of the integrand in each + hypercube and recalculate the samples per hypercube according + to VEGAS+ algorithm""" + + # neval_eff = int(self.n_events/2) + damped_arr_sdev = tf.pow(arr_var, BETA/2) + new_n_ev = tf.maximum(2, + damped_arr_sdev * self.n_events/tf.reduce_sum(damped_arr_sdev)/2) + self.n_ev = int_me(tf.math.floor(new_n_ev)) + self.n_events = int(tf.reduce_sum(self.n_ev)) + self.xjac = 1 / self.n_events + + def _run_event(self, integrand, ncalls=None): + + if ncalls is None: + n_events = self.n_events + else: + n_events = ncalls + + tech_cut = 1e-8 + # Generate all random number for this iteration + rnds = tf.random.uniform( + (self.n_dim, n_events), minval=tech_cut, maxval=1.0 - tech_cut, dtype=DTYPE + ) + + # Pass random numbers in hypercubes + + x, ind, w, segm = generate_samples_in_hypercubes(rnds, + self.n_strat, + self.n_ev, + self.hypercubes, + self.divisions) + + # compute integrand + xjac = self.xjac * w + if self.simplify_signature: + tmp = xjac * integrand(x) + else: + tmp = xjac * integrand(x, n_dim=self.n_dim, weight=xjac) + tmp2 = tf.square(tmp) + + # tensor containing resummed component for each hypercubes + ress = tf.math.segment_sum(tmp, segm) + ress2 = tf.math.segment_sum(tmp2, segm) + + Fn_ev = tf.cast(self.n_ev, DTYPE) + arr_var = (ress2 - tf.square(ress)/Fn_ev)/(Fn_ev - fone) + + arr_res2 = [] + if self.train: + # If the training is active, save the result of the integral sq + for j in range(self.n_dim): + arr_res2.append( + consume_array_into_indices(tmp2, ind[ :, j : j + 1], int_me(self.grid_bins - 1)) + ) + arr_res2 = tf.reshape(arr_res2, (self.n_dim, -1)) + + return ress, arr_var, arr_res2 + + def _iteration_content(self): + + # print("_iteration_content_vfp") + ress, arr_var, arr_res2 = self.run_event() + Fn_ev = tf.cast(self.n_ev, DTYPE) + # compute variance for each hypercube + sigmas2 = tf.maximum(arr_var, fzero) + + res = tf.reduce_sum(ress) + + sigma2 = tf.reduce_sum(sigmas2*Fn_ev) + sigma = tf.sqrt(sigma2) + + # If adaptive is active redistributes samples + if self.adaptive: + self.redistribute_samples(arr_var) + + if self.train: + self.refine_grid(arr_res2) + return res, sigma + + +def vegasflowplus_wrapper(integrand, n_dim, n_iter, total_n_events, **kwargs): + """Convenience wrapper + + Parameters + ---------- + `integrand`: tf.function + `n_dim`: number of dimensions + `n_iter`: number of iterations + `n_events`: number of events per iteration + + Returns + ------- + `final_result`: integral value + `sigma`: monte carlo error + """ + return wrapper(VegasFlowPlus, integrand, n_dim, n_iter, total_n_events, **kwargs) From 588d7ca1a943198ba6bb1da1cf02028288dab7c1 Mon Sep 17 00:00:00 2001 From: Andrea Pasquale Date: Mon, 8 Mar 2021 15:24:10 +0100 Subject: [PATCH 08/27] fixed problems in stratified.py and vflowplus.py --- src/vegasflow/stratified.py | 107 +++++++++++++++---------------- src/vegasflow/tests/test_algs.py | 20 +++++- src/vegasflow/vflowplus.py | 69 ++++++++++++-------- 3 files changed, 112 insertions(+), 84 deletions(-) diff --git a/src/vegasflow/stratified.py b/src/vegasflow/stratified.py index 65f5837..cdb84ff 100644 --- a/src/vegasflow/stratified.py +++ b/src/vegasflow/stratified.py @@ -1,13 +1,12 @@ """ - Implementation of the stratified sampling + Implementation of the adaptive stratified sampling """ -from itertools import chain, repeat, product +from itertools import product import tensorflow as tf import numpy as np -from vegasflow.configflow import DTYPE, fone, fzero, DTYPEINT +from vegasflow.configflow import DTYPE, fone, fzero, DTYPEINT, float_me, int_me from vegasflow.monte_carlo import MonteCarloFlow, wrapper -N_STRAT_MIN = 4 BETA = 0.75 @@ -32,60 +31,71 @@ def generate_random_array(rnds, n_strat, n_ev, hypercubes): """ # stratification width delta_y = tf.cast(1/n_strat, DTYPE) - random = rnds*delta_y points = tf.cast(tf.repeat(hypercubes, n_ev, axis=0), DTYPE) - n_dim = int(tf.shape(hypercubes)[1]) - x = points*delta_y + random - tf.cast([delta_y] * n_dim, DTYPE) + x = points*delta_y + random w = tf.cast(tf.repeat(1/n_ev, n_ev), DTYPE) segm = tf.cast(tf.repeat(tf.range(fzero, tf.shape(hypercubes)[0]), n_ev), DTYPEINT) + return x, w, segm class StratifiedFlow(MonteCarloFlow): """ - Simple Monte Carlo integrator with Stratified Sampling. + Monte Carlo integrator with Adaptive Stratified Sampling. """ def __init__(self, n_dim, n_events, adaptive=True, **kwargs): super().__init__(n_dim, n_events, **kwargs) + self.init_calls = n_events self.adaptive = adaptive # Initialize stratifications - self.n_strat = tf.math.floor(tf.math.pow(n_events/N_STRAT_MIN, - 1/n_dim)) - substratification_np = np.linspace(0, 1, int(self.n_strat)+1) - stratifications_np = substratification_np.repeat(n_dim).reshape(-1, - n_dim).T - self.stratifications = tf.Variable(stratifications_np, dtype=DTYPE) + if self.adaptive: + neval_eff = int(self.n_events/2) + self.n_strat = tf.math.floor(tf.math.pow(neval_eff/2, 1/n_dim)) + else: + neval_eff = self.n_events + self.n_strat = tf.math.floor(tf.math.pow(neval_eff/2, 1/n_dim)) + + self.n_strat = int_me(self.n_strat) + # Initialize hypercubes - hypercubes_one_dim = np.arange(1, int(self.n_strat)+1) + hypercubes_one_dim = np.arange(0, int(self.n_strat)) hypercubes = [list(p) for p in product(hypercubes_one_dim, repeat=int(n_dim))] - hypercubes_tf = tf.convert_to_tensor(hypercubes) - self.hypercubes = hypercubes_tf + self.hypercubes = tf.convert_to_tensor(hypercubes) + + if len(hypercubes) != int(tf.math.pow(self.n_strat, n_dim)): + raise ValueError("Hypercubes problem!") + + # Set min evaluations per hypercube + self.min_neval_hcube = int(neval_eff // len(hypercubes)) + if self.min_neval_hcube < 2: + self.min_neval_hcube = 2 # Initialize n_ev - n_ev = tf.math.floor(n_events/tf.math.pow(self.n_strat, n_dim)) - self.n_ev = tf.fill([1, len(hypercubes)], int(n_ev)) + self.n_ev = tf.fill([1, len(hypercubes)], self.min_neval_hcube) self.n_ev = tf.reshape(self.n_ev, [-1]) - - # Correction of self.n_events due to samples per hypercube - self.n_events = int(n_ev*tf.math.pow(self.n_strat, n_dim)) + self.n_events = int(tf.reduce_sum(self.n_ev)) + self.xjac = float_me(1/len(hypercubes)) def redistribute_samples(self, arr_var): """Receives an array with the variance of the integrand in each hypercube and recalculate the samples per hypercube according - to VEGAS+ algorithm""" - - damped_arr_sdev = tf.pow(arr_var, BETA/2) - new_n_ev = tf.maximum(2, - damped_arr_sdev * self.n_events/tf.reduce_sum(damped_arr_sdev)) - self.n_ev = tf.math.floor(new_n_ev) + to the VEGAS+ algorithm""" + + damped_arr_sdev = tf.pow(arr_var, BETA / 2) + new_n_ev = tf.maximum(self.min_neval_hcube, + damped_arr_sdev + * self.init_calls / 2 + / tf.reduce_sum(damped_arr_sdev)) + self.n_ev = int_me(tf.math.floor(new_n_ev)) + self.n_events = int(tf.reduce_sum(self.n_ev)) def _run_event(self, integrand, ncalls=None): @@ -99,11 +109,13 @@ def _run_event(self, integrand, ncalls=None): dtype=DTYPE) # Pass random numbers in hypercubes - x, w, segm = generate_random_array(rnds, self.n_strat, self.n_ev, + x, w, segm = generate_random_array(rnds, + self.n_strat, + self.n_ev, self.hypercubes) # compute integrand - xjac = w + xjac = self.xjac * w tmp = integrand(x, n_dim=self.n_dim, weight=xjac) * xjac tmp2 = tf.square(tmp) @@ -111,37 +123,22 @@ def _run_event(self, integrand, ncalls=None): ress = tf.math.segment_sum(tmp, segm) ress2 = tf.math.segment_sum(tmp2, segm) - # if adaptive save variance of each hypercube - arr_var = None - if self.adaptive: - hypercube_volume = tf.cast(tf.math.pow(1/self.n_strat, self.n_dim), - DTYPE) - arr_var = ((ress2 * - tf.cast(self.n_ev, DTYPE) * - tf.square(hypercube_volume)) - - tf.square(ress*hypercube_volume)) - - return ress, ress2, arr_var + Fn_ev = tf.cast(self.n_ev, DTYPE) + arr_var = ress2 * Fn_ev - tf.square(ress) - def _run_iteration(self): + return ress, arr_var - ress, raw_ress2, arr_var = self.run_event() + def _run_iteration(self): - # compute variance for each hypercube - ress2 = raw_ress2 * tf.cast(self.n_ev, DTYPE) - err_tmp2s = (ress2 - tf.square(ress))/(tf.cast(self.n_ev, DTYPE)-fone) - sigmas2 = tf.maximum(err_tmp2s, fzero) - - res = tf.reduce_sum(ress) / tf.cast(tf.math.pow(self.n_strat, - self.n_dim), - DTYPE) + ress, arr_var = self.run_event() - sigma2 = tf.reduce_sum(sigmas2)/tf.cast(tf.math.pow(self.n_strat, - self.n_dim), - DTYPE) / self.n_events + Fn_ev = tf.cast(self.n_ev, DTYPE) + sigmas2 = tf.maximum(arr_var, fzero) + res = tf.reduce_sum(ress) + sigma2 = tf.reduce_sum(sigmas2/(Fn_ev-fone)) sigma = tf.sqrt(sigma2) - # If adaptive is active redistributes samples + # If adaptive is True redistributes samples if self.adaptive: self.redistribute_samples(arr_var) diff --git a/src/vegasflow/tests/test_algs.py b/src/vegasflow/tests/test_algs.py index a2032d7..45e100d 100644 --- a/src/vegasflow/tests/test_algs.py +++ b/src/vegasflow/tests/test_algs.py @@ -14,6 +14,7 @@ from vegasflow.vflow import VegasFlow from vegasflow.plain import PlainFlow from vegasflow.stratified import StratifiedFlow +from vegasflow.vflowplus import VegasFlowPlus # Test setup dim = 2 @@ -111,7 +112,24 @@ def test_PlainFlow(): result = plain_instance.run_integration(n_iter) check_is_one(result) -def test_StratifiedFlow(): +def test_StratifiedFlow_ADAPTIVE_SAMPLING(): stratified_instance = instance_and_compile(StratifiedFlow) result = stratified_instance.run_integration(n_iter) check_is_one(result) + +def test_StratifiedFlow_NOT_ADAPTIVE_SAMPLING(): + stratified_instance = StratifiedFlow(dim, ncalls, adaptive=False) + stratified_instance.compile(example_integrand) + result = stratified_instance.run_integration(n_iter) + check_is_one(result) + +def test_VegasFlowPlus_ADAPTIVE_SAMPLING(): + vflowplus_instance = instance_and_compile(VegasFlowPlus) + result = vflowplus_instance.run_integration(n_iter) + check_is_one(result) + +def test_VegasFlowPlus_NOT_ADAPTIVE_SAMPLING(): + vflowplus_instance = VegasFlowPlus(dim, ncalls, adaptive=False) + vflowplus_instance.compile(example_integrand) + result = vflowplus_instance.run_integration(n_iter) + check_is_one(result) diff --git a/src/vegasflow/vflowplus.py b/src/vegasflow/vflowplus.py index 6a540a2..70f644f 100644 --- a/src/vegasflow/vflowplus.py +++ b/src/vegasflow/vflowplus.py @@ -1,10 +1,11 @@ #!/usr/env python """ Implementation of vegas+ algorithm: - adaptive importance sampling + adaptive stratified sampling - + + adaptive importance sampling + adaptive stratified sampling + """ -from itertools import chain, repeat, product +from itertools import product import tensorflow as tf import numpy as np from vegasflow.configflow import DTYPE, DTYPEINT, fone, fzero, float_me, ione, int_me, BINS_MAX @@ -12,14 +13,16 @@ from vegasflow.vflow import VegasFlow from vegasflow.utils import consume_array_into_indices - -N_STRAT_MIN = 4 BETA = 0.75 FBINS = float_me(BINS_MAX) @tf.function -def generate_samples_in_hypercubes(randoms, n_strat, n_ev, hypercubes, divisions): +def generate_samples_in_hypercubes(randoms, + n_strat, + n_ev, + hypercubes, + divisions): """Receives an array of random numbers 0 and 1 and distribute them in each hypercube according to the number of sample in each hypercube specified by n_ev @@ -39,7 +42,9 @@ def generate_samples_in_hypercubes(randoms, n_strat, n_ev, hypercubes, divisions div_index: division index in which each (n_dim) set of random numbers fall `segm` : segmentantion for later computations """ + points = tf.repeat(hypercubes, n_ev, axis=0) + n_evs = float_me(tf.repeat(n_ev, n_ev)) xn = tf.transpose((points+tf.transpose(randoms))*FBINS/float_me(n_strat)) segm = tf.cast(tf.repeat(tf.range(fzero, tf.shape(hypercubes)[0]), @@ -72,9 +77,10 @@ def compute_x(x_ini, xn, xdelta): # x = reg_i + rand_x * (reg_f - reg_i) # and the weight weights = tf.reduce_prod(xdelta * FBINS, axis=0) + final_weights = weights/n_evs x_t = tf.transpose(x) int_xn = tf.transpose(ind_i) - return x_t, int_xn, weights, segm + return x_t, int_xn, final_weights, segm class VegasFlowPlus(VegasFlow): @@ -85,10 +91,17 @@ class VegasFlowPlus(VegasFlow): def __init__(self, n_dim, n_events, train=True, adaptive=True, **kwargs): super().__init__(n_dim, n_events, train, **kwargs) + self.init_calls = n_events self.adaptive = adaptive - + # Initialize stratifications - self.n_strat = tf.math.floor(tf.math.pow(self.n_events/N_STRAT_MIN, 1/n_dim)) + if self.adaptive: + neval_eff = int(self.n_events/2) + self.n_strat = tf.math.floor(tf.math.pow(neval_eff/2, 1/n_dim)) + else: + neval_eff = self.n_events + self.n_strat = tf.math.floor(tf.math.pow(neval_eff/2, 1/n_dim)) + self.n_strat = int_me(self.n_strat) # Initialize hypercubes @@ -100,27 +113,28 @@ def __init__(self, n_dim, n_events, train=True, adaptive=True, **kwargs): if len(hypercubes) != int(tf.math.pow(self.n_strat, n_dim)): raise ValueError("Hypercubes problem!") - # Initialize n_ev - n_ev = int_me(tf.math.floordiv(self.n_events, tf.math.pow(self.n_strat, n_dim))) - n_ev = tf.math.maximum(n_ev, 2) - self.n_ev = tf.fill([1, len(hypercubes)], n_ev) + self.min_neval_hcube = int(neval_eff // len(hypercubes)) + if self.min_neval_hcube < 2: + self.min_neval_hcube = 2 + + self.n_ev = tf.fill([1, len(hypercubes)], self.min_neval_hcube) self.n_ev = tf.reshape(self.n_ev, [-1]) - self.n_events = int(tf.reduce_sum(self.n_ev)) - self.xjac = 1 / self.n_events + self.xjac = float_me(1/len(hypercubes)) def redistribute_samples(self, arr_var): """Receives an array with the variance of the integrand in each hypercube and recalculate the samples per hypercube according to VEGAS+ algorithm""" - - # neval_eff = int(self.n_events/2) + damped_arr_sdev = tf.pow(arr_var, BETA/2) - new_n_ev = tf.maximum(2, - damped_arr_sdev * self.n_events/tf.reduce_sum(damped_arr_sdev)/2) + new_n_ev = tf.maximum(self.min_neval_hcube, + damped_arr_sdev + * self.init_calls / 2 + / tf.reduce_sum(damped_arr_sdev)) + self.n_ev = int_me(tf.math.floor(new_n_ev)) self.n_events = int(tf.reduce_sum(self.n_ev)) - self.xjac = 1 / self.n_events def _run_event(self, integrand, ncalls=None): @@ -132,8 +146,10 @@ def _run_event(self, integrand, ncalls=None): tech_cut = 1e-8 # Generate all random number for this iteration rnds = tf.random.uniform( - (self.n_dim, n_events), minval=tech_cut, maxval=1.0 - tech_cut, dtype=DTYPE - ) + (self.n_dim, n_events), + minval=tech_cut, + maxval=1.0 - tech_cut, + dtype=DTYPE) # Pass random numbers in hypercubes @@ -156,7 +172,7 @@ def _run_event(self, integrand, ncalls=None): ress2 = tf.math.segment_sum(tmp2, segm) Fn_ev = tf.cast(self.n_ev, DTYPE) - arr_var = (ress2 - tf.square(ress)/Fn_ev)/(Fn_ev - fone) + arr_var = ress2*Fn_ev - tf.square(ress) arr_res2 = [] if self.train: @@ -171,15 +187,12 @@ def _run_event(self, integrand, ncalls=None): def _iteration_content(self): - # print("_iteration_content_vfp") ress, arr_var, arr_res2 = self.run_event() + Fn_ev = tf.cast(self.n_ev, DTYPE) - # compute variance for each hypercube sigmas2 = tf.maximum(arr_var, fzero) - res = tf.reduce_sum(ress) - - sigma2 = tf.reduce_sum(sigmas2*Fn_ev) + sigma2 = tf.reduce_sum(sigmas2/(Fn_ev-fone)) sigma = tf.sqrt(sigma2) # If adaptive is active redistributes samples From b4bee755c0feafa1738f58c9a9c8b017ab965958 Mon Sep 17 00:00:00 2001 From: Andrea Pasquale Date: Fri, 12 Mar 2021 09:29:09 +0100 Subject: [PATCH 09/27] first attempt at fixing retracing --- examples/retracing.py | 43 ++++++++++++++++++++++++++++++++++++++ src/vegasflow/vflowplus.py | 37 ++++++++++++++++++++++---------- 2 files changed, 69 insertions(+), 11 deletions(-) create mode 100644 examples/retracing.py diff --git a/examples/retracing.py b/examples/retracing.py new file mode 100644 index 0000000..a8bffb8 --- /dev/null +++ b/examples/retracing.py @@ -0,0 +1,43 @@ +""" + Retracing example in VegasFlowPlus +""" + +from vegasflow import VegasFlowPlus, VegasFlow, StratifiedFlow, PlainFlow +from vegasflow.configflow import DTYPE, DTYPEINT,run_eager, float_me +import time +import numpy as np +import tensorflow as tf +#import pineappl + +# MC integration setup +dim = 2 +ncalls = np.int32(1e3) +n_iter = 5 + +#@tf.function(input_signature=[ +# tf.TensorSpec(shape=[None,dim], dtype=DTYPE), +# tf.TensorSpec(shape=[], dtype=DTYPE), +# tf.TensorSpec(shape=[None], dtype=DTYPE) +# ] +# ) +def symgauss(xarr, n_dim=None,weight=None, **kwargs): + """symgauss test function""" + if n_dim is None: + n_dim = xarr.shape[-1] + a = tf.constant(0.1, dtype=DTYPE) + n100 = tf.cast(100 * n_dim, dtype=DTYPE) + pref = tf.pow(1.0 / a / np.sqrt(np.pi), n_dim) + coef = tf.reduce_sum(tf.range(n100 + 1)) + coef += tf.reduce_sum(tf.square((xarr - 1.0 / 2.0) / a), axis=1) + coef -= (n100 + 1) * n100 / 2.0 + return pref * tf.exp(-coef) + + + +if __name__ == "__main__": + """Testing several different integrations""" + + #run_eager() + vegas_instance = VegasFlowPlus(dim, ncalls,adaptive=True) + vegas_instance.compile(symgauss) + vegas_instance.run_integration(n_iter) diff --git a/src/vegasflow/vflowplus.py b/src/vegasflow/vflowplus.py index 70f644f..52ca530 100644 --- a/src/vegasflow/vflowplus.py +++ b/src/vegasflow/vflowplus.py @@ -17,8 +17,16 @@ FBINS = float_me(BINS_MAX) -@tf.function -def generate_samples_in_hypercubes(randoms, +@tf.function(input_signature=[ + tf.TensorSpec(shape=[],dtype=np.int32), + tf.TensorSpec(shape=[],dtype=np.int32), + tf.TensorSpec(shape=[],dtype=DTYPEINT), + tf.TensorSpec(shape=[None],dtype=DTYPEINT), + tf.TensorSpec(shape=[None,None],dtype=DTYPE), + tf.TensorSpec(shape=[None,None],dtype=DTYPE) +]) +def generate_samples_in_hypercubes(n_dim, + n_events, n_strat, n_ev, hypercubes, @@ -42,10 +50,17 @@ def generate_samples_in_hypercubes(randoms, div_index: division index in which each (n_dim) set of random numbers fall `segm` : segmentantion for later computations """ + tech_cut = 1e-8 + # Generate all random number for this iteration + rnds = tf.random.uniform( + (n_dim, n_events), + minval=tech_cut, + maxval=1.0 - tech_cut, + dtype=DTYPE) points = tf.repeat(hypercubes, n_ev, axis=0) n_evs = float_me(tf.repeat(n_ev, n_ev)) - xn = tf.transpose((points+tf.transpose(randoms))*FBINS/float_me(n_strat)) + xn = tf.transpose((points+tf.transpose(rnds))*FBINS/float_me(n_strat)) segm = tf.cast(tf.repeat(tf.range(fzero, tf.shape(hypercubes)[0]), n_ev), @@ -143,17 +158,17 @@ def _run_event(self, integrand, ncalls=None): else: n_events = ncalls - tech_cut = 1e-8 + #tech_cut = 1e-8 # Generate all random number for this iteration - rnds = tf.random.uniform( - (self.n_dim, n_events), - minval=tech_cut, - maxval=1.0 - tech_cut, - dtype=DTYPE) + #rnds = tf.random.uniform( + # (self.n_dim, n_events), + # minval=tech_cut, + # maxval=1.0 - tech_cut, + # dtype=DTYPE) # Pass random numbers in hypercubes - - x, ind, w, segm = generate_samples_in_hypercubes(rnds, + x, ind, w, segm = generate_samples_in_hypercubes(self.n_dim, + n_events, self.n_strat, self.n_ev, self.hypercubes, From d750c3c2c5674e287be425a75179e2c409286aa6 Mon Sep 17 00:00:00 2001 From: juacrumar Date: Mon, 15 Mar 2021 17:03:41 +0100 Subject: [PATCH 10/27] suggested changes --- examples/retracing.py | 18 ++-- src/vegasflow/monte_carlo.py | 9 +- src/vegasflow/vflowplus.py | 163 ++++++++++++++++------------------- 3 files changed, 90 insertions(+), 100 deletions(-) diff --git a/examples/retracing.py b/examples/retracing.py index a8bffb8..5747c1e 100644 --- a/examples/retracing.py +++ b/examples/retracing.py @@ -14,19 +14,19 @@ ncalls = np.int32(1e3) n_iter = 5 -#@tf.function(input_signature=[ -# tf.TensorSpec(shape=[None,dim], dtype=DTYPE), -# tf.TensorSpec(shape=[], dtype=DTYPE), -# tf.TensorSpec(shape=[None], dtype=DTYPE) -# ] -# ) -def symgauss(xarr, n_dim=None,weight=None, **kwargs): +@tf.function(input_signature=[ + tf.TensorSpec(shape=[None,dim], dtype=DTYPE), + tf.TensorSpec(shape=[], dtype=DTYPEINT), + tf.TensorSpec(shape=[None], dtype=DTYPE) + ] + ) +def symgauss(xarr, n_dim=None, weight=None, **kwargs): """symgauss test function""" if n_dim is None: n_dim = xarr.shape[-1] a = tf.constant(0.1, dtype=DTYPE) n100 = tf.cast(100 * n_dim, dtype=DTYPE) - pref = tf.pow(1.0 / a / np.sqrt(np.pi), n_dim) + pref = tf.pow(1.0 / a / np.sqrt(np.pi), float_me(n_dim)) coef = tf.reduce_sum(tf.range(n100 + 1)) coef += tf.reduce_sum(tf.square((xarr - 1.0 / 2.0) / a), axis=1) coef -= (n100 + 1) * n100 / 2.0 @@ -37,7 +37,7 @@ def symgauss(xarr, n_dim=None,weight=None, **kwargs): if __name__ == "__main__": """Testing several different integrations""" - #run_eager() +# run_eager() vegas_instance = VegasFlowPlus(dim, ncalls,adaptive=True) vegas_instance.compile(symgauss) vegas_instance.run_integration(n_iter) diff --git a/src/vegasflow/monte_carlo.py b/src/vegasflow/monte_carlo.py index c50837b..0ba36c7 100644 --- a/src/vegasflow/monte_carlo.py +++ b/src/vegasflow/monte_carlo.py @@ -39,7 +39,7 @@ import joblib import numpy as np import tensorflow as tf -from vegasflow.configflow import MAX_EVENTS_LIMIT, DEFAULT_ACTIVE_DEVICES, DTYPE, TECH_CUT, float_me +from vegasflow.configflow import MAX_EVENTS_LIMIT, DEFAULT_ACTIVE_DEVICES, DTYPE, TECH_CUT, float_me, int_me import logging @@ -269,7 +269,7 @@ def set_distribute(self, queue_object): self.cluster = queue_object self.distribute = True - def run_event(self, **kwargs): + def run_event(self, tensorize_events=False, **kwargs): """ Runs the Monte Carlo event. This corresponds to a number of calls decided by the `events_per_run` variable. The variable `acc` is exposed @@ -297,7 +297,10 @@ def run_event(self, **kwargs): ncalls = min(events_left, self.events_per_run) pc += ncalls / self.n_events * 100 percentages.append(pc) - events_to_do.append(ncalls) + if tensorize_events: + events_to_do.append(int_me(ncalls)) + else: + events_to_do.append(ncalls) events_left -= self.events_per_run if self.devices: diff --git a/src/vegasflow/vflowplus.py b/src/vegasflow/vflowplus.py index 52ca530..e886fb5 100644 --- a/src/vegasflow/vflowplus.py +++ b/src/vegasflow/vflowplus.py @@ -13,24 +13,42 @@ from vegasflow.vflow import VegasFlow from vegasflow.utils import consume_array_into_indices +import logging + +logger = logging.getLogger(__name__) + BETA = 0.75 FBINS = float_me(BINS_MAX) - -@tf.function(input_signature=[ - tf.TensorSpec(shape=[],dtype=np.int32), - tf.TensorSpec(shape=[],dtype=np.int32), - tf.TensorSpec(shape=[],dtype=DTYPEINT), - tf.TensorSpec(shape=[None],dtype=DTYPEINT), - tf.TensorSpec(shape=[None,None],dtype=DTYPE), - tf.TensorSpec(shape=[None,None],dtype=DTYPE) -]) -def generate_samples_in_hypercubes(n_dim, - n_events, - n_strat, - n_ev, - hypercubes, - divisions): +@tf.function(input_signature=3 * [tf.TensorSpec(shape=[None, None], dtype=DTYPE)]) +def _compute_x(x_ini, xn, xdelta): + """ Helper function for ``generate_samples_in_hypercubes`` """ + aux_rand = xn - tf.math.floor(xn) + return x_ini + xdelta * aux_rand + +# same as vegasflow generate_random_array +@tf.function(input_signature=2*[tf.TensorSpec(shape=[None, None], dtype=DTYPE)]) +def _digest(xn, divisions): + ind_i = tf.cast(xn, DTYPEINT) + # Get the value of the left and right sides of the bins + ind_f = ind_i + ione + x_ini = tf.gather(divisions, ind_i, batch_dims=1) + x_fin = tf.gather(divisions, ind_f, batch_dims=1) + # Compute the width of the bins + xdelta = x_fin - x_ini + return ind_i, x_ini, xdelta + + +@tf.function( + input_signature=[ + tf.TensorSpec(shape=[None, None], dtype=DTYPE), + tf.TensorSpec(shape=[], dtype=DTYPEINT), + tf.TensorSpec(shape=[None], dtype=DTYPEINT), + tf.TensorSpec(shape=[None, None], dtype=DTYPE), + tf.TensorSpec(shape=[None, None], dtype=DTYPE), + ] +) +def generate_samples_in_hypercubes(rnds, n_strat, n_ev, hypercubes, divisions): """Receives an array of random numbers 0 and 1 and distribute them in each hypercube according to the number of sample in each hypercube specified by n_ev @@ -50,49 +68,22 @@ def generate_samples_in_hypercubes(n_dim, div_index: division index in which each (n_dim) set of random numbers fall `segm` : segmentantion for later computations """ - tech_cut = 1e-8 - # Generate all random number for this iteration - rnds = tf.random.uniform( - (n_dim, n_events), - minval=tech_cut, - maxval=1.0 - tech_cut, - dtype=DTYPE) points = tf.repeat(hypercubes, n_ev, axis=0) n_evs = float_me(tf.repeat(n_ev, n_ev)) - xn = tf.transpose((points+tf.transpose(rnds))*FBINS/float_me(n_strat)) - segm = tf.cast(tf.repeat(tf.range(fzero, - tf.shape(hypercubes)[0]), - n_ev), - DTYPEINT) - - # same as vegasflow generate_random_array - @tf.function - def digest(xn): - ind_i = tf.cast(xn, DTYPEINT) - # Get the value of the left and right sides of the bins - ind_f = ind_i + ione - x_ini = tf.gather(divisions, ind_i, batch_dims=1) - x_fin = tf.gather(divisions, ind_f, batch_dims=1) - # Compute the width of the bins - xdelta = x_fin - x_ini - return ind_i, x_ini, xdelta - - ind_i, x_ini, xdelta = digest(xn) + xn = tf.transpose((points + tf.transpose(rnds)) * FBINS / float_me(n_strat)) + segm = tf.cast(tf.repeat(tf.range(fzero, tf.shape(hypercubes)[0]), n_ev), DTYPEINT) + + ind_i, x_ini, xdelta = _digest(xn, divisions) # Compute the random number between 0 and 1 # This is the heavy part of the calc - @tf.function - def compute_x(x_ini, xn, xdelta): - aux_rand = xn - tf.math.floor(xn) - return x_ini + xdelta * aux_rand - - x = compute_x(x_ini, xn, xdelta) + x = _compute_x(x_ini, xn, xdelta) # Compute the random number between the limits # x = reg_i + rand_x * (reg_f - reg_i) # and the weight weights = tf.reduce_prod(xdelta * FBINS, axis=0) - final_weights = weights/n_evs + final_weights = weights / n_evs x_t = tf.transpose(x) int_xn = tf.transpose(ind_i) return x_t, int_xn, final_weights, segm @@ -108,21 +99,20 @@ def __init__(self, n_dim, n_events, train=True, adaptive=True, **kwargs): self.init_calls = n_events self.adaptive = adaptive - + # Initialize stratifications if self.adaptive: - neval_eff = int(self.n_events/2) - self.n_strat = tf.math.floor(tf.math.pow(neval_eff/2, 1/n_dim)) + neval_eff = int(self.n_events / 2) + self.n_strat = tf.math.floor(tf.math.pow(neval_eff / 2, 1 / n_dim)) else: neval_eff = self.n_events - self.n_strat = tf.math.floor(tf.math.pow(neval_eff/2, 1/n_dim)) + self.n_strat = tf.math.floor(tf.math.pow(neval_eff / 2, 1 / n_dim)) self.n_strat = int_me(self.n_strat) # Initialize hypercubes hypercubes_one_dim = np.arange(0, int(self.n_strat)) - hypercubes = [list(p) for p in product(hypercubes_one_dim, - repeat=int(n_dim))] + hypercubes = [list(p) for p in product(hypercubes_one_dim, repeat=int(n_dim))] self.hypercubes = tf.convert_to_tensor(hypercubes, dtype=DTYPE) if len(hypercubes) != int(tf.math.pow(self.n_strat, n_dim)): @@ -131,83 +121,75 @@ def __init__(self, n_dim, n_events, train=True, adaptive=True, **kwargs): self.min_neval_hcube = int(neval_eff // len(hypercubes)) if self.min_neval_hcube < 2: self.min_neval_hcube = 2 - + self.n_ev = tf.fill([1, len(hypercubes)], self.min_neval_hcube) self.n_ev = tf.reshape(self.n_ev, [-1]) self.n_events = int(tf.reduce_sum(self.n_ev)) - self.xjac = float_me(1/len(hypercubes)) + self.xjac = float_me(1 / len(hypercubes)) def redistribute_samples(self, arr_var): """Receives an array with the variance of the integrand in each hypercube and recalculate the samples per hypercube according to VEGAS+ algorithm""" - damped_arr_sdev = tf.pow(arr_var, BETA/2) - new_n_ev = tf.maximum(self.min_neval_hcube, - damped_arr_sdev - * self.init_calls / 2 - / tf.reduce_sum(damped_arr_sdev)) + damped_arr_sdev = tf.pow(arr_var, BETA / 2) + # TODO: what if arr_var is negative? + new_n_ev = tf.maximum( + self.min_neval_hcube, + damped_arr_sdev * self.init_calls / 2 / tf.reduce_sum(damped_arr_sdev), + ) self.n_ev = int_me(tf.math.floor(new_n_ev)) self.n_events = int(tf.reduce_sum(self.n_ev)) - - def _run_event(self, integrand, ncalls=None): - if ncalls is None: - n_events = self.n_events - else: - n_events = ncalls + def _run_event(self, integrand, ncalls=None, n_ev=None): + # NOTE: needs to receive both ncalls and n_ev + n_events = ncalls - #tech_cut = 1e-8 + tech_cut = 1e-8 # Generate all random number for this iteration - #rnds = tf.random.uniform( - # (self.n_dim, n_events), - # minval=tech_cut, - # maxval=1.0 - tech_cut, - # dtype=DTYPE) + rnds = tf.random.uniform( + (self.n_dim, n_events), minval=tech_cut, maxval=1.0 - tech_cut, dtype=DTYPE + ) # Pass random numbers in hypercubes - x, ind, w, segm = generate_samples_in_hypercubes(self.n_dim, - n_events, - self.n_strat, - self.n_ev, - self.hypercubes, - self.divisions) - + x, ind, w, segm = generate_samples_in_hypercubes( + rnds, self.n_strat, n_ev, self.hypercubes, self.divisions + ) + # compute integrand xjac = self.xjac * w if self.simplify_signature: tmp = xjac * integrand(x) else: tmp = xjac * integrand(x, n_dim=self.n_dim, weight=xjac) - tmp2 = tf.square(tmp) + tmp2 = tf.square(tmp) # tensor containing resummed component for each hypercubes ress = tf.math.segment_sum(tmp, segm) ress2 = tf.math.segment_sum(tmp2, segm) Fn_ev = tf.cast(self.n_ev, DTYPE) - arr_var = ress2*Fn_ev - tf.square(ress) + arr_var = ress2 * Fn_ev - tf.square(ress) arr_res2 = [] if self.train: # If the training is active, save the result of the integral sq for j in range(self.n_dim): arr_res2.append( - consume_array_into_indices(tmp2, ind[ :, j : j + 1], int_me(self.grid_bins - 1)) + consume_array_into_indices(tmp2, ind[:, j : j + 1], int_me(self.grid_bins - 1)) ) arr_res2 = tf.reshape(arr_res2, (self.n_dim, -1)) - + return ress, arr_var, arr_res2 - + def _iteration_content(self): - - ress, arr_var, arr_res2 = self.run_event() + ress, arr_var, arr_res2 = self.run_event(n_ev=self.n_ev) Fn_ev = tf.cast(self.n_ev, DTYPE) sigmas2 = tf.maximum(arr_var, fzero) res = tf.reduce_sum(ress) - sigma2 = tf.reduce_sum(sigmas2/(Fn_ev-fone)) + sigma2 = tf.reduce_sum(sigmas2 / (Fn_ev - fone)) sigma = tf.sqrt(sigma2) # If adaptive is active redistributes samples @@ -218,6 +200,11 @@ def _iteration_content(self): self.refine_grid(arr_res2) return res, sigma + def run_event(self, tensorize_events=True, **kwargs): + """ Tensorizes the number of events so they are not python or numpy primitives """ + logger.warning("Variable number of events requires function signatures all across") + return super().run_event(tensorize_events=tensorize_events, **kwargs) + def vegasflowplus_wrapper(integrand, n_dim, n_iter, total_n_events, **kwargs): """Convenience wrapper From 26b5650a0aa031f4042b66d711f517864630d09e Mon Sep 17 00:00:00 2001 From: juacrumar Date: Mon, 15 Mar 2021 17:15:12 +0100 Subject: [PATCH 11/27] typo --- src/vegasflow/vflowplus.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/vegasflow/vflowplus.py b/src/vegasflow/vflowplus.py index e886fb5..fe8793e 100644 --- a/src/vegasflow/vflowplus.py +++ b/src/vegasflow/vflowplus.py @@ -169,7 +169,7 @@ def _run_event(self, integrand, ncalls=None, n_ev=None): ress = tf.math.segment_sum(tmp, segm) ress2 = tf.math.segment_sum(tmp2, segm) - Fn_ev = tf.cast(self.n_ev, DTYPE) + Fn_ev = tf.cast(n_ev, DTYPE) arr_var = ress2 * Fn_ev - tf.square(ress) arr_res2 = [] From dd7d54ca9aec807eedb8cf1f99df5cbcbb48faec Mon Sep 17 00:00:00 2001 From: Andrea Pasquale Date: Tue, 16 Mar 2021 10:53:22 +0100 Subject: [PATCH 12/27] fixed self.adaptive issue --- src/vegasflow/vflowplus.py | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/src/vegasflow/vflowplus.py b/src/vegasflow/vflowplus.py index fe8793e..b4beecf 100644 --- a/src/vegasflow/vflowplus.py +++ b/src/vegasflow/vflowplus.py @@ -127,13 +127,15 @@ def __init__(self, n_dim, n_events, train=True, adaptive=True, **kwargs): self.n_events = int(tf.reduce_sum(self.n_ev)) self.xjac = float_me(1 / len(hypercubes)) + if self.adaptive: + logger.warning("Variable number of events requires function signatures all across") + def redistribute_samples(self, arr_var): """Receives an array with the variance of the integrand in each hypercube and recalculate the samples per hypercube according to VEGAS+ algorithm""" damped_arr_sdev = tf.pow(arr_var, BETA / 2) - # TODO: what if arr_var is negative? new_n_ev = tf.maximum( self.min_neval_hcube, damped_arr_sdev * self.init_calls / 2 / tf.reduce_sum(damped_arr_sdev), @@ -200,10 +202,9 @@ def _iteration_content(self): self.refine_grid(arr_res2) return res, sigma - def run_event(self, tensorize_events=True, **kwargs): - """ Tensorizes the number of events so they are not python or numpy primitives """ - logger.warning("Variable number of events requires function signatures all across") - return super().run_event(tensorize_events=tensorize_events, **kwargs) + def run_event(self, tensorize_events=None, **kwargs): + """ Tensorizes the number of events so they are not python or numpy primitives if self.adaptive=True""" + return super().run_event(tensorize_events=self.adaptive, **kwargs) def vegasflowplus_wrapper(integrand, n_dim, n_iter, total_n_events, **kwargs): From ad897270d69a1a995ac3ff3f74d1562409c736ad Mon Sep 17 00:00:00 2001 From: Andrea Pasquale Date: Wed, 14 Apr 2021 15:17:56 +0200 Subject: [PATCH 13/27] fixed crash problem in vflowplus --- src/vegasflow/configflow.py | 3 ++- src/vegasflow/vflowplus.py | 47 ++++++++++++++++++++----------------- 2 files changed, 27 insertions(+), 23 deletions(-) diff --git a/src/vegasflow/configflow.py b/src/vegasflow/configflow.py index f5c0e3c..17ccd80 100644 --- a/src/vegasflow/configflow.py +++ b/src/vegasflow/configflow.py @@ -92,10 +92,11 @@ def float_me(i): BINS_MAX = 50 ALPHA = 1.5 TECH_CUT = float_me(1e-8) +BETA = 0.75 # Set up the logistics of the integration # Events Limit limits how many events are done in one single run of the event_loop # set it lower if hitting memory problems MAX_EVENTS_LIMIT = int(1e6) # Select the list of devices to look for -DEFAULT_ACTIVE_DEVICES = ["GPU"] # , 'CPU'] +DEFAULT_ACTIVE_DEVICES = ["GPU"] # , 'CPU'] diff --git a/src/vegasflow/vflowplus.py b/src/vegasflow/vflowplus.py index b4beecf..75498f5 100644 --- a/src/vegasflow/vflowplus.py +++ b/src/vegasflow/vflowplus.py @@ -1,31 +1,31 @@ #!/usr/env python """ Implementation of vegas+ algorithm: - + adaptive importance sampling + adaptive stratified sampling - """ from itertools import product +import logging import tensorflow as tf import numpy as np -from vegasflow.configflow import DTYPE, DTYPEINT, fone, fzero, float_me, ione, int_me, BINS_MAX +from vegasflow.configflow import DTYPE, DTYPEINT, fone, fzero, float_me, ione, int_me, BINS_MAX, BETA from vegasflow.monte_carlo import wrapper from vegasflow.vflow import VegasFlow from vegasflow.utils import consume_array_into_indices -import logging logger = logging.getLogger(__name__) -BETA = 0.75 FBINS = float_me(BINS_MAX) + @tf.function(input_signature=3 * [tf.TensorSpec(shape=[None, None], dtype=DTYPE)]) def _compute_x(x_ini, xn, xdelta): """ Helper function for ``generate_samples_in_hypercubes`` """ aux_rand = xn - tf.math.floor(xn) return x_ini + xdelta * aux_rand + # same as vegasflow generate_random_array @tf.function(input_signature=2*[tf.TensorSpec(shape=[None, None], dtype=DTYPE)]) def _digest(xn, divisions): @@ -44,20 +44,20 @@ def _digest(xn, divisions): tf.TensorSpec(shape=[None, None], dtype=DTYPE), tf.TensorSpec(shape=[], dtype=DTYPEINT), tf.TensorSpec(shape=[None], dtype=DTYPEINT), - tf.TensorSpec(shape=[None, None], dtype=DTYPE), + tf.TensorSpec(shape=[None, None], dtype=DTYPEINT), tf.TensorSpec(shape=[None, None], dtype=DTYPE), ] ) def generate_samples_in_hypercubes(rnds, n_strat, n_ev, hypercubes, divisions): """Receives an array of random numbers 0 and 1 and distribute them in each hypercube according to the - number of sample in each hypercube specified by n_ev + number of samples in each hypercube specified by n_ev Parameters ---------- - `rnds`: tensor of random number betweenn 0 and 1 + `rnds`: tensor of random number between 0 and 1 `n_strat`: tensor with number of stratifications in each dimension - `n_ev`: tensor containing number of sample per each hypercube + `n_ev`: tensor containing number of sample per hypercube `hypercubes`: tensor containing all different hypercube `divisions`: vegas grid @@ -65,14 +65,15 @@ def generate_samples_in_hypercubes(rnds, n_strat, n_ev, hypercubes, divisions): ------- `x` : random numbers collocated in hypercubes `w` : weight of each event - div_index: division index in which each (n_dim) set of random numbers fall + `ind`: division index in which each (n_dim) set of random numbers fall `segm` : segmentantion for later computations """ - points = tf.repeat(hypercubes, n_ev, axis=0) - n_evs = float_me(tf.repeat(n_ev, n_ev)) + indices = tf.repeat(tf.range(tf.shape(hypercubes, out_type=DTYPEINT)[0]), n_ev) + points = float_me(tf.gather(hypercubes, indices)) + n_evs = float_me(tf.gather(n_ev, indices)) xn = tf.transpose((points + tf.transpose(rnds)) * FBINS / float_me(n_strat)) - segm = tf.cast(tf.repeat(tf.range(fzero, tf.shape(hypercubes)[0]), n_ev), DTYPEINT) + segm = indices ind_i, x_ini, xdelta = _digest(xn, divisions) # Compute the random number between 0 and 1 @@ -86,6 +87,7 @@ def generate_samples_in_hypercubes(rnds, n_strat, n_ev, hypercubes, divisions): final_weights = weights / n_evs x_t = tf.transpose(x) int_xn = tf.transpose(ind_i) + return x_t, int_xn, final_weights, segm @@ -99,7 +101,6 @@ def __init__(self, n_dim, n_events, train=True, adaptive=True, **kwargs): self.init_calls = n_events self.adaptive = adaptive - # Initialize stratifications if self.adaptive: neval_eff = int(self.n_events / 2) @@ -108,22 +109,25 @@ def __init__(self, n_dim, n_events, train=True, adaptive=True, **kwargs): neval_eff = self.n_events self.n_strat = tf.math.floor(tf.math.pow(neval_eff / 2, 1 / n_dim)) + if tf.math.pow(self.n_strat, n_dim) > 1e4: + self.n_strat = tf.math.floor(tf.math.pow(1e4, 1/n_dim)) + self.n_strat = int_me(self.n_strat) # Initialize hypercubes hypercubes_one_dim = np.arange(0, int(self.n_strat)) hypercubes = [list(p) for p in product(hypercubes_one_dim, repeat=int(n_dim))] - self.hypercubes = tf.convert_to_tensor(hypercubes, dtype=DTYPE) + self.hypercubes = tf.convert_to_tensor(hypercubes, dtype=DTYPEINT) if len(hypercubes) != int(tf.math.pow(self.n_strat, n_dim)): - raise ValueError("Hypercubes problem!") + raise ValueError("Hypercubes are not equal to n_strat^n_dim") self.min_neval_hcube = int(neval_eff // len(hypercubes)) if self.min_neval_hcube < 2: self.min_neval_hcube = 2 self.n_ev = tf.fill([1, len(hypercubes)], self.min_neval_hcube) - self.n_ev = tf.reshape(self.n_ev, [-1]) + self.n_ev = tf.cast(tf.reshape(self.n_ev, [-1]), dtype=DTYPEINT) self.n_events = int(tf.reduce_sum(self.n_ev)) self.xjac = float_me(1 / len(hypercubes)) @@ -134,18 +138,17 @@ def redistribute_samples(self, arr_var): """Receives an array with the variance of the integrand in each hypercube and recalculate the samples per hypercube according to VEGAS+ algorithm""" - damped_arr_sdev = tf.pow(arr_var, BETA / 2) new_n_ev = tf.maximum( self.min_neval_hcube, damped_arr_sdev * self.init_calls / 2 / tf.reduce_sum(damped_arr_sdev), ) - - self.n_ev = int_me(tf.math.floor(new_n_ev)) + self.n_ev = int_me(new_n_ev) self.n_events = int(tf.reduce_sum(self.n_ev)) def _run_event(self, integrand, ncalls=None, n_ev=None): # NOTE: needs to receive both ncalls and n_ev + n_events = ncalls tech_cut = 1e-8 @@ -156,7 +159,7 @@ def _run_event(self, integrand, ncalls=None, n_ev=None): # Pass random numbers in hypercubes x, ind, w, segm = generate_samples_in_hypercubes( - rnds, self.n_strat, n_ev, self.hypercubes, self.divisions + rnds, self.n_strat, n_ev, self.hypercubes, self.divisions, ) # compute integrand @@ -179,7 +182,7 @@ def _run_event(self, integrand, ncalls=None, n_ev=None): # If the training is active, save the result of the integral sq for j in range(self.n_dim): arr_res2.append( - consume_array_into_indices(tmp2, ind[:, j : j + 1], int_me(self.grid_bins - 1)) + consume_array_into_indices(tmp2, ind[: , j : j + 1], int_me(self.grid_bins - 1)) ) arr_res2 = tf.reshape(arr_res2, (self.n_dim, -1)) From 016abe237fd731659643c0575fe6605a2a68bd7e Mon Sep 17 00:00:00 2001 From: Andrea Pasquale Date: Wed, 14 Apr 2021 15:25:57 +0200 Subject: [PATCH 14/27] set adaptive off is dim > 13 --- src/vegasflow/vflowplus.py | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/src/vegasflow/vflowplus.py b/src/vegasflow/vflowplus.py index 75498f5..5bdbf67 100644 --- a/src/vegasflow/vflowplus.py +++ b/src/vegasflow/vflowplus.py @@ -96,11 +96,17 @@ class VegasFlowPlus(VegasFlow): Implementation of the VEGAS+ algorithm """ - def __init__(self, n_dim, n_events, train=True, adaptive=True, **kwargs): + def __init__(self, n_dim, n_events, train=True, adaptive=None, **kwargs): super().__init__(n_dim, n_events, train, **kwargs) self.init_calls = n_events - self.adaptive = adaptive + + # naive check not to use adaptive if n_dim > 13 + if n_dim > 13 and self.adaptive == None: + self.adaptive = False + else: + self.adaptive = adaptive + # Initialize stratifications if self.adaptive: neval_eff = int(self.n_events / 2) From 38b87a0ac6ee783c734846aa82b563d147e742ea Mon Sep 17 00:00:00 2001 From: Andrea Pasquale Date: Wed, 14 Apr 2021 15:33:32 +0200 Subject: [PATCH 15/27] typo --- src/vegasflow/vflowplus.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/vegasflow/vflowplus.py b/src/vegasflow/vflowplus.py index 5bdbf67..08401ec 100644 --- a/src/vegasflow/vflowplus.py +++ b/src/vegasflow/vflowplus.py @@ -102,7 +102,7 @@ def __init__(self, n_dim, n_events, train=True, adaptive=None, **kwargs): self.init_calls = n_events # naive check not to use adaptive if n_dim > 13 - if n_dim > 13 and self.adaptive == None: + if n_dim > 13 and adaptive == None: self.adaptive = False else: self.adaptive = adaptive From 462373e8aa7220e0288544209e982d9d2c4a5376 Mon Sep 17 00:00:00 2001 From: Andrea Pasquale Date: Thu, 22 Apr 2021 08:51:44 +0200 Subject: [PATCH 16/27] removed StratifiedFlow --- examples/example_stratified.py | 42 --------- src/vegasflow/stratified.py | 150 --------------------------------- 2 files changed, 192 deletions(-) delete mode 100644 examples/example_stratified.py delete mode 100644 src/vegasflow/stratified.py diff --git a/examples/example_stratified.py b/examples/example_stratified.py deleted file mode 100644 index 894e05c..0000000 --- a/examples/example_stratified.py +++ /dev/null @@ -1,42 +0,0 @@ -""" - Example: integration with stratified sampling - - Basic example using the stratified_wrapper helper -""" - -from vegasflow.configflow import DTYPE -import time -import numpy as np -import tensorflow as tf -from vegasflow.vflow import vegas_wrapper -from vegasflow.plain import plain_wrapper -from vegasflow.stratified import stratified_wrapper - -# MC integration setup -dim = 4 -ncalls = np.int32(1e5) -n_iter = 5 - - -@tf.function -def symgauss(xarr, n_dim=None, **kwargs): - """symgauss test function""" - if n_dim is None: - n_dim = xarr.shape[-1] - a = tf.constant(0.1, dtype=DTYPE) - n100 = tf.cast(100 * n_dim, dtype=DTYPE) - pref = tf.pow(1.0 / a / np.sqrt(np.pi), n_dim) - coef = tf.reduce_sum(tf.range(n100 + 1)) - coef += tf.reduce_sum(tf.square((xarr - 1.0 / 2.0) / a), axis=1) - coef -= (n100 + 1) * n100 / 2.0 - return pref * tf.exp(-coef) - - -if __name__ == "__main__": - """Testing several different integrations""" - print(f"VEGAS MC, ncalls={ncalls}:") - start = time.time() - ncalls = 10*ncalls - r = stratified_wrapper(symgauss, dim, n_iter, ncalls) - end = time.time() - print(f"Vegas took: time (s): {end-start}") \ No newline at end of file diff --git a/src/vegasflow/stratified.py b/src/vegasflow/stratified.py deleted file mode 100644 index cdb84ff..0000000 --- a/src/vegasflow/stratified.py +++ /dev/null @@ -1,150 +0,0 @@ -""" - Implementation of the adaptive stratified sampling -""" -from itertools import product -import tensorflow as tf -import numpy as np -from vegasflow.configflow import DTYPE, fone, fzero, DTYPEINT, float_me, int_me -from vegasflow.monte_carlo import MonteCarloFlow, wrapper - -BETA = 0.75 - - -@tf.function -def generate_random_array(rnds, n_strat, n_ev, hypercubes): - """Receives an array of random numbers 0 and 1 and - distribute them in each hypercube according to the - number of sample in each hypercube specified by n_ev - - Parameters - ---------- - `rnds`: tensor of random number betweenn 0 and 1 - `n_strat`: tensor with number of stratifications in each dimension - `n_ev`: tensor containing number of sample per each hypercube - `hypercubes`: tensor containing all different hypercube - - Returns - ------- - `x` : random numbers collocated in hypercubes - `w` : weight of each event - `segm` : segmentantion for later computations - """ - # stratification width - delta_y = tf.cast(1/n_strat, DTYPE) - random = rnds*delta_y - points = tf.cast(tf.repeat(hypercubes, n_ev, axis=0), DTYPE) - x = points*delta_y + random - w = tf.cast(tf.repeat(1/n_ev, n_ev), DTYPE) - segm = tf.cast(tf.repeat(tf.range(fzero, - tf.shape(hypercubes)[0]), - n_ev), - DTYPEINT) - - return x, w, segm - - -class StratifiedFlow(MonteCarloFlow): - """ - Monte Carlo integrator with Adaptive Stratified Sampling. - """ - - def __init__(self, n_dim, n_events, adaptive=True, **kwargs): - super().__init__(n_dim, n_events, **kwargs) - - self.init_calls = n_events - self.adaptive = adaptive - - # Initialize stratifications - if self.adaptive: - neval_eff = int(self.n_events/2) - self.n_strat = tf.math.floor(tf.math.pow(neval_eff/2, 1/n_dim)) - else: - neval_eff = self.n_events - self.n_strat = tf.math.floor(tf.math.pow(neval_eff/2, 1/n_dim)) - - self.n_strat = int_me(self.n_strat) - - # Initialize hypercubes - hypercubes_one_dim = np.arange(0, int(self.n_strat)) - hypercubes = [list(p) for p in product(hypercubes_one_dim, - repeat=int(n_dim))] - self.hypercubes = tf.convert_to_tensor(hypercubes) - - if len(hypercubes) != int(tf.math.pow(self.n_strat, n_dim)): - raise ValueError("Hypercubes problem!") - - # Set min evaluations per hypercube - self.min_neval_hcube = int(neval_eff // len(hypercubes)) - if self.min_neval_hcube < 2: - self.min_neval_hcube = 2 - - # Initialize n_ev - self.n_ev = tf.fill([1, len(hypercubes)], self.min_neval_hcube) - self.n_ev = tf.reshape(self.n_ev, [-1]) - self.n_events = int(tf.reduce_sum(self.n_ev)) - self.xjac = float_me(1/len(hypercubes)) - - def redistribute_samples(self, arr_var): - """Receives an array with the variance of the integrand in each - hypercube and recalculate the samples per hypercube according - to the VEGAS+ algorithm""" - - damped_arr_sdev = tf.pow(arr_var, BETA / 2) - new_n_ev = tf.maximum(self.min_neval_hcube, - damped_arr_sdev - * self.init_calls / 2 - / tf.reduce_sum(damped_arr_sdev)) - self.n_ev = int_me(tf.math.floor(new_n_ev)) - self.n_events = int(tf.reduce_sum(self.n_ev)) - - def _run_event(self, integrand, ncalls=None): - - if ncalls is None: - n_events = self.n_events - else: - n_events = ncalls - - # Generate all random numbers - rnds = tf.random.uniform((n_events, self.n_dim), minval=0, maxval=1, - dtype=DTYPE) - - # Pass random numbers in hypercubes - x, w, segm = generate_random_array(rnds, - self.n_strat, - self.n_ev, - self.hypercubes) - - # compute integrand - xjac = self.xjac * w - tmp = integrand(x, n_dim=self.n_dim, weight=xjac) * xjac - tmp2 = tf.square(tmp) - - # tensor containing resummed component for each hypercubes - ress = tf.math.segment_sum(tmp, segm) - ress2 = tf.math.segment_sum(tmp2, segm) - - Fn_ev = tf.cast(self.n_ev, DTYPE) - arr_var = ress2 * Fn_ev - tf.square(ress) - - return ress, arr_var - - def _run_iteration(self): - - ress, arr_var = self.run_event() - - Fn_ev = tf.cast(self.n_ev, DTYPE) - sigmas2 = tf.maximum(arr_var, fzero) - res = tf.reduce_sum(ress) - sigma2 = tf.reduce_sum(sigmas2/(Fn_ev-fone)) - sigma = tf.sqrt(sigma2) - - # If adaptive is True redistributes samples - if self.adaptive: - self.redistribute_samples(arr_var) - - return res, sigma - - -def stratified_wrapper(*args): - """ Wrapper around PlainFlow """ - return wrapper(StratifiedFlow, *args) From 18ca36926e858396565be8ade951f353d4aa8267 Mon Sep 17 00:00:00 2001 From: Andrea Pasquale Date: Thu, 22 Apr 2021 08:56:52 +0200 Subject: [PATCH 17/27] file update from previous commit --- examples/retracing.py | 5 ++--- src/vegasflow/__init__.py | 1 - src/vegasflow/configflow.py | 2 ++ src/vegasflow/tests/test_algs.py | 12 ------------ src/vegasflow/vflowplus.py | 4 ++-- 5 files changed, 6 insertions(+), 18 deletions(-) diff --git a/examples/retracing.py b/examples/retracing.py index 5747c1e..7cff7d2 100644 --- a/examples/retracing.py +++ b/examples/retracing.py @@ -2,12 +2,11 @@ Retracing example in VegasFlowPlus """ -from vegasflow import VegasFlowPlus, VegasFlow, StratifiedFlow, PlainFlow -from vegasflow.configflow import DTYPE, DTYPEINT,run_eager, float_me +from vegasflow import VegasFlowPlus, VegasFlow, PlainFlow +from vegasflow.configflow import DTYPE, DTYPEINT, run_eager, float_me import time import numpy as np import tensorflow as tf -#import pineappl # MC integration setup dim = 2 diff --git a/src/vegasflow/__init__.py b/src/vegasflow/__init__.py index 7af16fb..c6596de 100644 --- a/src/vegasflow/__init__.py +++ b/src/vegasflow/__init__.py @@ -4,7 +4,6 @@ # Expose the main interfaces from vegasflow.vflow import VegasFlow, vegas_wrapper from vegasflow.plain import PlainFlow, plain_wrapper -from vegasflow.stratified import StratifiedFlow, stratified_wrapper from vegasflow.vflowplus import VegasFlowPlus, vegasflowplus_wrapper __version__ = "1.2.0" diff --git a/src/vegasflow/configflow.py b/src/vegasflow/configflow.py index 17ccd80..8bc71e9 100644 --- a/src/vegasflow/configflow.py +++ b/src/vegasflow/configflow.py @@ -98,5 +98,7 @@ def float_me(i): # Events Limit limits how many events are done in one single run of the event_loop # set it lower if hitting memory problems MAX_EVENTS_LIMIT = int(1e6) +# Maximum number of evaluation per hypercube for the class VegasFlowPlus +MAX_NEVAL_HCUBE = int(1e4) # Select the list of devices to look for DEFAULT_ACTIVE_DEVICES = ["GPU"] # , 'CPU'] diff --git a/src/vegasflow/tests/test_algs.py b/src/vegasflow/tests/test_algs.py index 385981c..d56941f 100644 --- a/src/vegasflow/tests/test_algs.py +++ b/src/vegasflow/tests/test_algs.py @@ -13,7 +13,6 @@ from vegasflow.configflow import DTYPE from vegasflow.vflow import VegasFlow from vegasflow.plain import PlainFlow -from vegasflow.stratified import StratifiedFlow from vegasflow.vflowplus import VegasFlowPlus # Test setup @@ -112,17 +111,6 @@ def test_PlainFlow(): result = plain_instance.run_integration(n_iter) check_is_one(result) -def test_StratifiedFlow_ADAPTIVE_SAMPLING(): - stratified_instance = instance_and_compile(StratifiedFlow) - result = stratified_instance.run_integration(n_iter) - check_is_one(result) - -def test_StratifiedFlow_NOT_ADAPTIVE_SAMPLING(): - stratified_instance = StratifiedFlow(dim, ncalls, adaptive=False) - stratified_instance.compile(example_integrand) - result = stratified_instance.run_integration(n_iter) - check_is_one(result) - def test_VegasFlowPlus_ADAPTIVE_SAMPLING(): vflowplus_instance = instance_and_compile(VegasFlowPlus) result = vflowplus_instance.run_integration(n_iter) diff --git a/src/vegasflow/vflowplus.py b/src/vegasflow/vflowplus.py index 08401ec..59b456c 100644 --- a/src/vegasflow/vflowplus.py +++ b/src/vegasflow/vflowplus.py @@ -8,7 +8,7 @@ import logging import tensorflow as tf import numpy as np -from vegasflow.configflow import DTYPE, DTYPEINT, fone, fzero, float_me, ione, int_me, BINS_MAX, BETA +from vegasflow.configflow import DTYPE, DTYPEINT, fone, fzero, float_me, ione, int_me, BINS_MAX, BETA, MAX_NEVAL_HCUBE from vegasflow.monte_carlo import wrapper from vegasflow.vflow import VegasFlow from vegasflow.utils import consume_array_into_indices @@ -115,7 +115,7 @@ def __init__(self, n_dim, n_events, train=True, adaptive=None, **kwargs): neval_eff = self.n_events self.n_strat = tf.math.floor(tf.math.pow(neval_eff / 2, 1 / n_dim)) - if tf.math.pow(self.n_strat, n_dim) > 1e4: + if tf.math.pow(self.n_strat, n_dim) > MAX_NEVAL_HCUBE: self.n_strat = tf.math.floor(tf.math.pow(1e4, 1/n_dim)) self.n_strat = int_me(self.n_strat) From 07cf119f1f13f20266cee462665637467ed4ef99 Mon Sep 17 00:00:00 2001 From: Andrea Pasquale Date: Tue, 22 Jun 2021 14:55:33 +0200 Subject: [PATCH 18/27] all events in same device for vegas+ Co-authored-by: Juan M. Cruz-Martinez --- src/vegasflow/vflowplus.py | 1 + 1 file changed, 1 insertion(+) diff --git a/src/vegasflow/vflowplus.py b/src/vegasflow/vflowplus.py index 59b456c..b2c1a75 100644 --- a/src/vegasflow/vflowplus.py +++ b/src/vegasflow/vflowplus.py @@ -97,6 +97,7 @@ class VegasFlowPlus(VegasFlow): """ def __init__(self, n_dim, n_events, train=True, adaptive=None, **kwargs): + _ = kwargs.setdefault("events_limit", n_events) super().__init__(n_dim, n_events, train, **kwargs) self.init_calls = n_events From f577231ed4cef9389c02fbad89bf2193b7330a11 Mon Sep 17 00:00:00 2001 From: juacrumar Date: Mon, 16 Aug 2021 21:53:47 +0200 Subject: [PATCH 19/27] fixed some todos --- src/vegasflow/__init__.py | 3 +-- src/vegasflow/configflow.py | 25 +++++++------------------ src/vegasflow/tests/test_algs.py | 11 ++++++++--- src/vegasflow/vflowplus.py | 16 ++++++++++++++++ 4 files changed, 32 insertions(+), 23 deletions(-) diff --git a/src/vegasflow/__init__.py b/src/vegasflow/__init__.py index 53362b7..0a1b3fe 100644 --- a/src/vegasflow/__init__.py +++ b/src/vegasflow/__init__.py @@ -4,7 +4,6 @@ # Expose the main interfaces from vegasflow.vflow import VegasFlow, vegas_wrapper, vegas_sampler from vegasflow.plain import PlainFlow, plain_wrapper, plain_sampler -from vegasflow.vflowplus import VegasFlowPlus, vegasflowplus_wrapper -# TODO: create a sampler for vegasflowplus +from vegasflow.vflowplus import VegasFlowPlus, vegasflowplus_wrapper, vegasflowplus_sampler __version__ = "1.2.2" diff --git a/src/vegasflow/configflow.py b/src/vegasflow/configflow.py index 28ba664..d823079 100644 --- a/src/vegasflow/configflow.py +++ b/src/vegasflow/configflow.py @@ -9,11 +9,17 @@ # Some global parameters BINS_MAX = 50 ALPHA = 1.5 +BETA = 0.75 # Vegas + TECH_CUT = 1e-8 + # Set up the logistics of the integration +# set the limits lower if hitting memory problems + # Events Limit limits how many events are done in one single run of the event_loop -# set it lower if hitting memory problems MAX_EVENTS_LIMIT = int(1e6) +# Maximum number of evaluation per hypercube for VegasFlowPlus +MAX_NEVAL_HCUBE = int(1e4) + # Select the list of devices to look for DEFAULT_ACTIVE_DEVICES = ["GPU"] # , 'CPU'] @@ -103,20 +109,3 @@ def float_me(i): izero = int_me(0) fone = float_me(1) fzero = float_me(0) - - -# TODO: everything below should be moved to the right classes -# Define some default parameters for the MC algorithms -BINS_MAX = 50 -ALPHA = 1.5 -TECH_CUT = float_me(1e-8) -BETA = 0.75 - -# Set up the logistics of the integration -# Events Limit limits how many events are done in one single run of the event_loop -# set it lower if hitting memory problems -MAX_EVENTS_LIMIT = int(1e6) -# Maximum number of evaluation per hypercube for the class VegasFlowPlus -MAX_NEVAL_HCUBE = int(1e4) -# Select the list of devices to look for -DEFAULT_ACTIVE_DEVICES = ["GPU"] # , 'CPU'] diff --git a/src/vegasflow/tests/test_algs.py b/src/vegasflow/tests/test_algs.py index e82f1f6..837e13b 100644 --- a/src/vegasflow/tests/test_algs.py +++ b/src/vegasflow/tests/test_algs.py @@ -67,6 +67,7 @@ def check_is_one(result, sigmas=3): def test_VegasFlow(): + """ Test VegasFlow class, importance sampling algorithm""" for mode in range(4): vegas_instance = instance_and_compile(VegasFlow, mode) _ = vegas_instance.run_integration(n_iter) @@ -91,6 +92,7 @@ def test_VegasFlow(): def test_VegasFlow_save_grid(): + """ Test the grid saving feature of vegasflow """ tmp_filename = tempfile.mktemp() vegas_instance = instance_and_compile(VegasFlow) # Run an iteration so the grid is not trivial @@ -176,11 +178,14 @@ def test_rng_generation(n_events=100): _ = helper_rng_tester(v, n_events) def test_VegasFlowPlus_ADAPTIVE_SAMPLING(): - vflowplus_instance = instance_and_compile(VegasFlowPlus) - result = vflowplus_instance.run_integration(n_iter) - check_is_one(result) + """ Test Vegasflow with Adaptive Sampling on (the default) """ + for mode in range(4): + vflowplus_instance = instance_and_compile(VegasFlowPlus, mode) + result = vflowplus_instance.run_integration(n_iter) + check_is_one(result) def test_VegasFlowPlus_NOT_ADAPTIVE_SAMPLING(): + """ Test Vegasflow with Adaptive Sampling off (non-default) """ vflowplus_instance = VegasFlowPlus(dim, ncalls, adaptive=False) vflowplus_instance.compile(example_integrand) result = vflowplus_instance.run_integration(n_iter) diff --git a/src/vegasflow/vflowplus.py b/src/vegasflow/vflowplus.py index 995c2f4..95fb207 100644 --- a/src/vegasflow/vflowplus.py +++ b/src/vegasflow/vflowplus.py @@ -230,3 +230,19 @@ def vegasflowplus_wrapper(integrand, n_dim, n_iter, total_n_events, **kwargs): `sigma`: monte carlo error """ return wrapper(VegasFlowPlus, integrand, n_dim, n_iter, total_n_events, **kwargs) + +def vegasflowplus_sampler(*args, **kwargs): + """Convenience wrapper for sampling random numbers + + Parameters + ---------- + `integrand`: tf.function + `n_dim`: number of dimensions + `n_events`: number of events per iteration + `training_steps`: number of training_iterations + + Returns + ------- + `sampler`: a reference to the generate_random_array method of the integrator class + """ + return sampler(VegasFlowPlus, *args, **kwargs) From 0b9ae0bf88eb7dfa4d97687706766b323877d5f5 Mon Sep 17 00:00:00 2001 From: juacrumar Date: Mon, 16 Aug 2021 23:15:16 +0200 Subject: [PATCH 20/27] stylize --- src/vegasflow/configflow.py | 2 +- src/vegasflow/monte_carlo.py | 2 +- src/vegasflow/vflow.py | 2 +- src/vegasflow/vflowplus.py | 134 +++++++++++++++++++++++------------ 4 files changed, 91 insertions(+), 49 deletions(-) diff --git a/src/vegasflow/configflow.py b/src/vegasflow/configflow.py index d823079..cc01659 100644 --- a/src/vegasflow/configflow.py +++ b/src/vegasflow/configflow.py @@ -9,7 +9,7 @@ # Some global parameters BINS_MAX = 50 ALPHA = 1.5 -BETA = 0.75 # Vegas + +BETA = 0.75 # Vegas+ TECH_CUT = 1e-8 # Set up the logistics of the integration diff --git a/src/vegasflow/monte_carlo.py b/src/vegasflow/monte_carlo.py index 4821d6e..13fbd10 100644 --- a/src/vegasflow/monte_carlo.py +++ b/src/vegasflow/monte_carlo.py @@ -48,7 +48,7 @@ DTYPEINT, TECH_CUT, float_me, - int_me + int_me, ) diff --git a/src/vegasflow/vflow.py b/src/vegasflow/vflow.py index 5cc7b14..5737fe2 100644 --- a/src/vegasflow/vflow.py +++ b/src/vegasflow/vflow.py @@ -1,4 +1,3 @@ -#!/usr/env python """ This module contains the VegasFlow class and all its auxuliary functions @@ -341,6 +340,7 @@ def _run_event(self, integrand, ncalls=None): return res, res2, arr_res2 def _iteration_content(self): + """Steps to follow per iteration""" # Compute the result res, res2, arr_res2 = self.run_event() # Compute the error diff --git a/src/vegasflow/vflowplus.py b/src/vegasflow/vflowplus.py index 95fb207..0b48f9f 100644 --- a/src/vegasflow/vflowplus.py +++ b/src/vegasflow/vflowplus.py @@ -1,18 +1,32 @@ -#!/usr/env python """ Implementation of vegas+ algorithm: adaptive importance sampling + adaptive stratified sampling + from https://arxiv.org/abs/2009.05112 + + The main interface is the `VegasFlowPlus` class. """ from itertools import product -import logging -import tensorflow as tf import numpy as np -from vegasflow.configflow import DTYPE, DTYPEINT, fone, fzero, float_me, ione, int_me, BINS_MAX, BETA, MAX_NEVAL_HCUBE -from vegasflow.monte_carlo import wrapper +import tensorflow as tf + +from vegasflow.configflow import ( + DTYPE, + DTYPEINT, + fone, + fzero, + float_me, + ione, + int_me, + BINS_MAX, + BETA, + MAX_NEVAL_HCUBE, +) +from vegasflow.monte_carlo import wrapper, sampler from vegasflow.vflow import VegasFlow from vegasflow.utils import consume_array_into_indices +import logging logger = logging.getLogger(__name__) @@ -21,13 +35,13 @@ @tf.function(input_signature=3 * [tf.TensorSpec(shape=[None, None], dtype=DTYPE)]) def _compute_x(x_ini, xn, xdelta): - """ Helper function for ``generate_samples_in_hypercubes`` """ + """Helper function for ``generate_samples_in_hypercubes``""" aux_rand = xn - tf.math.floor(xn) return x_ini + xdelta * aux_rand # same as vegasflow generate_random_array -@tf.function(input_signature=2*[tf.TensorSpec(shape=[None, None], dtype=DTYPE)]) +@tf.function(input_signature=2 * [tf.TensorSpec(shape=[None, None], dtype=DTYPE)]) def _digest(xn, divisions): ind_i = tf.cast(xn, DTYPEINT) # Get the value of the left and right sides of the bins @@ -48,7 +62,7 @@ def _digest(xn, divisions): tf.TensorSpec(shape=[None, None], dtype=DTYPE), ] ) -def generate_samples_in_hypercubes(rnds, n_strat, n_ev, hypercubes, divisions): +def generate_samples_in_hypercubes(rnds_t, n_strat, n_ev, hypercubes, divisions): """Receives an array of random numbers 0 and 1 and distribute them in each hypercube according to the number of samples in each hypercube specified by n_ev @@ -57,7 +71,7 @@ def generate_samples_in_hypercubes(rnds, n_strat, n_ev, hypercubes, divisions): ---------- `rnds`: tensor of random number between 0 and 1 `n_strat`: tensor with number of stratifications in each dimension - `n_ev`: tensor containing number of sample per hypercube + `n_ev`: tensor containing number of samples per hypercube `hypercubes`: tensor containing all different hypercube `divisions`: vegas grid @@ -68,7 +82,7 @@ def generate_samples_in_hypercubes(rnds, n_strat, n_ev, hypercubes, divisions): `ind`: division index in which each (n_dim) set of random numbers fall `segm` : segmentantion for later computations """ - + rnds = tf.transpose(rnds_t) indices = tf.repeat(tf.range(tf.shape(hypercubes, out_type=DTYPEINT)[0]), n_ev) points = float_me(tf.gather(hypercubes, indices)) n_evs = float_me(tf.gather(n_ev, indices)) @@ -100,33 +114,35 @@ def __init__(self, n_dim, n_events, train=True, adaptive=None, **kwargs): _ = kwargs.setdefault("events_limit", n_events) super().__init__(n_dim, n_events, train, **kwargs) - self.init_calls = n_events + # Save the initial number of events + self._init_calls = n_events - # naive check not to use adaptive if n_dim > 13 - if n_dim > 13 and adaptive == None: - self.adaptive = False + # Don't use adaptive if the number of dimension is too big + if n_dim > 13 and adaptive is None: + self._adaptive = False + logger.warning("Disabling adaptive mode from VegasFlowPlus, too many dimensions!") else: - self.adaptive = adaptive + self._adaptive = adaptive # Initialize stratifications - if self.adaptive: + if self._adaptive: neval_eff = int(self.n_events / 2) - self.n_strat = tf.math.floor(tf.math.pow(neval_eff / 2, 1 / n_dim)) + self._n_strat = tf.math.floor(tf.math.pow(neval_eff / 2, 1 / n_dim)) else: neval_eff = self.n_events - self.n_strat = tf.math.floor(tf.math.pow(neval_eff / 2, 1 / n_dim)) + self._n_strat = tf.math.floor(tf.math.pow(neval_eff / 2, 1 / n_dim)) - if tf.math.pow(self.n_strat, n_dim) > MAX_NEVAL_HCUBE: - self.n_strat = tf.math.floor(tf.math.pow(1e4, 1/n_dim)) + if tf.math.pow(self._n_strat, n_dim) > MAX_NEVAL_HCUBE: + self._n_strat = tf.math.floor(tf.math.pow(1e4, 1 / n_dim)) - self.n_strat = int_me(self.n_strat) + self._n_strat = int_me(self._n_strat) # Initialize hypercubes - hypercubes_one_dim = np.arange(0, int(self.n_strat)) + hypercubes_one_dim = np.arange(0, int(self._n_strat)) hypercubes = [list(p) for p in product(hypercubes_one_dim, repeat=int(n_dim))] - self.hypercubes = tf.convert_to_tensor(hypercubes, dtype=DTYPEINT) + self._hypercubes = tf.convert_to_tensor(hypercubes, dtype=DTYPEINT) - if len(hypercubes) != int(tf.math.pow(self.n_strat, n_dim)): + if len(hypercubes) != int(tf.math.pow(self._n_strat, n_dim)): raise ValueError("Hypercubes are not equal to n_strat^n_dim") self.min_neval_hcube = int(neval_eff // len(hypercubes)) @@ -138,7 +154,7 @@ def __init__(self, n_dim, n_events, train=True, adaptive=None, **kwargs): self._n_events = int(tf.reduce_sum(self.n_ev)) self.my_xjac = float_me(1 / len(hypercubes)) - if self.adaptive: + if self._adaptive: logger.warning("Variable number of events requires function signatures all across") def redistribute_samples(self, arr_var): @@ -148,29 +164,52 @@ def redistribute_samples(self, arr_var): damped_arr_sdev = tf.pow(arr_var, BETA / 2) new_n_ev = tf.maximum( self.min_neval_hcube, - damped_arr_sdev * self.init_calls / 2 / tf.reduce_sum(damped_arr_sdev), + damped_arr_sdev * self._init_calls / 2 / tf.reduce_sum(damped_arr_sdev), ) self.n_ev = int_me(new_n_ev) self.n_events = int(tf.reduce_sum(self.n_ev)) - def _run_event(self, integrand, ncalls=None, n_ev=None): - # NOTE: needs to receive both ncalls and n_ev - - n_events = ncalls - - tech_cut = 1e-8 - # Generate all random number for this iteration - rnds = tf.random.uniform( - (self.n_dim, n_events), minval=tech_cut, maxval=1.0 - tech_cut, dtype=DTYPE - ) + def _generate_random_array(self, n_events): + """Interface compatible with other algorithms dropping the segmentation in hypercubes""" + x, ind, w, _ = self._generate_random_array_plus(n_events, self.n_ev) + return x, ind, w - # Pass random numbers in hypercubes + def _generate_random_array_plus(self, n_events, n_ev): + """Generate a random array for a given number of events divided in hypercubes""" + rnds, _, _ = super(VegasFlow, self)._generate_random_array(n_events) + # Get random numbers from hypercubes x, ind, w, segm = generate_samples_in_hypercubes( - rnds, self.n_strat, n_ev, self.hypercubes, self.divisions, + rnds, + self._n_strat, + n_ev, + self._hypercubes, + self.divisions, ) + return x, ind, w * self.my_xjac, segm + + def _run_event(self, integrand, ncalls=None, n_ev=None): + """Run one step of VegasFlowPlus + Similar to the event step for importance sampling VegasFlow + adding the n_ev argument for the segmentation into hypercubes + n_ev is a tensor containing the number of samples per hypercube + + Parameters + ---------- + `integrand`: function to integrate + `ncalls`: how many events to run in this step + `n_ev`: number of samples per hypercube + + Returns + ------- + `res`: sum of the result of the integrand for all events + `res2`: sum of the result squared of the integrand for all events + `arr_res2`: result of the integrand squared per dimension and grid bin + """ + # NOTE: needs to receive both ncalls and n_ev + + x, ind, xjac, segm = self._generate_random_array_plus(ncalls, n_ev) # compute integrand - xjac = self.my_xjac * w tmp = xjac * integrand(x, weight=xjac) tmp2 = tf.square(tmp) @@ -186,32 +225,34 @@ def _run_event(self, integrand, ncalls=None, n_ev=None): # If the training is active, save the result of the integral sq for j in range(self.n_dim): arr_res2.append( - consume_array_into_indices(tmp2, ind[: , j : j + 1], int_me(self.grid_bins - 1)) + consume_array_into_indices(tmp2, ind[:, j : j + 1], int_me(self.grid_bins - 1)) ) arr_res2 = tf.reshape(arr_res2, (self.n_dim, -1)) return ress, arr_var, arr_res2 def _iteration_content(self): + """Steps to follow per iteration""" ress, arr_var, arr_res2 = self.run_event(n_ev=self.n_ev) - Fn_ev = tf.cast(self.n_ev, DTYPE) + # Compute the rror sigmas2 = tf.maximum(arr_var, fzero) res = tf.reduce_sum(ress) - sigma2 = tf.reduce_sum(sigmas2 / (Fn_ev - fone)) + sigma2 = tf.reduce_sum(sigmas2 / (float_me(self.n_ev) - fone)) sigma = tf.sqrt(sigma2) - # If adaptive is active redistributes samples - if self.adaptive: + # If adaptive is active redistribute the samples + if self._adaptive: self.redistribute_samples(arr_var) if self.train: self.refine_grid(arr_res2) + return res, sigma def run_event(self, tensorize_events=None, **kwargs): - """ Tensorizes the number of events so they are not python or numpy primitives if self.adaptive=True""" - return super().run_event(tensorize_events=self.adaptive, **kwargs) + """Tensorizes the number of events so they are not python or numpy primitives if self._adaptive=True""" + return super().run_event(tensorize_events=self._adaptive, **kwargs) def vegasflowplus_wrapper(integrand, n_dim, n_iter, total_n_events, **kwargs): @@ -231,6 +272,7 @@ def vegasflowplus_wrapper(integrand, n_dim, n_iter, total_n_events, **kwargs): """ return wrapper(VegasFlowPlus, integrand, n_dim, n_iter, total_n_events, **kwargs) + def vegasflowplus_sampler(*args, **kwargs): """Convenience wrapper for sampling random numbers From 9dadcf6505f842d1f9fe504223f98c43a946dfb4 Mon Sep 17 00:00:00 2001 From: juacrumar Date: Mon, 16 Aug 2021 23:15:48 +0200 Subject: [PATCH 21/27] bump version to 1.3.0 --- src/vegasflow/__init__.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/vegasflow/__init__.py b/src/vegasflow/__init__.py index 0a1b3fe..e3aed64 100644 --- a/src/vegasflow/__init__.py +++ b/src/vegasflow/__init__.py @@ -6,4 +6,4 @@ from vegasflow.plain import PlainFlow, plain_wrapper, plain_sampler from vegasflow.vflowplus import VegasFlowPlus, vegasflowplus_wrapper, vegasflowplus_sampler -__version__ = "1.2.2" +__version__ = "1.3.0" From b3b8c3614866eb9b1f32207cf22cda8adcdab7d9 Mon Sep 17 00:00:00 2001 From: juacrumar Date: Tue, 17 Aug 2021 14:42:34 +0200 Subject: [PATCH 22/27] more style --- src/vegasflow/plain.py | 2 +- src/vegasflow/utils.py | 39 ++++++++++++++++++++------------------ src/vegasflow/vflow.py | 1 - src/vegasflow/vflowplus.py | 11 ++++++----- 4 files changed, 28 insertions(+), 25 deletions(-) diff --git a/src/vegasflow/plain.py b/src/vegasflow/plain.py index f43e5b1..32d5433 100644 --- a/src/vegasflow/plain.py +++ b/src/vegasflow/plain.py @@ -2,7 +2,7 @@ Plain implementation of the plainest possible MonteCarlo """ -from vegasflow.configflow import DTYPE, fone, fzero +from vegasflow.configflow import fone, fzero from vegasflow.monte_carlo import MonteCarloFlow, wrapper, sampler import tensorflow as tf diff --git a/src/vegasflow/utils.py b/src/vegasflow/utils.py index c62b872..c748546 100644 --- a/src/vegasflow/utils.py +++ b/src/vegasflow/utils.py @@ -6,9 +6,13 @@ import tensorflow as tf -@tf.function(input_signature=[tf.TensorSpec(shape=[None], dtype=DTYPE), - tf.TensorSpec(shape=[None, None], dtype=DTYPEINT), - tf.TensorSpec(shape=[], dtype=DTYPEINT)]) +@tf.function( + input_signature=[ + tf.TensorSpec(shape=[None], dtype=DTYPE), + tf.TensorSpec(shape=[None, None], dtype=DTYPEINT), + tf.TensorSpec(shape=[], dtype=DTYPEINT), + ] +) def consume_array_into_indices(input_arr, indices, result_size): """ Accumulate the input tensor `input_arr` into an output tensor of @@ -38,6 +42,7 @@ def consume_array_into_indices(input_arr, indices, result_size): final_result = tf.reduce_sum(res_tmp, axis=1) return final_result + def py_consume_array_into_indices(input_arr, indices, result_size): """ Python interface wrapper for ``consume_array_into_indices``. @@ -46,8 +51,8 @@ def py_consume_array_into_indices(input_arr, indices, result_size): return consume_array_into_indices(float_me(input_arr), int_me(indices), int_me(result_size)) -def generate_condition_function(n_mask, condition='and'): - """ Generates a function that takes a number of masks +def generate_condition_function(n_mask, condition="and"): + """Generates a function that takes a number of masks and returns a combination of all n_masks for the given condition. It is possible to pass a list of allowed conditions, in that case @@ -76,7 +81,7 @@ def generate_condition_function(n_mask, condition='and'): full_mask= indices= - + Parameters ---------- @@ -90,16 +95,14 @@ def generate_condition_function(n_mask, condition='and'): `condition_to_idx`: function function(*masks) -> full_mask, true indices """ - allowed_conditions = { - 'and': tf.math.logical_and, - 'or' : tf.math.logical_or - } + allowed_conditions = {"and": tf.math.logical_and, "or": tf.math.logical_or} allo = list(allowed_conditions.keys()) # Check that the user is not asking for anything weird if n_mask < 2: - raise ValueError(f"At least two masks needed to generate a wrapper") - elif isinstance(condition, str): + raise ValueError("At least two masks needed to generate a wrapper") + + if isinstance(condition, str): if condition not in allowed_conditions: raise ValueError(f"Wrong condition {condition}, allowed values are {allo}") is_list = False @@ -112,7 +115,7 @@ def generate_condition_function(n_mask, condition='and'): is_list = True def py_condition(*masks): - """ Receives a list of conditions and returns a result mask + """Receives a list of conditions and returns a result mask and the list of indices in which the result mask is True Returns @@ -125,15 +128,15 @@ def py_condition(*masks): if is_list: res = masks[0] for i, cond in enumerate(condition): - res = allowed_conditions[cond](res, masks[i+1]) - elif condition == 'and': + res = allowed_conditions[cond](res, masks[i + 1]) + elif condition == "and": res = tf.math.reduce_all(masks, axis=0) - elif condition == 'or': + elif condition == "or": res = tf.math.reduce_any(masks, axis=0) indices = int_me(tf.where(res)) return res, indices - signature = n_mask*[tf.TensorSpec(shape=[None], dtype=tf.bool)] + signature = n_mask * [tf.TensorSpec(shape=[None], dtype=tf.bool)] condition_to_idx = tf.function(py_condition, input_signature=signature) - return condition_to_idx \ No newline at end of file + return condition_to_idx diff --git a/src/vegasflow/vflow.py b/src/vegasflow/vflow.py index 5737fe2..0162346 100644 --- a/src/vegasflow/vflow.py +++ b/src/vegasflow/vflow.py @@ -5,7 +5,6 @@ `vegas_wrapper` """ from vegasflow.configflow import DTYPE, DTYPEINT, fone, fzero, float_me, ione, int_me -import json import numpy as np import tensorflow as tf diff --git a/src/vegasflow/vflowplus.py b/src/vegasflow/vflowplus.py index 0b48f9f..a56bc64 100644 --- a/src/vegasflow/vflowplus.py +++ b/src/vegasflow/vflowplus.py @@ -22,7 +22,7 @@ BETA, MAX_NEVAL_HCUBE, ) -from vegasflow.monte_carlo import wrapper, sampler +from vegasflow.monte_carlo import wrapper, sampler, MonteCarloFlow from vegasflow.vflow import VegasFlow from vegasflow.utils import consume_array_into_indices @@ -146,8 +146,7 @@ def __init__(self, n_dim, n_events, train=True, adaptive=None, **kwargs): raise ValueError("Hypercubes are not equal to n_strat^n_dim") self.min_neval_hcube = int(neval_eff // len(hypercubes)) - if self.min_neval_hcube < 2: - self.min_neval_hcube = 2 + self.min_neval_hcube = max(self.min_neval_hcube, 2) self.n_ev = tf.fill([1, len(hypercubes)], self.min_neval_hcube) self.n_ev = tf.cast(tf.reshape(self.n_ev, [-1]), dtype=DTYPEINT) @@ -176,7 +175,8 @@ def _generate_random_array(self, n_events): def _generate_random_array_plus(self, n_events, n_ev): """Generate a random array for a given number of events divided in hypercubes""" - rnds, _, _ = super(VegasFlow, self)._generate_random_array(n_events) + # Needs to skip parent and go directly to the random array generation of MonteCarloFlow + rnds, _, _ = MonteCarloFlow._generate_random_array(self, n_events) # Get random numbers from hypercubes x, ind, w, segm = generate_samples_in_hypercubes( rnds, @@ -251,7 +251,8 @@ def _iteration_content(self): return res, sigma def run_event(self, tensorize_events=None, **kwargs): - """Tensorizes the number of events so they are not python or numpy primitives if self._adaptive=True""" + """Tensorizes the number of events + so they are not python or numpy primitives if self._adaptive=True""" return super().run_event(tensorize_events=self._adaptive, **kwargs) From d013f4fabc7154783e812defc802bff03f14d454 Mon Sep 17 00:00:00 2001 From: juacrumar Date: Tue, 17 Aug 2021 14:48:13 +0200 Subject: [PATCH 23/27] ups --- src/vegasflow/vflow.py | 1 + 1 file changed, 1 insertion(+) diff --git a/src/vegasflow/vflow.py b/src/vegasflow/vflow.py index 0162346..5737fe2 100644 --- a/src/vegasflow/vflow.py +++ b/src/vegasflow/vflow.py @@ -5,6 +5,7 @@ `vegas_wrapper` """ from vegasflow.configflow import DTYPE, DTYPEINT, fone, fzero, float_me, ione, int_me +import json import numpy as np import tensorflow as tf From 34def6fde73a6fa24b07d725cd4964bacb20ac76 Mon Sep 17 00:00:00 2001 From: juacrumar Date: Wed, 18 Aug 2021 10:47:48 +0200 Subject: [PATCH 24/27] remove duplicated code --- src/vegasflow/__init__.py | 1 + src/vegasflow/configflow.py | 2 +- src/vegasflow/plain.py | 5 +- src/vegasflow/vflow.py | 95 +++++++++++++++++++++++++++---------- src/vegasflow/vflowplus.py | 47 ++++-------------- 5 files changed, 84 insertions(+), 66 deletions(-) diff --git a/src/vegasflow/__init__.py b/src/vegasflow/__init__.py index e3aed64..9b0b057 100644 --- a/src/vegasflow/__init__.py +++ b/src/vegasflow/__init__.py @@ -1,6 +1,7 @@ """Monte Carlo integration with Tensorflow""" from vegasflow.configflow import int_me, float_me, run_eager + # Expose the main interfaces from vegasflow.vflow import VegasFlow, vegas_wrapper, vegas_sampler from vegasflow.plain import PlainFlow, plain_wrapper, plain_sampler diff --git a/src/vegasflow/configflow.py b/src/vegasflow/configflow.py index cc01659..4887d3f 100644 --- a/src/vegasflow/configflow.py +++ b/src/vegasflow/configflow.py @@ -9,7 +9,7 @@ # Some global parameters BINS_MAX = 50 ALPHA = 1.5 -BETA = 0.75 # Vegas+ +BETA = 0.75 # Vegas+ TECH_CUT = 1e-8 # Set up the logistics of the integration diff --git a/src/vegasflow/plain.py b/src/vegasflow/plain.py index 32d5433..86df754 100644 --- a/src/vegasflow/plain.py +++ b/src/vegasflow/plain.py @@ -38,9 +38,10 @@ def _run_iteration(self): def plain_wrapper(*args, **kwargs): - """ Wrapper around PlainFlow """ + """Wrapper around PlainFlow""" return wrapper(PlainFlow, *args, **kwargs) + def plain_sampler(*args, **kwargs): - """ Wrapper sampler around PlainFlow """ + """Wrapper sampler around PlainFlow""" return sampler(PlainFlow, *args, **kwargs) diff --git a/src/vegasflow/vflow.py b/src/vegasflow/vflow.py index 5737fe2..5e9656c 100644 --- a/src/vegasflow/vflow.py +++ b/src/vegasflow/vflow.py @@ -4,12 +4,21 @@ The main interfaces of this class are the class `VegasFlow` and the `vegas_wrapper` """ -from vegasflow.configflow import DTYPE, DTYPEINT, fone, fzero, float_me, ione, int_me import json import numpy as np import tensorflow as tf -from vegasflow.configflow import BINS_MAX, ALPHA +from vegasflow.configflow import ( + DTYPE, + DTYPEINT, + fone, + fzero, + float_me, + ione, + int_me, + BINS_MAX, + ALPHA, +) from vegasflow.monte_carlo import MonteCarloFlow, wrapper, sampler from vegasflow.utils import consume_array_into_indices @@ -19,6 +28,59 @@ FBINS = float_me(BINS_MAX) + +@tf.function( + input_signature=[ + tf.TensorSpec(shape=[None, None], dtype=DTYPE), + tf.TensorSpec(shape=[None, BINS_MAX + 1], dtype=DTYPE), + ] +) +def importance_sampling_digest(xn, divisions): + """Importance sampling algorithm: + receives a random array (number of dimensions, number of dim) + containing information about from which bins in the + grid (n_dims, BINS_MAX+1) the random points have to be sampled + + This algorithm is shared between the simplest form of Vegas + (VegasFlow: only importance sampling) + and Vegas+ (VegasFlowPlus: importance and stratified sampling) + and so it has been lifted to its own function + + Parameters: + ---------- + xn: float tensor (n_dim, n_events) + which bins to sample from + divisions: float tensor (n_dims, BINS_MAX+1) + grid of divisions for the importance sampling algorithm + + Returns + ------- + ind_i: integer tensor (n_events, n_dim) + index in the divisions grid from which the points should be sampled + x: float tensor (n_events, n_dim) + random values sampled in the divisions grid + xdelta: float tensor (n_events,) + weight of the random points + """ + ind_i = int_me(xn) + # Get the value of the left and right sides of the bins + ind_f = ind_i + ione + x_ini = tf.gather(divisions, ind_i, batch_dims=1) + x_fin = tf.gather(divisions, ind_f, batch_dims=1) + # Compute the width of the bins + xdelta = x_fin - x_ini + # Take the decimal part of bin (i.e., how deep within the bin) + aux_rand = xn - tf.math.floor(xn) + x = x_ini + xdelta * aux_rand + # Compute the weight of the points + weights = tf.reduce_prod(xdelta * FBINS, axis=0) + + # Tranpose the output to be what the external functions expect + x = tf.transpose(x) + ind_i = tf.transpose(ind_i) + return ind_i, x, weights + + # Auxiliary functions for Vegas @tf.function( input_signature=[ @@ -34,7 +96,7 @@ def _generate_random_array(rnds, divisions): ---------- rnds: array shaped (None, n_dim) Random numbers used as an input for Vegas - divisions: array shaped (n_dim, BINS_MAX) + divisions: array shaped (n_dim, BINS_MAX+1) vegas grid Returns @@ -52,33 +114,14 @@ def _generate_random_array(rnds, divisions): # Get the index of the division we are interested in xn = FBINS * (fone - tf.transpose(rnds)) - @tf.function( - input_signature=[ - tf.TensorSpec(shape=[None, None], dtype=DTYPE), - ] - ) - def digest(xn): - ind_i = tf.cast(xn, DTYPEINT) - # Get the value of the left and right sides of the bins - ind_f = ind_i + ione - x_ini = tf.gather(divisions, ind_i, batch_dims=1) - x_fin = tf.gather(divisions, ind_f, batch_dims=1) - # Compute the width of the bins - xdelta = x_fin - x_ini - return ind_i, x_ini, xdelta - - ind_i, x_ini, xdelta = digest(xn) # Compute the random number between 0 and 1 - aux_rand = xn - tf.math.floor(xn) - x = x_ini + xdelta * aux_rand + # and the index of the bin where it has been sampled from + ind_xn, x, weights = importance_sampling_digest(xn, divisions) # Compute the random number between the limits + # commented, for now only from 0 to 1 # x = reg_i + rand_x * (reg_f - reg_i) - # and the weight - weights = tf.reduce_prod(xdelta * FBINS, axis=0) - x_t = tf.transpose(x) - int_xn = tf.transpose(ind_i) - return x_t, int_xn, weights + return x, ind_xn, weights @tf.function( diff --git a/src/vegasflow/vflowplus.py b/src/vegasflow/vflowplus.py index a56bc64..b09d41f 100644 --- a/src/vegasflow/vflowplus.py +++ b/src/vegasflow/vflowplus.py @@ -16,14 +16,13 @@ fone, fzero, float_me, - ione, int_me, BINS_MAX, BETA, MAX_NEVAL_HCUBE, ) from vegasflow.monte_carlo import wrapper, sampler, MonteCarloFlow -from vegasflow.vflow import VegasFlow +from vegasflow.vflow import VegasFlow, importance_sampling_digest from vegasflow.utils import consume_array_into_indices import logging @@ -33,26 +32,6 @@ FBINS = float_me(BINS_MAX) -@tf.function(input_signature=3 * [tf.TensorSpec(shape=[None, None], dtype=DTYPE)]) -def _compute_x(x_ini, xn, xdelta): - """Helper function for ``generate_samples_in_hypercubes``""" - aux_rand = xn - tf.math.floor(xn) - return x_ini + xdelta * aux_rand - - -# same as vegasflow generate_random_array -@tf.function(input_signature=2 * [tf.TensorSpec(shape=[None, None], dtype=DTYPE)]) -def _digest(xn, divisions): - ind_i = tf.cast(xn, DTYPEINT) - # Get the value of the left and right sides of the bins - ind_f = ind_i + ione - x_ini = tf.gather(divisions, ind_i, batch_dims=1) - x_fin = tf.gather(divisions, ind_f, batch_dims=1) - # Compute the width of the bins - xdelta = x_fin - x_ini - return ind_i, x_ini, xdelta - - @tf.function( input_signature=[ tf.TensorSpec(shape=[None, None], dtype=DTYPE), @@ -62,7 +41,7 @@ def _digest(xn, divisions): tf.TensorSpec(shape=[None, None], dtype=DTYPE), ] ) -def generate_samples_in_hypercubes(rnds_t, n_strat, n_ev, hypercubes, divisions): +def generate_samples_in_hypercubes(rnds, n_strat, n_ev, hypercubes, divisions): """Receives an array of random numbers 0 and 1 and distribute them in each hypercube according to the number of samples in each hypercube specified by n_ev @@ -82,27 +61,21 @@ def generate_samples_in_hypercubes(rnds_t, n_strat, n_ev, hypercubes, divisions) `ind`: division index in which each (n_dim) set of random numbers fall `segm` : segmentantion for later computations """ - rnds = tf.transpose(rnds_t) + # Use the event-per-hypercube information to fix each random event to a hypercub indices = tf.repeat(tf.range(tf.shape(hypercubes, out_type=DTYPEINT)[0]), n_ev) points = float_me(tf.gather(hypercubes, indices)) n_evs = float_me(tf.gather(n_ev, indices)) - xn = tf.transpose((points + tf.transpose(rnds)) * FBINS / float_me(n_strat)) - segm = indices - ind_i, x_ini, xdelta = _digest(xn, divisions) - # Compute the random number between 0 and 1 - # This is the heavy part of the calc + # Compute in which division of the importance_sampling grid the points fall + xn = tf.transpose(points + rnds) * FBINS / float_me(n_strat) + + ind_xn, x, weights = importance_sampling_digest(xn, divisions) - x = _compute_x(x_ini, xn, xdelta) - # Compute the random number between the limits - # x = reg_i + rand_x * (reg_f - reg_i) - # and the weight - weights = tf.reduce_prod(xdelta * FBINS, axis=0) + # Reweight taking into account the number of events per hypercub final_weights = weights / n_evs - x_t = tf.transpose(x) - int_xn = tf.transpose(ind_i) - return x_t, int_xn, final_weights, segm + segm = indices + return x, ind_xn, final_weights, segm class VegasFlowPlus(VegasFlow): From b1a47138ac7870fa6e13fc05711b3a7d663414db Mon Sep 17 00:00:00 2001 From: juacrumar Date: Wed, 18 Aug 2021 11:13:26 +0200 Subject: [PATCH 25/27] remove more duplicates --- src/vegasflow/vflow.py | 28 ++++++++++++++++++++-------- src/vegasflow/vflowplus.py | 26 ++++++++++---------------- 2 files changed, 30 insertions(+), 24 deletions(-) diff --git a/src/vegasflow/vflow.py b/src/vegasflow/vflow.py index 5e9656c..837647f 100644 --- a/src/vegasflow/vflow.py +++ b/src/vegasflow/vflow.py @@ -341,6 +341,25 @@ def _generate_random_array(self, n_events): x, ind, w = _generate_random_array(rnds, self.divisions) return x, ind, w * xjac + def _importance_sampling_array_filling(self, results2, indices): + """Receives an array of results squared for every event + and an array of indices describing in which bin each result fall. + Fills a array with the total result in each bin to be used by + the importance sampling algorithm + """ + if not self.train: + return [] + + arr_res2 = [] + # If the training is active, save the result of the integral sq + for j in range(self.n_dim): + arr_res2.append( + consume_array_into_indices( + results2, indices[:, j : j + 1], int_me(self.grid_bins - 1) + ) + ) + return tf.reshape(arr_res2, (self.n_dim, -1)) + def _run_event(self, integrand, ncalls=None): """Runs one step of Vegas. @@ -371,14 +390,7 @@ def _run_event(self, integrand, ncalls=None): res = tf.reduce_sum(tmp) res2 = tf.reduce_sum(tmp2) - arr_res2 = [] - if self.train: - # If the training is active, save the result of the integral sq - for j in range(self.n_dim): - arr_res2.append( - consume_array_into_indices(tmp2, ind[:, j : j + 1], int_me(self.grid_bins - 1)) - ) - arr_res2 = tf.reshape(arr_res2, (self.n_dim, -1)) + arr_res2 = self._importance_sampling_array_filling(tmp2, ind) return res, res2, arr_res2 diff --git a/src/vegasflow/vflowplus.py b/src/vegasflow/vflowplus.py index b09d41f..7468b56 100644 --- a/src/vegasflow/vflowplus.py +++ b/src/vegasflow/vflowplus.py @@ -122,7 +122,7 @@ def __init__(self, n_dim, n_events, train=True, adaptive=None, **kwargs): self.min_neval_hcube = max(self.min_neval_hcube, 2) self.n_ev = tf.fill([1, len(hypercubes)], self.min_neval_hcube) - self.n_ev = tf.cast(tf.reshape(self.n_ev, [-1]), dtype=DTYPEINT) + self.n_ev = int_me(tf.reshape(self.n_ev, [-1])) self._n_events = int(tf.reduce_sum(self.n_ev)) self.my_xjac = float_me(1 / len(hypercubes)) @@ -174,12 +174,11 @@ def _run_event(self, integrand, ncalls=None, n_ev=None): Returns ------- - `res`: sum of the result of the integrand for all events - `res2`: sum of the result squared of the integrand for all events + `res`: sum of the result of the integrand for all events per segement + `res2`: sum of the result squared of the integrand for all events per segment `arr_res2`: result of the integrand squared per dimension and grid bin """ # NOTE: needs to receive both ncalls and n_ev - x, ind, xjac, segm = self._generate_random_array_plus(ncalls, n_ev) # compute integrand @@ -190,22 +189,17 @@ def _run_event(self, integrand, ncalls=None, n_ev=None): ress = tf.math.segment_sum(tmp, segm) ress2 = tf.math.segment_sum(tmp2, segm) - Fn_ev = tf.cast(n_ev, DTYPE) - arr_var = ress2 * Fn_ev - tf.square(ress) - - arr_res2 = [] - if self.train: - # If the training is active, save the result of the integral sq - for j in range(self.n_dim): - arr_res2.append( - consume_array_into_indices(tmp2, ind[:, j : j + 1], int_me(self.grid_bins - 1)) - ) - arr_res2 = tf.reshape(arr_res2, (self.n_dim, -1)) + fn_ev = float_me(n_ev) + arr_var = ress2 * fn_ev - tf.square(ress) + arr_res2 = self._importance_sampling_array_filling(tmp2, ind) return ress, arr_var, arr_res2 def _iteration_content(self): - """Steps to follow per iteration""" + """Steps to follow per iteration + Differently from importance-sampling Vegas, the result of the integration + is a result _per segment_ and thus the total result needs to be computed at this point + """ ress, arr_var, arr_res2 = self.run_event(n_ev=self.n_ev) # Compute the rror From 516610db30cf3ef667125891bf1c3090e04d2710 Mon Sep 17 00:00:00 2001 From: juacrumar Date: Wed, 25 Aug 2021 10:33:48 +0200 Subject: [PATCH 26/27] bugfix --- src/vegasflow/vflowplus.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/vegasflow/vflowplus.py b/src/vegasflow/vflowplus.py index 7468b56..ee81a63 100644 --- a/src/vegasflow/vflowplus.py +++ b/src/vegasflow/vflowplus.py @@ -83,7 +83,7 @@ class VegasFlowPlus(VegasFlow): Implementation of the VEGAS+ algorithm """ - def __init__(self, n_dim, n_events, train=True, adaptive=None, **kwargs): + def __init__(self, n_dim, n_events, train=True, adaptive=False, **kwargs): _ = kwargs.setdefault("events_limit", n_events) super().__init__(n_dim, n_events, train, **kwargs) @@ -91,7 +91,7 @@ def __init__(self, n_dim, n_events, train=True, adaptive=None, **kwargs): self._init_calls = n_events # Don't use adaptive if the number of dimension is too big - if n_dim > 13 and adaptive is None: + if n_dim > 13 and adaptive: self._adaptive = False logger.warning("Disabling adaptive mode from VegasFlowPlus, too many dimensions!") else: From 1e1d4e3fe51209c469110e2793de79321f7e483a Mon Sep 17 00:00:00 2001 From: juacrumar Date: Thu, 26 Aug 2021 15:54:30 +0200 Subject: [PATCH 27/27] add warning when changing event limit --- src/vegasflow/vflowplus.py | 15 ++++++++++++--- 1 file changed, 12 insertions(+), 3 deletions(-) diff --git a/src/vegasflow/vflowplus.py b/src/vegasflow/vflowplus.py index ee81a63..adc017e 100644 --- a/src/vegasflow/vflowplus.py +++ b/src/vegasflow/vflowplus.py @@ -83,9 +83,18 @@ class VegasFlowPlus(VegasFlow): Implementation of the VEGAS+ algorithm """ - def __init__(self, n_dim, n_events, train=True, adaptive=False, **kwargs): - _ = kwargs.setdefault("events_limit", n_events) - super().__init__(n_dim, n_events, train, **kwargs) + def __init__(self, n_dim, n_events, train=True, adaptive=False, events_limit=None, **kwargs): + # https://github.com/N3PDF/vegasflow/issues/78 + if events_limit is None: + logger.info("Events per device limit set to %d", n_events) + events_limit = n_events + elif events_limit < n_events: + logger.warning("VegasFlowPlus needs to hold all events in memory at once, " + "setting the `events_limit` to be equal to `n_events=%d`", n_events) + events_limit = n_events + super().__init__(n_dim, n_events, train, events_limit=events_limit, **kwargs) + + # Save the initial number of events self._init_calls = n_events