super-droplet attribute "alpha" sampling (incl. new example based on Fig 1 from Matsushima et al. 2023)

PySDM/initialisation/sampling/spectral_sampling.py

@@ -7,6 +7,7 @@
 
 import numpy as np
 from scipy import optimize
+from scipy.interpolate import interp1d
 
 default_cdf_range = (0.00001, 0.99999)
 
@@ -62,12 +63,6 @@ def _sample(self, grid, spectrum):
         return x, y_float
 
 
-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)
-
-
 class Logarithmic(
     DeterministicSpectralSampling
 ):  # pylint: disable=too-few-public-methods
@@ -86,36 +81,109 @@ def sample(self, n_sd, *, backend=None):  # pylint: disable=unused-argument
         return self._sample(grid, self.spectrum)
 
 
-class ConstantMultiplicity(
+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(
     DeterministicSpectralSampling
 ):  # pylint: disable=too-few-public-methods
-    def __init__(self, spectrum, size_range=None):
+    """as in [Matsushima et al. 2023](https://doi.org/10.5194/gmd-16-6211-2023)"""
+
+    def __init__(self, spectrum, alpha, size_range=None, dist_0=None, dist_1=None):
         super().__init__(spectrum, size_range)
+        self.alpha = alpha
+        if dist_0 is None:
+            dist_0 = self.spectrum
+        if dist_1 is None:
 
-        self.cdf_range = (
-            spectrum.cumulative(self.size_range[0]),
-            spectrum.cumulative(self.size_range[1]),
-        )
-        assert 0 < self.cdf_range[0] < self.cdf_range[1]
+            def dist_1_inv(y):
+                return (self.size_range[1] - self.size_range[0]) * y
 
-    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)
+        else:
+            dist_1_inv = dist_1.percentiles
+        self.dist_0_cdf = dist_0.cdf
+        self.dist_1_inv = dist_1_inv
+
+    def sample(
+        self, n_sd, *, backend=None, xprime=None
+    ):  # pylint: disable=unused-argument
+
+        if xprime is None:
+            even_spec = np.linspace(
+                default_cdf_range[0], default_cdf_range[1], num=2 * n_sd + 3
+            )
+            x_prime = self.spectrum.percentiles(even_spec)
+        sd_cdf = self.dist_0_cdf(x_prime)
 
-        assert np.isfinite(percentiles).all()
+        x_sd_cdf = (1 - self.alpha) * x_prime + self.alpha * self.dist_1_inv(sd_cdf)
+
+        inv_cdf = interp1d(sd_cdf, x_sd_cdf)
+
+        percent_values = self._find_percentiles(n_sd, backend)
+        percentiles = inv_cdf(percent_values)
 
         return self._sample(percentiles, self.spectrum)
 
+    def _find_percentiles(self, n_sd, backend):
+        percent_values = np.linspace(
+            default_cdf_range[0], default_cdf_range[1], num=2 * n_sd + 1
+        )
+        return percent_values
+
 
-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)
+class AlphaSamplingPseudoRandom(
+    AlphaSampling
+):  # pylint: disable=too-few-public-methods
+    """Alpha sampling with pseudo-random values within deterministic percentile bins"""
+
+    def _find_percentiles(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()
 
-        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)
+        percent_values = np.linspace(
+            default_cdf_range[0], default_cdf_range[1], num=2 * n_sd + 1
+        )
+
+        for i in range(1, len(percent_values) - 1, 2):
+            percent_values[i] = percent_values[i - 1] + u01[i // 2] * (
+                percent_values[i + 1] - percent_values[i - 1]
+            )
+
+        return percent_values
+
+
+class AlphaSamplingRandom(AlphaSampling):  # pylint: disable=too-few-public-methods
+    """Alpha sampling with uniform random percentile bins"""
+
+    def _find_percentiles(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()
+
+        percent_values = np.sort(
+            default_cdf_range[0] + u01 * (default_cdf_range[1] - default_cdf_range[0])
+        )
+        return percent_values
+
+
+class ConstantMultiplicity(AlphaSampling):  # pylint: disable=too-few-public-methods
+    def __init__(self, spectrum, size_range=None):
+        super().__init__(spectrum, 0, size_range)
+
+
+class Linear(AlphaSampling):  # pylint: disable=too-few-public-methods
+    def __init__(self, spectrum, size_range=None):
+        super().__init__(spectrum, 1, size_range)

PySDM/initialisation/spectra/sum.py

@@ -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

@@ -946,5 +946,14 @@
     ],
     "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"
+    ]
   }
 }

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