Merge pull request #692 from slayoo/freezing

no more cyclic imports (moving Formulae out of physics)

.github/workflows/pylint.yml

@@ -27,4 +27,4 @@ jobs:
     - name: Analysing the code with pylint
       run: |
         # TODO #682
-        pylint --disable=fixme,invalid-name,missing-function-docstring,missing-class-docstring,missing-module-docstring,too-many-arguments,too-many-locals,too-many-instance-attributes,too-few-public-methods,protected-access,too-many-statements,duplicate-code,too-many-branches,cyclic-import $(git ls-files '*.py')
+        pylint --disable=fixme,invalid-name,missing-function-docstring,missing-class-docstring,missing-module-docstring,too-many-arguments,too-many-locals,too-many-instance-attributes,too-few-public-methods,protected-access,too-many-statements,duplicate-code,too-many-branches $(git ls-files '*.py')

PySDM/__init__.py

@@ -25,6 +25,7 @@
 from pkg_resources import get_distribution, DistributionNotFound, VersionConflict
 from .builder import Builder
 from .particulator import Particulator
+from .formulae import Formulae
 
 try:
     __version__ = get_distribution(__name__).version

PySDM/backends/numba.py

@@ -13,7 +13,7 @@
 from PySDM.backends.impl_numba.methods.displacement_methods import DisplacementMethods
 from PySDM.backends.impl_numba.random import Random as ImportedRandom
 from PySDM.backends.impl_numba.storage import Storage as ImportedStorage
-from PySDM.physics import Formulae
+from PySDM.formulae import Formulae
 
 
 class Numba(  # pylint: disable=too-many-ancestors,duplicate-code

PySDM/backends/thrust_rtc.py

@@ -14,7 +14,7 @@
 from PySDM.backends.impl_thrust_rtc.methods.displacement_methods import DisplacementMethods
 from PySDM.backends.impl_thrust_rtc.storage import make_storage_class
 from PySDM.backends.impl_thrust_rtc.random import Random as ImportedRandom
-from PySDM.physics import Formulae
+from PySDM.formulae import Formulae
 from PySDM.backends.impl_thrust_rtc.conf import trtc
 
 

PySDM/builder.py

@@ -4,7 +4,7 @@
 import inspect
 import numpy as np
 from PySDM.particulator import Particulator
-from PySDM.initialisation.multiplicities import discretise_n  # TODO #324
+from PySDM.initialisation.discretise_multiplicities import discretise_multiplicities  # TODO #324
 from PySDM.impl.particle_attributes_factory import ParticlesFactory
 from PySDM.impl.wall_timer import WallTimer
 from PySDM.attributes.impl.mapper import get_class as attr_class
@@ -55,7 +55,7 @@ def request_attribute(self, attribute, variant=None):
         if variant is not None:
             assert variant == self.req_attr[attribute]
 
-    def build(self, attributes: dict, products: tuple = (), int_caster=discretise_n):
+    def build(self, attributes: dict, products: tuple = (), int_caster=discretise_multiplicities):
         assert self.particulator.environment is not None
 
         for dynamic in self.particulator.dynamics.values():

PySDM/dynamics/freezing.py

@@ -3,6 +3,7 @@
     SingularAttributes, TimeDependentAttributes
 )
 
+
 class Freezing:
     def __init__(self, *, singular=True):
         self.singular = singular

PySDM/environments/kinematic_1d.py

@@ -5,8 +5,8 @@
 
 import numpy as np
 from ._moist import _Moist
-from ..initialisation.r_wet_init import r_wet_init
-from ..initialisation.multiplicities import discretise_n
+from ..initialisation.equilibrate_wet_radii import equilibrate_wet_radii
+from ..initialisation.discretise_multiplicities import discretise_multiplicities
 from ..impl import arakawa_c
 
 
@@ -47,14 +47,16 @@ def init_attributes(self, *,
             r_dry, n_per_kg = spectral_discretisation.sample(self.particulator.n_sd)
             attributes['dry volume'] = self.formulae.trivia.volume(radius=r_dry)
             attributes['kappa times dry volume'] = attributes['dry volume'] * kappa
-            r_wet = r_wet_init(r_dry, self, cell_id=attributes['cell id'],
-                               kappa_times_dry_volume=attributes['kappa times dry volume'])
+            r_wet = equilibrate_wet_radii(
+                r_dry, self, cell_id=attributes['cell id'],
+                kappa_times_dry_volume=attributes['kappa times dry volume']
+            )
 
             rhod = self['rhod'].to_ndarray()
             cell_id = attributes['cell id']
             domain_volume = np.prod(np.array(self.mesh.size))
 
-        attributes['n'] = discretise_n(n_per_kg * rhod[cell_id] * domain_volume)
+        attributes['n'] = discretise_multiplicities(n_per_kg * rhod[cell_id] * domain_volume)
         attributes['volume'] = self.formulae.trivia.volume(radius=r_wet)
 
         return attributes

PySDM/environments/kinematic_2d.py

@@ -5,8 +5,8 @@
 
 import numpy as np
 from PySDM.impl.mesh import Mesh
-from PySDM.initialisation.r_wet_init import r_wet_init, default_rtol
-from PySDM.initialisation.multiplicities import discretise_n
+from PySDM.initialisation.equilibrate_wet_radii import equilibrate_wet_radii, default_rtol
+from PySDM.initialisation.discretise_multiplicities import discretise_multiplicities
 from ._moist import _Moist
 from ..impl import arakawa_c
 
@@ -60,7 +60,7 @@ def init_attributes(self, *,
             if kappa == 0:
                 r_wet = r_dry
             else:
-                r_wet = r_wet_init(
+                r_wet = equilibrate_wet_radii(
                     r_dry=r_dry,
                     environment=self,
                     kappa_times_dry_volume=attributes['kappa times dry volume'],
@@ -71,7 +71,7 @@ def init_attributes(self, *,
             cell_id = attributes['cell id']
             domain_volume = np.prod(np.array(self.mesh.size))
 
-        attributes['n'] = discretise_n(n_per_kg * rhod[cell_id] * domain_volume)
+        attributes['n'] = discretise_multiplicities(n_per_kg * rhod[cell_id] * domain_volume)
         attributes['volume'] = self.formulae.trivia.volume(radius=r_wet)
 
         return attributes

PySDM/environments/parcel.py

@@ -4,8 +4,8 @@
 
 import numpy as np
 from PySDM.impl.mesh import Mesh
-from PySDM.initialisation.r_wet_init import r_wet_init, default_rtol
-from PySDM.initialisation.multiplicities import discretise_n
+from PySDM.initialisation.equilibrate_wet_radii import equilibrate_wet_radii, default_rtol
+from PySDM.initialisation.discretise_multiplicities import discretise_multiplicities
 from ..physics import constants as const
 from ._moist import _Moist
 
@@ -73,8 +73,8 @@ def init_attributes(
         attributes = {}
         attributes['dry volume'] = self.formulae.trivia.volume(radius=r_dry)
         attributes['kappa times dry volume'] = attributes['dry volume'] * kappa
-        attributes['n'] = discretise_n(n_in_dv)
-        r_wet = r_wet_init(
+        attributes['n'] = discretise_multiplicities(n_in_dv)
+        r_wet = equilibrate_wet_radii(
             r_dry=r_dry,
             environment=self,
             kappa_times_dry_volume=attributes['kappa times dry volume'],

PySDM/initialisation/__init__.py

@@ -1,10 +1,6 @@
 """
-Initialisation logic incl.
-`PySDM.initialisation.spectral_sampling`,
-`PySDM.initialisation.spatial_sampling`,
-`PySDM.initialisation.spectra`
+Initialisation logic particle size spectra, sampling methods and
+wet radii equilibration
 """
-from .r_wet_init import r_wet_init
-from . import spectral_sampling
-from . import spatial_sampling
-from . import multiplicities
+from .discretise_multiplicities import discretise_multiplicities
+from .equilibrate_wet_radii import equilibrate_wet_radii

PySDM/initialisation/multiplicities.py → PySDM/initialisation/discretise_multiplicities.py

@@ -4,7 +4,7 @@
 import numpy as np
 
 
-def discretise_n(y_float):
+def discretise_multiplicities(y_float):
     y_int = y_float.round().astype(np.int64)
 
     percent_diff = 100 * abs(1 - np.sum(y_float) / np.sum(y_int.astype(float)))

PySDM/initialisation/r_wet_init.py → PySDM/initialisation/equilibrate_wet_radii.py

@@ -12,11 +12,11 @@
 default_max_iters = 64
 
 
-def r_wet_init(r_dry: np.ndarray, environment,
-               kappa_times_dry_volume: np.ndarray,
-               f_org: np.ndarray = None,
-               cell_id: np.ndarray = None,
-               rtol=default_rtol, max_iters=default_max_iters):
+def equilibrate_wet_radii(r_dry: np.ndarray, environment,
+                          kappa_times_dry_volume: np.ndarray,
+                          f_org: np.ndarray = None,
+                          cell_id: np.ndarray = None,
+                          rtol=default_rtol, max_iters=default_max_iters):
     if cell_id is None:
         cell_id = np.zeros_like(r_dry, dtype=int)
     if f_org is None:

PySDM/initialisation/impl/spectrum.py

@@ -0,0 +1,34 @@
+import numpy as np
+from PySDM.initialisation.sampling.spectral_sampling import default_cdf_range
+
+default_interpolation_grid = tuple(np.linspace(*default_cdf_range, 999))
+
+
+class Spectrum:
+
+    def __init__(self, distribution, distribution_params, norm_factor):
+        self.distribution_params = distribution_params  # (loc, scale)
+        self.norm_factor = norm_factor
+        self.distribution = distribution
+
+    def size_distribution(self, arg):
+        result = self.norm_factor * self.distribution.pdf(arg, *self.distribution_params)
+        return result
+
+    def pdf(self, arg):
+        return self.size_distribution(arg) / self.norm_factor
+
+    def cdf(self, arg):
+        return self.distribution.cdf(arg, *self.distribution_params)
+
+    def stats(self, moments):
+        result = self.distribution.stats(*self.distribution_params, moments)
+        return result
+
+    def cumulative(self, arg):
+        result = self.norm_factor * self.distribution.cdf(arg, *self.distribution_params)
+        return result
+
+    def percentiles(self, cdf_values):
+        result = self.distribution.ppf(cdf_values, *self.distribution_params)
+        return result

PySDM/initialisation/spectral_sampling.py → PySDM/initialisation/sampling/spectral_sampling.py

@@ -1,4 +1,4 @@
-from typing import Tuple
+from typing import Tuple, Optional
 import numpy as np
 from scipy import optimize
 from PySDM.physics import constants as const
@@ -7,8 +7,13 @@
 
 
 class SpectralSampling:
-    def __init__(self, spectrum, size_range: [None, Tuple[float, float]] = None):
+    def __init__(self,
+                 spectrum,
+                 size_range: [None, Tuple[float, float]] = None,
+                 error_threshold: Optional[float] = None
+     ):
         self.spectrum = spectrum
+        self.error_threshold = error_threshold or .01
 
         if size_range is None:
             if hasattr(spectrum, 'percentiles'):
@@ -28,15 +33,14 @@ def __init__(self, spectrum, size_range: [None, Tuple[float, float]] = None):
             assert size_range[1] > size_range[0]
             self.size_range = size_range
 
-    @staticmethod
-    def _sample(grid, spectrum):
+    def _sample(self, grid, spectrum):
         x = grid[1: -1: 2]
         cdf = spectrum.cumulative(grid[0::2])
         y_float = cdf[1:] - cdf[0:-1]
 
-        percent_diff = 100 * abs(1 - np.sum(y_float) / spectrum.norm_factor)
-        if percent_diff > 1:
-            raise Exception(f"{percent_diff}% error in total real-droplet number due to sampling")
+        diff = abs(1 - np.sum(y_float) / spectrum.norm_factor)
+        if diff > self.error_threshold:
+            raise Exception(f"{100*diff}% error in total real-droplet number due to sampling")
 
         return x, y_float
 
@@ -48,8 +52,12 @@ def sample(self, n_sd):
 
 
 class Logarithmic(SpectralSampling):
-    def __init__(self, spectrum, size_range: [None, Tuple[float, float]] = None):
-        super().__init__(spectrum, size_range)
+    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])
 

PySDM/initialisation/spectro_glacial.py → PySDM/initialisation/sampling/spectro_glacial_sampling.py

@@ -1,5 +1,5 @@
 import numpy as np
-from PySDM.initialisation.spectral_sampling import default_cdf_range
+from PySDM.initialisation.sampling.spectral_sampling import default_cdf_range
 from PySDM.physics import constants as const
 
 # DIM_SIZE = 0

PySDM/initialisation/spectra/__init__.py

@@ -0,0 +1,8 @@
+"""
+Classes representing logic around size spectra and more generally
+probability density functions (based on SciPy.stats logic)
+"""
+from .exponential import Exponential
+from .lognormal import Lognormal
+from .top_hat import TopHat
+from .sum import Sum

PySDM/initialisation/spectra/exponential.py

@@ -0,0 +1,10 @@
+from scipy.stats import expon
+from PySDM.initialisation.impl.spectrum import Spectrum
+
+class Exponential(Spectrum):
+
+    def __init__(self, norm_factor, scale):
+        super().__init__(expon, (
+            0,     # loc
+            scale  # scale = 1/lambda
+        ), norm_factor)

PySDM/initialisation/spectra/lognormal.py

@@ -0,0 +1,31 @@
+import math
+from scipy.stats import lognorm
+from PySDM.initialisation.impl.spectrum import Spectrum
+
+
+class Lognormal(Spectrum):
+
+    def __init__(self, norm_factor: float, m_mode: float, s_geom: float):
+        super().__init__(lognorm, (math.log(s_geom), 0, m_mode), norm_factor)
+
+    @property
+    def s_geom(self):
+        return math.exp(self.distribution_params[0])
+
+    @property
+    def m_mode(self):
+        return self.distribution_params[2]
+
+    @property
+    def median(self):
+        return self.m_mode
+
+    @property
+    def geometric_mean(self):
+        return self.s_geom
+
+    def __str__(self):
+        return f"{self.__class__.__name__}:"\
+               f" (N={self.norm_factor:.3g},"\
+               f" m_mode={self.m_mode:.3g},"\
+               f" s_geom={self.s_geom:.3g})"

PySDM/initialisation/spectra/sum.py

@@ -0,0 +1,30 @@
+import numpy as np
+from scipy.interpolate import interp1d
+from PySDM.initialisation.impl.spectrum import default_interpolation_grid
+
+
+class Sum:
+
+    def __init__(self, spectra: tuple, interpolation_grid=default_interpolation_grid):
+        self.spectra = spectra
+        self.norm_factor = sum((s.norm_factor for s in self.spectra))
+        percentiles = [s.percentiles(interpolation_grid) for s in self.spectra]
+        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)
+
+    def size_distribution(self, arg):
+        result = 0.
+        for spectrum in self.spectra:
+            result += spectrum.size_distribution(arg)
+        return result
+
+    def cumulative(self, arg):
+        result = 0.
+        for spectrum in self.spectra:
+            result += spectrum.cumulative(arg)
+        return result
+
+    def percentiles(self, cdf_values):
+        return self.inverse_cdf(cdf_values)

PySDM/initialisation/spectra/top_hat.py

@@ -0,0 +1,16 @@
+import numpy as np
+
+
+class TopHat:
+    def __init__(self, norm_factor, endpoints):
+        self.norm_factor = norm_factor
+        self.endpoints = endpoints
+        self._mn = endpoints[0]
+        self._mx = endpoints[1]
+
+    def cumulative(self, arg):
+        cdf = np.minimum(1, np.maximum(0, (arg - self._mn) / (self._mx - self._mn)))
+        return self.norm_factor * cdf
+
+    def percentiles(self, cdf_values):
+        return (self._mx - self._mn) * (np.asarray(cdf_values) + self._mn / (self._mx - self._mn))

PySDM/physics/__init__.py

@@ -6,4 +6,3 @@
     ventilation, state_variable_triplet, trivia, particle_advection, hydrostatics,
     freezing_temperature_spectrum, heterogeneous_ice_nucleation_rate)
 from .constants import si
-from .formulae import Formulae

PySDM/physics/dimensional_analysis.py

@@ -3,7 +3,8 @@
 for use in unit tests which disables Numba and enables Pint
 """
 from importlib import reload
-from . import constants, formulae
+from PySDM import formulae
+from . import constants
 from .impl import flag