From 62b89e4137273a4b9d44fe4f2def52d90be770ef Mon Sep 17 00:00:00 2001 From: juacrumar Date: Fri, 5 Feb 2021 16:25:14 +0100 Subject: [PATCH 01/19] initial rtbm commit --- examples/simgauss_tf.py | 12 +++- examples/singletop_lo_tf.py | 8 +++ src/vegasflow/monte_carlo.py | 60 +++++++++------- src/vegasflow/rtbm.py | 133 +++++++++++++++++++++++++++++++++++ src/vegasflow/vflow.py | 6 +- 5 files changed, 189 insertions(+), 30 deletions(-) create mode 100644 src/vegasflow/rtbm.py diff --git a/examples/simgauss_tf.py b/examples/simgauss_tf.py index 3605a4e..a17261b 100644 --- a/examples/simgauss_tf.py +++ b/examples/simgauss_tf.py @@ -10,10 +10,11 @@ import tensorflow as tf from vegasflow.vflow import vegas_wrapper from vegasflow.plain import plain_wrapper +from vegasflow.rtbm import rtbm_wrapper # MC integration setup -dim = 4 +dim = 3 ncalls = np.int32(1e5) n_iter = 5 @@ -34,9 +35,16 @@ def symgauss(xarr, n_dim=None, **kwargs): if __name__ == "__main__": """Testing several different integrations""" + print(f"RTBM MC, ncalls={ncalls}:") + start = time.time() + ncalls = ncalls + r = rtbm_wrapper(symgauss, dim, n_iter, ncalls) + end = time.time() + print(f"RTBM took: time (s): {end-start}") + print(f"VEGAS MC, ncalls={ncalls}:") start = time.time() - ncalls = 10*ncalls + ncalls = ncalls r = vegas_wrapper(symgauss, dim, n_iter, ncalls) end = time.time() print(f"Vegas took: time (s): {end-start}") diff --git a/examples/singletop_lo_tf.py b/examples/singletop_lo_tf.py index 7f687d4..0ed90af 100644 --- a/examples/singletop_lo_tf.py +++ b/examples/singletop_lo_tf.py @@ -8,6 +8,7 @@ from vegasflow.configflow import DTYPE import tensorflow as tf from vegasflow.vflow import vegas_wrapper +from vegasflow.rtbm import rtbm_wrapper # MC integration setup dim = 3 @@ -272,6 +273,13 @@ def singletop(xarr, n_dim=None, **kwargs): if __name__ == "__main__": """Testing a basic integration""" + + print(f"RTBM MC, ncalls={ncalls}:") + start = time.time() + r = rtbm_wrapper(singletop, dim, n_iter, ncalls) + end = time.time() + print(f"RTBM took: time (s): {end-start}") + print(f"VEGAS MC, ncalls={ncalls}:") start = time.time() r = vegas_wrapper(singletop, dim, n_iter, ncalls) diff --git a/src/vegasflow/monte_carlo.py b/src/vegasflow/monte_carlo.py index c50837b..d26c88d 100644 --- a/src/vegasflow/monte_carlo.py +++ b/src/vegasflow/monte_carlo.py @@ -57,27 +57,7 @@ def print_iteration(it, res, error, extra="", threshold=0.1): logger.info(f"Result for iteration {it}: {res:.4f} +/- {error:.4f}" + extra) -def _accumulate(accumulators): - """Accumulate all the quantities in accumulators - The default accumulation is implemented for tensorflow tensors - as a sum of all partial results. - Parameters - ---------- - `accumulators`: list of tensorflow tensors - - Returns - ------- - `results`: `sum` for each element of the accumulators - - Function not compiled - """ - results = [] - len_acc = len(accumulators[0]) - for i in range(len_acc): - total = tf.reduce_sum([acc[i] for acc in accumulators], axis=0) - results.append(total) - return results class MonteCarloFlow(ABC): @@ -99,13 +79,15 @@ class MonteCarloFlow(ABC): def __init__( self, - n_dim, - n_events, + n_dim=None, + n_events=None, events_limit=MAX_EVENTS_LIMIT, list_devices=DEFAULT_ACTIVE_DEVICES, verbose=True, simplify_signature=False, ): + if not ( isinstance(n_dim, int) and n_dim > 0 ): + raise ValueError(f"Can't integrate an integrand with non-positive or non-integer dimensionality, received: {n_dim}") # Save some parameters self.n_dim = n_dim self.integrand = None @@ -132,6 +114,11 @@ def __init__( self.pool = joblib.Parallel(n_jobs=len(devices), prefer="threads") else: self.devices = None + self._train = False + + @property + def train(self): + return self._train @property def events_per_run(self): @@ -162,6 +149,29 @@ def history(self): """ return self._history + @staticmethod + def _accumulate(accumulators): + """Accumulate all the quantities in accumulators + The default accumulation is implemented for tensorflow tensors + as a sum of all partial results. + + Parameters + ---------- + `accumulators`: list of tensorflow tensors + + Returns + ------- + `results`: `sum` for each element of the accumulators + + Function not compiled + """ + results = [] + len_acc = len(accumulators[0]) + for i in range(len_acc): + total = tf.reduce_sum([acc[i] for acc in accumulators], axis=0) + results.append(total) + return results + def generate_random_array(self, n_events): """Generate a 2D array of (n_events, n_dim) points Parameters @@ -314,7 +324,7 @@ def run_event(self, **kwargs): # Generate the client to control the distribution using the cluster variable client = Client(cluster) accumulators_future = client.map(self.device_run, events_to_do, percentages) - result_future = client.submit(_accumulate, accumulators_future) + result_future = client.submit(self._accumulate, accumulators_future) result = result_future.result() # Liberate the client client.close() @@ -325,7 +335,7 @@ def run_event(self, **kwargs): for ncalls, pc in zip(events_to_do, percentages): res = self.device_run(ncalls, sent_pc=pc, **kwargs) accumulators.append(res) - return _accumulate(accumulators) + return self._accumulate(accumulators) def compile(self, integrand, compilable=True): """Receives an integrand, prepares it for integration @@ -482,6 +492,6 @@ def wrapper(integrator_class, integrand, n_dim, n_iter, total_n_events, compilab `final_result`: integral value `sigma`: monte carlo error """ - mc_instance = integrator_class(n_dim, total_n_events) + mc_instance = integrator_class(n_dim = n_dim, n_events=total_n_events) mc_instance.compile(integrand, compilable=compilable) return mc_instance.run_integration(n_iter) diff --git a/src/vegasflow/rtbm.py b/src/vegasflow/rtbm.py new file mode 100644 index 0000000..356f8f4 --- /dev/null +++ b/src/vegasflow/rtbm.py @@ -0,0 +1,133 @@ +""" + Plain implementation of the plainest possible MonteCarlo +""" + +from vegasflow.configflow import DTYPE, fone, fzero, float_me +from vegasflow.monte_carlo import MonteCarloFlow, wrapper +import numpy as np +import tensorflow as tf + +from theta.rtbm import RTBM +from theta.minimizer import CMA +from theta import costfunctions + +import logging + +logger = logging.getLogger(__name__) + + +class RTBMFlow(MonteCarloFlow): + """ + RTBM based Monte Carlo integrator + """ + + def __init__(self, n_hidden=2, rtbm = None, train = True, *args, **kwargs): + super().__init__(*args, **kwargs) + self._train = train + self._first_run = True + if rtbm is None: + logger.info("Generating a RTBM with %d visible nodes and %d hidden" % (self.n_dim, n_hidden)) + self._rtbm = RTBM(self.n_dim, n_hidden, init_max_param_bound=50, random_bound=1, diagonal_T=True) + else: + # Check whether it is a valid rtbm model + if not hasattr(rtbm, "make_sample"): + raise TypeError(f"{rtbm} is not a valid boltzman machine") + self._rtbm = rtbm + + def freeze(self): + self.train=False + + def unfreeze(self): + self.train=True + + def compile(self, integrand, compilable=False, **kwargs): + if compilable: + logger.warning("RTBMFlow is still WIP and not compilable") + self._trainer = CMA(ncores=1, verbose=True) + super().compile(integrand, compilable=False, **kwargs) + + def generate_random_array(self, n_events): + if self._first_run: + return super().generate_random_array(n_events) + xrand, _ = self._rtbm.make_sample(n_events) + xjac = 1.0/self._rtbm(xrand.T)/n_events + return float_me(xrand), None, xjac + + def _train_machine(self, x, yraw): + options = { + "popsize": 50, + "tolfun": 1e-1, + "maxiter" : 3, + } + # The output is a probability, therefore: + y = yraw / tf.reduce_sum(yraw) + # We don't want to train tf variables because it would be slow here... + xnp = x.numpy().T + ynp = y.numpy() + self._trainer.train(costfunctions.kullbackLeibler, self._rtbm, xnp, ynp, **options) + self._first_run = False + + @staticmethod + def _accumulate(accumulators): + """ For the RTBM accumulation strategy we need to keep track + of all results and who produced it""" + # In the accumulators we should receive a number of items with + # (res, xrands) which have shape ( (n_events,), (n_events, n_dim) ) + all_res = [] + all_rnd = [] + for res, rnds, in accumulators: + all_res.append(res) + all_rnd.append(rnds) + return tf.concat(all_res, axis=0), tf.concat(all_rnd, axis=0) + + def _run_event(self, integrand, ncalls=None): + if ncalls is None: + n_events = self.n_events + else: + n_events = ncalls + + # Generate all random number for this iteration + rnds, _, xjac = self.generate_random_array(n_events) + # Compute the integrand + res = integrand(rnds, n_dim=self.n_dim, weight=xjac) * xjac + + # Clean up the array from numbers outside the 0-1 range + if not self._first_run: + # Accept for now only random number between 0 and 1 + condition = tf.reduce_all(rnds>=0.0, axis=1) & tf.reduce_all(rnds<=1.0, axis=1) + res = tf.where(condition, res, fzero)[0] + if np.count_nonzero(res) != n_events: + logger.warning(f"Found only {np.count_nonzero(res)} of {n_events} valid values\n") + + return res, rnds + + def _run_iteration(self): + all_res, rnds = self.run_event() + if self.train: + self._train_machine(rnds, all_res) + + res = tf.reduce_sum(all_res) + all_res2 = all_res**2 + res2 = tf.reduce_sum(all_res2)*self.n_events + + # Compute the error + err_tmp2 = (res2 - tf.square(res)) / (self.n_events - fone) + sigma = tf.sqrt(tf.maximum(err_tmp2, fzero)) + return res, sigma + +def rtbm_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(RTBMFlow, integrand, n_dim, n_iter, total_n_events, **kwargs) diff --git a/src/vegasflow/vflow.py b/src/vegasflow/vflow.py index 11d3b5f..78ac772 100644 --- a/src/vegasflow/vflow.py +++ b/src/vegasflow/vflow.py @@ -168,12 +168,12 @@ class VegasFlow(MonteCarloFlow): Implementation of the important sampling algorithm from Vegas """ - def __init__(self, n_dim, n_events, train=True, **kwargs): - super().__init__(n_dim, n_events, **kwargs) + def __init__(self, n_dim=None, n_events=None, train=True, **kwargs): + super().__init__(n_dim=n_dim, n_events=n_events, **kwargs) # If training is True, the grid will be changed after every iteration # otherwise it will be frozen - self.train = train + self._train = train self.iteration_content = None self.compile_args = None From 5ab22d5274964df0c2bf4561f1f41f211a1404d2 Mon Sep 17 00:00:00 2001 From: juacrumar Date: Wed, 10 Feb 2021 12:55:19 +0100 Subject: [PATCH 02/19] some changes --- examples/simgauss_tf.py | 14 +++++++++----- src/vegasflow/rtbm.py | 26 +++++++++++++++++++------- 2 files changed, 28 insertions(+), 12 deletions(-) diff --git a/examples/simgauss_tf.py b/examples/simgauss_tf.py index a17261b..52ef294 100644 --- a/examples/simgauss_tf.py +++ b/examples/simgauss_tf.py @@ -4,19 +4,23 @@ Basic example using the vegas_wrapper helper """ -from vegasflow.configflow import DTYPE +import tensorflow as tf +tf.config.threading.set_inter_op_parallelism_threads(1) +tf.config.threading.set_intra_op_parallelism_threads(1) +from vegasflow.configflow import DTYPE, run_eager +run_eager(True) + import time import numpy as np -import tensorflow as tf from vegasflow.vflow import vegas_wrapper from vegasflow.plain import plain_wrapper from vegasflow.rtbm import rtbm_wrapper # MC integration setup -dim = 3 -ncalls = np.int32(1e5) -n_iter = 5 +dim = 4 +ncalls = np.int32(1e6) +n_iter = 3 @tf.function diff --git a/src/vegasflow/rtbm.py b/src/vegasflow/rtbm.py index 356f8f4..fbdb8bb 100644 --- a/src/vegasflow/rtbm.py +++ b/src/vegasflow/rtbm.py @@ -27,7 +27,7 @@ def __init__(self, n_hidden=2, rtbm = None, train = True, *args, **kwargs): self._first_run = True if rtbm is None: logger.info("Generating a RTBM with %d visible nodes and %d hidden" % (self.n_dim, n_hidden)) - self._rtbm = RTBM(self.n_dim, n_hidden, init_max_param_bound=50, random_bound=1, diagonal_T=True) + self._rtbm = RTBM(self.n_dim, n_hidden, init_max_param_bound=50, random_bound=2, diagonal_T=False, positive_T=True, positive_Q=True) else: # Check whether it is a valid rtbm model if not hasattr(rtbm, "make_sample"): @@ -43,21 +43,33 @@ def unfreeze(self): def compile(self, integrand, compilable=False, **kwargs): if compilable: logger.warning("RTBMFlow is still WIP and not compilable") - self._trainer = CMA(ncores=1, verbose=True) + self._trainer = CMA(parallel=False, ncores=1, verbose=True) super().compile(integrand, compilable=False, **kwargs) def generate_random_array(self, n_events): if self._first_run: return super().generate_random_array(n_events) xrand, _ = self._rtbm.make_sample(n_events) - xjac = 1.0/self._rtbm(xrand.T)/n_events - return float_me(xrand), None, xjac + xjac_raw = 1.0/self._rtbm(xrand.T)/n_events + # Now re-scale the values to the 0-1 range per dimension. + # NOTE: this is only valid while the points in_device are equal to the points being utilized + # so when the RTBM can run in the GPU it will need to be modified! + # mainly to carry the first limits to the rest of the calculation + epsilon = np.abs(np.max(xrand)/4.0) + max_per_d = np.max(xrand, axis=0) + epsilon + min_per_d = np.min(xrand, axis=0) - epsilon + delta = max_per_d - min_per_d + new_rand = (xrand - min_per_d)/delta + xjac = xjac_raw/np.prod(delta) + print(delta) + print("") + return float_me(new_rand), None, xjac def _train_machine(self, x, yraw): options = { - "popsize": 50, - "tolfun": 1e-1, - "maxiter" : 3, + "popsize": 75, + "tolfun": 1e-4, + "maxiter" : 4, } # The output is a probability, therefore: y = yraw / tf.reduce_sum(yraw) From cc74726ea4cb564a1454182659047a17bb239167 Mon Sep 17 00:00:00 2001 From: juacrumar Date: Wed, 10 Feb 2021 13:01:51 +0100 Subject: [PATCH 03/19] fix --- examples/simgauss_tf.py | 2 +- src/vegasflow/monte_carlo.py | 4 +++ src/vegasflow/rtbm.py | 54 ++++++++++++++++++++++-------------- src/vegasflow/vflow.py | 4 +-- 4 files changed, 40 insertions(+), 24 deletions(-) diff --git a/examples/simgauss_tf.py b/examples/simgauss_tf.py index 52ef294..b2e680c 100644 --- a/examples/simgauss_tf.py +++ b/examples/simgauss_tf.py @@ -18,7 +18,7 @@ # MC integration setup -dim = 4 +dim = 2 ncalls = np.int32(1e6) n_iter = 3 diff --git a/src/vegasflow/monte_carlo.py b/src/vegasflow/monte_carlo.py index d26c88d..10bf832 100644 --- a/src/vegasflow/monte_carlo.py +++ b/src/vegasflow/monte_carlo.py @@ -120,6 +120,10 @@ def __init__( def train(self): return self._train + @train.setter + def train(self, new_val): + self._train = new_val + @property def events_per_run(self): """Number of events to run in a single step. diff --git a/src/vegasflow/rtbm.py b/src/vegasflow/rtbm.py index fbdb8bb..19e46f8 100644 --- a/src/vegasflow/rtbm.py +++ b/src/vegasflow/rtbm.py @@ -21,13 +21,23 @@ class RTBMFlow(MonteCarloFlow): RTBM based Monte Carlo integrator """ - def __init__(self, n_hidden=2, rtbm = None, train = True, *args, **kwargs): + def __init__(self, n_hidden=2, rtbm=None, train=True, *args, **kwargs): super().__init__(*args, **kwargs) self._train = train self._first_run = True if rtbm is None: - logger.info("Generating a RTBM with %d visible nodes and %d hidden" % (self.n_dim, n_hidden)) - self._rtbm = RTBM(self.n_dim, n_hidden, init_max_param_bound=50, random_bound=2, diagonal_T=False, positive_T=True, positive_Q=True) + logger.info( + "Generating a RTBM with %d visible nodes and %d hidden" % (self.n_dim, n_hidden) + ) + self._rtbm = RTBM( + self.n_dim, + n_hidden, + minimization_bound=50, + gaussian_init=True, + diagonal_T=False, + positive_T=True, + positive_Q=True, + ) else: # Check whether it is a valid rtbm model if not hasattr(rtbm, "make_sample"): @@ -35,10 +45,10 @@ def __init__(self, n_hidden=2, rtbm = None, train = True, *args, **kwargs): self._rtbm = rtbm def freeze(self): - self.train=False - + self.train = False + def unfreeze(self): - self.train=True + self.train = True def compile(self, integrand, compilable=False, **kwargs): if compilable: @@ -50,27 +60,25 @@ def generate_random_array(self, n_events): if self._first_run: return super().generate_random_array(n_events) xrand, _ = self._rtbm.make_sample(n_events) - xjac_raw = 1.0/self._rtbm(xrand.T)/n_events + xjac_raw = 1.0 / self._rtbm(xrand.T) / n_events # Now re-scale the values to the 0-1 range per dimension. # NOTE: this is only valid while the points in_device are equal to the points being utilized # so when the RTBM can run in the GPU it will need to be modified! # mainly to carry the first limits to the rest of the calculation - epsilon = np.abs(np.max(xrand)/4.0) + epsilon = np.abs(np.max(xrand) / 4.0) max_per_d = np.max(xrand, axis=0) + epsilon min_per_d = np.min(xrand, axis=0) - epsilon delta = max_per_d - min_per_d - new_rand = (xrand - min_per_d)/delta - xjac = xjac_raw/np.prod(delta) - print(delta) - print("") + new_rand = (xrand - min_per_d) / delta + xjac = xjac_raw / np.prod(delta) return float_me(new_rand), None, xjac def _train_machine(self, x, yraw): options = { - "popsize": 75, - "tolfun": 1e-4, - "maxiter" : 4, - } + "popsize": 75, + "tolfun": 1e-4, + "maxiter": 4, + } # The output is a probability, therefore: y = yraw / tf.reduce_sum(yraw) # We don't want to train tf variables because it would be slow here... @@ -81,13 +89,16 @@ def _train_machine(self, x, yraw): @staticmethod def _accumulate(accumulators): - """ For the RTBM accumulation strategy we need to keep track + """For the RTBM accumulation strategy we need to keep track of all results and who produced it""" # In the accumulators we should receive a number of items with # (res, xrands) which have shape ( (n_events,), (n_events, n_dim) ) all_res = [] all_rnd = [] - for res, rnds, in accumulators: + for ( + res, + rnds, + ) in accumulators: all_res.append(res) all_rnd.append(rnds) return tf.concat(all_res, axis=0), tf.concat(all_rnd, axis=0) @@ -106,7 +117,7 @@ def _run_event(self, integrand, ncalls=None): # Clean up the array from numbers outside the 0-1 range if not self._first_run: # Accept for now only random number between 0 and 1 - condition = tf.reduce_all(rnds>=0.0, axis=1) & tf.reduce_all(rnds<=1.0, axis=1) + condition = tf.reduce_all(rnds >= 0.0, axis=1) & tf.reduce_all(rnds <= 1.0, axis=1) res = tf.where(condition, res, fzero)[0] if np.count_nonzero(res) != n_events: logger.warning(f"Found only {np.count_nonzero(res)} of {n_events} valid values\n") @@ -119,14 +130,15 @@ def _run_iteration(self): self._train_machine(rnds, all_res) res = tf.reduce_sum(all_res) - all_res2 = all_res**2 - res2 = tf.reduce_sum(all_res2)*self.n_events + all_res2 = all_res ** 2 + res2 = tf.reduce_sum(all_res2) * self.n_events # Compute the error err_tmp2 = (res2 - tf.square(res)) / (self.n_events - fone) sigma = tf.sqrt(tf.maximum(err_tmp2, fzero)) return res, sigma + def rtbm_wrapper(integrand, n_dim, n_iter, total_n_events, **kwargs): """Convenience wrapper diff --git a/src/vegasflow/vflow.py b/src/vegasflow/vflow.py index 78ac772..648c2c7 100644 --- a/src/vegasflow/vflow.py +++ b/src/vegasflow/vflow.py @@ -185,12 +185,12 @@ def __init__(self, n_dim=None, n_events=None, train=True, **kwargs): def freeze_grid(self): """ Stops the grid from refining any more """ - self.train = False + self._train = False self.recompile() def unfreeze_grid(self): """ Enable the refining of the grid """ - self.train = True + self._train = True self.recompile() def save_grid(self, file_name): From a234db47fee0f8e81da3baaebd0f19b7eeef573f Mon Sep 17 00:00:00 2001 From: juacrumar Date: Wed, 10 Feb 2021 17:07:49 +0100 Subject: [PATCH 04/19] use directly cma --- src/vegasflow/rtbm.py | 51 +++++++++++++++++++++++++++++++------------ 1 file changed, 37 insertions(+), 14 deletions(-) diff --git a/src/vegasflow/rtbm.py b/src/vegasflow/rtbm.py index 19e46f8..4396968 100644 --- a/src/vegasflow/rtbm.py +++ b/src/vegasflow/rtbm.py @@ -1,15 +1,15 @@ """ Plain implementation of the plainest possible MonteCarlo """ - +import copy +import numpy as np from vegasflow.configflow import DTYPE, fone, fzero, float_me from vegasflow.monte_carlo import MonteCarloFlow, wrapper -import numpy as np import tensorflow as tf from theta.rtbm import RTBM -from theta.minimizer import CMA from theta import costfunctions +from cma import CMAEvolutionStrategy import logging @@ -32,7 +32,7 @@ def __init__(self, n_hidden=2, rtbm=None, train=True, *args, **kwargs): self._rtbm = RTBM( self.n_dim, n_hidden, - minimization_bound=50, + minimization_bound=80, gaussian_init=True, diagonal_T=False, positive_T=True, @@ -53,7 +53,6 @@ def unfreeze(self): def compile(self, integrand, compilable=False, **kwargs): if compilable: logger.warning("RTBMFlow is still WIP and not compilable") - self._trainer = CMA(parallel=False, ncores=1, verbose=True) super().compile(integrand, compilable=False, **kwargs) def generate_random_array(self, n_events): @@ -71,20 +70,44 @@ def generate_random_array(self, n_events): delta = max_per_d - min_per_d new_rand = (xrand - min_per_d) / delta xjac = xjac_raw / np.prod(delta) + print(delta) + print("") return float_me(new_rand), None, xjac def _train_machine(self, x, yraw): + # Get a reference to the initial solution of the CMA + x0 = copy.deepcopy(self._rtbm.get_parameters()) + bounds = self._rtbm.get_bounds() + options = { - "popsize": 75, - "tolfun": 1e-4, - "maxiter": 4, - } - # The output is a probability, therefore: - y = yraw / tf.reduce_sum(yraw) - # We don't want to train tf variables because it would be slow here... + "bounds": bounds, + "maxiter": 250, + } + xnp = x.numpy().T - ynp = y.numpy() - self._trainer.train(costfunctions.kullbackLeibler, self._rtbm, xnp, ynp, **options) + ynp = yraw.numpy() + + sol_found = False + def optimization(n=1): + sigma = np.min(bounds[1])/(4.0*n) + es = CMAEvolutionStrategy(x0, sigma, options) + + def target(params): + if not self._rtbm.set_parameters(params): + return np.NaN + prob = self._rtbm(xnp) + return costfunctions.kullbackLeibler.cost(prob, ynp) + + es.optimize(target) + return es.result + + n = 1 + while not sol_found: + res = optimization(n) + sol_found = self._rtbm.set_parameters(res.xbest) + logger.warning("Optimization failed, trying again!") + n+=1 + self._first_run = False @staticmethod From aa4d4966916d322ff33161a05eeb35f20ea9e08b Mon Sep 17 00:00:00 2001 From: juacrumar Date: Thu, 11 Feb 2021 12:20:10 +0100 Subject: [PATCH 05/19] only show the warning on failure --- examples/simgauss_tf.py | 4 ++-- examples/singletop_lo_tf.py | 4 ++-- src/vegasflow/rtbm.py | 3 ++- 3 files changed, 6 insertions(+), 5 deletions(-) diff --git a/examples/simgauss_tf.py b/examples/simgauss_tf.py index b2e680c..8066379 100644 --- a/examples/simgauss_tf.py +++ b/examples/simgauss_tf.py @@ -18,8 +18,8 @@ # MC integration setup -dim = 2 -ncalls = np.int32(1e6) +dim = 4 +ncalls = np.int32(1e5) n_iter = 3 diff --git a/examples/singletop_lo_tf.py b/examples/singletop_lo_tf.py index 0ed90af..6f48da9 100644 --- a/examples/singletop_lo_tf.py +++ b/examples/singletop_lo_tf.py @@ -12,8 +12,8 @@ # MC integration setup dim = 3 -ncalls = np.int32(1e5) -n_iter = 5 +ncalls = np.int32(1e6) +n_iter = 4 # Physics setup # top mass diff --git a/src/vegasflow/rtbm.py b/src/vegasflow/rtbm.py index 4396968..f717063 100644 --- a/src/vegasflow/rtbm.py +++ b/src/vegasflow/rtbm.py @@ -105,7 +105,8 @@ def target(params): while not sol_found: res = optimization(n) sol_found = self._rtbm.set_parameters(res.xbest) - logger.warning("Optimization failed, trying again!") + if not sol_found: + logger.warning("Optimization failed, trying again!") n+=1 self._first_run = False From fa66ae4a9f36035247618d7a60ea52f5440b07a2 Mon Sep 17 00:00:00 2001 From: juacrumar Date: Thu, 11 Feb 2021 16:58:16 +0100 Subject: [PATCH 06/19] use the rtbm own system --- src/vegasflow/rtbm.py | 36 +++++++++++++++--------------------- 1 file changed, 15 insertions(+), 21 deletions(-) diff --git a/src/vegasflow/rtbm.py b/src/vegasflow/rtbm.py index f717063..e93ab01 100644 --- a/src/vegasflow/rtbm.py +++ b/src/vegasflow/rtbm.py @@ -32,7 +32,7 @@ def __init__(self, n_hidden=2, rtbm=None, train=True, *args, **kwargs): self._rtbm = RTBM( self.n_dim, n_hidden, - minimization_bound=80, + minimization_bound=50, gaussian_init=True, diagonal_T=False, positive_T=True, @@ -56,25 +56,19 @@ def compile(self, integrand, compilable=False, **kwargs): super().compile(integrand, compilable=False, **kwargs) def generate_random_array(self, n_events): - if self._first_run: - return super().generate_random_array(n_events) - xrand, _ = self._rtbm.make_sample(n_events) - xjac_raw = 1.0 / self._rtbm(xrand.T) / n_events - # Now re-scale the values to the 0-1 range per dimension. - # NOTE: this is only valid while the points in_device are equal to the points being utilized - # so when the RTBM can run in the GPU it will need to be modified! - # mainly to carry the first limits to the rest of the calculation - epsilon = np.abs(np.max(xrand) / 4.0) - max_per_d = np.max(xrand, axis=0) + epsilon - min_per_d = np.min(xrand, axis=0) - epsilon - delta = max_per_d - min_per_d - new_rand = (xrand - min_per_d) / delta - xjac = xjac_raw / np.prod(delta) - print(delta) - print("") - return float_me(new_rand), None, xjac + """ + Returns (xrand, original_r, xjac) + where xrand is the integration variable between 0 and 1 + and xjac the correspondent jacobian factor. + original_r is the original random values sampled by the RTBM to be used at training + """ + xrand, original_r, px = self._rtbm.make_sample_rho(n_events) + xjac = float_me(1.0/px) + return float_me(xrand), original_r, xjac def _train_machine(self, x, yraw): + """ Takes as input the random sample from the previous step + and the values of the calculation""" # Get a reference to the initial solution of the CMA x0 = copy.deepcopy(self._rtbm.get_parameters()) bounds = self._rtbm.get_bounds() @@ -134,19 +128,19 @@ def _run_event(self, integrand, ncalls=None): n_events = ncalls # Generate all random number for this iteration - rnds, _, xjac = self.generate_random_array(n_events) + rnds, original_r, xjac = self.generate_random_array(n_events) # Compute the integrand res = integrand(rnds, n_dim=self.n_dim, weight=xjac) * xjac # Clean up the array from numbers outside the 0-1 range - if not self._first_run: + if True: # They should always be between 0 and 1 anyway # Accept for now only random number between 0 and 1 condition = tf.reduce_all(rnds >= 0.0, axis=1) & tf.reduce_all(rnds <= 1.0, axis=1) res = tf.where(condition, res, fzero)[0] if np.count_nonzero(res) != n_events: logger.warning(f"Found only {np.count_nonzero(res)} of {n_events} valid values\n") - return res, rnds + return res, original_r def _run_iteration(self): all_res, rnds = self.run_event() From 53c8e081b1d41541019ee507f2891ebe7efbae5f Mon Sep 17 00:00:00 2001 From: juacrumar Date: Thu, 11 Feb 2021 17:20:36 +0100 Subject: [PATCH 07/19] fix --- src/vegasflow/rtbm.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/vegasflow/rtbm.py b/src/vegasflow/rtbm.py index e93ab01..e1bbfc2 100644 --- a/src/vegasflow/rtbm.py +++ b/src/vegasflow/rtbm.py @@ -62,8 +62,8 @@ def generate_random_array(self, n_events): and xjac the correspondent jacobian factor. original_r is the original random values sampled by the RTBM to be used at training """ - xrand, original_r, px = self._rtbm.make_sample_rho(n_events) - xjac = float_me(1.0/px) + xrand, px, original_r = self._rtbm.make_sample_rho(n_events) + xjac = float_me(1.0/px/n_events) return float_me(xrand), original_r, xjac def _train_machine(self, x, yraw): @@ -136,7 +136,7 @@ def _run_event(self, integrand, ncalls=None): if True: # They should always be between 0 and 1 anyway # Accept for now only random number between 0 and 1 condition = tf.reduce_all(rnds >= 0.0, axis=1) & tf.reduce_all(rnds <= 1.0, axis=1) - res = tf.where(condition, res, fzero)[0] + res = tf.where(condition, res, fzero) if np.count_nonzero(res) != n_events: logger.warning(f"Found only {np.count_nonzero(res)} of {n_events} valid values\n") From aab04e12139073e1e9f49d92e7cc738bd4276ed4 Mon Sep 17 00:00:00 2001 From: juacrumar Date: Thu, 11 Feb 2021 17:27:07 +0100 Subject: [PATCH 08/19] pylint --- src/vegasflow/rtbm.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/vegasflow/rtbm.py b/src/vegasflow/rtbm.py index e1bbfc2..9e074a9 100644 --- a/src/vegasflow/rtbm.py +++ b/src/vegasflow/rtbm.py @@ -7,9 +7,9 @@ from vegasflow.monte_carlo import MonteCarloFlow, wrapper import tensorflow as tf -from theta.rtbm import RTBM -from theta import costfunctions -from cma import CMAEvolutionStrategy +from theta.rtbm import RTBM # pylint:disable=import-error +from theta import costfunctions # pylint:disable=import-error +from cma import CMAEvolutionStrategy # pylint:disable=import-error import logging From 520d5f7fcdaf04e98c626286672929a654c40b9e Mon Sep 17 00:00:00 2001 From: juacrumar Date: Fri, 12 Feb 2021 17:35:33 +0100 Subject: [PATCH 09/19] add a simpler integrand --- examples/simgauss_tf.py | 1 - examples/sin_tf.py | 49 +++++++++++++++++++++++++++++++++++++++++ 2 files changed, 49 insertions(+), 1 deletion(-) create mode 100644 examples/sin_tf.py diff --git a/examples/simgauss_tf.py b/examples/simgauss_tf.py index 8066379..a179102 100644 --- a/examples/simgauss_tf.py +++ b/examples/simgauss_tf.py @@ -8,7 +8,6 @@ tf.config.threading.set_inter_op_parallelism_threads(1) tf.config.threading.set_intra_op_parallelism_threads(1) from vegasflow.configflow import DTYPE, run_eager -run_eager(True) import time import numpy as np diff --git a/examples/sin_tf.py b/examples/sin_tf.py new file mode 100644 index 0000000..871ee38 --- /dev/null +++ b/examples/sin_tf.py @@ -0,0 +1,49 @@ +""" + Example: basic integration + + Basic example using the vegas_wrapper helper +""" + +import tensorflow as tf +tf.config.threading.set_inter_op_parallelism_threads(1) +tf.config.threading.set_intra_op_parallelism_threads(1) +from vegasflow import run_eager, float_me +run_eager(True) + +import time +import numpy as np +from vegasflow.vflow import vegas_wrapper +from vegasflow.rtbm import rtbm_wrapper + + +# MC integration setup +dim = 2 +ncalls = np.int32(1e5) +n_iter = 3 +tf_pi = float_me(np.pi) + + +@tf.function +def sin_fun(xarr, **kwargs): + """symgauss test function""" + res = tf.pow(tf.sin(2.0*xarr*tf_pi),2) + return tf.reduce_prod(res, axis=1) + + +if __name__ == "__main__": + """Testing several different integrations""" + print(f"VEGAS MC, ncalls={ncalls}:") + start = time.time() + ncalls = ncalls + r = vegas_wrapper(sin_fun, dim, n_iter, ncalls) + end = time.time() + print(f"Vegas took: time (s): {end-start}") + + print(f"RTBM MC, ncalls={ncalls}:") + start = time.time() + ncalls = ncalls + r = rtbm_wrapper(sin_fun, dim, n_iter, ncalls) + end = time.time() + print(f"RTBM took: time (s): {end-start}") + + From 5ce4b677c71940a0a828a8b9244c960d5246d0b7 Mon Sep 17 00:00:00 2001 From: juacrumar Date: Tue, 16 Feb 2021 19:16:38 +0100 Subject: [PATCH 10/19] fixes --- examples/simgauss_tf.py | 2 +- examples/sin_tf.py | 15 +++++++++------ src/vegasflow/rtbm.py | 10 ++++++++-- 3 files changed, 18 insertions(+), 9 deletions(-) diff --git a/examples/simgauss_tf.py b/examples/simgauss_tf.py index a179102..909792b 100644 --- a/examples/simgauss_tf.py +++ b/examples/simgauss_tf.py @@ -18,7 +18,7 @@ # MC integration setup dim = 4 -ncalls = np.int32(1e5) +ncalls = np.int32(1e3) n_iter = 3 diff --git a/examples/sin_tf.py b/examples/sin_tf.py index 871ee38..75793da 100644 --- a/examples/sin_tf.py +++ b/examples/sin_tf.py @@ -5,8 +5,8 @@ """ import tensorflow as tf -tf.config.threading.set_inter_op_parallelism_threads(1) -tf.config.threading.set_intra_op_parallelism_threads(1) +tf.config.threading.set_inter_op_parallelism_threads(2) +tf.config.threading.set_intra_op_parallelism_threads(2) from vegasflow import run_eager, float_me run_eager(True) @@ -17,9 +17,9 @@ # MC integration setup -dim = 2 -ncalls = np.int32(1e5) -n_iter = 3 +dim = 3 +ncalls = np.int32(1e4) +n_iter = 12 tf_pi = float_me(np.pi) @@ -42,8 +42,11 @@ def sin_fun(xarr, **kwargs): print(f"RTBM MC, ncalls={ncalls}:") start = time.time() ncalls = ncalls - r = rtbm_wrapper(sin_fun, dim, n_iter, ncalls) + rt = rtbm_wrapper(sin_fun, dim, n_iter, ncalls) end = time.time() print(f"RTBM took: time (s): {end-start}") + print(f"Result computed by Vegas: {r}") + print(f"Result computed by RTBM: {rt}") + diff --git a/src/vegasflow/rtbm.py b/src/vegasflow/rtbm.py index 9e074a9..5f62ae5 100644 --- a/src/vegasflow/rtbm.py +++ b/src/vegasflow/rtbm.py @@ -21,7 +21,7 @@ class RTBMFlow(MonteCarloFlow): RTBM based Monte Carlo integrator """ - def __init__(self, n_hidden=2, rtbm=None, train=True, *args, **kwargs): + def __init__(self, n_hidden=3, rtbm=None, train=True, *args, **kwargs): super().__init__(*args, **kwargs) self._train = train self._first_run = True @@ -76,6 +76,7 @@ def _train_machine(self, x, yraw): options = { "bounds": bounds, "maxiter": 250, + "popsize" : 50, } xnp = x.numpy().T @@ -90,6 +91,8 @@ def target(params): if not self._rtbm.set_parameters(params): return np.NaN prob = self._rtbm(xnp) + if (prob <= 0.0).any(): + return np.NaN return costfunctions.kullbackLeibler.cost(prob, ynp) es.optimize(target) @@ -98,7 +101,10 @@ def target(params): n = 1 while not sol_found: res = optimization(n) - sol_found = self._rtbm.set_parameters(res.xbest) + try: + sol_found = self._rtbm.set_parameters(res.xbest) + except TypeError: + sol_found = False if not sol_found: logger.warning("Optimization failed, trying again!") n+=1 From ccb7d632ce607faf9fcbfca94d918e6bf69189d3 Mon Sep 17 00:00:00 2001 From: juacrumar Date: Thu, 18 Mar 2021 14:26:58 +0100 Subject: [PATCH 11/19] update the rtbm integration --- examples/sin_tf.py | 6 +- src/vegasflow/rtbm.py | 163 ++++++++++++++++++++++++++---------------- 2 files changed, 106 insertions(+), 63 deletions(-) diff --git a/examples/sin_tf.py b/examples/sin_tf.py index 75793da..c1c2e44 100644 --- a/examples/sin_tf.py +++ b/examples/sin_tf.py @@ -17,9 +17,9 @@ # MC integration setup -dim = 3 -ncalls = np.int32(1e4) -n_iter = 12 +dim = 2 +ncalls = np.int32(1e3) +n_iter = 5 tf_pi = float_me(np.pi) diff --git a/src/vegasflow/rtbm.py b/src/vegasflow/rtbm.py index 5f62ae5..9974d7f 100644 --- a/src/vegasflow/rtbm.py +++ b/src/vegasflow/rtbm.py @@ -12,9 +12,94 @@ from cma import CMAEvolutionStrategy # pylint:disable=import-error import logging +from joblib import Parallel, delayed +from time import time logger = logging.getLogger(__name__) +# Cost functions for training +def _kl(x, ytarget): + return ytarget*np.log(ytarget/x) + +_loss = _kl + +def _train_machine(rtbm, target_tf, original_r_tf): + + logger.info("Training RTBM") + + target = target_tf.numpy() + original_r = original_r_tf.numpy() + + n_jobs = 8 + max_iterations = 250 + + + def target_loss(params): + if not rtbm.set_parameters(params): + return np.NaN + _, prob = rtbm.get_transformation(original_r) + return np.sum(_loss(prob, target)) + + best_parameters = copy.deepcopy(rtbm.get_parameters()) + min_bound, max_bound = rtbm.get_bounds() + + with Parallel(n_jobs=n_jobs) as parallel: + prev = time() + n_parameters = len(best_parameters) + + # random hyperparameters + pop_per_rate = 32 + mutation_rates = np.array([0.2, 0.4, 0.6, 0.8]) + rates = np.concatenate([np.ones(pop_per_rate)*mr for mr in mutation_rates]) + original_sigma = 0.25 + sigma = original_sigma + repeats = 3 + + for it in range(max_iterations): + + # Get the best parameters from the previous iteration + p0 = copy.deepcopy(best_parameters) + loss_val = target_loss(p0) + + def compute_mutant(mutation_rate): + number_of_mut = int(mutation_rate*n_parameters) + mut_idx = np.random.choice(n_parameters, number_of_mut, replace=False) + r1, r2 = np.random.rand(2, number_of_mut)*sigma + + mutant = copy.deepcopy(p0) + var_plus = max_bound - p0 + var_minus = min_bound - p0 + mutant[mut_idx] += var_plus[mut_idx]*r1 + var_minus[mut_idx]*r2 + + return target_loss(mutant), mutant + + parallel_runs = [delayed(compute_mutant)(rate) for rate in rates] + result = parallel(parallel_runs) + losses, mutants = zip(*result) + + best_loss = np.nanmin(losses) + if best_loss < loss_val: + loss_val = best_loss + best_parameters = mutants[losses.index(best_loss)] + else: + sigma /= 2 + + if it % 10 == 0: + current = time() + print(f"Iteration {it}, best loss: {loss_val:.2f}, time; {current-prev:.2f}s") + prev = current + + if sigma < 1e-3: + sigma = original_sigma + print(f"Resetting sigma with loss: {loss_val:.2f}") + repeats -= 1 + break + + if not repeats: + print(f"No more repeats allowed, iteration: {it}, loss: {loss_val:2.f}") + + rtbm.set_parameters(best_parameters) + return rtbm class RTBMFlow(MonteCarloFlow): """ @@ -34,21 +119,27 @@ def __init__(self, n_hidden=3, rtbm=None, train=True, *args, **kwargs): n_hidden, minimization_bound=50, gaussian_init=True, - diagonal_T=False, positive_T=True, positive_Q=True, + gaussian_parameters={"mean":0.0, "std": 0.55}, + sampling_activation="tanh" ) else: # Check whether it is a valid rtbm model if not hasattr(rtbm, "make_sample"): raise TypeError(f"{rtbm} is not a valid boltzman machine") self._rtbm = rtbm + self._p0 = self._rtbm.get_parameters() def freeze(self): + """ Stop the training """ self.train = False - def unfreeze(self): + def unfreeze(self, reset_p0 = False): + """ Restart the training """ self.train = True + if reset_p0: + self._p0 = self._rtbm.get_parameters() def compile(self, integrand, compilable=False, **kwargs): if compilable: @@ -61,55 +152,13 @@ def generate_random_array(self, n_events): where xrand is the integration variable between 0 and 1 and xjac the correspondent jacobian factor. original_r is the original random values sampled by the RTBM to be used at training + + The return dimension of the random variables is (n_events,n_dim) and the jacobian (n_events) """ xrand, px, original_r = self._rtbm.make_sample_rho(n_events) - xjac = float_me(1.0/px/n_events) - return float_me(xrand), original_r, xjac - - def _train_machine(self, x, yraw): - """ Takes as input the random sample from the previous step - and the values of the calculation""" - # Get a reference to the initial solution of the CMA - x0 = copy.deepcopy(self._rtbm.get_parameters()) - bounds = self._rtbm.get_bounds() - - options = { - "bounds": bounds, - "maxiter": 250, - "popsize" : 50, - } - - xnp = x.numpy().T - ynp = yraw.numpy() - - sol_found = False - def optimization(n=1): - sigma = np.min(bounds[1])/(4.0*n) - es = CMAEvolutionStrategy(x0, sigma, options) - - def target(params): - if not self._rtbm.set_parameters(params): - return np.NaN - prob = self._rtbm(xnp) - if (prob <= 0.0).any(): - return np.NaN - return costfunctions.kullbackLeibler.cost(prob, ynp) - - es.optimize(target) - return es.result - - n = 1 - while not sol_found: - res = optimization(n) - try: - sol_found = self._rtbm.set_parameters(res.xbest) - except TypeError: - sol_found = False - if not sol_found: - logger.warning("Optimization failed, trying again!") - n+=1 - - self._first_run = False + # Since we are using the tanh function, the integration limits are (-1,1), move: + xjac = 1.0/px/n_events + return float_me(xrand), original_r, float_me(xjac) @staticmethod def _accumulate(accumulators): @@ -136,22 +185,16 @@ def _run_event(self, integrand, ncalls=None): # Generate all random number for this iteration rnds, original_r, xjac = self.generate_random_array(n_events) # Compute the integrand - res = integrand(rnds, n_dim=self.n_dim, weight=xjac) * xjac - - # Clean up the array from numbers outside the 0-1 range - if True: # They should always be between 0 and 1 anyway - # Accept for now only random number between 0 and 1 - condition = tf.reduce_all(rnds >= 0.0, axis=1) & tf.reduce_all(rnds <= 1.0, axis=1) - res = tf.where(condition, res, fzero) - if np.count_nonzero(res) != n_events: - logger.warning(f"Found only {np.count_nonzero(res)} of {n_events} valid values\n") + tmp = integrand(rnds, n_dim=self.n_dim, weight=xjac) + res = tmp*xjac return res, original_r def _run_iteration(self): - all_res, rnds = self.run_event() + all_res, original_r = self.run_event() + if self.train: - self._train_machine(rnds, all_res) + _train_machine(self._rtbm, all_res, original_r) res = tf.reduce_sum(all_res) all_res2 = all_res ** 2 From cfef2300970bff7cfe6b96c6c5ca8d53ab8d8515 Mon Sep 17 00:00:00 2001 From: juacrumar Date: Thu, 18 Mar 2021 14:29:00 +0100 Subject: [PATCH 12/19] modify examples --- examples/sin_tf.py | 21 +++++++++++---------- examples/singletop_lo_tf.py | 2 +- 2 files changed, 12 insertions(+), 11 deletions(-) diff --git a/examples/sin_tf.py b/examples/sin_tf.py index 75793da..cedb711 100644 --- a/examples/sin_tf.py +++ b/examples/sin_tf.py @@ -5,28 +5,28 @@ """ import tensorflow as tf -tf.config.threading.set_inter_op_parallelism_threads(2) -tf.config.threading.set_intra_op_parallelism_threads(2) +tf.config.threading.set_inter_op_parallelism_threads(4) +tf.config.threading.set_intra_op_parallelism_threads(4) from vegasflow import run_eager, float_me run_eager(True) import time import numpy as np from vegasflow.vflow import vegas_wrapper -from vegasflow.rtbm import rtbm_wrapper +from vegasflow.rtbm import RTBMFlow # MC integration setup dim = 3 -ncalls = np.int32(1e4) -n_iter = 12 +ncalls = np.int32(1e5) +n_iter = 5 tf_pi = float_me(np.pi) @tf.function def sin_fun(xarr, **kwargs): """symgauss test function""" - res = tf.pow(tf.sin(2.0*xarr*tf_pi),2) + res = tf.pow(tf.sin(xarr*tf_pi*4.0),2) return tf.reduce_prod(res, axis=1) @@ -41,12 +41,13 @@ def sin_fun(xarr, **kwargs): print(f"RTBM MC, ncalls={ncalls}:") start = time.time() - ncalls = ncalls - rt = rtbm_wrapper(sin_fun, dim, n_iter, ncalls) + rtbm = RTBMFlow(n_dim=dim, n_events=ncalls, gaussian=True, train=True, n_hidden=2) + rtbm.compile(sin_fun) + rtbm.run_integration(4) end = time.time() print(f"RTBM took: time (s): {end-start}") - print(f"Result computed by Vegas: {r}") - print(f"Result computed by RTBM: {rt}") +# print(f"Result computed by Vegas: {r}") +# print(f"Result computed by RTBM: {rt}") diff --git a/examples/singletop_lo_tf.py b/examples/singletop_lo_tf.py index 6f48da9..028ea5a 100644 --- a/examples/singletop_lo_tf.py +++ b/examples/singletop_lo_tf.py @@ -12,7 +12,7 @@ # MC integration setup dim = 3 -ncalls = np.int32(1e6) +ncalls = np.int32(1e5) n_iter = 4 # Physics setup From fe04a841db3bc413027ea62275694f3052fc135b Mon Sep 17 00:00:00 2001 From: juacrumar Date: Thu, 18 Mar 2021 18:30:45 +0100 Subject: [PATCH 13/19] outperforms Vegas! --- examples/sin_tf.py | 18 ++++---- src/vegasflow/configflow.py | 2 +- src/vegasflow/rtbm.py | 82 ++++++++++++++++++++++--------------- 3 files changed, 59 insertions(+), 43 deletions(-) diff --git a/examples/sin_tf.py b/examples/sin_tf.py index 1af7cbd..6987c01 100644 --- a/examples/sin_tf.py +++ b/examples/sin_tf.py @@ -4,10 +4,8 @@ Basic example using the vegas_wrapper helper """ -import tensorflow as tf -tf.config.threading.set_inter_op_parallelism_threads(4) -tf.config.threading.set_intra_op_parallelism_threads(4) from vegasflow import run_eager, float_me +import tensorflow as tf run_eager(True) import time @@ -17,8 +15,8 @@ # MC integration setup -dim = 3 -ncalls = np.int32(1e3) +dim = 1 +ncalls = np.int32(1e4) n_iter = 5 tf_pi = float_me(np.pi) @@ -26,7 +24,7 @@ @tf.function def sin_fun(xarr, **kwargs): """symgauss test function""" - res = tf.pow(tf.sin(xarr*tf_pi*4.0),2) + res = tf.pow(tf.sin(xarr*tf_pi),2) return tf.reduce_prod(res, axis=1) @@ -41,13 +39,13 @@ def sin_fun(xarr, **kwargs): print(f"RTBM MC, ncalls={ncalls}:") start = time.time() - rtbm = RTBMFlow(n_dim=dim, n_events=ncalls, gaussian=True, train=True, n_hidden=2) + rtbm = RTBMFlow(n_dim=dim, n_events=ncalls, train=True, n_hidden=1) rtbm.compile(sin_fun) - rtbm.run_integration(4) + rt = rtbm.run_integration(4) end = time.time() print(f"RTBM took: time (s): {end-start}") -# print(f"Result computed by Vegas: {r}") -# print(f"Result computed by RTBM: {rt}") + print(f"Result computed by Vegas: {r}") + print(f"Result computed by RTBM: {rt}") diff --git a/src/vegasflow/configflow.py b/src/vegasflow/configflow.py index f5c0e3c..bd91759 100644 --- a/src/vegasflow/configflow.py +++ b/src/vegasflow/configflow.py @@ -7,7 +7,7 @@ import numpy as np # Set TF to only log errors -os.environ.setdefault("TF_CPP_MIN_LOG_LEVEL", "1") +os.environ.setdefault("TF_CPP_MIN_LOG_LEVEL", "3") import tensorflow as tf # uncomment this line for debugging to avoid compiling any tf.function diff --git a/src/vegasflow/rtbm.py b/src/vegasflow/rtbm.py index 9974d7f..5eeacbf 100644 --- a/src/vegasflow/rtbm.py +++ b/src/vegasflow/rtbm.py @@ -7,9 +7,9 @@ from vegasflow.monte_carlo import MonteCarloFlow, wrapper import tensorflow as tf -from theta.rtbm import RTBM # pylint:disable=import-error -from theta import costfunctions # pylint:disable=import-error -from cma import CMAEvolutionStrategy # pylint:disable=import-error +from theta.rtbm import RTBM # pylint:disable=import-error +from theta import costfunctions # pylint:disable=import-error +from cma import CMAEvolutionStrategy # pylint:disable=import-error import logging from joblib import Parallel, delayed @@ -19,10 +19,12 @@ # Cost functions for training def _kl(x, ytarget): - return ytarget*np.log(ytarget/x) + return ytarget * np.log(ytarget / x) + _loss = _kl + def _train_machine(rtbm, target_tf, original_r_tf): logger.info("Training RTBM") @@ -30,9 +32,8 @@ def _train_machine(rtbm, target_tf, original_r_tf): target = target_tf.numpy() original_r = original_r_tf.numpy() - n_jobs = 8 - max_iterations = 250 - + n_jobs = 32 + max_iterations = 3000 def target_loss(params): if not rtbm.set_parameters(params): @@ -42,34 +43,35 @@ def target_loss(params): best_parameters = copy.deepcopy(rtbm.get_parameters()) min_bound, max_bound = rtbm.get_bounds() + loss_val = target_loss(best_parameters) with Parallel(n_jobs=n_jobs) as parallel: prev = time() n_parameters = len(best_parameters) - # random hyperparameters - pop_per_rate = 32 + # training hyperparameters + pop_per_rate = 256 mutation_rates = np.array([0.2, 0.4, 0.6, 0.8]) - rates = np.concatenate([np.ones(pop_per_rate)*mr for mr in mutation_rates]) + rates = np.concatenate([np.ones(pop_per_rate) * mr for mr in mutation_rates]) original_sigma = 0.25 sigma = original_sigma repeats = 3 + ####### for it in range(max_iterations): # Get the best parameters from the previous iteration p0 = copy.deepcopy(best_parameters) - loss_val = target_loss(p0) def compute_mutant(mutation_rate): - number_of_mut = int(mutation_rate*n_parameters) + number_of_mut = int(mutation_rate * n_parameters) mut_idx = np.random.choice(n_parameters, number_of_mut, replace=False) - r1, r2 = np.random.rand(2, number_of_mut)*sigma + r1, r2 = np.random.rand(2, number_of_mut) * sigma mutant = copy.deepcopy(p0) var_plus = max_bound - p0 var_minus = min_bound - p0 - mutant[mut_idx] += var_plus[mut_idx]*r1 + var_minus[mut_idx]*r2 + mutant[mut_idx] += var_plus[mut_idx] * r1 + var_minus[mut_idx] * r2 return target_loss(mutant), mutant @@ -84,14 +86,25 @@ def compute_mutant(mutation_rate): else: sigma /= 2 - if it % 10 == 0: + if it % 50 == 0: current = time() - print(f"Iteration {it}, best loss: {loss_val:.2f}, time; {current-prev:.2f}s") + logger.info( + "Iteration %d, best_loss: %.2f, (%2.fs)", + it, + loss_val, + current - prev, + ) + logger.debug( + "NaNs in last iteration: %d/%d, sigma=%.4f", + np.count_nonzero(np.isnan(losses)), + len(losses), + sigma, + ) prev = current if sigma < 1e-3: sigma = original_sigma - print(f"Resetting sigma with loss: {loss_val:.2f}") + logger.debug("Resetting sigma with loss: %.2f", loss_val) repeats -= 1 break @@ -101,6 +114,7 @@ def compute_mutant(mutation_rate): rtbm.set_parameters(best_parameters) return rtbm + class RTBMFlow(MonteCarloFlow): """ RTBM based Monte Carlo integrator @@ -112,7 +126,8 @@ def __init__(self, n_hidden=3, rtbm=None, train=True, *args, **kwargs): self._first_run = True if rtbm is None: logger.info( - "Generating a RTBM with %d visible nodes and %d hidden" % (self.n_dim, n_hidden) + "Generating a RTBM with %d visible nodes and %d hidden" + % (self.n_dim, n_hidden) ) self._rtbm = RTBM( self.n_dim, @@ -121,8 +136,8 @@ def __init__(self, n_hidden=3, rtbm=None, train=True, *args, **kwargs): gaussian_init=True, positive_T=True, positive_Q=True, - gaussian_parameters={"mean":0.0, "std": 0.55}, - sampling_activation="tanh" + gaussian_parameters={"mean": 0.0, "std": 0.55}, + sampling_activation="tanh", ) else: # Check whether it is a valid rtbm model @@ -135,7 +150,7 @@ def freeze(self): """ Stop the training """ self.train = False - def unfreeze(self, reset_p0 = False): + def unfreeze(self, reset_p0=False): """ Restart the training """ self.train = True if reset_p0: @@ -157,7 +172,7 @@ def generate_random_array(self, n_events): """ xrand, px, original_r = self._rtbm.make_sample_rho(n_events) # Since we are using the tanh function, the integration limits are (-1,1), move: - xjac = 1.0/px/n_events + xjac = 1.0 / px / n_events return float_me(xrand), original_r, float_me(xjac) @staticmethod @@ -167,14 +182,17 @@ def _accumulate(accumulators): # In the accumulators we should receive a number of items with # (res, xrands) which have shape ( (n_events,), (n_events, n_dim) ) all_res = [] + all_unw = [] all_rnd = [] - for ( - res, - rnds, - ) in accumulators: + for (res, unw, rnds) in accumulators: all_res.append(res) + all_unw.append(unw) all_rnd.append(rnds) - return tf.concat(all_res, axis=0), tf.concat(all_rnd, axis=0) + return ( + tf.concat(all_res, axis=0), + tf.concat(all_unw, axis=0), + tf.concat(all_rnd, axis=0), + ) def _run_event(self, integrand, ncalls=None): if ncalls is None: @@ -185,16 +203,16 @@ def _run_event(self, integrand, ncalls=None): # Generate all random number for this iteration rnds, original_r, xjac = self.generate_random_array(n_events) # Compute the integrand - tmp = integrand(rnds, n_dim=self.n_dim, weight=xjac) - res = tmp*xjac + unw = integrand(rnds, n_dim=self.n_dim, weight=xjac) + res = unw * xjac - return res, original_r + return res, unw, original_r def _run_iteration(self): - all_res, original_r = self.run_event() + all_res, unw, original_r = self.run_event() if self.train: - _train_machine(self._rtbm, all_res, original_r) + _train_machine(self._rtbm, unw, original_r) res = tf.reduce_sum(all_res) all_res2 = all_res ** 2 From 71592d6db47e322a5a94b9aaa698060c196ec9f3 Mon Sep 17 00:00:00 2001 From: juacrumar Date: Fri, 19 Mar 2021 10:44:52 +0100 Subject: [PATCH 14/19] still better for nh=2 --- examples/sin_tf.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/sin_tf.py b/examples/sin_tf.py index 6987c01..88f6779 100644 --- a/examples/sin_tf.py +++ b/examples/sin_tf.py @@ -15,7 +15,7 @@ # MC integration setup -dim = 1 +dim = 2 ncalls = np.int32(1e4) n_iter = 5 tf_pi = float_me(np.pi) From e3f165b74c90c9ffa892faae8ad31fc1cf1924e1 Mon Sep 17 00:00:00 2001 From: juacrumar Date: Fri, 19 Mar 2021 20:12:58 +0100 Subject: [PATCH 15/19] trying dim=5 --- examples/simgauss_tf.py | 2 -- examples/sin_tf.py | 20 ++++++++++++-------- src/vegasflow/rtbm.py | 4 ++-- 3 files changed, 14 insertions(+), 12 deletions(-) diff --git a/examples/simgauss_tf.py b/examples/simgauss_tf.py index 909792b..42c6235 100644 --- a/examples/simgauss_tf.py +++ b/examples/simgauss_tf.py @@ -5,8 +5,6 @@ """ import tensorflow as tf -tf.config.threading.set_inter_op_parallelism_threads(1) -tf.config.threading.set_intra_op_parallelism_threads(1) from vegasflow.configflow import DTYPE, run_eager import time diff --git a/examples/sin_tf.py b/examples/sin_tf.py index 88f6779..d9ceb61 100644 --- a/examples/sin_tf.py +++ b/examples/sin_tf.py @@ -15,7 +15,7 @@ # MC integration setup -dim = 2 +dim = 5 ncalls = np.int32(1e4) n_iter = 5 tf_pi = float_me(np.pi) @@ -27,25 +27,29 @@ def sin_fun(xarr, **kwargs): res = tf.pow(tf.sin(xarr*tf_pi),2) return tf.reduce_prod(res, axis=1) +# integrand = sin_fun +from simgauss_tf import symgauss as integrand if __name__ == "__main__": """Testing several different integrations""" print(f"VEGAS MC, ncalls={ncalls}:") start = time.time() ncalls = ncalls - r = vegas_wrapper(sin_fun, dim, n_iter, ncalls) + r = vegas_wrapper(integrand, dim, n_iter, ncalls) end = time.time() print(f"Vegas took: time (s): {end-start}") print(f"RTBM MC, ncalls={ncalls}:") start = time.time() - rtbm = RTBMFlow(n_dim=dim, n_events=ncalls, train=True, n_hidden=1) - rtbm.compile(sin_fun) - rt = rtbm.run_integration(4) + rtbm = RTBMFlow(n_dim=dim, n_events=ncalls, train=True, n_hidden=2) + rtbm.compile(integrand) + rt = rtbm.run_integration(n_iter) end = time.time() print(f"RTBM took: time (s): {end-start}") - print(f"Result computed by Vegas: {r}") - print(f"Result computed by RTBM: {rt}") - + print(f"Result computed by Vegas: {r[0]:.5f} +- {r[1]:.6f}") + print(f"Result computed by RTBM: {rt[0]:.5f} +- {rt[1]:.6f}") +# Notes +# For 1 and 2 dimensions is enough with n_hidden=1 to get a better per event error +# For 3 dimensions we need to go to n_hidden=2 diff --git a/src/vegasflow/rtbm.py b/src/vegasflow/rtbm.py index 5eeacbf..a837ad0 100644 --- a/src/vegasflow/rtbm.py +++ b/src/vegasflow/rtbm.py @@ -88,7 +88,7 @@ def compute_mutant(mutation_rate): if it % 50 == 0: current = time() - logger.info( + logger.debug( "Iteration %d, best_loss: %.2f, (%2.fs)", it, loss_val, @@ -136,7 +136,7 @@ def __init__(self, n_hidden=3, rtbm=None, train=True, *args, **kwargs): gaussian_init=True, positive_T=True, positive_Q=True, - gaussian_parameters={"mean": 0.0, "std": 0.55}, + gaussian_parameters={"mean": 0.0, "std": 0.75}, sampling_activation="tanh", ) else: From c118c1cc86edb8a33852ef6fafed1a43be060b55 Mon Sep 17 00:00:00 2001 From: juacrumar Date: Tue, 13 Apr 2021 10:20:18 +0200 Subject: [PATCH 16/19] add initialization --- examples/sin_tf.py | 10 ++- src/vegasflow/monte_carlo.py | 8 +- src/vegasflow/rtbm.py | 167 ++++++++++++++++++++++++++++------- 3 files changed, 146 insertions(+), 39 deletions(-) diff --git a/examples/sin_tf.py b/examples/sin_tf.py index d9ceb61..e5b728f 100644 --- a/examples/sin_tf.py +++ b/examples/sin_tf.py @@ -6,7 +6,6 @@ from vegasflow import run_eager, float_me import tensorflow as tf -run_eager(True) import time import numpy as np @@ -35,7 +34,8 @@ def sin_fun(xarr, **kwargs): print(f"VEGAS MC, ncalls={ncalls}:") start = time.time() ncalls = ncalls - r = vegas_wrapper(integrand, dim, n_iter, ncalls) + _, vegas_instance = vegas_wrapper(integrand, dim, n_iter, ncalls, return_instance=True) + vegas_instance.freeze_grid() end = time.time() print(f"Vegas took: time (s): {end-start}") @@ -43,10 +43,14 @@ def sin_fun(xarr, **kwargs): start = time.time() rtbm = RTBMFlow(n_dim=dim, n_events=ncalls, train=True, n_hidden=2) rtbm.compile(integrand) - rt = rtbm.run_integration(n_iter) + _ = rtbm.run_integration(n_iter) + rtbm.freeze() end = time.time() print(f"RTBM took: time (s): {end-start}") + print("Results with frozen grids") + vegas_instance.run_integration(5) + rt = rtbm.run_integration(5) print(f"Result computed by Vegas: {r[0]:.5f} +- {r[1]:.6f}") print(f"Result computed by RTBM: {rt[0]:.5f} +- {rt[1]:.6f}") diff --git a/src/vegasflow/monte_carlo.py b/src/vegasflow/monte_carlo.py index 10bf832..de76f0d 100644 --- a/src/vegasflow/monte_carlo.py +++ b/src/vegasflow/monte_carlo.py @@ -480,7 +480,7 @@ def run_integration(self, n_iter, log_time=True, histograms=None): return final_result.numpy(), sigma -def wrapper(integrator_class, integrand, n_dim, n_iter, total_n_events, compilable=True): +def wrapper(integrator_class, integrand, n_dim, n_iter, total_n_events, compilable=True, return_instance=False): """Convenience wrapper Parameters @@ -498,4 +498,8 @@ def wrapper(integrator_class, integrand, n_dim, n_iter, total_n_events, compilab """ mc_instance = integrator_class(n_dim = n_dim, n_events=total_n_events) mc_instance.compile(integrand, compilable=compilable) - return mc_instance.run_integration(n_iter) + if return_instance: + r = mc_instance.run_integration(n_iter) + return r, mc_instance + else: + return mc_instance.run_integration(n_iter) diff --git a/src/vegasflow/rtbm.py b/src/vegasflow/rtbm.py index a837ad0..1979b8d 100644 --- a/src/vegasflow/rtbm.py +++ b/src/vegasflow/rtbm.py @@ -2,14 +2,21 @@ Plain implementation of the plainest possible MonteCarlo """ import copy +import itertools import numpy as np -from vegasflow.configflow import DTYPE, fone, fzero, float_me +from vegasflow.configflow import DTYPE, fone, fzero, float_me, run_eager from vegasflow.monte_carlo import MonteCarloFlow, wrapper import tensorflow as tf -from theta.rtbm import RTBM # pylint:disable=import-error -from theta import costfunctions # pylint:disable=import-error -from cma import CMAEvolutionStrategy # pylint:disable=import-error +run_eager(True) + +try: + from theta.rtbm import RTBM # pylint:disable=import-error + from theta import costfunctions # pylint:disable=import-error +except ImportError as e: + raise ValueError( + "Cannot use the RTBM based integrator without installing the appropiate libraries" + ) from e import logging from joblib import Parallel, delayed @@ -25,22 +32,41 @@ def _kl(x, ytarget): _loss = _kl -def _train_machine(rtbm, target_tf, original_r_tf): - - logger.info("Training RTBM") - - target = target_tf.numpy() - original_r = original_r_tf.numpy() - - n_jobs = 32 - max_iterations = 3000 - +def _generate_target_loss(rtbm, original_r, target): def target_loss(params): if not rtbm.set_parameters(params): return np.NaN _, prob = rtbm.get_transformation(original_r) return np.sum(_loss(prob, target)) + return target_loss + + +def _train_machine( + rtbm, + target_tf, + original_r_tf, + n_jobs=32, + max_iterations=3000, + pop_per_rate=256, + verbose=True, +): + + if verbose: + logger.info("Training RTBM") + + if isinstance(target_tf, np.ndarray): + target = target_tf + else: + target = target_tf.numpy() + + if isinstance(original_r_tf, np.ndarray): + original_r = original_r_tf + else: + original_r = original_r_tf.numpy() + + target_loss = _generate_target_loss(rtbm, original_r, target) + best_parameters = copy.deepcopy(rtbm.get_parameters()) min_bound, max_bound = rtbm.get_bounds() loss_val = target_loss(best_parameters) @@ -50,7 +76,6 @@ def target_loss(params): n_parameters = len(best_parameters) # training hyperparameters - pop_per_rate = 256 mutation_rates = np.array([0.2, 0.4, 0.6, 0.8]) rates = np.concatenate([np.ones(pop_per_rate) * mr for mr in mutation_rates]) original_sigma = 0.25 @@ -86,7 +111,7 @@ def compute_mutant(mutation_rate): else: sigma /= 2 - if it % 50 == 0: + if it % 50 == 0 and verbose: current = time() logger.debug( "Iteration %d, best_loss: %.2f, (%2.fs)", @@ -112,7 +137,7 @@ def compute_mutant(mutation_rate): print(f"No more repeats allowed, iteration: {it}, loss: {loss_val:2.f}") rtbm.set_parameters(best_parameters) - return rtbm + return loss_val class RTBMFlow(MonteCarloFlow): @@ -124,37 +149,36 @@ def __init__(self, n_hidden=3, rtbm=None, train=True, *args, **kwargs): super().__init__(*args, **kwargs) self._train = train self._first_run = True + self._n_hidden = n_hidden + self._ga_generations = 3000 if rtbm is None: logger.info( - "Generating a RTBM with %d visible nodes and %d hidden" - % (self.n_dim, n_hidden) - ) - self._rtbm = RTBM( - self.n_dim, - n_hidden, - minimization_bound=50, - gaussian_init=True, - positive_T=True, - positive_Q=True, - gaussian_parameters={"mean": 0.0, "std": 0.75}, - sampling_activation="tanh", + "Generating a RTBM with %d visible nodes and %d hidden" % (self.n_dim, n_hidden) ) + # self._rtbm = RTBM( + # self.n_dim, + # n_hidden, + # minimization_bound=50, + # gaussian_init=True, + # positive_T=True, + # positive_Q=True, + # gaussian_parameters={"mean": 0.0, "std": 0.75}, + # sampling_activations="tanh", + # ) else: # Check whether it is a valid rtbm model if not hasattr(rtbm, "make_sample"): raise TypeError(f"{rtbm} is not a valid boltzman machine") self._rtbm = rtbm - self._p0 = self._rtbm.get_parameters() + self._first_run = False def freeze(self): """ Stop the training """ self.train = False - def unfreeze(self, reset_p0=False): + def unfreeze(self): """ Restart the training """ self.train = True - if reset_p0: - self._p0 = self._rtbm.get_parameters() def compile(self, integrand, compilable=False, **kwargs): if compilable: @@ -170,6 +194,10 @@ def generate_random_array(self, n_events): The return dimension of the random variables is (n_events,n_dim) and the jacobian (n_events) """ + if self._first_run: + rnds, _, xjac = super().generate_random_array(n_events) + return rnds, rnds, xjac + xrand, px, original_r = self._rtbm.make_sample_rho(n_events) # Since we are using the tanh function, the integration limits are (-1,1), move: xjac = 1.0 / px / n_events @@ -210,9 +238,15 @@ def _run_event(self, integrand, ncalls=None): def _run_iteration(self): all_res, unw, original_r = self.run_event() + generations = self._ga_generations + + if self.train and self._first_run: + original_r = self._run_first_run(unw, original_r) + self._first_run = False + generations = 250 if self.train: - _train_machine(self._rtbm, unw, original_r) + _train_machine(self._rtbm, unw, original_r, max_iterations=generations) res = tf.reduce_sum(all_res) all_res2 = all_res ** 2 @@ -223,6 +257,71 @@ def _run_iteration(self): sigma = tf.sqrt(tf.maximum(err_tmp2, fzero)) return res, sigma + def _run_first_run(self, unweighted_events, tf_rnds): + """ + Run the first iteration of the RTBM + """ + rnds = tf_rnds.numpy() + configurations_raw = [ + ("tanh", 0.0, 0.75), + ("sigmoid", 0.0, 1.5), +# ("softmax", 0.0, 0.75), + ] + names = ("name", "mean", "std") + + configuration_per_d = [] + for rnd in rnds.T: + vals, lims = np.histogram(rnd, bins=3, weights=unweighted_events, density=True) + mean = None + if vals[0] > 2 * vals[-1]: + mean = -0.5 + elif vals[-1] > 2 * vals[0]: + mean = 0.5 + tmp_list = copy.copy(configurations_raw) + if mean is not None: + tmp_list.append(("tanh", mean, 0.75)) +# tmp_list.append(("sigmoid", mean, 1.5)) + configuration_per_d.append([dict(zip(names, tup)) for tup in tmp_list]) + + # TODO: + # The number of combinations will intractably grow with the number of dimensions + # put some limits to this + all_configs = list(itertools.product(*configuration_per_d)) + logger.info(f"Generating {len(all_configs)} configurations") + + best_loss = 1e9 + best_params = None + winner_config = None + for configuration in all_configs: + acts, means, stds = zip(*[i.values() for i in configuration]) + rtbm = RTBM( + self.n_dim, + self._n_hidden, + minimization_bound=50, + gaussian_init=True, + positive_T=True, + positive_Q=True, + gaussian_parameters={"mean": means, "std": stds}, + sampling_activations=acts, + ) + guessed_original_r = rtbm.undo_transformation(rnds) + loss = _train_machine( + rtbm, + unweighted_events, + guessed_original_r, + max_iterations=10, + pop_per_rate=64, + verbose=False, + ) + if loss < best_loss: + best_loss = loss + best_params = rtbm.get_parameters() + winner_config = configuration + self._rtbm = rtbm + rtbm.set_parameters(best_params) + logger.info(f"And the winner is: {configuration}") + return rtbm.undo_transformation(rnds) + def rtbm_wrapper(integrand, n_dim, n_iter, total_n_events, **kwargs): """Convenience wrapper From 62b0e4c0ea8c738731c0011cd0a1ae2494ad134f Mon Sep 17 00:00:00 2001 From: juacrumar Date: Mon, 19 Apr 2021 12:28:50 +0200 Subject: [PATCH 17/19] add timeouts --- examples/sin_tf.py | 4 ++-- src/vegasflow/rtbm.py | 50 ++++++++++++++++++++++++++++++++----------- 2 files changed, 39 insertions(+), 15 deletions(-) diff --git a/examples/sin_tf.py b/examples/sin_tf.py index e5b728f..bafecd9 100644 --- a/examples/sin_tf.py +++ b/examples/sin_tf.py @@ -41,7 +41,7 @@ def sin_fun(xarr, **kwargs): print(f"RTBM MC, ncalls={ncalls}:") start = time.time() - rtbm = RTBMFlow(n_dim=dim, n_events=ncalls, train=True, n_hidden=2) + rtbm = RTBMFlow(n_dim=dim, n_events=ncalls, train=True, n_hidden=3) rtbm.compile(integrand) _ = rtbm.run_integration(n_iter) rtbm.freeze() @@ -49,7 +49,7 @@ def sin_fun(xarr, **kwargs): print(f"RTBM took: time (s): {end-start}") print("Results with frozen grids") - vegas_instance.run_integration(5) + r = vegas_instance.run_integration(5) rt = rtbm.run_integration(5) print(f"Result computed by Vegas: {r[0]:.5f} +- {r[1]:.6f}") print(f"Result computed by RTBM: {rt[0]:.5f} +- {rt[1]:.6f}") diff --git a/src/vegasflow/rtbm.py b/src/vegasflow/rtbm.py index 1979b8d..222529e 100644 --- a/src/vegasflow/rtbm.py +++ b/src/vegasflow/rtbm.py @@ -23,10 +23,11 @@ from time import time logger = logging.getLogger(__name__) +TOL = 1e-9 # Cost functions for training def _kl(x, ytarget): - return ytarget * np.log(ytarget / x) + return ytarget * np.log(ytarget / (x+TOL)) _loss = _kl @@ -50,6 +51,10 @@ def _train_machine( max_iterations=3000, pop_per_rate=256, verbose=True, + resets=True, + timeout=5*60, # dont wait more than 5 minutes per iteration + fail_ontimeout=False, + skip_on_nan=True, ): if verbose: @@ -71,7 +76,7 @@ def _train_machine( min_bound, max_bound = rtbm.get_bounds() loss_val = target_loss(best_parameters) - with Parallel(n_jobs=n_jobs) as parallel: + with Parallel(n_jobs=n_jobs, timeout=timeout) as parallel: prev = time() n_parameters = len(best_parameters) @@ -101,10 +106,18 @@ def compute_mutant(mutation_rate): return target_loss(mutant), mutant parallel_runs = [delayed(compute_mutant)(rate) for rate in rates] - result = parallel(parallel_runs) + try: + result = parallel(parallel_runs) + except: # TODO control better + logger.debug("Time'd out, skip me") + if fail_ontimeout: + raise ValueError + result = [(np.nan, None)] losses, mutants = zip(*result) - best_loss = np.nanmin(losses) + + # Insert the last best to avoid runtime errors + best_loss = np.nanmin(list(losses) + [loss_val+1.0]) if best_loss < loss_val: loss_val = best_loss best_parameters = mutants[losses.index(best_loss)] @@ -130,7 +143,10 @@ def compute_mutant(mutation_rate): if sigma < 1e-3: sigma = original_sigma logger.debug("Resetting sigma with loss: %.2f", loss_val) - repeats -= 1 + if resets: + repeats -= 1 + else: + repeats = 0 break if not repeats: @@ -292,7 +308,9 @@ def _run_first_run(self, unweighted_events, tf_rnds): best_loss = 1e9 best_params = None winner_config = None + for configuration in all_configs: + logger.info(f"Testing {configuration}") acts, means, stds = zip(*[i.values() for i in configuration]) rtbm = RTBM( self.n_dim, @@ -305,14 +323,20 @@ def _run_first_run(self, unweighted_events, tf_rnds): sampling_activations=acts, ) guessed_original_r = rtbm.undo_transformation(rnds) - loss = _train_machine( - rtbm, - unweighted_events, - guessed_original_r, - max_iterations=10, - pop_per_rate=64, - verbose=False, - ) + try: + loss = _train_machine( + rtbm, + unweighted_events, + guessed_original_r, + max_iterations=32, + pop_per_rate=32, + resets=False, + verbose=False, + fail_ontimeout=True, + timeout=10 # TODO if any initial iteration takes more than 10 seconds, get out + ) + except ValueError: + loss = 1e9 if loss < best_loss: best_loss = loss best_params = rtbm.get_parameters() From 4daa5edb38568fd4ebe54f4357e404b0dca89181 Mon Sep 17 00:00:00 2001 From: juacrumar Date: Fri, 2 Jul 2021 11:51:13 +0200 Subject: [PATCH 18/19] push my last version --- examples/sin_tf.py | 20 +++++++++---- src/vegasflow/rtbm.py | 69 ++++++++++++++++++++++++++++--------------- 2 files changed, 60 insertions(+), 29 deletions(-) diff --git a/examples/sin_tf.py b/examples/sin_tf.py index bafecd9..411d034 100644 --- a/examples/sin_tf.py +++ b/examples/sin_tf.py @@ -15,19 +15,21 @@ # MC integration setup dim = 5 -ncalls = np.int32(1e4) +ncalls = np.int32(5e4) n_iter = 5 +n_hidden = 2 tf_pi = float_me(np.pi) +npeaks = 2.0 @tf.function def sin_fun(xarr, **kwargs): """symgauss test function""" - res = tf.pow(tf.sin(xarr*tf_pi),2) + res = tf.pow(tf.sin(npeaks*xarr*tf_pi),2) return tf.reduce_prod(res, axis=1) -# integrand = sin_fun -from simgauss_tf import symgauss as integrand +integrand = sin_fun +#from simgauss_tf import symgauss as integrand if __name__ == "__main__": """Testing several different integrations""" @@ -41,13 +43,18 @@ def sin_fun(xarr, **kwargs): print(f"RTBM MC, ncalls={ncalls}:") start = time.time() - rtbm = RTBMFlow(n_dim=dim, n_events=ncalls, train=True, n_hidden=3) + rtbm = RTBMFlow(n_dim=dim, n_events=ncalls, train=True, n_hidden=n_hidden) rtbm.compile(integrand) _ = rtbm.run_integration(n_iter) rtbm.freeze() end = time.time() print(f"RTBM took: time (s): {end-start}") + parameters = rtbm._rtbm.get_parameters() + param_file = f"PARAMS_for_ndim={dim}_{npeaks}xpeaks_hidden={n_hidden}.npy" + np.save(param_file, parameters) + print(f" > Saved parameters to {param_file}") + print("Results with frozen grids") r = vegas_instance.run_integration(5) rt = rtbm.run_integration(5) @@ -57,3 +64,6 @@ def sin_fun(xarr, **kwargs): # Notes # For 1 and 2 dimensions is enough with n_hidden=1 to get a better per event error # For 3 dimensions we need to go to n_hidden=2 +# We might be actually overfitting big time because after a few iterations the integration stops being so good +# a good way of testing this would be plotting the distribution of points one gets after each training! +# or fitting wrongly anyway diff --git a/src/vegasflow/rtbm.py b/src/vegasflow/rtbm.py index 222529e..5236a71 100644 --- a/src/vegasflow/rtbm.py +++ b/src/vegasflow/rtbm.py @@ -23,21 +23,32 @@ from time import time logger = logging.getLogger(__name__) -TOL = 1e-9 +TOL = 1e-6 # Cost functions for training -def _kl(x, ytarget): - return ytarget * np.log(ytarget / (x+TOL)) +def _kl(x, ytarget_raw): + ytarget = ( ytarget_raw + TOL ) + ytarget /= np.sum(ytarget) + x /= np.sum(x) -_loss = _kl + return ytarget * np.log(ytarget / x) + +def _mse(x, y): + integral = np.sum(y) + x /= np.sum(x) + y /= integral + return integral*pow(x-y,2) +_loss = _kl def _generate_target_loss(rtbm, original_r, target): def target_loss(params): if not rtbm.set_parameters(params): return np.NaN _, prob = rtbm.get_transformation(original_r) + if (prob <= 0.0).any(): + return np.NaN return np.sum(_loss(prob, target)) return target_loss @@ -49,14 +60,19 @@ def _train_machine( original_r_tf, n_jobs=32, max_iterations=3000, - pop_per_rate=256, + pop_per_rate=512, verbose=True, resets=True, timeout=5*60, # dont wait more than 5 minutes per iteration + timeout_repeat=3, fail_ontimeout=False, skip_on_nan=True, ): + # note that if between 2 timeout there are severaliterations + # but they do not produce a better esult it doesn't count + timeout_original = timeout_repeat + if verbose: logger.info("Training RTBM") @@ -81,9 +97,9 @@ def _train_machine( n_parameters = len(best_parameters) # training hyperparameters - mutation_rates = np.array([0.2, 0.4, 0.6, 0.8]) + mutation_rates = np.array([0.1, 0.2, 0.4, 0.6, 0.8, 0.96]) rates = np.concatenate([np.ones(pop_per_rate) * mr for mr in mutation_rates]) - original_sigma = 0.25 + original_sigma = 0.55 sigma = original_sigma repeats = 3 ####### @@ -110,8 +126,12 @@ def compute_mutant(mutation_rate): result = parallel(parallel_runs) except: # TODO control better logger.debug("Time'd out, skip me") + timeout_repeat -= 1 if fail_ontimeout: raise ValueError + if timeout_repeat == 0: + logger.debug("Full timeout, out") + break result = [(np.nan, None)] losses, mutants = zip(*result) @@ -119,15 +139,16 @@ def compute_mutant(mutation_rate): # Insert the last best to avoid runtime errors best_loss = np.nanmin(list(losses) + [loss_val+1.0]) if best_loss < loss_val: + timeout_repeat = timeout_original # reset the timeouts loss_val = best_loss best_parameters = mutants[losses.index(best_loss)] else: - sigma /= 2 + sigma *= 0.9 if it % 50 == 0 and verbose: current = time() logger.debug( - "Iteration %d, best_loss: %.2f, (%2.fs)", + "Iteration %d, best_loss: %.4f, (%2.fs)", it, loss_val, current - prev, @@ -140,17 +161,17 @@ def compute_mutant(mutation_rate): ) prev = current - if sigma < 1e-3: + if sigma < 1e-2: sigma = original_sigma - logger.debug("Resetting sigma with loss: %.2f", loss_val) + logger.debug("Resetting sigma with loss: %.4f after %d iterations", loss_val, it) if resets: repeats -= 1 else: repeats = 0 - break if not repeats: - print(f"No more repeats allowed, iteration: {it}, loss: {loss_val:2.f}") + print(f"No more repeats allowed, iteration: {it}, loss: {loss_val:.4f}") + break rtbm.set_parameters(best_parameters) return loss_val @@ -164,23 +185,23 @@ class RTBMFlow(MonteCarloFlow): def __init__(self, n_hidden=3, rtbm=None, train=True, *args, **kwargs): super().__init__(*args, **kwargs) self._train = train - self._first_run = True + self._first_run = False self._n_hidden = n_hidden self._ga_generations = 3000 if rtbm is None: logger.info( "Generating a RTBM with %d visible nodes and %d hidden" % (self.n_dim, n_hidden) ) - # self._rtbm = RTBM( - # self.n_dim, - # n_hidden, - # minimization_bound=50, - # gaussian_init=True, - # positive_T=True, - # positive_Q=True, - # gaussian_parameters={"mean": 0.0, "std": 0.75}, - # sampling_activations="tanh", - # ) + self._rtbm = RTBM( + self.n_dim, + n_hidden, + minimization_bound=50, + gaussian_init=True, + positive_T=True, + positive_Q=True, + gaussian_parameters={"mean": 0.0, "std": 0.75}, + sampling_activations="tanh", + ) else: # Check whether it is a valid rtbm model if not hasattr(rtbm, "make_sample"): From 3206cf480eebcda76c1f1adc29076a93f1b8af81 Mon Sep 17 00:00:00 2001 From: juacrumar Date: Fri, 2 Jul 2021 13:17:47 +0200 Subject: [PATCH 19/19] seems to work --- examples/sin_tf.py | 15 ++++++++------- src/vegasflow/rtbm.py | 30 ++++++++++++++++-------------- 2 files changed, 24 insertions(+), 21 deletions(-) diff --git a/examples/sin_tf.py b/examples/sin_tf.py index 411d034..9e589fc 100644 --- a/examples/sin_tf.py +++ b/examples/sin_tf.py @@ -14,10 +14,10 @@ # MC integration setup -dim = 5 -ncalls = np.int32(5e4) +dim = 2 +ncalls = np.int32(1e2) n_iter = 5 -n_hidden = 2 +n_hidden = 1 tf_pi = float_me(np.pi) npeaks = 2.0 @@ -25,14 +25,15 @@ @tf.function def sin_fun(xarr, **kwargs): """symgauss test function""" - res = tf.pow(tf.sin(npeaks*xarr*tf_pi),2) + res = tf.pow(tf.sin(npeaks * xarr * tf_pi), 2) return tf.reduce_prod(res, axis=1) + integrand = sin_fun -#from simgauss_tf import symgauss as integrand +# from simgauss_tf import symgauss as integrand if __name__ == "__main__": - """Testing several different integrations""" + # Testing several different integrations print(f"VEGAS MC, ncalls={ncalls}:") start = time.time() ncalls = ncalls @@ -43,7 +44,7 @@ def sin_fun(xarr, **kwargs): print(f"RTBM MC, ncalls={ncalls}:") start = time.time() - rtbm = RTBMFlow(n_dim=dim, n_events=ncalls, train=True, n_hidden=n_hidden) + rtbm = RTBMFlow(n_dim=dim, n_events=ncalls, train=True, n_hidden=n_hidden, generations=500) rtbm.compile(integrand) _ = rtbm.run_integration(n_iter) rtbm.freeze() diff --git a/src/vegasflow/rtbm.py b/src/vegasflow/rtbm.py index 5236a71..7e4e6b4 100644 --- a/src/vegasflow/rtbm.py +++ b/src/vegasflow/rtbm.py @@ -27,21 +27,24 @@ # Cost functions for training def _kl(x, ytarget_raw): - ytarget = ( ytarget_raw + TOL ) + ytarget = ytarget_raw + TOL ytarget /= np.sum(ytarget) x /= np.sum(x) return ytarget * np.log(ytarget / x) + def _mse(x, y): integral = np.sum(y) x /= np.sum(x) y /= integral - return integral*pow(x-y,2) + return integral * pow(x - y, 2) + _loss = _kl + def _generate_target_loss(rtbm, original_r, target): def target_loss(params): if not rtbm.set_parameters(params): @@ -63,7 +66,7 @@ def _train_machine( pop_per_rate=512, verbose=True, resets=True, - timeout=5*60, # dont wait more than 5 minutes per iteration + timeout=5 * 60, # dont wait more than 5 minutes per iteration timeout_repeat=3, fail_ontimeout=False, skip_on_nan=True, @@ -124,7 +127,7 @@ def compute_mutant(mutation_rate): parallel_runs = [delayed(compute_mutant)(rate) for rate in rates] try: result = parallel(parallel_runs) - except: # TODO control better + except: # TODO control better logger.debug("Time'd out, skip me") timeout_repeat -= 1 if fail_ontimeout: @@ -135,11 +138,10 @@ def compute_mutant(mutation_rate): result = [(np.nan, None)] losses, mutants = zip(*result) - # Insert the last best to avoid runtime errors - best_loss = np.nanmin(list(losses) + [loss_val+1.0]) + best_loss = np.nanmin(list(losses) + [loss_val + 1.0]) if best_loss < loss_val: - timeout_repeat = timeout_original # reset the timeouts + timeout_repeat = timeout_original # reset the timeouts loss_val = best_loss best_parameters = mutants[losses.index(best_loss)] else: @@ -182,12 +184,12 @@ class RTBMFlow(MonteCarloFlow): RTBM based Monte Carlo integrator """ - def __init__(self, n_hidden=3, rtbm=None, train=True, *args, **kwargs): + def __init__(self, n_hidden=3, rtbm=None, train=True, generations=3000, *args, **kwargs): super().__init__(*args, **kwargs) self._train = train self._first_run = False self._n_hidden = n_hidden - self._ga_generations = 3000 + self._ga_generations = generations if rtbm is None: logger.info( "Generating a RTBM with %d visible nodes and %d hidden" % (self.n_dim, n_hidden) @@ -210,11 +212,11 @@ def __init__(self, n_hidden=3, rtbm=None, train=True, *args, **kwargs): self._first_run = False def freeze(self): - """ Stop the training """ + """Stop the training""" self.train = False def unfreeze(self): - """ Restart the training """ + """Restart the training""" self.train = True def compile(self, integrand, compilable=False, **kwargs): @@ -302,7 +304,7 @@ def _run_first_run(self, unweighted_events, tf_rnds): configurations_raw = [ ("tanh", 0.0, 0.75), ("sigmoid", 0.0, 1.5), -# ("softmax", 0.0, 0.75), + # ("softmax", 0.0, 0.75), ] names = ("name", "mean", "std") @@ -317,7 +319,7 @@ def _run_first_run(self, unweighted_events, tf_rnds): tmp_list = copy.copy(configurations_raw) if mean is not None: tmp_list.append(("tanh", mean, 0.75)) -# tmp_list.append(("sigmoid", mean, 1.5)) + # tmp_list.append(("sigmoid", mean, 1.5)) configuration_per_d.append([dict(zip(names, tup)) for tup in tmp_list]) # TODO: @@ -354,7 +356,7 @@ def _run_first_run(self, unweighted_events, tf_rnds): resets=False, verbose=False, fail_ontimeout=True, - timeout=10 # TODO if any initial iteration takes more than 10 seconds, get out + timeout=10, # TODO if any initial iteration takes more than 10 seconds, get out ) except ValueError: loss = 1e9