super-droplet attribute "alpha" sampling (incl. new example based on Fig 1 from Matsushima et al. 2023) + refactor of the spectral sampling logic to enable deterministic/pseudorandom/quasirandom choice for each sampling type

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)

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"]

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)

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

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",

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)
 ```
 </details>
 <details>
@@ -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}));
 ```
 </details>
@@ -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)
 ```
 </details>
 
@@ -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)
 

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
         )

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(

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,

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",