diff --git a/PySDM/initialisation/sampling/spectral_sampling.py b/PySDM/initialisation/sampling/spectral_sampling.py
index 554795a575..c8ce4ce463 100644
--- a/PySDM/initialisation/sampling/spectral_sampling.py
+++ b/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)
diff --git a/PySDM/initialisation/spectra/sum.py b/PySDM/initialisation/spectra/sum.py
index e5a9c6d52f..bff539484b 100644
--- a/PySDM/initialisation/spectra/sum.py
+++ b/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
diff --git a/docs/bibliography.json b/docs/bibliography.json
index deb122ac9e..ef2bb62a2b 100644
--- a/docs/bibliography.json
+++ b/docs/bibliography.json
@@ -950,6 +950,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"
+ ]
"https://doi.org/10.48550/arXiv.2509.05536": {
"usages": [
"examples/PySDM_examples/Ware_et_al_2025/__init__.py",
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..6c7a6104db
--- /dev/null
+++ b/examples/PySDM_examples/Matsushima_et_al_2023/figure_1.ipynb
@@ -0,0 +1,5155 @@
+{
+ "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": [
+ "[Matsushima 2023 (GMD)](https://doi.org/10.5194/gmd-16-6211-2023)"
+ ]
+ },
+ {
+ "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": 10,
+ "id": "10da5008-8e20-43a1-8d11-0ade29c9a492",
+ "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": "ec089f1977eb4599af0d102028f290db",
+ "version_major": 2,
+ "version_minor": 0
+ },
+ "text/plain": [
+ "HBox(children=(HTML(value=\"./fig_1.pdf
\"), HTML(value=\"