From c4db40aa21736aa9fda0b13562ef5e06b36efc84 Mon Sep 17 00:00:00 2001 From: Piotr Date: Mon, 30 Oct 2023 12:37:43 +0100 Subject: [PATCH 01/10] Backend as strict class + minor changes `formulae` or `Storage` should never be overridden --- PySDM/backends/impl_common/backend_methods.py | 2 ++ .../impl_numba/methods/collisions_methods.py | 19 ++++++++++------- .../impl_numba/methods/physics_methods.py | 21 ++++++++++++------- 3 files changed, 26 insertions(+), 16 deletions(-) diff --git a/PySDM/backends/impl_common/backend_methods.py b/PySDM/backends/impl_common/backend_methods.py index 6ab0824f3..d3a1cd2c4 100644 --- a/PySDM/backends/impl_common/backend_methods.py +++ b/PySDM/backends/impl_common/backend_methods.py @@ -3,8 +3,10 @@ class for all backend methods classes """ +from pystrict import strict # pylint: disable=too-few-public-methods +@strict class BackendMethods: def __init__(self): if not hasattr(self, "formulae"): diff --git a/PySDM/backends/impl_numba/methods/collisions_methods.py b/PySDM/backends/impl_numba/methods/collisions_methods.py index aba8bd698..3f3805514 100644 --- a/PySDM/backends/impl_numba/methods/collisions_methods.py +++ b/PySDM/backends/impl_numba/methods/collisions_methods.py @@ -2,6 +2,8 @@ CPU implementation of backend methods for particle collisions """ # pylint: disable=too-many-lines +from functools import cached_property + import numba import numpy as np @@ -374,14 +376,6 @@ def __collision_coalescence_breakup_body( self.__collision_coalescence_breakup_body = __collision_coalescence_breakup_body - @numba.njit(**{**conf.JIT_FLAGS, "fastmath": self.formulae.fastmath}) - def __ll82_coalescence_check_body(*, Ec, dl): - for i in numba.prange(len(Ec)): # pylint: disable=not-an-iterable - if dl[i] < 0.4e-3: - Ec[i] = 1.0 - - self.__ll82_coalescence_check_body = __ll82_coalescence_check_body - if self.formulae.fragmentation_function.__name__ == "Straub2010Nf": straub_sigma1 = self.formulae.fragmentation_function.params_sigma1 straub_mu1 = self.formulae.fragmentation_function.params_mu1 @@ -969,6 +963,15 @@ def ll82_fragmentation( nfmax=nfmax, ) + @cached_property + def __ll82_coalescence_check_body(self): + @numba.njit(**{**conf.JIT_FLAGS, "fastmath": self.formulae.fastmath}) + def __ll82_coalescence_check_body(*, Ec, dl): + for i in numba.prange(len(Ec)): # pylint: disable=not-an-iterable + if dl[i] < 0.4e-3: + Ec[i] = 1.0 + return __ll82_coalescence_check_body + def ll82_coalescence_check(self, *, Ec, dl): self.__ll82_coalescence_check_body( Ec=Ec.data, diff --git a/PySDM/backends/impl_numba/methods/physics_methods.py b/PySDM/backends/impl_numba/methods/physics_methods.py index 8e886cc3b..e408a1365 100644 --- a/PySDM/backends/impl_numba/methods/physics_methods.py +++ b/PySDM/backends/impl_numba/methods/physics_methods.py @@ -1,6 +1,8 @@ """ CPU implementation of backend methods wrapping basic physics formulae """ +from functools import cached_property + import numba from numba import prange @@ -21,7 +23,6 @@ def __init__(self): # pylint: disable=too-many-locals phys_volume = self.formulae.trivia.volume phys_r_cr = self.formulae.hygroscopicity.r_cr phys_mass_to_volume = self.formulae.particle_shape_and_density.mass_to_volume - phys_volume_to_mass = self.formulae.particle_shape_and_density.volume_to_mass const = self.formulae.constants @numba.njit(**{**conf.JIT_FLAGS, "fastmath": self.formulae.fastmath}) @@ -77,13 +78,6 @@ def volume_of_mass(volume, mass): self.volume_of_mass = volume_of_mass - @numba.njit(**{**conf.JIT_FLAGS, "fastmath": self.formulae.fastmath}) - def mass_of_volume(mass, volume): - for i in prange(volume.shape[0]): # pylint: disable=not-an-iterable - mass[i] = phys_volume_to_mass(volume[i]) - - self.mass_of_volume = mass_of_volume - def temperature_pressure_RH( self, *, rhod, thd, water_vapour_mixing_ratio, T, p, RH ): @@ -122,5 +116,16 @@ def a_w_ice(self, *, T, p, RH, water_vapour_mixing_ratio, a_w_ice): def volume_of_water_mass(self, volume, mass): self.volume_of_mass(volume.data, mass.data) + @cached_property + def mass_of_volume(self): + phys_volume_to_mass = self.formulae.particle_shape_and_density.volume_to_mass + + @numba.njit(**{**conf.JIT_FLAGS, "fastmath": self.formulae.fastmath}) + def mass_of_volume(mass, volume): + for i in prange(volume.shape[0]): # pylint: disable=not-an-iterable + mass[i] = phys_volume_to_mass(volume[i]) + + return mass_of_volume + def mass_of_water_volume(self, mass, volume): self.mass_of_volume(mass.data, volume.data) From d08651526681b7389f96b7c8c712bd5730ce4868 Mon Sep 17 00:00:00 2001 From: Piotr Date: Mon, 6 Nov 2023 12:13:09 +0100 Subject: [PATCH 02/10] CPU Time instead of wall time --- examples/PySDM_examples/Shima_et_al_2009/example.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/examples/PySDM_examples/Shima_et_al_2009/example.py b/examples/PySDM_examples/Shima_et_al_2009/example.py index 315e675ce..e486b8def 100644 --- a/examples/PySDM_examples/Shima_et_al_2009/example.py +++ b/examples/PySDM_examples/Shima_et_al_2009/example.py @@ -10,7 +10,7 @@ from PySDM.dynamics import Coalescence from PySDM.environments import Box from PySDM.initialisation.sampling.spectral_sampling import ConstantMultiplicity -from PySDM.products import ParticleVolumeVersusRadiusLogarithmSpectrum, WallTime +from PySDM.products import ParticleVolumeVersusRadiusLogarithmSpectrum, CPUTime def run(settings, backend=CPU, observers=()): @@ -27,7 +27,7 @@ def run(settings, backend=CPU, observers=()): ParticleVolumeVersusRadiusLogarithmSpectrum( settings.radius_bins_edges, name="dv/dlnr" ), - WallTime(), + CPUTime(), ) particulator = builder.build(attributes, products) if hasattr(settings, "u_term") and "terminal velocity" in particulator.attributes: @@ -39,13 +39,13 @@ def run(settings, backend=CPU, observers=()): particulator.observers.append(observer) vals = {} - particulator.products["wall time"].reset() + particulator.products["CPU Time"].reset() for step in settings.output_steps: particulator.run(step - particulator.n_steps) vals[step] = particulator.products["dv/dlnr"].get()[0] vals[step][:] *= settings.rho - exec_time = particulator.products["wall time"].get() + exec_time = particulator.products["CPU Time"].get() return vals, exec_time From 523dc5a337ffed98afb806e5dafb6c25f8448658 Mon Sep 17 00:00:00 2001 From: Piotr Date: Mon, 6 Nov 2023 12:13:52 +0100 Subject: [PATCH 03/10] pass seed to Formulae --- examples/PySDM_examples/Shima_et_al_2009/settings.py | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/examples/PySDM_examples/Shima_et_al_2009/settings.py b/examples/PySDM_examples/Shima_et_al_2009/settings.py index 5e0d524ba..3764e797e 100644 --- a/examples/PySDM_examples/Shima_et_al_2009/settings.py +++ b/examples/PySDM_examples/Shima_et_al_2009/settings.py @@ -11,9 +11,9 @@ @strict class Settings: - def __init__(self, steps: Optional[list] = None): + def __init__(self, seed: int, steps: Optional[list] = None): steps = steps or [0, 1200, 2400, 3600] - self.formulae = Formulae() + self.formulae = Formulae(seed=seed) self.n_sd = 2**13 self.n_part = 2**23 / si.metre**3 self.X0 = self.formulae.trivia.volume(radius=30.531 * si.micrometres) @@ -22,7 +22,6 @@ def __init__(self, steps: Optional[list] = None): self.rho = 1000 * si.kilogram / si.metre**3 self.dt = 1 * si.seconds self.adaptive = False - self.seed = 44 self.steps = steps self.kernel = Golovin(b=1.5e3 / si.second) self.spectrum = spectra.Exponential(norm_factor=self.norm_factor, scale=self.X0) From 633a3dd1441dd0465c728c298b90862783423786 Mon Sep 17 00:00:00 2001 From: Piotr Date: Mon, 6 Nov 2023 12:14:48 +0100 Subject: [PATCH 04/10] simpler title --- examples/PySDM_examples/Shima_et_al_2009/spectrum_plotter.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/examples/PySDM_examples/Shima_et_al_2009/spectrum_plotter.py b/examples/PySDM_examples/Shima_et_al_2009/spectrum_plotter.py index 20285f0a8..9a2710d7d 100644 --- a/examples/PySDM_examples/Shima_et_al_2009/spectrum_plotter.py +++ b/examples/PySDM_examples/Shima_et_al_2009/spectrum_plotter.py @@ -123,7 +123,8 @@ def analytic_solution(x): if spectrum is not None: y = spectrum * si.kilograms / si.grams error = error_measure(y, y_true, x) - self.title = f"error measure: {error:.2f}" # TODO #327 relative error + self.title = f"error: {error:.2f}" # TODO #327 relative error + self.title = f"error: {error:.2f}" # TODO #327 relative error return error return None From d29a20f2e8c1effdf9507ad5c5d154cc92551605 Mon Sep 17 00:00:00 2001 From: Piotr Date: Mon, 6 Nov 2023 12:15:08 +0100 Subject: [PATCH 05/10] WIP: simpler error_measure --- examples/PySDM_examples/Shima_et_al_2009/error_measure.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/examples/PySDM_examples/Shima_et_al_2009/error_measure.py b/examples/PySDM_examples/Shima_et_al_2009/error_measure.py index 5097f09be..e230bf0a4 100644 --- a/examples/PySDM_examples/Shima_et_al_2009/error_measure.py +++ b/examples/PySDM_examples/Shima_et_al_2009/error_measure.py @@ -1,7 +1,7 @@ import numpy as np -def error_measure(y, y_true, x): +def error_measure_old(y, y_true, x): errors = y_true - y errors = errors[0:-1] + errors[1:] dx = np.diff(x) @@ -9,3 +9,6 @@ def error_measure(y, y_true, x): errors /= 2 error = np.sum(np.abs(errors)) return error + +def error_measure(y, y_true, _): + return np.sqrt(np.mean(np.square(y - y_true))) \ No newline at end of file From 4bfb0d14f1f1e42f721de5b722d574d6a626308d Mon Sep 17 00:00:00 2001 From: Piotr Date: Mon, 6 Nov 2023 12:15:48 +0100 Subject: [PATCH 06/10] WIP --- .../fig_4_adaptive_sdm.py | 32 ++++++++++++------- 1 file changed, 20 insertions(+), 12 deletions(-) diff --git a/examples/PySDM_examples/Bartman_2020_MasterThesis/fig_4_adaptive_sdm.py b/examples/PySDM_examples/Bartman_2020_MasterThesis/fig_4_adaptive_sdm.py index a0fef1492..071fbb889 100644 --- a/examples/PySDM_examples/Bartman_2020_MasterThesis/fig_4_adaptive_sdm.py +++ b/examples/PySDM_examples/Bartman_2020_MasterThesis/fig_4_adaptive_sdm.py @@ -1,4 +1,5 @@ import os +import numpy as np from matplotlib import pyplot as plt from PySDM_examples.Shima_et_al_2009.example import run @@ -7,10 +8,11 @@ def main(plot: bool = True, save: str = None): - n_sds = [13, 15, 17] - dts = [10, 5, 1, "adaptive"] - iters = 10 + n_sds = [13, 15, 17, 19] + dts = [10, 1, "adaptive"] + iters_without_warmup = 5 base_time = None + base_error = None plt.ioff() fig, axs = plt.subplots( @@ -20,22 +22,25 @@ def main(plot: bool = True, save: str = None): for i, dt in enumerate(dts): for j, n_sd in enumerate(n_sds): outputs = [] - exec_time = 0 - for _ in range(iters): - settings = Settings() + exec_times = [] + one_for_warmup = 1 + for it in range(iters_without_warmup + one_for_warmup): + settings = Settings(seed=it) settings.n_sd = 2**n_sd - settings.dt = dt if dt != "adaptive" else 10 + settings.dt = dt if dt != "adaptive" else max(dts[:-1]) settings.adaptive = dt == "adaptive" states, exec_time = run(settings) + print(f"{dt=}, {n_sd=}, {exec_time=}, {it=}") + exec_times.append(exec_time) outputs.append(states) - mean_time = exec_time / iters + mean_time = np.mean(exec_times[one_for_warmup:]) if base_time is None: base_time = mean_time norm_time = mean_time / base_time mean_output = {} - for key in outputs[0].keys(): + for key in outputs[0].keys (): mean_output[key] = sum((output[key] for output in outputs)) / len( outputs ) @@ -45,10 +50,13 @@ def main(plot: bool = True, save: str = None): plotter.ax = axs[i, j] plotter.smooth = True for step, vals in mean_output.items(): - plotter.plot(vals, step * settings.dt) + error = plotter.plot(vals, step * settings.dt) + last_step_error = error + if base_error is None: + base_error = last_step_error plotter.ylabel = ( - r"$\bf{dt: " + str(dt) + "}$\ndm/dlnr [g/m^3/(unit dr/r)]" + r"$\bf{dt: " + str(settings.dt) + ("+ adaptivity" if settings.adaptive else "") + "}$\ndm/dlnr [g/m^3/(unit dr/r)]" if j == 0 else None ) @@ -57,7 +65,7 @@ def main(plot: bool = True, save: str = None): if i == len(dts) - 1 else None ) - plotter.title = f"norm. time: {norm_time:.2f}; " + plotter.title + plotter.title = f"norm. time: {norm_time:.2f}; error: {last_step_error / base_error:.2f}" plotter.finished = False plotter.finish() if save is not None: From 81a92f4b838170c59cd53e3837d424454d0922a3 Mon Sep 17 00:00:00 2001 From: emmacware Date: Thu, 10 Oct 2024 17:51:17 +0200 Subject: [PATCH 07/10] test --- examples/PySDM_examples/Shima_et_al_2009/settings.py | 1 - 1 file changed, 1 deletion(-) diff --git a/examples/PySDM_examples/Shima_et_al_2009/settings.py b/examples/PySDM_examples/Shima_et_al_2009/settings.py index 3764e797e..ffc0fc4bb 100644 --- a/examples/PySDM_examples/Shima_et_al_2009/settings.py +++ b/examples/PySDM_examples/Shima_et_al_2009/settings.py @@ -8,7 +8,6 @@ from PySDM.initialisation import spectra from PySDM.physics import si - @strict class Settings: def __init__(self, seed: int, steps: Optional[list] = None): From 9d164a3c2d3fea467689d3bc58caa8cd97931b88 Mon Sep 17 00:00:00 2001 From: Sylwester Arabas Date: Thu, 10 Oct 2024 22:39:01 +0200 Subject: [PATCH 08/10] address pylint hints --- PySDM/backends/impl_numba/methods/collisions_methods.py | 1 - .../Bartman_2020_MasterThesis/fig_4_adaptive_sdm.py | 9 ++++++--- .../PySDM_examples/Shima_et_al_2009/error_measure.py | 3 ++- examples/PySDM_examples/Shima_et_al_2009/settings.py | 3 ++- 4 files changed, 10 insertions(+), 6 deletions(-) diff --git a/PySDM/backends/impl_numba/methods/collisions_methods.py b/PySDM/backends/impl_numba/methods/collisions_methods.py index f3d36008d..06573bfc6 100644 --- a/PySDM/backends/impl_numba/methods/collisions_methods.py +++ b/PySDM/backends/impl_numba/methods/collisions_methods.py @@ -515,7 +515,6 @@ def collision_coalescence_breakup( def _compute_gamma_body(self): @numba.njit(**self.default_jit_flags) # pylint: disable=too-many-arguments,too-many-locals - def body( prob, rand, diff --git a/examples/PySDM_examples/Bartman_2020_MasterThesis/fig_4_adaptive_sdm.py b/examples/PySDM_examples/Bartman_2020_MasterThesis/fig_4_adaptive_sdm.py index 071fbb889..88c07a57f 100644 --- a/examples/PySDM_examples/Bartman_2020_MasterThesis/fig_4_adaptive_sdm.py +++ b/examples/PySDM_examples/Bartman_2020_MasterThesis/fig_4_adaptive_sdm.py @@ -40,7 +40,7 @@ def main(plot: bool = True, save: str = None): base_time = mean_time norm_time = mean_time / base_time mean_output = {} - for key in outputs[0].keys (): + for key in outputs[0].keys(): mean_output[key] = sum((output[key] for output in outputs)) / len( outputs ) @@ -56,7 +56,10 @@ def main(plot: bool = True, save: str = None): base_error = last_step_error plotter.ylabel = ( - r"$\bf{dt: " + str(settings.dt) + ("+ adaptivity" if settings.adaptive else "") + "}$\ndm/dlnr [g/m^3/(unit dr/r)]" + r"$\bf{dt: " + + str(settings.dt) + + ("+ adaptivity" if settings.adaptive else "") + + "}$\ndm/dlnr [g/m^3/(unit dr/r)]" if j == 0 else None ) @@ -65,7 +68,7 @@ def main(plot: bool = True, save: str = None): if i == len(dts) - 1 else None ) - plotter.title = f"norm. time: {norm_time:.2f}; error: {last_step_error / base_error:.2f}" + plotter.title = f"time: {norm_time:.2f} error: {last_step_error / base_error:.2f} (normalised)" plotter.finished = False plotter.finish() if save is not None: diff --git a/examples/PySDM_examples/Shima_et_al_2009/error_measure.py b/examples/PySDM_examples/Shima_et_al_2009/error_measure.py index e230bf0a4..e297d7d35 100644 --- a/examples/PySDM_examples/Shima_et_al_2009/error_measure.py +++ b/examples/PySDM_examples/Shima_et_al_2009/error_measure.py @@ -10,5 +10,6 @@ def error_measure_old(y, y_true, x): error = np.sum(np.abs(errors)) return error + def error_measure(y, y_true, _): - return np.sqrt(np.mean(np.square(y - y_true))) \ No newline at end of file + return np.sqrt(np.mean(np.square(y - y_true))) diff --git a/examples/PySDM_examples/Shima_et_al_2009/settings.py b/examples/PySDM_examples/Shima_et_al_2009/settings.py index ffc0fc4bb..dd586ae60 100644 --- a/examples/PySDM_examples/Shima_et_al_2009/settings.py +++ b/examples/PySDM_examples/Shima_et_al_2009/settings.py @@ -8,9 +8,10 @@ from PySDM.initialisation import spectra from PySDM.physics import si + @strict class Settings: - def __init__(self, seed: int, steps: Optional[list] = None): + def __init__(self, seed: Optional[int] = None, steps: Optional[list] = None): steps = steps or [0, 1200, 2400, 3600] self.formulae = Formulae(seed=seed) self.n_sd = 2**13 From d04064c0889efa4daf7a1bbf848a7cff5d0cad2a Mon Sep 17 00:00:00 2001 From: Sylwester Arabas Date: Thu, 10 Oct 2024 22:40:49 +0200 Subject: [PATCH 09/10] address pylint hints --- .../Bartman_2020_MasterThesis/fig_4_adaptive_sdm.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/examples/PySDM_examples/Bartman_2020_MasterThesis/fig_4_adaptive_sdm.py b/examples/PySDM_examples/Bartman_2020_MasterThesis/fig_4_adaptive_sdm.py index 88c07a57f..76926f599 100644 --- a/examples/PySDM_examples/Bartman_2020_MasterThesis/fig_4_adaptive_sdm.py +++ b/examples/PySDM_examples/Bartman_2020_MasterThesis/fig_4_adaptive_sdm.py @@ -54,6 +54,7 @@ def main(plot: bool = True, save: str = None): last_step_error = error if base_error is None: base_error = last_step_error + norm_error = last_step_error / base_error plotter.ylabel = ( r"$\bf{dt: " @@ -68,7 +69,9 @@ def main(plot: bool = True, save: str = None): if i == len(dts) - 1 else None ) - plotter.title = f"time: {norm_time:.2f} error: {last_step_error / base_error:.2f} (normalised)" + plotter.title = ( + f"time: {norm_time:.2f} error: {norm_error:.2f} (normalised)" + ) plotter.finished = False plotter.finish() if save is not None: From dc17c24045c8b7b6d6ed6c76f184ac31be12e249 Mon Sep 17 00:00:00 2001 From: Sylwester Arabas Date: Fri, 11 Oct 2024 01:14:20 +0200 Subject: [PATCH 10/10] add missing dependency on pystrict; address pylint find --- PySDM/backends/impl_numba/methods/physics_methods.py | 2 +- setup.py | 1 + 2 files changed, 2 insertions(+), 1 deletion(-) diff --git a/PySDM/backends/impl_numba/methods/physics_methods.py b/PySDM/backends/impl_numba/methods/physics_methods.py index 450d35493..f73b7a466 100644 --- a/PySDM/backends/impl_numba/methods/physics_methods.py +++ b/PySDM/backends/impl_numba/methods/physics_methods.py @@ -129,7 +129,7 @@ def body(mass, volume): def mass_of_volume(self): phys_volume_to_mass = self.formulae.particle_shape_and_density.volume_to_mass - @numba.njit(**{**conf.JIT_FLAGS, "fastmath": self.formulae.fastmath}) + @numba.njit(**self.default_jit_flags) def mass_of_volume(mass, volume): for i in prange(volume.shape[0]): # pylint: disable=not-an-iterable mass[i] = phys_volume_to_mass(volume[i]) diff --git a/setup.py b/setup.py index bf46134d4..4e3cdc79a 100644 --- a/setup.py +++ b/setup.py @@ -69,6 +69,7 @@ def get_long_description(): else "" ), "pyevtk" + ("==1.2.0" if CI else ""), + "pystrict", ], extras_require={ "tests": [