diff --git a/PySDM/environments/kinematic_1d.py b/PySDM/environments/kinematic_1d.py
index 4eaa690843..81590173e7 100644
--- a/PySDM/environments/kinematic_1d.py
+++ b/PySDM/environments/kinematic_1d.py
@@ -60,12 +60,12 @@ def init_attributes(
) = self.mesh.cellular_attributes(positions)
if collisions_only:
- v_wet, n_per_kg = spectral_discretisation.sample(
+ v_wet, n_per_kg = spectral_discretisation.sample_deterministic(
backend=self.particulator.backend, n_sd=self.particulator.n_sd
)
attributes["volume"] = v_wet
else:
- r_dry, n_per_kg = spectral_discretisation.sample(
+ r_dry, n_per_kg = spectral_discretisation.sample_deterministic(
backend=self.particulator.backend, n_sd=self.particulator.n_sd
)
attributes["dry volume"] = self.formulae.trivia.volume(radius=r_dry)
diff --git a/PySDM/environments/kinematic_2d.py b/PySDM/environments/kinematic_2d.py
index 8e7087097d..5034a7bbc0 100644
--- a/PySDM/environments/kinematic_2d.py
+++ b/PySDM/environments/kinematic_2d.py
@@ -60,9 +60,9 @@ def init_attributes(
attributes["position in cell"],
) = self.mesh.cellular_attributes(positions)
- r_dry, n_per_kg = spectral_sampling(spectrum=dry_radius_spectrum).sample(
- n_sd=n_sd, backend=self.particulator.backend
- )
+ r_dry, n_per_kg = spectral_sampling(
+ spectrum=dry_radius_spectrum
+ ).sample_deterministic(n_sd=n_sd, backend=self.particulator.backend)
attributes["dry volume"] = self.formulae.trivia.volume(radius=r_dry)
attributes["kappa times dry volume"] = kappa * attributes["dry volume"]
diff --git a/PySDM/initialisation/sampling/spectral_sampling.py b/PySDM/initialisation/sampling/spectral_sampling.py
index 554795a575..c1198ec307 100644
--- a/PySDM/initialisation/sampling/spectral_sampling.py
+++ b/PySDM/initialisation/sampling/spectral_sampling.py
@@ -1,58 +1,51 @@
"""
-spectral sampling logic incl. linear, logarithmic, uniform-random and constant-multiplicity
- sampling classes
+spectral discretisation logic incl. linear, logarithmic, and constant-multiplicity
+ layouts with deterministic, pseudorandom and quasirandom sampling
"""
+import abc
from typing import Optional, Tuple
import numpy as np
-from scipy import optimize
+from scipy.interpolate import interp1d
default_cdf_range = (0.00001, 0.99999)
-class SpectralSampling: # pylint: disable=too-few-public-methods
- def __init__(self, spectrum, size_range: Optional[Tuple[float, float]] = None):
+class SpectralSampling:
+ def __init__(
+ self,
+ spectrum,
+ *,
+ size_range: Optional[Tuple[float, float]] = None,
+ error_threshold: Optional[float] = None,
+ ):
self.spectrum = spectrum
+ self.error_threshold = error_threshold or 0.01
if size_range is None:
- if hasattr(spectrum, "percentiles"):
- self.size_range = spectrum.percentiles(default_cdf_range)
- else:
- self.size_range = [np.nan, np.nan]
- for i in (0, 1):
- result = optimize.root(
- lambda x, value=default_cdf_range[i]: spectrum.cdf(x) - value,
- x0=spectrum.median(),
- )
- assert result.success
- self.size_range[i] = result.x
+ self.cdf_range = default_cdf_range
+ self.size_range = spectrum.percentiles(self.cdf_range)
else:
assert len(size_range) == 2
assert size_range[0] > 0
assert size_range[1] > size_range[0]
self.size_range = size_range
+ self.cdf_range = (
+ spectrum.cdf(size_range[0]),
+ spectrum.cdf(size_range[1]),
+ )
+ @abc.abstractmethod
+ def _sample(self, frac_values):
+ pass
-class DeterministicSpectralSampling(
- SpectralSampling
-): # pylint: disable=too-few-public-methods
- # TODO #1031 - error_threshold will be also used in non-deterministic sampling
- def __init__(
- self,
- spectrum,
- size_range: Optional[Tuple[float, float]] = None,
- error_threshold: Optional[float] = None,
- ):
- super().__init__(spectrum, size_range)
- self.error_threshold = error_threshold or 0.01
-
- def _sample(self, grid, spectrum):
+ def _sample_with_grid(self, grid):
x = grid[1:-1:2]
- cdf = spectrum.cumulative(grid[0::2])
+ cdf = self.spectrum.cumulative(grid[0::2])
y_float = cdf[1:] - cdf[0:-1]
- diff = abs(1 - np.sum(y_float) / spectrum.norm_factor)
+ diff = abs(1 - np.sum(y_float) / self.spectrum.norm_factor)
if diff > self.error_threshold:
raise ValueError(
f"{diff * 100:.3g}% error in total real-droplet number due to sampling "
@@ -61,61 +54,100 @@ def _sample(self, grid, spectrum):
return x, y_float
+ def sample_deterministic(
+ self, n_sd, *, backend=None
+ ): # pylint: disable=unused-argument
+ return self._sample(
+ frac_values=np.linspace(
+ self.cdf_range[0], self.cdf_range[1], num=2 * n_sd + 1
+ )
+ )
-class Linear(DeterministicSpectralSampling): # pylint: disable=too-few-public-methods
- def sample(self, n_sd, *, backend=None): # pylint: disable=unused-argument
- grid = np.linspace(*self.size_range, num=2 * n_sd + 1)
- return self._sample(grid, self.spectrum)
+ def sample_quasirandom(self, n_sd, *, backend):
+ num_elements = n_sd
+ storage = backend.Storage.empty(num_elements, dtype=float)
+ backend.Random(seed=backend.formulae.seed, size=num_elements)(storage)
+ u01 = storage.to_ndarray()
+ frac_values = np.linspace(
+ self.cdf_range[0], self.cdf_range[1], num=2 * n_sd + 1
+ )
-class Logarithmic(
- DeterministicSpectralSampling
-): # pylint: disable=too-few-public-methods
- def __init__(
- self,
- spectrum,
- size_range: [None, Tuple[float, float]] = None,
- error_threshold: Optional[float] = None,
- ):
- super().__init__(spectrum, size_range, error_threshold)
- self.start = np.log10(self.size_range[0])
- self.stop = np.log10(self.size_range[1])
+ for i in range(1, len(frac_values) - 1, 2):
+ frac_values[i] = frac_values[i - 1] + u01[i // 2] * (
+ frac_values[i + 1] - frac_values[i - 1]
+ )
- def sample(self, n_sd, *, backend=None): # pylint: disable=unused-argument
- grid = np.logspace(self.start, self.stop, num=2 * n_sd + 1)
- return self._sample(grid, self.spectrum)
+ return self._sample(frac_values=frac_values)
+ def sample_pseudorandom(self, n_sd, *, backend):
+ num_elements = 2 * n_sd + 1
+ storage = backend.Storage.empty(num_elements, dtype=float)
+ backend.Random(seed=backend.formulae.seed, size=num_elements)(storage)
+ u01 = storage.to_ndarray()
+
+ frac_values = np.sort(
+ self.cdf_range[0] + u01 * (self.cdf_range[1] - self.cdf_range[0])
+ )
+ return self._sample(frac_values=frac_values)
-class ConstantMultiplicity(
- DeterministicSpectralSampling
-): # pylint: disable=too-few-public-methods
- def __init__(self, spectrum, size_range=None):
- super().__init__(spectrum, size_range)
- self.cdf_range = (
- spectrum.cumulative(self.size_range[0]),
- spectrum.cumulative(self.size_range[1]),
+class Logarithmic(SpectralSampling):
+ def _sample(self, frac_values):
+ grid = np.exp(
+ (np.log(self.size_range[1]) - np.log(self.size_range[0])) * frac_values
+ + np.log(self.size_range[0])
)
- assert 0 < self.cdf_range[0] < self.cdf_range[1]
+ return self._sample_with_grid(grid)
- def sample(self, n_sd, *, backend=None): # pylint: disable=unused-argument
- cdf_arg = np.linspace(self.cdf_range[0], self.cdf_range[1], num=2 * n_sd + 1)
- cdf_arg /= self.spectrum.norm_factor
- percentiles = self.spectrum.percentiles(cdf_arg)
- assert np.isfinite(percentiles).all()
+class Linear(SpectralSampling):
+ def _sample(self, frac_values):
+ grid = self.size_range[0] + frac_values * (
+ self.size_range[1] - self.size_range[0]
+ )
+ return self._sample_with_grid(grid)
- return self._sample(percentiles, self.spectrum)
+class ConstantMultiplicity(SpectralSampling):
+ def _sample(self, frac_values):
+ grid = self.spectrum.percentiles(frac_values)
+ assert np.isfinite(grid).all()
+ return self._sample_with_grid(grid)
-class UniformRandom(SpectralSampling): # pylint: disable=too-few-public-methods
- def sample(self, n_sd, *, backend):
- n_elements = n_sd
- storage = backend.Storage.empty(n_elements, dtype=float)
- backend.Random(seed=backend.formulae.seed, size=n_elements)(storage)
- u01 = storage.to_ndarray()
- pdf_arg = self.size_range[0] + u01 * (self.size_range[1] - self.size_range[0])
- dr = abs(self.size_range[1] - self.size_range[0]) / n_sd
- # TODO #1031 - should also handle error_threshold check
- return pdf_arg, dr * self.spectrum.size_distribution(pdf_arg)
+class AlphaSampling(SpectralSampling):
+ """as in [Matsushima et al. 2023](https://doi.org/10.5194/gmd-16-6211-2023)"""
+
+ def __init__(
+ self,
+ spectrum,
+ *,
+ alpha,
+ dist_1_inv,
+ interp_points,
+ size_range=None,
+ error_threshold: Optional[float] = None,
+ ):
+ super().__init__(
+ spectrum, size_range=size_range, error_threshold=error_threshold
+ )
+ self.alpha = alpha
+ self.dist_0_cdf = self.spectrum.cdf
+ self.dist_1_inv = dist_1_inv
+ self.x_prime = np.linspace(
+ self.size_range[0], self.size_range[1], num=interp_points
+ )
+
+ def _sample(self, frac_values):
+ if self.alpha == 0:
+ frac_values = self.spectrum.percentiles(frac_values)
+ elif self.alpha == 1:
+ frac_values = self.dist_1_inv(frac_values, self.size_range)
+ else:
+ sd_cdf = self.dist_0_cdf(self.x_prime)
+ x_sd_cdf = (1 - self.alpha) * self.x_prime + self.alpha * self.dist_1_inv(
+ sd_cdf, self.size_range
+ )
+ frac_values = interp1d(sd_cdf, x_sd_cdf)(frac_values)
+ return self._sample_with_grid(frac_values)
diff --git a/PySDM/initialisation/spectra/sum.py b/PySDM/initialisation/spectra/sum.py
index e5a9c6d52f..25f626a259 100644
--- a/PySDM/initialisation/spectra/sum.py
+++ b/PySDM/initialisation/spectra/sum.py
@@ -20,7 +20,7 @@ def __init__(self, spectra: tuple, interpolation_grid=None):
cdf_arg = np.zeros(len(interpolation_grid) * len(self.spectra) + 1)
cdf_arg[1:] = np.concatenate(percentiles)
cdf = self.cumulative(cdf_arg) / self.norm_factor
- self.inverse_cdf = interp1d(cdf, cdf_arg)
+ self.inverse_cdf = interp1d(cdf, cdf_arg, bounds_error=False)
def size_distribution(self, arg):
result = 0.0
@@ -36,3 +36,9 @@ def cumulative(self, arg):
def percentiles(self, cdf_values):
return self.inverse_cdf(cdf_values)
+
+ def pdf(self, arg):
+ return self.size_distribution(arg) / self.norm_factor
+
+ def cdf(self, arg):
+ return self.cumulative(arg) / self.norm_factor
diff --git a/docs/bibliography.json b/docs/bibliography.json
index f8d14ba18f..c455491e5b 100644
--- a/docs/bibliography.json
+++ b/docs/bibliography.json
@@ -951,6 +951,15 @@
"title": "A physically constrained classical description of the homogeneous nucleation of ice in water",
"label": "Koop and Murray 2016 (J. Chem. Phys. 145)"
},
+ "https://doi.org/10.5194/gmd-16-6211-2023": {
+ "title": "Overcoming computational challenges to realize meter- to submeter-scale resolution in cloud simulations using the super-droplet method",
+ "label": "Matsushima et al. 2023 (Geosci. Model Dev. 16)",
+ "usages": [
+ "PySDM/initialisation/sampling/spectral_sampling.py",
+ "examples/PySDM_examples/Matsushima_et_al_2023/figure_1.ipynb",
+ "examples/PySDM_examples/Matsushima_et_al_2023/__init__.py"
+ ]
+ },
"https://doi.org/10.48550/arXiv.2509.05536": {
"usages": [
"examples/PySDM_examples/Ware_et_al_2025/__init__.py",
diff --git a/docs/markdown/pysdm_landing.md b/docs/markdown/pysdm_landing.md
index 30623d05b4..48b8a2a2b8 100644
--- a/docs/markdown/pysdm_landing.md
+++ b/docs/markdown/pysdm_landing.md
@@ -103,7 +103,7 @@ Exponential = pyimport("PySDM.initialisation.spectra").Exponential
n_sd = 2^15
initial_spectrum = Exponential(norm_factor=8.39e12, scale=1.19e5 * si.um^3)
attributes = Dict()
-attributes["volume"], attributes["multiplicity"] = ConstantMultiplicity(spectrum=initial_spectrum).sample(n_sd)
+attributes["volume"], attributes["multiplicity"] = ConstantMultiplicity(spectrum=initial_spectrum).sample_deterministic(n_sd)
```
@@ -119,7 +119,7 @@ initial_spectrum = Exponential(pyargs(...
'norm_factor', 8.39e12, ...
'scale', 1.19e5 * si.um ^ 3 ...
));
-tmp = ConstantMultiplicity(initial_spectrum).sample(int32(n_sd));
+tmp = ConstantMultiplicity(initial_spectrum).sample_deterministic(int32(n_sd));
attributes = py.dict(pyargs('volume', tmp{1}, 'multiplicity', tmp{2}));
```
@@ -133,7 +133,7 @@ from PySDM.initialisation.spectra.exponential import Exponential
n_sd = 2 ** 15
initial_spectrum = Exponential(norm_factor=8.39e12, scale=1.19e5 * si.um ** 3)
attributes = {}
-attributes['volume'], attributes['multiplicity'] = ConstantMultiplicity(initial_spectrum).sample(n_sd)
+attributes['volume'], attributes['multiplicity'] = ConstantMultiplicity(initial_spectrum).sample_deterministic(n_sd)
```
@@ -323,7 +323,7 @@ The component submodules used to create this simulation are visualized below:
IS["initial_spectrum :Exponential"] -->|passed as arg to| CM_INIT
CM_INIT(["ConstantMultiplicity.__init__()"]) -->|instantiates| CM_INSTANCE
CM_INSTANCE[":ConstantMultiplicity"] -.-|has a method| SAMPLE
- SAMPLE(["ConstantMultiplicity.sample()"]) -->|returns| n
+ SAMPLE(["ConstantMultiplicity.sample_deterministic()"]) -->|returns| n
SAMPLE -->|returns| volume
n -->|added as element of| ATTRIBUTES
PARTICULATOR_INSTANCE -.-|has a method| PARTICULATOR_RUN(["Particulator.run()"])
@@ -414,7 +414,7 @@ builder = Builder(backend=CPU(formulae), n_sd=n_sd, environment=env)
builder.add_dynamic(AmbientThermodynamics())
builder.add_dynamic(Condensation())
-r_dry, specific_concentration = spectral_sampling.Logarithmic(spectrum).sample(n_sd)
+r_dry, specific_concentration = spectral_sampling.Logarithmic(spectrum).sample_deterministic(n_sd)
v_dry = formulae.trivia.volume(radius=r_dry)
r_wet = equilibrate_wet_radii(r_dry=r_dry, environment=builder.particulator.environment, kappa_times_dry_volume=kappa * v_dry)
@@ -495,7 +495,7 @@ builder = Builder(pyargs('backend', CPU(formulae), 'n_sd', int32(n_sd), 'environ
builder.add_dynamic(AmbientThermodynamics());
builder.add_dynamic(Condensation());
-tmp = spectral_sampling.Logarithmic(spectrum).sample(int32(n_sd));
+tmp = spectral_sampling.Logarithmic(spectrum).sample_deterministic(int32(n_sd));
r_dry = tmp{1};
v_dry = formulae.trivia.volume(pyargs('radius', r_dry));
specific_concentration = tmp{2};
@@ -596,7 +596,7 @@ builder = Builder(backend=CPU(formulae), n_sd=n_sd, environment=env)
builder.add_dynamic(AmbientThermodynamics())
builder.add_dynamic(Condensation())
-r_dry, specific_concentration = spectral_sampling.Logarithmic(spectrum).sample(n_sd)
+r_dry, specific_concentration = spectral_sampling.Logarithmic(spectrum).sample_deterministic(n_sd)
v_dry = formulae.trivia.volume(radius=r_dry)
r_wet = equilibrate_wet_radii(r_dry=r_dry, environment=builder.particulator.environment, kappa_times_dry_volume=kappa * v_dry)
diff --git a/examples/PySDM_examples/Abade_and_Albuquerque_2024/simulation.py b/examples/PySDM_examples/Abade_and_Albuquerque_2024/simulation.py
index 3838804104..9f790a1242 100644
--- a/examples/PySDM_examples/Abade_and_Albuquerque_2024/simulation.py
+++ b/examples/PySDM_examples/Abade_and_Albuquerque_2024/simulation.py
@@ -49,9 +49,9 @@ def __init__(self, settings):
if settings.enable_vapour_deposition_on_ice:
builder.add_dynamic(VapourDepositionOnIce(adaptive=True))
- r_dry, n_in_dv = ConstantMultiplicity(settings.soluble_aerosol).sample(
- n_sd=settings.n_sd
- )
+ r_dry, n_in_dv = ConstantMultiplicity(
+ settings.soluble_aerosol
+ ).sample_deterministic(n_sd=settings.n_sd)
attributes = builder.particulator.environment.init_attributes(
n_in_dv=n_in_dv, kappa=settings.kappa, r_dry=r_dry
)
diff --git a/examples/PySDM_examples/Abdul_Razzak_Ghan_2000/run_ARG_parcel.py b/examples/PySDM_examples/Abdul_Razzak_Ghan_2000/run_ARG_parcel.py
index 1321db40e6..ccd8e8ed1e 100644
--- a/examples/PySDM_examples/Abdul_Razzak_Ghan_2000/run_ARG_parcel.py
+++ b/examples/PySDM_examples/Abdul_Razzak_Ghan_2000/run_ARG_parcel.py
@@ -61,7 +61,9 @@ def run_parcel(
}
for i, mode in enumerate(aerosol.modes):
kappa, spectrum = mode["kappa"]["CompressedFilmOvadnevaite"], mode["spectrum"]
- r_dry, concentration = ConstantMultiplicity(spectrum).sample(n_sd_per_mode)
+ r_dry, concentration = ConstantMultiplicity(spectrum).sample_deterministic(
+ n_sd_per_mode
+ )
v_dry = builder.formulae.trivia.volume(radius=r_dry)
specific_concentration = concentration / builder.formulae.constants.rho_STP
attributes["multiplicity"] = np.append(
diff --git a/examples/PySDM_examples/Alpert_and_Knopf_2016/simulation.py b/examples/PySDM_examples/Alpert_and_Knopf_2016/simulation.py
index 1c3c9d32e4..c121210dab 100644
--- a/examples/PySDM_examples/Alpert_and_Knopf_2016/simulation.py
+++ b/examples/PySDM_examples/Alpert_and_Knopf_2016/simulation.py
@@ -227,7 +227,9 @@ def simulation(
n_sd, multiplicity / volume
)
else:
- _isa, _conc = spectral_sampling.ConstantMultiplicity(spectrum).sample(n_sd)
+ _isa, _conc = spectral_sampling.ConstantMultiplicity(
+ spectrum
+ ).sample_deterministic(n_sd)
attributes = {
"multiplicity": discretise_multiplicities(_conc * volume),
"immersed surface area": _isa,
diff --git a/examples/PySDM_examples/Arabas_et_al_2025/aida.ipynb b/examples/PySDM_examples/Arabas_et_al_2025/aida.ipynb
index c224b4daaf..393dbd82f3 100644
--- a/examples/PySDM_examples/Arabas_et_al_2025/aida.ipynb
+++ b/examples/PySDM_examples/Arabas_et_al_2025/aida.ipynb
@@ -301,7 +301,7 @@
}
},
"source": [
- "d_dry, conc = spectral_sampling.ConstantMultiplicity(dNdD).sample(common['n_sd'])\n",
+ "d_dry, conc = spectral_sampling.ConstantMultiplicity(dNdD).sample_deterministic(common['n_sd'])\n",
"\n",
"ATTRIBUTES = {\n",
" 'n': conc * common['volume'],\n",
diff --git a/examples/PySDM_examples/Bieli_et_al_2022/simulation.py b/examples/PySDM_examples/Bieli_et_al_2022/simulation.py
index ec108c476a..cadac6d489 100644
--- a/examples/PySDM_examples/Bieli_et_al_2022/simulation.py
+++ b/examples/PySDM_examples/Bieli_et_al_2022/simulation.py
@@ -18,7 +18,7 @@ def make_core(settings, coal_eff):
attributes = {}
attributes["volume"], attributes["multiplicity"] = ConstantMultiplicity(
settings.spectrum
- ).sample(settings.n_sd)
+ ).sample_deterministic(settings.n_sd)
collision = Collision(
collision_kernel=settings.kernel,
coalescence_efficiency=coal_eff,
diff --git a/examples/PySDM_examples/Ervens_and_Feingold_2012/settings.py b/examples/PySDM_examples/Ervens_and_Feingold_2012/settings.py
index 5e588639c7..9872afb777 100644
--- a/examples/PySDM_examples/Ervens_and_Feingold_2012/settings.py
+++ b/examples/PySDM_examples/Ervens_and_Feingold_2012/settings.py
@@ -9,4 +9,4 @@ def sampled_ccn_diameter_number_concentration_spectrum(
return Logarithmic(
spectrum=Lognormal(s_geom=1.4, m_mode=0.04 * si.um, norm_factor=100 / si.cm**3),
size_range=size_range,
- ).sample(n_sd)
+ ).sample_deterministic(n_sd)
diff --git a/examples/PySDM_examples/Grabowski_and_Pawlowska_2023/simulation.py b/examples/PySDM_examples/Grabowski_and_Pawlowska_2023/simulation.py
index d0d46799a1..93348f3149 100644
--- a/examples/PySDM_examples/Grabowski_and_Pawlowska_2023/simulation.py
+++ b/examples/PySDM_examples/Grabowski_and_Pawlowska_2023/simulation.py
@@ -17,7 +17,6 @@ def __init__(
settings,
products=None,
scipy_solver=False,
- sampling_class=ConstantMultiplicity,
):
builder = Builder(
n_sd=settings.n_sd,
@@ -55,7 +54,9 @@ def __init__(
kappa = tuple(settings.aerosol_modes_by_kappa.keys())[0]
spectrum = settings.aerosol_modes_by_kappa[kappa]
- r_dry, n_per_volume = sampling_class(spectrum).sample(settings.n_sd)
+ r_dry, n_per_volume = ConstantMultiplicity(spectrum).sample_deterministic(
+ settings.n_sd
+ )
v_dry = settings.formulae.trivia.volume(radius=r_dry)
attributes["multiplicity"] = np.append(
attributes["multiplicity"], n_per_volume * volume
diff --git a/examples/PySDM_examples/Graf_et_al_2019/figure_4.ipynb b/examples/PySDM_examples/Graf_et_al_2019/figure_4.ipynb
index 0303e5c9ea..bf1466d2ed 100644
--- a/examples/PySDM_examples/Graf_et_al_2019/figure_4.ipynb
+++ b/examples/PySDM_examples/Graf_et_al_2019/figure_4.ipynb
@@ -176,7 +176,7 @@
},
"outputs": [],
"source": [
- "r_dry, specific_concentration = spectral_sampling.Logarithmic(spectrum).sample(n_sd)"
+ "r_dry, specific_concentration = spectral_sampling.Logarithmic(spectrum).sample_deterministic(n_sd)"
]
},
{
diff --git a/examples/PySDM_examples/Jensen_and_Nugent_2017/simulation.py b/examples/PySDM_examples/Jensen_and_Nugent_2017/simulation.py
index 289e8f27cc..827d833fff 100644
--- a/examples/PySDM_examples/Jensen_and_Nugent_2017/simulation.py
+++ b/examples/PySDM_examples/Jensen_and_Nugent_2017/simulation.py
@@ -65,7 +65,7 @@ def __init__(
self.r_dry, n_in_unit_volume = Logarithmic(
spectrum=settings.dry_radii_spectrum,
- ).sample(builder.particulator.n_sd - n_gccn)
+ ).sample_deterministic(builder.particulator.n_sd - n_gccn)
if gccn:
nonzero_concentration_mask = np.nonzero(table_3.NA)
diff --git a/examples/PySDM_examples/Kreidenweis_et_al_2003/settings.py b/examples/PySDM_examples/Kreidenweis_et_al_2003/settings.py
index eba95a6120..a80eed0fdf 100644
--- a/examples/PySDM_examples/Kreidenweis_et_al_2003/settings.py
+++ b/examples/PySDM_examples/Kreidenweis_et_al_2003/settings.py
@@ -17,7 +17,8 @@ def __init__(
dt: float,
n_sd: int,
n_substep: int,
- spectral_sampling: spec_sampling.SpectralSampling = spec_sampling.Logarithmic,
+ spectral_sampling_class: spec_sampling.SpectralSampling = spec_sampling.Logarithmic,
+ spectral_sampling_method: str = "sample_deterministic",
):
self.formulae = Formulae(
saturation_vapour_pressure="AugustRocheMagnus",
@@ -51,13 +52,16 @@ def __init__(
# note: rho is not specified in the paper
rho0 = 1
- self.r_dry, self.n_in_dv = spectral_sampling(
- spectrum=spectra.Lognormal(
- norm_factor=566 / si.cm**3 / rho0 * self.mass_of_dry_air,
- m_mode=0.08 * si.um / 2,
- s_geom=2,
- )
- ).sample(n_sd)
+ self.r_dry, self.n_in_dv = getattr(
+ spectral_sampling_class(
+ spectrum=spectra.Lognormal(
+ norm_factor=566 / si.cm**3 / rho0 * self.mass_of_dry_air,
+ m_mode=0.08 * si.um / 2,
+ s_geom=2,
+ )
+ ),
+ spectral_sampling_method,
+ )(n_sd)
self.ENVIRONMENT_MOLE_FRACTIONS = {
"SO2": 0.2 * PPB,
diff --git a/examples/PySDM_examples/Lowe_et_al_2019/fig_2.ipynb b/examples/PySDM_examples/Lowe_et_al_2019/fig_2.ipynb
index 0bce96f9b5..81918cb04f 100644
--- a/examples/PySDM_examples/Lowe_et_al_2019/fig_2.ipynb
+++ b/examples/PySDM_examples/Lowe_et_al_2019/fig_2.ipynb
@@ -50,7 +50,6 @@
"from open_atmos_jupyter_utils import show_plot\n",
"\n",
"from PySDM import Formulae\n",
- "from PySDM.initialisation.sampling import spectral_sampling as spec_sampling\n",
"from PySDM.initialisation.spectra import Sum\n",
"from PySDM.physics import si\n",
"\n",
@@ -85,7 +84,6 @@
" dz=1*si.m, n_sd_per_mode=100, \n",
" model=model,\n",
" aerosol=aerosol,\n",
- " spectral_sampling=spec_sampling.ConstantMultiplicity,\n",
" )\n",
" simulation = Simulation(settings)\n",
" output[key] = simulation.run()\n",
diff --git a/examples/PySDM_examples/Lowe_et_al_2019/fig_3.ipynb b/examples/PySDM_examples/Lowe_et_al_2019/fig_3.ipynb
index 9c9c3a41fb..5bccaca919 100644
--- a/examples/PySDM_examples/Lowe_et_al_2019/fig_3.ipynb
+++ b/examples/PySDM_examples/Lowe_et_al_2019/fig_3.ipynb
@@ -46,7 +46,6 @@
"from PySDM_examples.Lowe_et_al_2019.constants_def import LOWE_CONSTS\n",
"\n",
"from PySDM import Formulae\n",
- "from PySDM.initialisation.sampling import spectral_sampling as spec_sampling\n",
"from PySDM.physics import si\n",
"\n",
"import numpy as np\n",
@@ -123,7 +122,6 @@
" \"d\": AerosolBoreal(water_molar_volume=WATER_MOLAR_VOLUME, Forg=Forg, Acc_N2=Acc[\"d\"])\n",
" }[subplot],\n",
" w = w * si.m / si.s,\n",
- " spectral_sampling = spec_sampling.ConstantMultiplicity,\n",
" ))\n",
" for w in updraft_list\n",
" for Forg in forg_list\n",
diff --git a/examples/PySDM_examples/Lowe_et_al_2019/fig_s2.ipynb b/examples/PySDM_examples/Lowe_et_al_2019/fig_s2.ipynb
index b32a46e213..e6851bec5d 100644
--- a/examples/PySDM_examples/Lowe_et_al_2019/fig_s2.ipynb
+++ b/examples/PySDM_examples/Lowe_et_al_2019/fig_s2.ipynb
@@ -55,7 +55,6 @@
"from open_atmos_jupyter_utils import show_plot\n",
"\n",
"from PySDM import Formulae\n",
- "from PySDM.initialisation.sampling import spectral_sampling as spec_sampling\n",
"from PySDM.initialisation.spectra import Sum\n",
"from PySDM.physics import si, in_unit\n",
"\n",
@@ -147,7 +146,6 @@
" model = model,\n",
" aerosol = aerosol,\n",
" w = w * si.m / si.s,\n",
- " spectral_sampling = spec_sampling.ConstantMultiplicity,\n",
" ))\n",
" for w in updraft_list\n",
" for model in models\n",
diff --git a/examples/PySDM_examples/Lowe_et_al_2019/settings.py b/examples/PySDM_examples/Lowe_et_al_2019/settings.py
index c8d5b0f44e..847788e15e 100644
--- a/examples/PySDM_examples/Lowe_et_al_2019/settings.py
+++ b/examples/PySDM_examples/Lowe_et_al_2019/settings.py
@@ -7,7 +7,6 @@
from PySDM import Formulae
from PySDM.backends import CPU
from PySDM.initialisation.aerosol_composition import DryAerosolMixture
-from PySDM.initialisation.sampling import spectral_sampling as spec_sampling
from PySDM.physics import si
@@ -36,7 +35,6 @@ def __init__(
n_sd_per_mode: int,
aerosol: DryAerosolMixture,
model: str,
- spectral_sampling: type(spec_sampling.SpectralSampling),
w: float = 0.32 * si.m / si.s,
):
assert model in ("Constant", "CompressedFilmOvadnevaite")
@@ -46,7 +44,6 @@ def __init__(
self.formulae = self.backend.formulae
const = self.formulae.constants
self.aerosol = aerosol
- self.spectral_sampling = spectral_sampling
max_altitude = 200 * si.m
self.w = w
diff --git a/examples/PySDM_examples/Lowe_et_al_2019/simulation.py b/examples/PySDM_examples/Lowe_et_al_2019/simulation.py
index 850430181d..0fd85de952 100644
--- a/examples/PySDM_examples/Lowe_et_al_2019/simulation.py
+++ b/examples/PySDM_examples/Lowe_et_al_2019/simulation.py
@@ -7,6 +7,7 @@
from PySDM.environments import Parcel
from PySDM.initialisation.hygroscopic_equilibrium import equilibrate_wet_radii
from PySDM.initialisation.spectra import Sum
+from PySDM.initialisation.sampling.spectral_sampling import ConstantMultiplicity
class Simulation(BasicSimulation):
@@ -33,9 +34,9 @@ def __init__(self, settings, products=None):
}
initial_volume = settings.mass_of_dry_air / settings.rho0
for mode in settings.aerosol.modes:
- r_dry, n_in_dv = settings.spectral_sampling(
+ r_dry, n_in_dv = ConstantMultiplicity(
spectrum=mode["spectrum"]
- ).sample(settings.n_sd_per_mode)
+ ).sample_deterministic(settings.n_sd_per_mode)
v_dry = settings.formulae.trivia.volume(radius=r_dry)
attributes["multiplicity"] = np.append(
attributes["multiplicity"], n_in_dv * initial_volume
diff --git a/examples/PySDM_examples/Matsushima_et_al_2023/__init__.py b/examples/PySDM_examples/Matsushima_et_al_2023/__init__.py
new file mode 100644
index 0000000000..dd4ffe2934
--- /dev/null
+++ b/examples/PySDM_examples/Matsushima_et_al_2023/__init__.py
@@ -0,0 +1,7 @@
+# pylint: disable=invalid-name
+"""
+Fig. 1 from [Matsushima et al. 2023 (GMD)](https://doi.org/10.5194/gmd-16-6211-2023)
+
+figure_1.ipynb:
+.. include:: ./figure_1.ipynb.badges.md
+"""
diff --git a/examples/PySDM_examples/Matsushima_et_al_2023/figure_1.ipynb b/examples/PySDM_examples/Matsushima_et_al_2023/figure_1.ipynb
new file mode 100644
index 0000000000..b2a866e9eb
--- /dev/null
+++ b/examples/PySDM_examples/Matsushima_et_al_2023/figure_1.ipynb
@@ -0,0 +1,16027 @@
+{
+ "cells": [
+ {
+ "cell_type": "markdown",
+ "id": "1790321c-2424-4ebd-9e26-14637e1ec29f",
+ "metadata": {},
+ "source": [
+ "[](https://github.com/open-atmos/PySDM/blob/main/examples/PySDM_examples/Matsushima_et_al_2023/figure_1.ipynb)\n",
+ "[](https://mybinder.org/v2/gh/open-atmos/PySDM.git/main?urlpath=lab/tree/examples/PySDM_examples/Matsushima_et_al_2023/figure_1.ipynb)\n",
+ "[](https://colab.research.google.com/github/open-atmos/PySDM/blob/main/examples/PySDM_examples/Matsushima_et_al_2023/figure_1.ipynb)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "6bc29eca-b273-4418-bedc-4787635bc19d",
+ "metadata": {},
+ "source": [
+ "Fig. 1 from [Matsushima 2023 (GMD)](https://doi.org/10.5194/gmd-16-6211-2023) depicting alpha-sampling in PySDM"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 1,
+ "id": "3642a633-e611-4001-8bb0-cbe95e12f232",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "import sys\n",
+ "if 'google.colab' in sys.modules:\n",
+ " !pip --quiet install open-atmos-jupyter-utils\n",
+ " from open_atmos_jupyter_utils import pip_install_on_colab\n",
+ " pip_install_on_colab('PySDM-examples')"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 2,
+ "id": "c87ff404-3a80-4d0c-a955-3bc7d0a4d31e",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "import numpy as np\n",
+ "\n",
+ "from PySDM.physics import si, in_unit\n",
+ "from PySDM.initialisation.sampling.spectral_sampling import AlphaSampling\n",
+ "from PySDM.initialisation import spectra\n",
+ "\n",
+ "from matplotlib import pyplot, colors\n",
+ "from open_atmos_jupyter_utils import show_plot"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 3,
+ "id": "10da5008-8e20-43a1-8d11-0ade29c9a492",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "n_sd = 2**16\n",
+ "spectrum = spectra.Sum((\n",
+ " spectra.Lognormal(\n",
+ " norm_factor=90 / si.cm**3,\n",
+ " m_mode=0.03 * si.um,\n",
+ " s_geom=1.28,\n",
+ " ),\n",
+ " spectra.Lognormal(\n",
+ " norm_factor=15 / si.cm**3,\n",
+ " m_mode=0.14 * si.um,\n",
+ " s_geom=1.75,\n",
+ " ),\n",
+ "))\n",
+ "\n",
+ "samplings = {\n",
+ " alpha: AlphaSampling(\n",
+ " spectrum,\n",
+ " alpha=alpha,\n",
+ " interp_points=100000,\n",
+ " size_range=(10 * si.nm, 5000 * si.nm),\n",
+ " dist_1_inv=lambda y, size_range: np.exp((np.log(size_range[1]) - np.log(size_range[0])) * y + np.log(size_range[0]))\n",
+ " ) for alpha in np.linspace(0, 1, 11) \n",
+ "}\n",
+ "\n",
+ "xas, yas = {}, {}\n",
+ "for alpha, sampling in samplings.items():\n",
+ " xas[alpha], yas[alpha] = sampling.sample_deterministic(n_sd)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 4,
+ "id": "858d3450-6107-4d33-b88b-daa5402e4ac7",
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "image/svg+xml": [
+ "\n",
+ "\n",
+ "\n"
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ },
+ {
+ "data": {
+ "application/vnd.jupyter.widget-view+json": {
+ "model_id": "daeba983b9ce4c04b04221c099f18313",
+ "version_major": 2,
+ "version_minor": 0
+ },
+ "text/plain": [
+ "HBox(children=(HTML(value=\"./fig_1.pdf
\"), HTML(value=\"= m_range[0]
+ assert np.max(m) <= m_range[1]
+ actual = np.sum(n)
+ desired = spectrum.cumulative(m_range[1]) - spectrum.cumulative(m_range[0])
+ np.testing.assert_approx_equal(
+ actual=actual / desired,
+ desired=1.0,
+ significant=2 if method == "pseudorandom" else 4,
+ )
- # Assert
- assert m.shape == n.shape
- assert n.shape == (n_sd,)
- assert np.min(m) >= m_range[0]
- assert np.max(m) <= m_range[1]
- actual = np.sum(n)
- desired = spectrum.cumulative(m_range[1]) - spectrum.cumulative(m_range[0])
- quotient = actual / desired
- np.testing.assert_approx_equal(actual=quotient, desired=1.0, significant=2)
+ @staticmethod
+ @pytest.mark.parametrize(
+ "sampling_class, error_threshold",
+ (
+ (spectral_sampling.ConstantMultiplicity, None),
+ (spectral_sampling.Logarithmic, None),
+ (spectral_sampling.Linear, None),
+ pytest.param(
+ spectral_sampling.ConstantMultiplicity,
+ 1e-5,
+ marks=pytest.mark.xfail(strict=True, raises=ValueError),
+ ),
+ pytest.param(
+ spectral_sampling.Logarithmic,
+ 1e-5,
+ marks=pytest.mark.xfail(strict=True, raises=ValueError),
+ ),
+ pytest.param(
+ spectral_sampling.Linear,
+ 1e-5,
+ marks=pytest.mark.xfail(strict=True, raises=ValueError),
+ ),
+ ),
+ )
+ def test_error_threshold_with_deterministic_sampling(
+ sampling_class, error_threshold
+ ):
+ sampling_class(
+ spectra.Lognormal(n_part, m_mode, s_geom), error_threshold=error_threshold
+ ).sample_deterministic(n_sd=10)
diff --git a/tests/unit_tests/products/test_activation_criteria.py b/tests/unit_tests/products/test_activation_criteria.py
index ef37d6931e..37b4dde0f9 100644
--- a/tests/unit_tests/products/test_activation_criteria.py
+++ b/tests/unit_tests/products/test_activation_criteria.py
@@ -45,7 +45,7 @@ def test_activation_criteria(backend, plot=False):
r_dry, specific_concentration = spectral_sampling.ConstantMultiplicity(
Lognormal(norm_factor=1e4 / si.mg, m_mode=50 * si.nm, s_geom=1.5)
- ).sample(builder.particulator.n_sd)
+ ).sample_deterministic(builder.particulator.n_sd)
particulator = builder.build(
attributes=builder.particulator.environment.init_attributes(
diff --git a/tests/unit_tests/test_imports.py b/tests/unit_tests/test_imports.py
index 941543bc00..9dd92194ed 100644
--- a/tests/unit_tests/test_imports.py
+++ b/tests/unit_tests/test_imports.py
@@ -24,7 +24,7 @@
"exporters.VTKExporter",
"initialisation.aerosol_composition.DryAerosolMixture",
"initialisation.init_fall_momenta",
- "initialisation.sampling.spectral_sampling.DeterministicSpectralSampling",
+ "initialisation.sampling.spectral_sampling.AlphaSampling",
"initialisation.spectra.Lognormal",
"physics.constants_defaults",
"physics.diffusion_thermics.LoweEtAl2019",