From f2147d0d2d88e24e0b8382bf23cf37d5876a5324 Mon Sep 17 00:00:00 2001 From: Konrad Bodzioch Date: Tue, 21 Oct 2025 12:32:54 +0200 Subject: [PATCH 1/4] backends: add JAX backend class (stripped) New draft backend class consisting of stripped Numba backend code sufficient to run the tutorial collision example. --- .github/workflows/readme_snippets.yml | 2 + PySDM/backends/__init__.py | 5 + PySDM/backends/impl_jax/__init__.py | 0 PySDM/backends/impl_jax/methods/__init__.py | 7 + .../impl_jax/methods/collisions_methods.py | 249 ++++++++++++++++++ .../impl_jax/methods/index_methods.py | 26 ++ .../impl_jax/methods/moments_methods.py | 182 +++++++++++++ .../backends/impl_jax/methods/pair_methods.py | 91 +++++++ .../impl_jax/methods/physics_methods.py | 29 ++ PySDM/backends/impl_jax/storage.py | 93 +++++++ PySDM/backends/impl_jax/storage_impl.py | 11 + PySDM/backends/jax.py | 46 ++++ PySDM/dynamics/collisions/collision.py | 4 +- .../impl/random_generator_optimizer.py | 5 +- PySDM/impl/particle_attributes.py | 2 +- PySDM/impl/particle_attributes_factory.py | 2 +- .../products/impl/spectrum_moment_product.py | 4 +- docs/markdown/pysdm_landing.md | 6 +- 18 files changed, 754 insertions(+), 10 deletions(-) create mode 100644 PySDM/backends/impl_jax/__init__.py create mode 100644 PySDM/backends/impl_jax/methods/__init__.py create mode 100644 PySDM/backends/impl_jax/methods/collisions_methods.py create mode 100644 PySDM/backends/impl_jax/methods/index_methods.py create mode 100644 PySDM/backends/impl_jax/methods/moments_methods.py create mode 100644 PySDM/backends/impl_jax/methods/pair_methods.py create mode 100644 PySDM/backends/impl_jax/methods/physics_methods.py create mode 100644 PySDM/backends/impl_jax/storage.py create mode 100644 PySDM/backends/impl_jax/storage_impl.py create mode 100644 PySDM/backends/jax.py diff --git a/.github/workflows/readme_snippets.yml b/.github/workflows/readme_snippets.yml index 21524cb1b3..95fd801f5f 100644 --- a/.github/workflows/readme_snippets.yml +++ b/.github/workflows/readme_snippets.yml @@ -34,6 +34,8 @@ jobs: python -We readme.py sed -i -e 's/CPU/GPU/g' readme.py python -We readme.py + sed -i -e 's/GPU/JAX/g' readme.py + python -We readme.py - name: artefacts if: github.ref == 'refs/heads/main' && matrix.platform == 'ubuntu-latest' diff --git a/PySDM/backends/__init__.py b/PySDM/backends/__init__.py index 3ec1b6e492..226bf1a60e 100644 --- a/PySDM/backends/__init__.py +++ b/PySDM/backends/__init__.py @@ -11,11 +11,14 @@ from numba import cuda from . import numba as _numba +from . import jax as _jax # for pdoc CPU = None GPU = None +JAX = None Numba = _numba.Numba +Jax = _jax.Jax ThrustRTC = None @@ -93,3 +96,5 @@ def _cached_backend(formulae=None, backend_class=None, **kwargs): GPU = partial(_cached_backend, backend_class=ThrustRTC) """ returns a cached instance of the ThrustRTC backend (cache key including formulae parameters) """ + +JAX = partial(_cached_backend, backend_class=Jax) diff --git a/PySDM/backends/impl_jax/__init__.py b/PySDM/backends/impl_jax/__init__.py new file mode 100644 index 0000000000..e69de29bb2 diff --git a/PySDM/backends/impl_jax/methods/__init__.py b/PySDM/backends/impl_jax/methods/__init__.py new file mode 100644 index 0000000000..0e09071787 --- /dev/null +++ b/PySDM/backends/impl_jax/methods/__init__.py @@ -0,0 +1,7 @@ +"""method classes of the JAX backend""" + +from .collisions_methods import CollisionsMethods +from .index_methods import IndexMethods +from .moments_methods import MomentsMethods +from .pair_methods import PairMethods +from .physics_methods import PhysicsMethods diff --git a/PySDM/backends/impl_jax/methods/collisions_methods.py b/PySDM/backends/impl_jax/methods/collisions_methods.py new file mode 100644 index 0000000000..48c3a1078d --- /dev/null +++ b/PySDM/backends/impl_jax/methods/collisions_methods.py @@ -0,0 +1,249 @@ +""" +CPU implementation of backend methods for particle collisions +""" + +from functools import cached_property +import numba +import numpy as np + +from PySDM.backends.impl_common.backend_methods import BackendMethods +from PySDM.backends.impl_numba import conf +from PySDM.backends.impl_jax.storage import Storage + + +@numba.njit(**{**conf.JIT_FLAGS, **{"parallel": False}}) +def pair_indices(i, idx, is_first_in_pair, prob_like): + """given permutation array `idx` and `is_first_in_pair` flag array, + returns indices `j` and `k` of droplets within pair `i` and a `skip_pair` flag, + `j` points to the droplet that is first in pair (higher or equal multiplicity) + output is valid only if `2*i` or `2*i+1` points to a valid pair start index (within one cell) + otherwise the `skip_pair` flag is set to True and returned `j` & `k` indices are set to -1. + In addition, the `prob_like` array is checked for zeros at position `i`, in which case + the `skip_pair` is also set to `True` + """ + skip_pair = False + + if prob_like[i] == 0: + skip_pair = True + j, k = -1, -1 + else: + offset = 1 - is_first_in_pair[2 * i] + j = idx[2 * i + offset] + k = idx[2 * i + 1 + offset] + return j, k, skip_pair + + +@numba.njit(**{**conf.JIT_FLAGS, **{"parallel": False}}) +def flag_zero_multiplicity(j, k, multiplicity, healthy): + if multiplicity[k] == 0 or multiplicity[j] == 0: + healthy[0] = 0 + + +@numba.njit(**{**conf.JIT_FLAGS, **{"parallel": False}}) +def coalesce( # pylint: disable=too-many-arguments + i, j, k, cid, multiplicity, gamma, attributes, coalescence_rate +): + new_n = multiplicity[j] - gamma[i] * multiplicity[k] + if new_n > 0: + multiplicity[j] = new_n + for a in range(len(attributes)): + attributes[a, k] += gamma[i] * attributes[a, j] + else: # new_n == 0 + multiplicity[j] = multiplicity[k] // 2 + multiplicity[k] = multiplicity[k] - multiplicity[j] + for a in range(len(attributes)): + attributes[a, j] = gamma[i] * attributes[a, j] + attributes[a, k] + attributes[a, k] = attributes[a, j] + + +class CollisionsMethods(BackendMethods): + @cached_property + def _collision_coalescence_body(self): + @numba.njit(**self.default_jit_flags) + def body( + *, + multiplicity, + idx, + length, + attributes, + gamma, + healthy, + cell_id, + coalescence_rate, + is_first_in_pair, + ): + for ( + i + ) in numba.prange( # pylint: disable=not-an-iterable,too-many-nested-blocks + length // 2 + ): + j, k, skip_pair = pair_indices(i, idx, is_first_in_pair, gamma) + if skip_pair: + continue + coalesce( + i, + j, + k, + cell_id[j], + multiplicity, + gamma, + attributes, + coalescence_rate, + ) + flag_zero_multiplicity(j, k, multiplicity, healthy) + + return body + + def collision_coalescence( + self, + *, + multiplicity, + idx, + attributes, + gamma, + healthy, + cell_id, + coalescence_rate, + is_first_in_pair, + ): + self._collision_coalescence_body( + multiplicity=multiplicity.data, + idx=idx.data, + length=len(idx), + attributes=attributes.data, + gamma=gamma.data, + healthy=healthy.data, + cell_id=cell_id.data, + coalescence_rate=coalescence_rate.data, + is_first_in_pair=is_first_in_pair.indicator.data, + ) + + @cached_property + def _compute_gamma_body(self): + @numba.njit(**self.default_jit_flags) + # pylint: disable=too-many-arguments,too-many-locals + def body( + prob, + rand, + idx, + length, + multiplicity, + cell_id, + collision_rate_deficit, + collision_rate, + is_first_in_pair, + out, + ): + # TODO #1731 - shared docstring for all backends + for i in numba.prange(length // 2): # pylint: disable=not-an-iterable + out[i] = np.ceil(prob[i] - rand[i]) + j, k, skip_pair = pair_indices(i, idx, is_first_in_pair, out) + if skip_pair: + continue + prop = multiplicity[j] // multiplicity[k] + g = min(int(out[i]), prop) + out[i] = g + + return body + + def compute_gamma( + self, + *, + prob, + rand, + multiplicity, + cell_id, + collision_rate_deficit, + collision_rate, + is_first_in_pair, + out, + ): + return self._compute_gamma_body( + prob.data, + rand.data, + multiplicity.idx.data, + len(multiplicity), + multiplicity.data, + cell_id.data, + collision_rate_deficit.data, + collision_rate.data, + is_first_in_pair.indicator.data, + out.data, + ) + + @staticmethod + def make_cell_caretaker(idx_shape, idx_dtype, cell_start_len, scheme="default"): + class CellCaretaker: # pylint: disable=too-few-public-methods + def __init__(self, idx_shape, idx_dtype, cell_start_len, scheme): + if scheme == "default": + if conf.JIT_FLAGS["parallel"]: + scheme = "counting_sort_parallel" + else: + scheme = "counting_sort" + self.scheme = scheme + if scheme in ("counting_sort", "counting_sort_parallel"): + self.tmp_idx = Storage.empty(idx_shape, idx_dtype) + + def __call__(self, cell_id, cell_idx, cell_start, idx): + length = len(idx) + if self.scheme == "counting_sort": + CollisionsMethods._counting_sort_by_cell_id_and_update_cell_start( + self.tmp_idx.data, + idx.data, + cell_id.data, + cell_idx.data, + length, + cell_start.data, + ) + idx.data, self.tmp_idx.data = self.tmp_idx.data, idx.data + + return CellCaretaker(idx_shape, idx_dtype, cell_start_len, scheme) + + @cached_property + def _normalize_body(self): + @numba.njit(**{**self.default_jit_flags, **{"parallel": False}}) + # pylint: disable=too-many-arguments + def body(prob, cell_id, cell_idx, cell_start, norm_factor, timestep, dv): + n_cell = cell_start.shape[0] - 1 + for i in range(n_cell): + sd_num = cell_start[i + 1] - cell_start[i] + if sd_num < 2: + norm_factor[i] = 0 + else: + norm_factor[i] = ( + timestep / dv * sd_num * (sd_num - 1) / 2 / (sd_num // 2) + ) + for d in numba.prange(prob.shape[0]): # pylint: disable=not-an-iterable + prob[d] *= norm_factor[cell_idx[cell_id[d]]] + + return body + + # pylint: disable=too-many-arguments + def normalize(self, prob, cell_id, cell_idx, cell_start, norm_factor, timestep, dv): + return self._normalize_body( + prob.data, + cell_id.data, + cell_idx.data, + cell_start.data, + norm_factor.data, + timestep, + dv, + ) + + + @staticmethod + @numba.njit(**conf.JIT_FLAGS) + # pylint: disable=too-many-arguments + def _counting_sort_by_cell_id_and_update_cell_start( + new_idx, idx, cell_id, cell_idx, length, cell_start + ): + cell_end = cell_start + # Warning: Assuming len(cell_end) == n_cell+1 + cell_end[:] = 0 + for i in range(length): + cell_end[cell_idx[cell_id[idx[i]]]] += 1 + for i in range(1, len(cell_end)): + cell_end[i] += cell_end[i - 1] + for i in range(length - 1, -1, -1): + cell_end[cell_idx[cell_id[idx[i]]]] -= 1 + new_idx[cell_end[cell_idx[cell_id[idx[i]]]]] = idx[i] diff --git a/PySDM/backends/impl_jax/methods/index_methods.py b/PySDM/backends/impl_jax/methods/index_methods.py new file mode 100644 index 0000000000..2cbabfa458 --- /dev/null +++ b/PySDM/backends/impl_jax/methods/index_methods.py @@ -0,0 +1,26 @@ +""" +CPU implementation of shuffling and sorting backend methods +""" + +from functools import cached_property + +import numba + +from PySDM.backends.impl_common.backend_methods import BackendMethods + + +class IndexMethods(BackendMethods): + + @cached_property + def shuffle_local(self): + @numba.njit(**self.default_jit_flags) + def body(idx, u01, cell_start): + # pylint: disable=not-an-iterable + for c in numba.prange(len(cell_start) - 1): + for i in range(cell_start[c + 1] - 1, cell_start[c], -1): + j = int( + cell_start[c] + u01[i] * (cell_start[c + 1] - cell_start[c]) + ) + idx[i], idx[j] = idx[j], idx[i] + + return body diff --git a/PySDM/backends/impl_jax/methods/moments_methods.py b/PySDM/backends/impl_jax/methods/moments_methods.py new file mode 100644 index 0000000000..29cfba19ac --- /dev/null +++ b/PySDM/backends/impl_jax/methods/moments_methods.py @@ -0,0 +1,182 @@ +""" +CPU implementation of moment calculation backend methods +""" + +from functools import cached_property + +import numba + +from PySDM.backends.impl_common.backend_methods import BackendMethods +from PySDM.backends.impl_numba.atomic_operations import atomic_add + + +class MomentsMethods(BackendMethods): + @cached_property + def _moments_body(self): + @numba.njit(**self.default_jit_flags) + def body( + *, + moment_0, + moments, + multiplicity, + attr_data, + cell_id, + idx, + length, + ranks, + min_x, + max_x, + x_attr, + weighting_attribute, + weighting_rank, + skip_division_by_m0, + ): + # pylint: disable=too-many-locals + moment_0[:] = 0 + moments[:, :] = 0 + for idx_i in numba.prange(length): # pylint: disable=not-an-iterable + i = idx[idx_i] + if min_x <= x_attr[i] < max_x: + atomic_add( + moment_0, + cell_id[i], + multiplicity[i] * weighting_attribute[i] ** weighting_rank, + ) + for k in range(ranks.shape[0]): + atomic_add( + moments, + (k, cell_id[i]), + ( + multiplicity[i] + * weighting_attribute[i] ** weighting_rank + * attr_data[i] ** ranks[k] + ), + ) + if not skip_division_by_m0: + for c_id in range(moment_0.shape[0]): + for k in range(ranks.shape[0]): + moments[k, c_id] = ( + moments[k, c_id] / moment_0[c_id] + if moment_0[c_id] != 0 + else 0 + ) + + return body + + def moments( + self, + *, + moment_0, + moments, + multiplicity, + attr_data, + cell_id, + idx, + length, + ranks, + min_x, + max_x, + x_attr, + weighting_attribute, + weighting_rank, + skip_division_by_m0, + ): + return self._moments_body( + moment_0=moment_0.data, + moments=moments.data, + multiplicity=multiplicity.data, + attr_data=attr_data.data, + cell_id=cell_id.data, + idx=idx.data, + length=length, + ranks=ranks.data, + min_x=min_x, + max_x=max_x, + x_attr=x_attr.data, + weighting_attribute=weighting_attribute.data, + weighting_rank=weighting_rank, + skip_division_by_m0=skip_division_by_m0, + ) + + @cached_property + def _spectrum_moments_body(self): + @numba.njit(**self.default_jit_flags) + def body( + *, + moment_0, + moments, + multiplicity, + attr_data, + cell_id, + idx, + length, + rank, + x_bins, + x_attr, + weighting_attribute, + weighting_rank, + ): + # pylint: disable=too-many-locals + moment_0[:, :] = 0 + moments[:, :] = 0 + for idx_i in numba.prange(length): # pylint: disable=not-an-iterable + i = idx[idx_i] + for k in range(x_bins.shape[0] - 1): + if x_bins[k] <= x_attr[i] < x_bins[k + 1]: + atomic_add( + moment_0, + (k, cell_id[i]), + multiplicity[i] * weighting_attribute[i] ** weighting_rank, + ) + atomic_add( + moments, + (k, cell_id[i]), + ( + multiplicity[i] + * weighting_attribute[i] ** weighting_rank + * attr_data[i] ** rank + ), + ) + break + for c_id in range(moment_0.shape[1]): + for k in range(x_bins.shape[0] - 1): + moments[k, c_id] = ( + moments[k, c_id] / moment_0[k, c_id] + if moment_0[k, c_id] != 0 + else 0 + ) + + return body + + def spectrum_moments( + self, + *, + moment_0, + moments, + multiplicity, + attr_data, + cell_id, + idx, + length, + rank, + x_bins, + x_attr, + weighting_attribute, + weighting_rank, + ): + assert moments.shape[0] == x_bins.shape[0] - 1 + assert moment_0.shape == moments.shape + return self._spectrum_moments_body( + moment_0=moment_0.data, + moments=moments.data, + multiplicity=multiplicity.data, + attr_data=attr_data.data, + cell_id=cell_id.data, + idx=idx.data, + length=length, + rank=rank, + x_bins=x_bins.data, + x_attr=x_attr.data, + weighting_attribute=weighting_attribute.data, + weighting_rank=weighting_rank, + ) diff --git a/PySDM/backends/impl_jax/methods/pair_methods.py b/PySDM/backends/impl_jax/methods/pair_methods.py new file mode 100644 index 0000000000..cc187804cd --- /dev/null +++ b/PySDM/backends/impl_jax/methods/pair_methods.py @@ -0,0 +1,91 @@ +""" +CPU implementation of pairwise operations backend methods +""" + +from functools import cached_property + +import numba + +from PySDM.backends.impl_common.backend_methods import BackendMethods + + +class PairMethods(BackendMethods): + + @cached_property + def _find_pairs_body(self): + @numba.njit(**self.default_jit_flags) + def body(*, cell_start, is_first_in_pair, cell_id, cell_idx, idx, length): + for i in numba.prange(length - 1): # pylint: disable=not-an-iterable + is_in_same_cell = cell_id[idx[i]] == cell_id[idx[i + 1]] + is_even_index = (i - cell_start[cell_idx[cell_id[idx[i]]]]) % 2 == 0 + is_first_in_pair[i] = is_in_same_cell and is_even_index + is_first_in_pair[length - 1] = False + + return body + + # pylint: disable=too-many-arguments + def find_pairs(self, cell_start, is_first_in_pair, cell_id, cell_idx, idx): + return self._find_pairs_body( + cell_start=cell_start.data, + is_first_in_pair=is_first_in_pair.indicator.data, + cell_id=cell_id.data, + cell_idx=cell_idx.data, + idx=idx.data, + length=len(idx), + ) + + @cached_property + def _max_pair_body(self): + @numba.njit(**self.default_jit_flags) + def body(data_out, data_in, is_first_in_pair, idx, length): + data_out[:] = 0 + for i in numba.prange(length - 1): # pylint: disable=not-an-iterable + if is_first_in_pair[i]: + data_out[i // 2] = max(data_in[idx[i]], data_in[idx[i + 1]]) + + return body + + def max_pair(self, data_out, data_in, is_first_in_pair, idx): + return self._max_pair_body( + data_out.data, + data_in.data, + is_first_in_pair.indicator.data, + idx.data, + len(idx), + ) + + @cached_property + def _sort_within_pair_by_attr_body(self): + @numba.njit(**self.default_jit_flags) + def body(idx, length, is_first_in_pair, attr): + for i in numba.prange(length - 1): # pylint: disable=not-an-iterable + if is_first_in_pair[i]: + if attr[idx[i]] < attr[idx[i + 1]]: + idx[i], idx[i + 1] = idx[i + 1], idx[i] + + return body + + def sort_within_pair_by_attr(self, idx, is_first_in_pair, attr): + self._sort_within_pair_by_attr_body( + idx.data, len(idx), is_first_in_pair.indicator.data, attr.data + ) + + @cached_property + def _sum_pair_body(self): + @numba.njit(**self.default_jit_flags) + def body(data_out, data_in, is_first_in_pair, idx, length): + data_out[:] = 0 + for i in numba.prange(length): # pylint: disable=not-an-iterable + if is_first_in_pair[i]: + data_out[i // 2] = data_in[idx[i]] + data_in[idx[i + 1]] + + return body + + def sum_pair(self, data_out, data_in, is_first_in_pair, idx): + return self._sum_pair_body( + data_out.data, + data_in.data, + is_first_in_pair.indicator.data, + idx.data, + len(idx), + ) diff --git a/PySDM/backends/impl_jax/methods/physics_methods.py b/PySDM/backends/impl_jax/methods/physics_methods.py new file mode 100644 index 0000000000..de80965434 --- /dev/null +++ b/PySDM/backends/impl_jax/methods/physics_methods.py @@ -0,0 +1,29 @@ +""" +CPU implementation of backend methods wrapping basic physics formulae +""" + +from functools import cached_property + +import numba +from numba import prange + +from PySDM.backends.impl_common.backend_methods import BackendMethods + + +class PhysicsMethods(BackendMethods): + def __init__(self): + BackendMethods.__init__(self) + + @cached_property + def _volume_of_mass_body(self): + ff = self.formulae_flattened + + @numba.njit(**self.default_jit_flags) + def body(volume, mass): + for i in prange(volume.shape[0]): # pylint: disable=not-an-iterable + volume[i] = ff.particle_shape_and_density__mass_to_volume(mass[i]) + + return body + + def volume_of_water_mass(self, volume, mass): + self._volume_of_mass_body(volume.data, mass.data) diff --git a/PySDM/backends/impl_jax/storage.py b/PySDM/backends/impl_jax/storage.py new file mode 100644 index 0000000000..340a41ca6b --- /dev/null +++ b/PySDM/backends/impl_jax/storage.py @@ -0,0 +1,93 @@ +""" +CPU Numpy-based implementation of Storage class +""" + +import numpy as np + +from PySDM.backends.impl_common.storage_utils import ( + StorageBase, + StorageSignature, + empty, + get_data_from_ndarray, +) +from PySDM.backends.impl_jax import storage_impl as impl + + +class Storage(StorageBase): + FLOAT = np.float64 + INT = np.int64 + BOOL = np.bool_ + + def row_view(self, i): + return Storage( + StorageSignature(self.data[i], (*self.shape[1:],), self.dtype) + ) + + def at(self, index): + assert self.shape == (1,), "Cannot call at() on Storage of shape other than (1,)" + return self.data[index] + + def __imul__(self, other): + if hasattr(other, "data"): + impl.multiply(self.data, other.data) + else: + impl.multiply(self.data, other) + return self + + def __itruediv__(self, other): + if hasattr(other, "data"): + self.data[:] /= other.data[:] + else: + self.data[:] /= other + return self + + def download(self, target, reshape=False): + if reshape: + data = self.data.reshape(target.shape) + else: + data = self.data + np.copyto(target, data, casting="safe") + + @staticmethod + def _get_empty_data(shape, dtype): + if dtype in (float, Storage.FLOAT): + data = np.full(shape, np.nan, dtype=Storage.FLOAT) + dtype = Storage.FLOAT + elif dtype in (int, Storage.INT): + data = np.full(shape, -1, dtype=Storage.INT) + dtype = Storage.INT + elif dtype in (bool, Storage.BOOL): + data = np.full(shape, -1, dtype=Storage.BOOL) + dtype = Storage.BOOL + else: + raise NotImplementedError() + + return StorageSignature(data, shape, dtype) + + @staticmethod + def empty(shape, dtype): + return empty(shape, dtype, Storage) + + @staticmethod + def _get_data_from_ndarray(array): + return get_data_from_ndarray( + array=array, + storage_class=Storage, + copy_fun=lambda array_astype: array_astype.copy(), + ) + + @staticmethod + def from_ndarray(array): + result = Storage(Storage._get_data_from_ndarray(array)) + return result + + def urand(self, generator): + generator(self) + + def upload(self, data): + np.copyto(self.data, data, casting="safe") + def fill(self, other): + if isinstance(other, Storage): + self.data[:] = other.data + else: + self.data[:] = other diff --git a/PySDM/backends/impl_jax/storage_impl.py b/PySDM/backends/impl_jax/storage_impl.py new file mode 100644 index 0000000000..fe9d492706 --- /dev/null +++ b/PySDM/backends/impl_jax/storage_impl.py @@ -0,0 +1,11 @@ +""" +Numba njit-ted basic arithmetics routines for CPU backend +""" + +import numba +from PySDM.backends.impl_numba import conf + + +@numba.njit(**{**conf.JIT_FLAGS, **{"parallel": False}}) +def multiply(output, multiplier): + output *= multiplier diff --git a/PySDM/backends/jax.py b/PySDM/backends/jax.py new file mode 100644 index 0000000000..22f25dbea1 --- /dev/null +++ b/PySDM/backends/jax.py @@ -0,0 +1,46 @@ +""" +Multi-threaded CPU backend using LLVM-powered just-in-time compilation +""" + +from PySDM.backends.impl_jax import methods +from PySDM.backends.impl_numba.random import Random as ImportedRandom +from PySDM.backends.impl_jax.storage import Storage as ImportedStorage +from PySDM.formulae import Formulae + + +class Jax( + methods.CollisionsMethods, + methods.PairMethods, + methods.IndexMethods, + methods.PhysicsMethods, + methods.MomentsMethods, +): + Storage = ImportedStorage + Random = ImportedRandom + + default_croupier = "local" + + def __init__( + self, formulae=None, *, double_precision=True, override_jit_flags=None + ): + if not double_precision: + raise NotImplementedError() + + self.formulae = formulae or Formulae() + self.formulae_flattened = self.formulae.flatten + + # assert "fastmath" not in (override_jit_flags or {}) + # self.default_jit_flags = { + # **JIT_FLAGS, # here parallel=False (for out-of-backend code) + # **{"fastmath": self.formulae.fastmath, "parallel": parallel_default}, + # **(override_jit_flags or {}), + # } + self.default_jit_flags = { + "parallel": False + } + + methods.CollisionsMethods.__init__(self) + methods.PairMethods.__init__(self) + methods.IndexMethods.__init__(self) + methods.PhysicsMethods.__init__(self) + methods.MomentsMethods.__init__(self) diff --git a/PySDM/dynamics/collisions/collision.py b/PySDM/dynamics/collisions/collision.py index 18b24d8637..97af236a37 100644 --- a/PySDM/dynamics/collisions/collision.py +++ b/PySDM/dynamics/collisions/collision.py @@ -131,9 +131,9 @@ def register(self, builder): self.stats_n_substep = self.particulator.Storage.empty( self.particulator.mesh.n_cell, dtype=int ) - self.stats_n_substep[:] = 0 if self.adaptive else self.__substeps + self.stats_n_substep.fill(0 if self.adaptive else self.__substeps) self.stats_dt_min = self.particulator.Storage.empty(**empty_args_cellwise) - self.stats_dt_min[:] = np.nan + self.stats_dt_min.fill(np.nan) self.rnd_opt_coll.register(builder) self.collision_kernel.register(builder) diff --git a/PySDM/dynamics/impl/random_generator_optimizer.py b/PySDM/dynamics/impl/random_generator_optimizer.py index 4e2bb54094..793897d219 100644 --- a/PySDM/dynamics/impl/random_generator_optimizer.py +++ b/PySDM/dynamics/impl/random_generator_optimizer.py @@ -45,4 +45,7 @@ def get_random_arrays(self): self.pairs_rand.urand(self.rnd) self.rand.urand(self.rnd) self.substep += 1 - return self.pairs_rand[shift : self.particulator.n_sd + shift], self.rand + if self.optimized_random: + return self.pairs_rand[shift : self.particulator.n_sd + shift], self.rand + else: + return self.pairs_rand, self.rand diff --git a/PySDM/impl/particle_attributes.py b/PySDM/impl/particle_attributes.py index 5e486db29f..d91038ffd3 100644 --- a/PySDM/impl/particle_attributes.py +++ b/PySDM/impl/particle_attributes.py @@ -41,7 +41,7 @@ def __init__( @property def healthy(self) -> bool: - return bool(self.__healthy_memory[0]) + return bool(self.__healthy_memory.at(0)) @healthy.setter def healthy(self, value: bool): diff --git a/PySDM/impl/particle_attributes_factory.py b/PySDM/impl/particle_attributes_factory.py index f8f1b4546b..e5dad088c2 100644 --- a/PySDM/impl/particle_attributes_factory.py +++ b/PySDM/impl/particle_attributes_factory.py @@ -61,7 +61,7 @@ def attributes(particulator, req_attr, attributes): def helper(req_attr, all_attr, names, data, keys): for i, attr in enumerate(names): keys[attr] = i - req_attr[attr].set_data(data[i, :]) + req_attr[attr].set_data(data.row_view(i)) try: req_attr[attr].init(all_attr[attr]) except KeyError as err: diff --git a/PySDM/products/impl/spectrum_moment_product.py b/PySDM/products/impl/spectrum_moment_product.py index b02fe781d5..3762a9a507 100644 --- a/PySDM/products/impl/spectrum_moment_product.py +++ b/PySDM/products/impl/spectrum_moment_product.py @@ -48,6 +48,6 @@ def _recalculate_spectrum_moment( def _download_spectrum_moment_to_buffer(self, rank, bin_number): if rank == 0: # TODO #217 - self._download_to_buffer(self.moment_0[bin_number, :]) + self._download_to_buffer(self.moment_0.row_view(bin_number)) else: - self._download_to_buffer(self.moments[bin_number, :]) + self._download_to_buffer(self.moments.row_view(bin_number)) diff --git a/docs/markdown/pysdm_landing.md b/docs/markdown/pysdm_landing.md index d373d487f7..ef5fa29c28 100644 --- a/docs/markdown/pysdm_landing.md +++ b/docs/markdown/pysdm_landing.md @@ -157,7 +157,7 @@ radius_bins_edges = 10 .^ range(log10(10*si.um), log10(5e3*si.um), length=32) env = Box(dt=1 * si.s, dv=1e6 * si.m^3) builder = Builder(n_sd=n_sd, backend=CPU(), environment=env) -builder.add_dynamic(Coalescence(collision_kernel=Golovin(b=1.5e3 / si.s))) +builder.add_dynamic(Coalescence(collision_kernel=Golovin(b=1.5e3 / si.s), adaptive=false)) products = [ParticleVolumeVersusRadiusLogarithmSpectrum(radius_bins_edges=radius_bins_edges, name="dv/dlnr")] particulator = builder.build(attributes, products) ``` @@ -177,7 +177,7 @@ radius_bins_edges = logspace(log10(10 * si.um), log10(5e3 * si.um), 32); env = Box(pyargs('dt', 1 * si.s, 'dv', 1e6 * si.m ^ 3)); builder = Builder(pyargs('n_sd', int32(n_sd), 'backend', CPU(), 'environment', env)); -builder.add_dynamic(Coalescence(pyargs('collision_kernel', Golovin(1.5e3 / si.s)))); +builder.add_dynamic(Coalescence(pyargs('collision_kernel', Golovin(1.5e3 / si.s), 'adaptive', false))); products = py.list({ ParticleVolumeVersusRadiusLogarithmSpectrum(pyargs( ... 'radius_bins_edges', py.numpy.array(radius_bins_edges), ... 'name', 'dv/dlnr' ... @@ -201,7 +201,7 @@ radius_bins_edges = np.logspace(np.log10(10 * si.um), np.log10(5e3 * si.um), num env = Box(dt=1 * si.s, dv=1e6 * si.m ** 3) builder = Builder(n_sd=n_sd, backend=CPU(), environment=env) -builder.add_dynamic(Coalescence(collision_kernel=Golovin(b=1.5e3 / si.s))) +builder.add_dynamic(Coalescence(collision_kernel=Golovin(b=1.5e3 / si.s), adaptive=False)) products = [ParticleVolumeVersusRadiusLogarithmSpectrum(radius_bins_edges=radius_bins_edges, name='dv/dlnr')] particulator = builder.build(attributes, products) ``` From 58e7f519434d600e0632268e732115b070eb21c5 Mon Sep 17 00:00:00 2001 From: Konrad Bodzioch Date: Tue, 21 Oct 2025 12:53:43 +0200 Subject: [PATCH 2/4] readme: remove unimplemented readme examples for JAX --- .github/workflows/readme_snippets.yml | 1 + 1 file changed, 1 insertion(+) diff --git a/.github/workflows/readme_snippets.yml b/.github/workflows/readme_snippets.yml index 95fd801f5f..87b2db2257 100644 --- a/.github/workflows/readme_snippets.yml +++ b/.github/workflows/readme_snippets.yml @@ -35,6 +35,7 @@ jobs: sed -i -e 's/CPU/GPU/g' readme.py python -We readme.py sed -i -e 's/GPU/JAX/g' readme.py + sed -i '/pyplot\.legend()/,$d' readme.py python -We readme.py - name: artefacts From 3e539f393f8bddf72125af27370780916fb9f485 Mon Sep 17 00:00:00 2001 From: Konrad Bodzioch Date: Tue, 16 Dec 2025 19:15:13 +0100 Subject: [PATCH 3/4] WIP --- PySDM/backends/impl_common/storage_utils.py | 4 + .../impl_jax/methods/moments_methods.py | 32 ++- .../impl_jax/methods/physics_methods.py | 26 +- PySDM/backends/impl_jax/storage.py | 49 +++- .../impl_numba/methods/moments_methods.py | 6 + PySDM/backends/impl_numba/storage.py | 10 + PySDM/backends/jax.py | 2 + PySDM/particulator.py | 2 + jax_test.ipynb | 226 ++++++++++++++++++ readme.png | Bin 0 -> 14636 bytes 10 files changed, 332 insertions(+), 25 deletions(-) create mode 100644 jax_test.ipynb create mode 100644 readme.png diff --git a/PySDM/backends/impl_common/storage_utils.py b/PySDM/backends/impl_common/storage_utils.py index 3b24e6be7d..83fde0c2c0 100644 --- a/PySDM/backends/impl_common/storage_utils.py +++ b/PySDM/backends/impl_common/storage_utils.py @@ -58,6 +58,10 @@ def upload(self, data): @abstractmethod def fill(self, other): raise NotImplementedError() + + @abstractmethod + def rowview(self, i): + raise NotImplementedError() def get_data_from_ndarray(array, storage_class: Type[StorageBase], copy_fun): diff --git a/PySDM/backends/impl_jax/methods/moments_methods.py b/PySDM/backends/impl_jax/methods/moments_methods.py index 29cfba19ac..147191db2e 100644 --- a/PySDM/backends/impl_jax/methods/moments_methods.py +++ b/PySDM/backends/impl_jax/methods/moments_methods.py @@ -3,8 +3,10 @@ """ from functools import cached_property +import time +import jax -import numba +# import numba from PySDM.backends.impl_common.backend_methods import BackendMethods from PySDM.backends.impl_numba.atomic_operations import atomic_add @@ -13,7 +15,7 @@ class MomentsMethods(BackendMethods): @cached_property def _moments_body(self): - @numba.njit(**self.default_jit_flags) + # @numba.njit(**self.default_jit_flags) def body( *, moment_0, @@ -34,7 +36,7 @@ def body( # pylint: disable=too-many-locals moment_0[:] = 0 moments[:, :] = 0 - for idx_i in numba.prange(length): # pylint: disable=not-an-iterable + for idx_i in range(length): # pylint: disable=not-an-iterable i = idx[idx_i] if min_x <= x_attr[i] < max_x: atomic_add( @@ -99,8 +101,9 @@ def moments( ) @cached_property + # @jax.jit def _spectrum_moments_body(self): - @numba.njit(**self.default_jit_flags) + # @numba.njit(**self.default_jit_flags) def body( *, moment_0, @@ -116,13 +119,22 @@ def body( weighting_attribute, weighting_rank, ): + t1 = time.time() # pylint: disable=too-many-locals - moment_0[:, :] = 0 - moments[:, :] = 0 - for idx_i in numba.prange(length): # pylint: disable=not-an-iterable + moment_0.at[:, :].set(0) + moments.at[:, :].set(0) + # print(f"initial moments setter: {time.time() - t1}") + # print(length) + for idx_i in range(length): # pylint: disable=not-an-iterable i = idx[idx_i] + t4 = time.time() + # print(f"x_bins type: {type(x_bins)}") + # print(f"x_attr type: {type(x_attr)}") for k in range(x_bins.shape[0] - 1): + # print(f"x_attr[i]: {x_attr[i]}") if x_bins[k] <= x_attr[i] < x_bins[k + 1]: + # if x_bins[k] != x_bins[k + 1] and x_bins[k] == x_bins[k+1]: + t2 = time.time() atomic_add( moment_0, (k, cell_id[i]), @@ -137,14 +149,18 @@ def body( * attr_data[i] ** rank ), ) + # print(f"adding... : {time.time() - t2}") break + # print(f"oop: {time.time() - t4}") + # print(f"After loop 1 checkpoint: {time.time() - t1}") for c_id in range(moment_0.shape[1]): for k in range(x_bins.shape[0] - 1): - moments[k, c_id] = ( + moments.at[k, c_id].set( moments[k, c_id] / moment_0[k, c_id] if moment_0[k, c_id] != 0 else 0 ) + # print(f"spectrum moments body: {time.time() - t1}") return body diff --git a/PySDM/backends/impl_jax/methods/physics_methods.py b/PySDM/backends/impl_jax/methods/physics_methods.py index de80965434..75b98c9af8 100644 --- a/PySDM/backends/impl_jax/methods/physics_methods.py +++ b/PySDM/backends/impl_jax/methods/physics_methods.py @@ -2,10 +2,12 @@ CPU implementation of backend methods wrapping basic physics formulae """ -from functools import cached_property +from functools import cached_property, partial +import time -import numba -from numba import prange +import jax +# import numba +# from numba import prange from PySDM.backends.impl_common.backend_methods import BackendMethods @@ -18,12 +20,22 @@ def __init__(self): def _volume_of_mass_body(self): ff = self.formulae_flattened - @numba.njit(**self.default_jit_flags) + # @numba.njit(**self.default_jit_flags) + # @jax.jit def body(volume, mass): - for i in prange(volume.shape[0]): # pylint: disable=not-an-iterable - volume[i] = ff.particle_shape_and_density__mass_to_volume(mass[i]) + t1 = time.time() + for i in range(volume.shape[0]): # pylint: disable=not-an-iterable + # volume.at[i].set(ff.particle_shape_and_density__mass_to_volume.py_func(mass[i])) + volume = volume.at[i].set(ff.particle_shape_and_density__mass_to_volume(mass[i])) + # print(f"Mass[i]: {mass[i]}, Volume: {volume[i]}") + print("Post loop: ", volume) + + return volume + # print(f"volume of water mass run time: {time.time() - t1}") return body def volume_of_water_mass(self, volume, mass): - self._volume_of_mass_body(volume.data, mass.data) + print("Pre volume: ", volume.data) + volume.data = self._volume_of_mass_body(volume.data, mass.data) + print("Post volume: ", volume.data) diff --git a/PySDM/backends/impl_jax/storage.py b/PySDM/backends/impl_jax/storage.py index 340a41ca6b..3f8b6ca6f8 100644 --- a/PySDM/backends/impl_jax/storage.py +++ b/PySDM/backends/impl_jax/storage.py @@ -2,7 +2,10 @@ CPU Numpy-based implementation of Storage class """ +import jax.numpy as jnp import numpy as np +import jax +import time from PySDM.backends.impl_common.storage_utils import ( StorageBase, @@ -14,20 +17,23 @@ class Storage(StorageBase): - FLOAT = np.float64 - INT = np.int64 - BOOL = np.bool_ + FLOAT = jnp.float64 + INT = jnp.int64 + BOOL = jnp.bool_ def row_view(self, i): + print('rowview') return Storage( StorageSignature(self.data[i], (*self.shape[1:],), self.dtype) ) def at(self, index): + print('at') assert self.shape == (1,), "Cannot call at() on Storage of shape other than (1,)" return self.data[index] def __imul__(self, other): + print('imul') if hasattr(other, "data"): impl.multiply(self.data, other.data) else: @@ -35,6 +41,7 @@ def __imul__(self, other): return self def __itruediv__(self, other): + print('itruediv') if hasattr(other, "data"): self.data[:] /= other.data[:] else: @@ -42,34 +49,44 @@ def __itruediv__(self, other): return self def download(self, target, reshape=False): + # print('download') + t1 = time.time() if reshape: data = self.data.reshape(target.shape) else: data = self.data - np.copyto(target, data, casting="safe") + target = jnp.asarray(data) + print(f'download: {time.time() - t1}') + # np.copyto(target, np.asarray(data), casting="safe") @staticmethod def _get_empty_data(shape, dtype): + # print('get_empty_data') + t1 = time.time() if dtype in (float, Storage.FLOAT): - data = np.full(shape, np.nan, dtype=Storage.FLOAT) + data = jnp.full(shape, jnp.nan, dtype=Storage.FLOAT) dtype = Storage.FLOAT elif dtype in (int, Storage.INT): - data = np.full(shape, -1, dtype=Storage.INT) + data = jnp.full(shape, -1, dtype=Storage.INT) dtype = Storage.INT elif dtype in (bool, Storage.BOOL): - data = np.full(shape, -1, dtype=Storage.BOOL) + data = jnp.full(shape, -1, dtype=Storage.BOOL) dtype = Storage.BOOL else: raise NotImplementedError() + print(f'get_empty_data: {time.time() - t1}') return StorageSignature(data, shape, dtype) @staticmethod def empty(shape, dtype): + print('empty') return empty(shape, dtype, Storage) @staticmethod def _get_data_from_ndarray(array): + print('get_data_from_ndarray') + return get_data_from_ndarray( array=array, storage_class=Storage, @@ -78,16 +95,28 @@ def _get_data_from_ndarray(array): @staticmethod def from_ndarray(array): + print('from_ndarray') result = Storage(Storage._get_data_from_ndarray(array)) return result def urand(self, generator): + print('urand') generator(self) + def upload(self, data): - np.copyto(self.data, data, casting="safe") + print('upload') + self.fill(data) + def fill(self, other): + t1 = time.time() + print('fill') if isinstance(other, Storage): - self.data[:] = other.data + # self.data[:] = other.data + # self.data.at[:].set(other.data) + self.data = other.data else: - self.data[:] = other + # self.data[:] = other + # self.data.at[:].set(other) + self.data = other + print(f'fill: {time.time() - t1}') diff --git a/PySDM/backends/impl_numba/methods/moments_methods.py b/PySDM/backends/impl_numba/methods/moments_methods.py index 29cfba19ac..b10368c009 100644 --- a/PySDM/backends/impl_numba/methods/moments_methods.py +++ b/PySDM/backends/impl_numba/methods/moments_methods.py @@ -5,6 +5,7 @@ from functools import cached_property import numba +import time from PySDM.backends.impl_common.backend_methods import BackendMethods from PySDM.backends.impl_numba.atomic_operations import atomic_add @@ -98,6 +99,7 @@ def moments( skip_division_by_m0=skip_division_by_m0, ) + @cached_property def _spectrum_moments_body(self): @numba.njit(**self.default_jit_flags) @@ -119,10 +121,13 @@ def body( # pylint: disable=too-many-locals moment_0[:, :] = 0 moments[:, :] = 0 + print(length) for idx_i in numba.prange(length): # pylint: disable=not-an-iterable i = idx[idx_i] + t4 = time.time() for k in range(x_bins.shape[0] - 1): if x_bins[k] <= x_attr[i] < x_bins[k + 1]: + # print("adagio") atomic_add( moment_0, (k, cell_id[i]), @@ -138,6 +143,7 @@ def body( ), ) break + print(f"oop: {time.time() - t4}") for c_id in range(moment_0.shape[1]): for k in range(x_bins.shape[0] - 1): moments[k, c_id] = ( diff --git a/PySDM/backends/impl_numba/storage.py b/PySDM/backends/impl_numba/storage.py index 64b2291a14..c48ceb8480 100644 --- a/PySDM/backends/impl_numba/storage.py +++ b/PySDM/backends/impl_numba/storage.py @@ -211,3 +211,13 @@ def exp(self): def abs(self): self.data[:] = np.abs(self.data) + + + def row_view(self, i): + return Storage( + StorageSignature(self.data[i], (*self.shape[1:],), self.dtype) + ) + + def at(self, index): + assert self.shape == (1,), "Cannot call at() on Storage of shape other than (1,)" + return self.data[index] \ No newline at end of file diff --git a/PySDM/backends/jax.py b/PySDM/backends/jax.py index 22f25dbea1..aef3dbf84f 100644 --- a/PySDM/backends/jax.py +++ b/PySDM/backends/jax.py @@ -1,6 +1,7 @@ """ Multi-threaded CPU backend using LLVM-powered just-in-time compilation """ +import jax from PySDM.backends.impl_jax import methods from PySDM.backends.impl_numba.random import Random as ImportedRandom @@ -23,6 +24,7 @@ class Jax( def __init__( self, formulae=None, *, double_precision=True, override_jit_flags=None ): + jax.config.update('jax_enable_x64', True) if not double_precision: raise NotImplementedError() diff --git a/PySDM/particulator.py b/PySDM/particulator.py index 5919feb805..a7948912c3 100644 --- a/PySDM/particulator.py +++ b/PySDM/particulator.py @@ -351,6 +351,7 @@ def moments( ranks = self.backend.Storage.from_ndarray(np.array(ranks, dtype=float)) + # print("Attr name: " + attr_name) self.backend.moments( moment_0=moment_0, moments=moments, @@ -380,6 +381,7 @@ def spectrum_moments( weighting_attribute="water mass", weighting_rank=0, ): + print("Attr name: " + attr_name) attr_data = self.attributes[attr] self.backend.spectrum_moments( moment_0=moment_0, diff --git a/jax_test.ipynb b/jax_test.ipynb new file mode 100644 index 0000000000..c2a078f3cb --- /dev/null +++ b/jax_test.ipynb @@ -0,0 +1,226 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 1, + "id": "4c6293f6", + "metadata": {}, + "outputs": [], + "source": [ + "import os\n", + "os.environ['NUMBA_DISABLE_JIT'] = \"1\"\n", + "\n", + "from PySDM.physics import si\n", + "from PySDM.initialisation.sampling.spectral_sampling import ConstantMultiplicity\n", + "from PySDM.initialisation.spectra.exponential import Exponential\n", + "\n", + "n_sd = 2 ** 15\n", + "initial_spectrum = Exponential(norm_factor=8.39e12, scale=1.19e5 * si.um ** 3)\n", + "attributes = {}\n", + "attributes['volume'], attributes['multiplicity'] = ConstantMultiplicity(initial_spectrum).sample(n_sd)" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "fbb45435", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "empty\n", + "get_empty_data: 0.051486968994140625\n", + "get_empty_data: 0.01377105712890625\n", + "empty\n", + "get_empty_data: 0.004712820053100586\n", + "get_empty_data: 7.104873657226562e-05\n", + "empty\n", + "get_empty_data: 0.014365196228027344\n", + "empty\n", + "get_empty_data: 7.319450378417969e-05\n", + "empty\n", + "get_empty_data: 0.007919073104858398\n", + "fill\n", + "fill: 5.0067901611328125e-06\n", + "empty\n", + "get_empty_data: 7.295608520507812e-05\n", + "fill\n", + "fill: 5.0067901611328125e-06\n", + "empty\n", + "get_empty_data: 0.011842012405395508\n", + "empty\n", + "get_empty_data: 6.699562072753906e-05\n", + "from_ndarray\n", + "get_data_from_ndarray\n", + "from_ndarray\n", + "get_data_from_ndarray\n", + "from_ndarray\n", + "get_data_from_ndarray\n", + "from_ndarray\n", + "get_data_from_ndarray\n", + "empty\n", + "get_empty_data: 0.012073993682861328\n", + "empty\n", + "get_empty_data: 5.817413330078125e-05\n", + "get_data_from_ndarray\n", + "empty\n", + "get_empty_data: 0.010919809341430664\n", + "empty\n", + "get_empty_data: 0.004091978073120117\n", + "empty\n", + "get_empty_data: 6.985664367675781e-05\n", + "rowview\n", + "upload\n", + "fill\n", + "fill: 5.9604644775390625e-06\n", + "empty\n", + "get_empty_data: 0.011930227279663086\n", + "upload\n", + "fill\n", + "fill: 5.0067901611328125e-06\n", + "empty\n", + "get_empty_data: 8.20159912109375e-05\n", + "upload\n", + "fill\n", + "fill: 3.814697265625e-06\n", + "from_ndarray\n", + "get_data_from_ndarray\n", + "get_data_from_ndarray\n", + "from_ndarray\n", + "get_data_from_ndarray\n", + "empty\n", + "get_empty_data: 4.8160552978515625e-05\n", + "dict_keys(['multiplicity', 'signed water mass', 'water mass', 'cell id', 'volume'])\n", + "Pre volume: [nan nan nan ... nan nan nan]\n", + "Post loop: [3.00579754e-18 6.63746385e-18 1.02692410e-17 ... 1.11357449e-12\n", + " 1.16550703e-12 1.25977547e-12]\n", + "Post volume: [3.00579754e-18 6.63746385e-18 1.02692410e-17 ... 1.11357449e-12\n", + " 1.16550703e-12 1.25977547e-12]\n", + "[3.00579754e-18 6.63746385e-18 1.02692410e-17 ... 1.11357449e-12\n", + " 1.16550703e-12 1.25977547e-12]\n" + ] + } + ], + "source": [ + "import numpy as np\n", + "from PySDM import Builder\n", + "from PySDM.environments import Box\n", + "from PySDM.dynamics import Coalescence\n", + "from PySDM.dynamics.collisions.collision_kernels import Golovin\n", + "from PySDM.backends import JAX, CPU\n", + "from PySDM.products import ParticleVolumeVersusRadiusLogarithmSpectrum\n", + "\n", + "radius_bins_edges = np.logspace(np.log10(10 * si.um), np.log10(5e3 * si.um), num=32)\n", + "\n", + "env = Box(dt=1 * si.s, dv=1e6 * si.m ** 3)\n", + "builder = Builder(n_sd=n_sd, backend=JAX(), environment=env)\n", + "builder.add_dynamic(Coalescence(collision_kernel=Golovin(b=1.5e3 / si.s), adaptive=False))\n", + "products = [ParticleVolumeVersusRadiusLogarithmSpectrum(radius_bins_edges=radius_bins_edges, name='dv/dlnr')]\n", + "particulator = builder.build(attributes, products)\n", + "\n", + "print(particulator.attributes.keys())\n", + "print(particulator.attributes['volume'].data)" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "id": "3c3564c8", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Attr name: volume\n", + "at\n" + ] + }, + { + "ename": "TypeError", + "evalue": "JAX arrays are immutable and do not support in-place item assignment. Instead of x[idx] = y, use x = x.at[idx].set(y) or another .at[] method: https://docs.jax.dev/en/latest/_autosummary/jax.numpy.ndarray.at.html", + "output_type": "error", + "traceback": [ + "\u001b[31m---------------------------------------------------------------------------\u001b[39m", + "\u001b[31mTypeError\u001b[39m Traceback (most recent call last)", + "\u001b[36mCell\u001b[39m\u001b[36m \u001b[39m\u001b[32mIn[3]\u001b[39m\u001b[32m, line 8\u001b[39m\n\u001b[32m 3\u001b[39m \u001b[38;5;28;01mfor\u001b[39;00m step \u001b[38;5;129;01min\u001b[39;00m [\u001b[32m0\u001b[39m]: \u001b[38;5;66;03m#, 1200, 2400, 3600]:\u001b[39;00m\n\u001b[32m 4\u001b[39m particulator.run(step - particulator.n_steps)\n\u001b[32m 5\u001b[39m pyplot.step(\n\u001b[32m 6\u001b[39m x=radius_bins_edges[:-\u001b[32m1\u001b[39m] / si.um,\n\u001b[32m 7\u001b[39m y=particulator.formulae.particle_shape_and_density.volume_to_mass(\n\u001b[32m----> \u001b[39m\u001b[32m8\u001b[39m \u001b[43mparticulator\u001b[49m\u001b[43m.\u001b[49m\u001b[43mproducts\u001b[49m\u001b[43m[\u001b[49m\u001b[33;43m'\u001b[39;49m\u001b[33;43mdv/dlnr\u001b[39;49m\u001b[33;43m'\u001b[39;49m\u001b[43m]\u001b[49m\u001b[43m.\u001b[49m\u001b[43mget\u001b[49m\u001b[43m(\u001b[49m\u001b[43m)\u001b[49m[\u001b[32m0\u001b[39m]\n\u001b[32m 9\u001b[39m ) / si.g,\n\u001b[32m 10\u001b[39m where=\u001b[33m'\u001b[39m\u001b[33mpost\u001b[39m\u001b[33m'\u001b[39m, label=\u001b[33mf\u001b[39m\u001b[33m\"\u001b[39m\u001b[33mt = \u001b[39m\u001b[38;5;132;01m{\u001b[39;00mstep\u001b[38;5;132;01m}\u001b[39;00m\u001b[33ms\u001b[39m\u001b[33m\"\u001b[39m\n\u001b[32m 11\u001b[39m )\n\u001b[32m 13\u001b[39m pyplot.xscale(\u001b[33m'\u001b[39m\u001b[33mlog\u001b[39m\u001b[33m'\u001b[39m)\n\u001b[32m 14\u001b[39m pyplot.xlabel(\u001b[33m'\u001b[39m\u001b[33mparticle radius [µm]\u001b[39m\u001b[33m'\u001b[39m)\n", + "\u001b[36mFile \u001b[39m\u001b[32m~/Desktop/mgr/PySDM/PySDM/products/impl/product.py:99\u001b[39m, in \u001b[36mProduct.get\u001b[39m\u001b[34m(self, **kwargs)\u001b[39m\n\u001b[32m 98\u001b[39m \u001b[38;5;28;01mdef\u001b[39;00m\u001b[38;5;250m \u001b[39m\u001b[34mget\u001b[39m(\u001b[38;5;28mself\u001b[39m, **kwargs):\n\u001b[32m---> \u001b[39m\u001b[32m99\u001b[39m result = \u001b[38;5;28;43mself\u001b[39;49m\u001b[43m.\u001b[49m\u001b[43m_impl\u001b[49m\u001b[43m(\u001b[49m\u001b[43m*\u001b[49m\u001b[43m*\u001b[49m\u001b[43mkwargs\u001b[49m\u001b[43m)\u001b[49m\n\u001b[32m 100\u001b[39m result /= \u001b[38;5;28mself\u001b[39m.unit_magnitude_in_base_units\n\u001b[32m 101\u001b[39m \u001b[38;5;28;01mreturn\u001b[39;00m result\n", + "\u001b[36mFile \u001b[39m\u001b[32m~/Desktop/mgr/PySDM/PySDM/products/size_spectral/particle_volume_versus_radius_logarithm_spectrum.py:37\u001b[39m, in \u001b[36mParticleVolumeVersusRadiusLogarithmSpectrum._impl\u001b[39m\u001b[34m(self, **kwargs)\u001b[39m\n\u001b[32m 35\u001b[39m \u001b[38;5;28;01mdef\u001b[39;00m\u001b[38;5;250m \u001b[39m\u001b[34m_impl\u001b[39m(\u001b[38;5;28mself\u001b[39m, **kwargs):\n\u001b[32m 36\u001b[39m vals = np.empty([\u001b[38;5;28mself\u001b[39m.particulator.mesh.n_cell, \u001b[38;5;28mlen\u001b[39m(\u001b[38;5;28mself\u001b[39m.attr_bins_edges) - \u001b[32m1\u001b[39m])\n\u001b[32m---> \u001b[39m\u001b[32m37\u001b[39m \u001b[38;5;28;43mself\u001b[39;49m\u001b[43m.\u001b[49m\u001b[43m_recalculate_spectrum_moment\u001b[49m\u001b[43m(\u001b[49m\u001b[43mattr\u001b[49m\u001b[43m=\u001b[49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[43m.\u001b[49m\u001b[43mattr\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mrank\u001b[49m\u001b[43m=\u001b[49m\u001b[32;43m1\u001b[39;49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mfilter_attr\u001b[49m\u001b[43m=\u001b[49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[43m.\u001b[49m\u001b[43mattr\u001b[49m\u001b[43m)\u001b[49m\n\u001b[32m 39\u001b[39m \u001b[38;5;28;01mfor\u001b[39;00m i \u001b[38;5;129;01min\u001b[39;00m \u001b[38;5;28mrange\u001b[39m(vals.shape[\u001b[32m1\u001b[39m]):\n\u001b[32m 40\u001b[39m \u001b[38;5;28mself\u001b[39m._download_spectrum_moment_to_buffer(rank=\u001b[32m1\u001b[39m, bin_number=i)\n", + "\u001b[36mFile \u001b[39m\u001b[32m~/Desktop/mgr/PySDM/PySDM/products/impl/spectrum_moment_product.py:38\u001b[39m, in \u001b[36mSpectrumMomentProduct._recalculate_spectrum_moment\u001b[39m\u001b[34m(self, attr, rank, filter_attr, weighting_attribute, weighting_rank)\u001b[39m\n\u001b[32m 29\u001b[39m \u001b[38;5;28;01mdef\u001b[39;00m\u001b[38;5;250m \u001b[39m\u001b[34m_recalculate_spectrum_moment\u001b[39m(\n\u001b[32m 30\u001b[39m \u001b[38;5;28mself\u001b[39m,\n\u001b[32m 31\u001b[39m *,\n\u001b[32m (...)\u001b[39m\u001b[32m 36\u001b[39m weighting_rank=\u001b[32m0\u001b[39m,\n\u001b[32m 37\u001b[39m ):\n\u001b[32m---> \u001b[39m\u001b[32m38\u001b[39m \u001b[38;5;28;43mself\u001b[39;49m\u001b[43m.\u001b[49m\u001b[43mparticulator\u001b[49m\u001b[43m.\u001b[49m\u001b[43mspectrum_moments\u001b[49m\u001b[43m(\u001b[49m\n\u001b[32m 39\u001b[39m \u001b[43m \u001b[49m\u001b[43mmoment_0\u001b[49m\u001b[43m=\u001b[49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[43m.\u001b[49m\u001b[43mmoment_0\u001b[49m\u001b[43m,\u001b[49m\n\u001b[32m 40\u001b[39m \u001b[43m \u001b[49m\u001b[43mmoments\u001b[49m\u001b[43m=\u001b[49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[43m.\u001b[49m\u001b[43mmoments\u001b[49m\u001b[43m,\u001b[49m\n\u001b[32m 41\u001b[39m \u001b[43m \u001b[49m\u001b[43mattr\u001b[49m\u001b[43m=\u001b[49m\u001b[43mattr\u001b[49m\u001b[43m,\u001b[49m\n\u001b[32m 42\u001b[39m \u001b[43m \u001b[49m\u001b[43mrank\u001b[49m\u001b[43m=\u001b[49m\u001b[43mrank\u001b[49m\u001b[43m,\u001b[49m\n\u001b[32m 43\u001b[39m \u001b[43m \u001b[49m\u001b[43mattr_bins\u001b[49m\u001b[43m=\u001b[49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[43m.\u001b[49m\u001b[43mattr_bins_edges\u001b[49m\u001b[43m,\u001b[49m\n\u001b[32m 44\u001b[39m \u001b[43m \u001b[49m\u001b[43mattr_name\u001b[49m\u001b[43m=\u001b[49m\u001b[43mfilter_attr\u001b[49m\u001b[43m,\u001b[49m\n\u001b[32m 45\u001b[39m \u001b[43m \u001b[49m\u001b[43mweighting_attribute\u001b[49m\u001b[43m=\u001b[49m\u001b[43mweighting_attribute\u001b[49m\u001b[43m,\u001b[49m\n\u001b[32m 46\u001b[39m \u001b[43m \u001b[49m\u001b[43mweighting_rank\u001b[49m\u001b[43m=\u001b[49m\u001b[43mweighting_rank\u001b[49m\u001b[43m,\u001b[49m\n\u001b[32m 47\u001b[39m \u001b[43m \u001b[49m\u001b[43m)\u001b[49m\n", + "\u001b[36mFile \u001b[39m\u001b[32m~/Desktop/mgr/PySDM/PySDM/particulator.py:386\u001b[39m, in \u001b[36mParticulator.spectrum_moments\u001b[39m\u001b[34m(self, moment_0, moments, attr, rank, attr_bins, attr_name, weighting_attribute, weighting_rank)\u001b[39m\n\u001b[32m 384\u001b[39m \u001b[38;5;28mprint\u001b[39m(\u001b[33m\"\u001b[39m\u001b[33mAttr name: \u001b[39m\u001b[33m\"\u001b[39m + attr_name)\n\u001b[32m 385\u001b[39m attr_data = \u001b[38;5;28mself\u001b[39m.attributes[attr]\n\u001b[32m--> \u001b[39m\u001b[32m386\u001b[39m \u001b[38;5;28;43mself\u001b[39;49m\u001b[43m.\u001b[49m\u001b[43mbackend\u001b[49m\u001b[43m.\u001b[49m\u001b[43mspectrum_moments\u001b[49m\u001b[43m(\u001b[49m\n\u001b[32m 387\u001b[39m \u001b[43m \u001b[49m\u001b[43mmoment_0\u001b[49m\u001b[43m=\u001b[49m\u001b[43mmoment_0\u001b[49m\u001b[43m,\u001b[49m\n\u001b[32m 388\u001b[39m \u001b[43m \u001b[49m\u001b[43mmoments\u001b[49m\u001b[43m=\u001b[49m\u001b[43mmoments\u001b[49m\u001b[43m,\u001b[49m\n\u001b[32m 389\u001b[39m \u001b[43m \u001b[49m\u001b[43mmultiplicity\u001b[49m\u001b[43m=\u001b[49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[43m.\u001b[49m\u001b[43mattributes\u001b[49m\u001b[43m[\u001b[49m\u001b[33;43m\"\u001b[39;49m\u001b[33;43mmultiplicity\u001b[39;49m\u001b[33;43m\"\u001b[39;49m\u001b[43m]\u001b[49m\u001b[43m,\u001b[49m\n\u001b[32m 390\u001b[39m \u001b[43m \u001b[49m\u001b[43mattr_data\u001b[49m\u001b[43m=\u001b[49m\u001b[43mattr_data\u001b[49m\u001b[43m,\u001b[49m\n\u001b[32m 391\u001b[39m \u001b[43m \u001b[49m\u001b[43mcell_id\u001b[49m\u001b[43m=\u001b[49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[43m.\u001b[49m\u001b[43mattributes\u001b[49m\u001b[43m[\u001b[49m\u001b[33;43m\"\u001b[39;49m\u001b[33;43mcell id\u001b[39;49m\u001b[33;43m\"\u001b[39;49m\u001b[43m]\u001b[49m\u001b[43m,\u001b[49m\n\u001b[32m 392\u001b[39m \u001b[43m \u001b[49m\u001b[43midx\u001b[49m\u001b[43m=\u001b[49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[43m.\u001b[49m\u001b[43mattributes\u001b[49m\u001b[43m.\u001b[49m\u001b[43m_ParticleAttributes__idx\u001b[49m\u001b[43m,\u001b[49m\n\u001b[32m 393\u001b[39m \u001b[43m \u001b[49m\u001b[43mlength\u001b[49m\u001b[43m=\u001b[49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[43m.\u001b[49m\u001b[43mattributes\u001b[49m\u001b[43m.\u001b[49m\u001b[43msuper_droplet_count\u001b[49m\u001b[43m,\u001b[49m\n\u001b[32m 394\u001b[39m \u001b[43m \u001b[49m\u001b[43mrank\u001b[49m\u001b[43m=\u001b[49m\u001b[43mrank\u001b[49m\u001b[43m,\u001b[49m\n\u001b[32m 395\u001b[39m \u001b[43m \u001b[49m\u001b[43mx_bins\u001b[49m\u001b[43m=\u001b[49m\u001b[43mattr_bins\u001b[49m\u001b[43m,\u001b[49m\n\u001b[32m 396\u001b[39m \u001b[43m \u001b[49m\u001b[43mx_attr\u001b[49m\u001b[43m=\u001b[49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[43m.\u001b[49m\u001b[43mattributes\u001b[49m\u001b[43m[\u001b[49m\u001b[43mattr_name\u001b[49m\u001b[43m]\u001b[49m\u001b[43m,\u001b[49m\n\u001b[32m 397\u001b[39m \u001b[43m \u001b[49m\u001b[43mweighting_attribute\u001b[49m\u001b[43m=\u001b[49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[43m.\u001b[49m\u001b[43mattributes\u001b[49m\u001b[43m[\u001b[49m\u001b[43mweighting_attribute\u001b[49m\u001b[43m]\u001b[49m\u001b[43m,\u001b[49m\n\u001b[32m 398\u001b[39m \u001b[43m \u001b[49m\u001b[43mweighting_rank\u001b[49m\u001b[43m=\u001b[49m\u001b[43mweighting_rank\u001b[49m\u001b[43m,\u001b[49m\n\u001b[32m 399\u001b[39m \u001b[43m\u001b[49m\u001b[43m)\u001b[49m\n", + "\u001b[36mFile \u001b[39m\u001b[32m~/Desktop/mgr/PySDM/PySDM/backends/impl_jax/methods/moments_methods.py:185\u001b[39m, in \u001b[36mMomentsMethods.spectrum_moments\u001b[39m\u001b[34m(self, moment_0, moments, multiplicity, attr_data, cell_id, idx, length, rank, x_bins, x_attr, weighting_attribute, weighting_rank)\u001b[39m\n\u001b[32m 183\u001b[39m \u001b[38;5;28;01massert\u001b[39;00m moments.shape[\u001b[32m0\u001b[39m] == x_bins.shape[\u001b[32m0\u001b[39m] - \u001b[32m1\u001b[39m\n\u001b[32m 184\u001b[39m \u001b[38;5;28;01massert\u001b[39;00m moment_0.shape == moments.shape\n\u001b[32m--> \u001b[39m\u001b[32m185\u001b[39m \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[38;5;28;43mself\u001b[39;49m\u001b[43m.\u001b[49m\u001b[43m_spectrum_moments_body\u001b[49m\u001b[43m(\u001b[49m\n\u001b[32m 186\u001b[39m \u001b[43m \u001b[49m\u001b[43mmoment_0\u001b[49m\u001b[43m=\u001b[49m\u001b[43mmoment_0\u001b[49m\u001b[43m.\u001b[49m\u001b[43mdata\u001b[49m\u001b[43m,\u001b[49m\n\u001b[32m 187\u001b[39m \u001b[43m \u001b[49m\u001b[43mmoments\u001b[49m\u001b[43m=\u001b[49m\u001b[43mmoments\u001b[49m\u001b[43m.\u001b[49m\u001b[43mdata\u001b[49m\u001b[43m,\u001b[49m\n\u001b[32m 188\u001b[39m \u001b[43m \u001b[49m\u001b[43mmultiplicity\u001b[49m\u001b[43m=\u001b[49m\u001b[43mmultiplicity\u001b[49m\u001b[43m.\u001b[49m\u001b[43mdata\u001b[49m\u001b[43m,\u001b[49m\n\u001b[32m 189\u001b[39m \u001b[43m \u001b[49m\u001b[43mattr_data\u001b[49m\u001b[43m=\u001b[49m\u001b[43mattr_data\u001b[49m\u001b[43m.\u001b[49m\u001b[43mdata\u001b[49m\u001b[43m,\u001b[49m\n\u001b[32m 190\u001b[39m \u001b[43m \u001b[49m\u001b[43mcell_id\u001b[49m\u001b[43m=\u001b[49m\u001b[43mcell_id\u001b[49m\u001b[43m.\u001b[49m\u001b[43mdata\u001b[49m\u001b[43m,\u001b[49m\n\u001b[32m 191\u001b[39m \u001b[43m \u001b[49m\u001b[43midx\u001b[49m\u001b[43m=\u001b[49m\u001b[43midx\u001b[49m\u001b[43m.\u001b[49m\u001b[43mdata\u001b[49m\u001b[43m,\u001b[49m\n\u001b[32m 192\u001b[39m \u001b[43m \u001b[49m\u001b[43mlength\u001b[49m\u001b[43m=\u001b[49m\u001b[43mlength\u001b[49m\u001b[43m,\u001b[49m\n\u001b[32m 193\u001b[39m \u001b[43m \u001b[49m\u001b[43mrank\u001b[49m\u001b[43m=\u001b[49m\u001b[43mrank\u001b[49m\u001b[43m,\u001b[49m\n\u001b[32m 194\u001b[39m \u001b[43m \u001b[49m\u001b[43mx_bins\u001b[49m\u001b[43m=\u001b[49m\u001b[43mx_bins\u001b[49m\u001b[43m.\u001b[49m\u001b[43mdata\u001b[49m\u001b[43m,\u001b[49m\n\u001b[32m 195\u001b[39m \u001b[43m \u001b[49m\u001b[43mx_attr\u001b[49m\u001b[43m=\u001b[49m\u001b[43mx_attr\u001b[49m\u001b[43m.\u001b[49m\u001b[43mdata\u001b[49m\u001b[43m,\u001b[49m\n\u001b[32m 196\u001b[39m \u001b[43m \u001b[49m\u001b[43mweighting_attribute\u001b[49m\u001b[43m=\u001b[49m\u001b[43mweighting_attribute\u001b[49m\u001b[43m.\u001b[49m\u001b[43mdata\u001b[49m\u001b[43m,\u001b[49m\n\u001b[32m 197\u001b[39m \u001b[43m \u001b[49m\u001b[43mweighting_rank\u001b[49m\u001b[43m=\u001b[49m\u001b[43mweighting_rank\u001b[49m\u001b[43m,\u001b[49m\n\u001b[32m 198\u001b[39m \u001b[43m\u001b[49m\u001b[43m)\u001b[49m\n", + "\u001b[36mFile \u001b[39m\u001b[32m~/Desktop/mgr/PySDM/PySDM/backends/impl_jax/methods/moments_methods.py:138\u001b[39m, in \u001b[36mMomentsMethods._spectrum_moments_body..body\u001b[39m\u001b[34m(moment_0, moments, multiplicity, attr_data, cell_id, idx, length, rank, x_bins, x_attr, weighting_attribute, weighting_rank)\u001b[39m\n\u001b[32m 135\u001b[39m \u001b[38;5;28;01mif\u001b[39;00m x_bins[k] <= x_attr[i] < x_bins[k + \u001b[32m1\u001b[39m]:\n\u001b[32m 136\u001b[39m \u001b[38;5;66;03m# if x_bins[k] != x_bins[k + 1] and x_bins[k] == x_bins[k+1]:\u001b[39;00m\n\u001b[32m 137\u001b[39m t2 = time.time()\n\u001b[32m--> \u001b[39m\u001b[32m138\u001b[39m \u001b[43matomic_add\u001b[49m\u001b[43m(\u001b[49m\n\u001b[32m 139\u001b[39m \u001b[43m \u001b[49m\u001b[43mmoment_0\u001b[49m\u001b[43m,\u001b[49m\n\u001b[32m 140\u001b[39m \u001b[43m \u001b[49m\u001b[43m(\u001b[49m\u001b[43mk\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mcell_id\u001b[49m\u001b[43m[\u001b[49m\u001b[43mi\u001b[49m\u001b[43m]\u001b[49m\u001b[43m)\u001b[49m\u001b[43m,\u001b[49m\n\u001b[32m 141\u001b[39m \u001b[43m \u001b[49m\u001b[43mmultiplicity\u001b[49m\u001b[43m[\u001b[49m\u001b[43mi\u001b[49m\u001b[43m]\u001b[49m\u001b[43m \u001b[49m\u001b[43m*\u001b[49m\u001b[43m \u001b[49m\u001b[43mweighting_attribute\u001b[49m\u001b[43m[\u001b[49m\u001b[43mi\u001b[49m\u001b[43m]\u001b[49m\u001b[43m \u001b[49m\u001b[43m*\u001b[49m\u001b[43m*\u001b[49m\u001b[43m \u001b[49m\u001b[43mweighting_rank\u001b[49m\u001b[43m,\u001b[49m\n\u001b[32m 142\u001b[39m \u001b[43m \u001b[49m\u001b[43m)\u001b[49m\n\u001b[32m 143\u001b[39m atomic_add(\n\u001b[32m 144\u001b[39m moments,\n\u001b[32m 145\u001b[39m (k, cell_id[i]),\n\u001b[32m (...)\u001b[39m\u001b[32m 150\u001b[39m ),\n\u001b[32m 151\u001b[39m )\n\u001b[32m 152\u001b[39m \u001b[38;5;66;03m# print(f\"adding... : {time.time() - t2}\")\u001b[39;00m\n", + "\u001b[36mFile \u001b[39m\u001b[32m~/Desktop/mgr/PySDM/PySDM/backends/impl_numba/atomic_operations.py:102\u001b[39m, in \u001b[36matomic_add\u001b[39m\u001b[34m(ary, index, value)\u001b[39m\n\u001b[32m 93\u001b[39m \u001b[38;5;250m\u001b[39m\u001b[33;03m\"\"\"\u001b[39;00m\n\u001b[32m 94\u001b[39m \u001b[33;03mAtomically, perform `ary[i] += v` and return the previous value of `ary[i]`.\u001b[39;00m\n\u001b[32m 95\u001b[39m \n\u001b[32m (...)\u001b[39m\u001b[32m 99\u001b[39m \u001b[33;03mThis should be used from numba compiled code.\u001b[39;00m\n\u001b[32m 100\u001b[39m \u001b[33;03m\"\"\"\u001b[39;00m\n\u001b[32m 101\u001b[39m orig = ary[index]\n\u001b[32m--> \u001b[39m\u001b[32m102\u001b[39m \u001b[43mary\u001b[49m\u001b[43m[\u001b[49m\u001b[43mindex\u001b[49m\u001b[43m]\u001b[49m += value\n\u001b[32m 103\u001b[39m \u001b[38;5;28;01mreturn\u001b[39;00m orig\n", + "\u001b[36mFile \u001b[39m\u001b[32m~/Desktop/mgr/PySDM/.venv/lib/python3.11/site-packages/jax/_src/numpy/array_methods.py:617\u001b[39m, in \u001b[36m_unimplemented_setitem\u001b[39m\u001b[34m(self, i, x)\u001b[39m\n\u001b[32m 613\u001b[39m \u001b[38;5;28;01mdef\u001b[39;00m\u001b[38;5;250m \u001b[39m\u001b[34m_unimplemented_setitem\u001b[39m(\u001b[38;5;28mself\u001b[39m, i, x):\n\u001b[32m 614\u001b[39m msg = (\u001b[33m\"\u001b[39m\u001b[33mJAX arrays are immutable and do not support in-place item assignment.\u001b[39m\u001b[33m\"\u001b[39m\n\u001b[32m 615\u001b[39m \u001b[33m\"\u001b[39m\u001b[33m Instead of x[idx] = y, use x = x.at[idx].set(y) or another .at[] method:\u001b[39m\u001b[33m\"\u001b[39m\n\u001b[32m 616\u001b[39m \u001b[33m\"\u001b[39m\u001b[33m https://docs.jax.dev/en/latest/_autosummary/jax.numpy.ndarray.at.html\u001b[39m\u001b[33m\"\u001b[39m)\n\u001b[32m--> \u001b[39m\u001b[32m617\u001b[39m \u001b[38;5;28;01mraise\u001b[39;00m \u001b[38;5;167;01mTypeError\u001b[39;00m(msg.format(\u001b[38;5;28mtype\u001b[39m(\u001b[38;5;28mself\u001b[39m)))\n", + "\u001b[31mTypeError\u001b[39m: JAX arrays are immutable and do not support in-place item assignment. Instead of x[idx] = y, use x = x.at[idx].set(y) or another .at[] method: https://docs.jax.dev/en/latest/_autosummary/jax.numpy.ndarray.at.html" + ] + } + ], + "source": [ + "from matplotlib import pyplot\n", + "\n", + "for step in [0]: #, 1200, 2400, 3600]:\n", + " particulator.run(step - particulator.n_steps)\n", + " pyplot.step(\n", + " x=radius_bins_edges[:-1] / si.um,\n", + " y=particulator.formulae.particle_shape_and_density.volume_to_mass(\n", + " particulator.products['dv/dlnr'].get()[0]\n", + " ) / si.g,\n", + " where='post', label=f\"t = {step}s\"\n", + " )\n", + "\n", + "pyplot.xscale('log')\n", + "pyplot.xlabel('particle radius [µm]')\n", + "pyplot.ylabel(\"dm/dlnr [g/m$^3$/(unit dr/r)]\")\n", + "pyplot.legend()\n", + "pyplot.savefig('readme.png')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d5e22881", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "31" + ] + }, + "execution_count": 4, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "32768\n", + "\n", + "31" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": ".venv", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.11.0" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/readme.png b/readme.png new file mode 100644 index 0000000000000000000000000000000000000000..1ec42f09932c62c40e69927327a733526661c1d3 GIT binary patch literal 14636 zcmeHu2T)Yq*5xIOiU}1#K$I*93N$&1l7l2A2NhbN2}*8atAGKOoTFqUCkc{_T2>-2k8zX2tuW(a9tBY zhzSTnG<<*@KEZShjKYhkqpYr@mL1m7<+i;!qI%o$zO|jBwdEZqXLEZ8OFLVBo-3Dm zgt(b393AgFi1PB<{CxzEo&7ytl>%mMxC-Tc1w97@IeHuYN0cFzZiyhmk&4%^YP-hH z{-E+`Tix4V?6rOyWjOxkisShQpRPM=D`6B<4Z|k$4Oi2(`-AL{P>am^q{#T6P~*(Z z7|P`M6&TLvqm_H%+V}JJVJ~hx_-z31Tm#0W`Zxw|Nj^NJ5v&cn9wAO zH#Ikh4tZ^>zR}LS;N|62&2t1nIDE3Rxkd0xDL#IFO1=l3dwh>0%F-WAbaf;1{D~0c z`R$=<&(6fFBn8Kb@=1%w+m^}R=$E|7^Cv-&gXdXUKWh)WdIfy1#}M6L{r!~>zhtJU zpb+?&;%MWq&;F?frM6_o4s!z{QI~H&s$6ar3l0rsi7;}G`t<2jeeFY%Lrynt6y^Cz zp=&(mpYZudsUbC@rKOd6GwT*xE_Zr))FtV~J z9IGA0sY3U4+f2QpG?AXN2X}eHk0D5+b-u%-|5J~xsr1gbHU?1w zmciY5CMLPnTUN9p_9?L<4qtj@;ZYVY_I4wu!+Lf+ms^BR(9>IQ9b%Ml(`NKs(b(Kw z9~zr29cnFYOSBac3nE6;yEzqF6Qp=)X=$r@p5r#Rx0R)(312;T_MVmysc&yh(Tx0X z849!d`X=hGR=Qf4ag|#`a}>sKv_5Zl1;aO_r0f~04Euk3@W zz`$etW&D>fxA)3&hGbI_wo=>EOm1vwwM1k1y`yG_YD}U9EbTim$?^f=larJFenBLN zdgSvb67Fk}{hVpZXEkk(WY^NwWT<-I-eU z4XhE7cg37%E}I1Ju=yT7c#!Wbg2<-6zJ4w77QMyn$x3;H#f!JHIo#LhFbVIR*JW&T zpM9bm6}0I$6HLuCImgb%hV30(ogb3g-l_K7ZcL3ZzDPoTI55F$tD%{BW2NgxjHTI< z)gBp=D4>3J?6Bm-qw?`Cy&A8oV%s16)=Du}s|9K^*~M+X?1wMrtZx6#${YRkj2bJj zaP6oBi_tshFdY-mulMF*k!vY)Es5v86(n1(x3vwFJMfA+w_tiBkD(db=LxuWSeBVLJJG`G<#%N z|ESbXzMkmKQXOMFt9zokC2zu;eYP7hf~sF$&@!+8a7hpm^xpHxG$_sG)+@Zu=&^8> zKS^hM!+WJe-hOA}!SX=+>uW92H%?L!9qSPdiPuP4Nbi}Kur+393EU4uPV%ly)|{_a z&rI@uz32M1J-sWX#G*rnO(oLw)z@4ju{Hk?#JhMvbeFk4gsmQJrmsm1*`%8Yw-Uc-2-V#fNl@W^1R*lD5_KIFQUU4O)b zlamv6Q2y6}YP%CpkzFSwht@2?W*1ci~iN@egPglRaxIH9_Aj4XszjPjd%d6S5 zy4Bg4AJ(ghgAw&y&UTFqH^_*vNDBEB4}PyDoic8@V+VA&1xJ zulYq%vC5kP1IRy&&Ot37sWgHhO~()j4(@9R>A!tJ0jYGe;=rw82c9$NffU-=x*a>> zo?FFv{&&dQCT2DVg3d0(ZHAG!GF_c)VULaczN3_sSC$7~ui*9scBUq>6H2}2%Kn`F zWU|n_k(&6Y1l^u^f)I}!+^6=u{zy8m@3&j|)YMcp&nui|5iiD@j2XS;Wh!>klN%=< zqeJ`WLH)O*{ja`wfg>LN{WU8U*PDBSj?St_->}e3@pT@+K@KS?se-F(nNAagz53}g zHY6=AO)2YEjv2(AY{eX4M`sO8=V#D-mYf%?aW}6!Izh3C#I8N*uq~h7=mnx!-N`st?!;WKyKQ;$?aG%msV>= z*jdKh=i$08x_Ts$K|!~pJ?jF_@d8LWz0UK0`4UBXbGzHOZyR4?hogjaX^|W8;+)gH zc@aN9K8fj<`T;S_B(=AZjNe)9HJQNwfW*csA|m4cz0@vY#GC?##mN;zfG4_W@YcG) z$$GpwifzMl^*!k-wL;7BWghuF<=~_4S!P<`_YnVA$G9AN{CQ zXMa2-p)lLq+0=Y3LuT2LB%3=u-JQt{8;^zoZgKx?FR8RmXVq$tLD710a&k*sMg~We zm`g6bxU2dJYHCki8FcF2V~Fza*1ZZ23k#cHrXxPI1aTZi1v_sxo1PTl)|4Mj*8h(Z zGRg%fBEr&~qrO~0ZkBe3$QWZec&T(QUbr~CW9^OHxmsa&juV*PAB1m0%8l@|+ zn$8lO4qySMMm>UCPjh(locB(FBOsksZ|)@^Dq(B0-#o^TOW%cL0&LyswJPbf)_E9mMq)Tb8e34tOm!urEpNFH4b1pA;qa zBdr+sK93M&H)dD&&9KwA^h?`Z15U{2y|*)i9fHk?9}ClA zf#en(9NfI96m~Ic&lAqK9a}{pj+<2_3k!@0AYzy0ZW_NPZVp;3jAn=5- z;rEY&BO@bsce3@0Mt3%sm%iN%h-&%hAPKWzfR{j=r+aeh9kONz%C&OKem%cn)LU>@ zzM-+Paeed9wI6S8W<6ioLT~m{;heOM@mKuX!6IKDA7Kp^_z(x@TZLF@L_6#+Mbt!= z-~o%()nW)vFMlMt?ymM$-8v0N=}J-^hMT>lPatQLkod7eiGp!AL~d<~%byc)C`wpH zJKL3yo$jj1doc2)Gh2_J8aa!xV!7sbIOv?$#`|X|HtDua6Wxat%eG)?B5h=4(-6#x zeKgJ>=B(i1QTgJmP-}CY3nbe5IF~`oGy+^(TEv#PC%OtUc~ipPO1gN0U2ksOVR~Wp zU%!5-C>R!DUkX@ue)|3y)-XT3mn1Q_XnJZp9pgc_vUJe(>eZ`!j+5$kG=LekZm(Vv z=js+DSr1nBWqyZCCvvNhB!`IC0E(y$o`Ats$CArYxC6bF59c-EvFR@kY7FBt5ikX; zMR*PppaKY+&VAP^-lI_~vbnh_qIVL$z=_ZyI=O~ZqqFd4I7^an<;s=dz`)^8Iq+Wn zH1Y7>?sh5wI9hnn3R*{cEPkd54i6vylmkN;_D3>Eq#Som;;8dsSMc5WbRhF)R`|+H zpVoL=0xx8Hhom!bo%{B5kq1v*RCscP!K^byvCwndSr4M@`t|GeXYDG_zkdDN9>{)E zD5usaY&*46P>I$my=oPxw@Y9M z=g*(#a~M}a>#dea(;u%G#IUc7$}TEsWEq97w{ zxT#~r3Jm^|1}zM{kwl&(@rhqm_$7nroRs%%%oxl~{Sqryerz*G_on;WOwhaeY6C!v zTYbfrAy7BSkRN5_ZVJD|6}AMYhu+MM!j%4D4oVIyV-Y6t5?f_@McAblp`rRF-WqOL z^R+vQ@X*}ZMBImfb+hvE@%3Jz@^TZGnrMyVXO42Ph3kyHEFrqg$*K6eCnsh)at@av z1`ihBJ4x^NH_~3bgvKcq;n}likZYYf9!R9Mg*l<6-RS7(j1A0FH4khh>ApFV49~eE zjN%+lzq=cp4CRTZ9iZ@h6_yGl6zbyFhINiLKtMDzHB-g&Iyyd7!01Bnz3M@%R1?%X z>|X<4P`rOXKSVYN%~Q`$aoIDV6<6~E7^jd!3i)MG(f*z81tB23{O1y%^larAl&XFB z_)*#MC{@rq=h;R~i@(~uW)W7^@>~9E%W8@H!-VDpsUdMGK#3?x&x#jx7$x!}^c*JM zyCohA?wv}NXTK!=!|G` zQc^Q%TmMRIQNn9M5!^TT-`aZG9o3 z!2(7RK~@XJg)$Ied?l2uhHY1v;G|k#qQGZ~a5ayBKQG)P7F3C?M7n{Ik!R`}A&5=6 zV1yvTILxM}r#Hn*xU(=ZN#|BBDgEJZ--(9-A$Uc@A)wK+)^UFD$3H@K!V`h$vI>8C zR?Im&&mX-MSptQ0SZu5{kZYtYn{wl;j!8>NLRc6p#Zksk@7-yW&QxU@_CJSSg1bVb z#NujXd=F~AbGCW*{CRLkvb@;e!0+^l7cOHsNOAefl>&u{bHev)apv$K1QMc!t?9go z{bldHfSrYpM>pe>Nu;rD@s}>?7Y7DNZq!kM112qXA9yYb? zN_z!`7*Mzgs5aEd7PhZ!@)Z{XqHTq4k$cXZv z6xDLQjfzen5U5K7jfuy3OG`_&0Xq1$`T)UUyg4W(g#`};DCJ)#sY|9ttI5}%9Koi{ zLx7t9sBk)Q{P=mh%0)lgQy%|FP>~d9?RSZ2K?Vww5RGKHhW-+qIbZ}Y zyM}fl8OuBPXqNzbkqanCjar_UsbE}4y!v=F^Zk~w_pZ~_VAZ>9-2xd@ z+*mfWG$_0O`ntS)Q`v(NGpIp2-`wPy8YnNKm-G;BTc+XGWonb$&KB99Nbn9X9dv7| z{4%t+R;tMhdK2N)GE%+uJ1gz&JdlUZ{M*<1itYtuPIZ|N*#aQ~N_jtX@CJJfBcIt* zG%?N8eg52VFToQ_o@mkjuDx`~i%V$GHF;(Bd+!XcrJ+uOu)xy2?-{dgUacl9e8qPW@ZN}(GE|e+bq071NR@7f9X=G1nqEhXSq%9 zC2WVxT4rlPR_48$3Up}!H*#b{P$7LTerevrM!*(jxgw(QZJM{Of-)N1qW+BpNCzW7~6~RfUfcA(W64#uDDu)z+DAhDFFCKvu#$2D0yXm8&3wFr2ME=$b=*`QBtxrdpgRmP=JM zv$#hX zVhBFcf39{38JCao6|x(&mYcw0vEy!YgYNT}dT`|AtS1w~a?uznwdJPEx!L zSIoBbAu#>Cre(-|<6^lktt{Y5>c&9gK!3dm317%_o+pLmTJbrP2bL|dYHT@Ot;M+B z*5<`@j9a|Md$%DNS`60kD-=@v_GJKK*w$Jo?4$IrO3YC{f>pB}wLm;}7VF0Lw&x>& zxBHV@hicCcR<)RQZS_~o-3RcajE&!BVPhMAvv`&t`+_oVODiIwOQVr)Z*MWcJJhNi zuF^oZmH(Zn^U?sBTts`q%%0hN~{&BqsSP>U)nB#uJIe{@t)^MLU*L6hqAV3e7vAl=LTmIQgj=4!y5|3>2FC>*252{HVb6 zI*3sCz^ApWME60sR2zh3hE)G1$!IEKWnmG1zy)tZeUUQneXl7f7X>WZr9lZ~Lp3$p z<7i@wAg>j17+28L430Wb@vmWt?o8Zv9tuk*7ybh*0nRo8C~#|~sRWTw91*_B0$1kR zw3Zu(lCvR0BiRf%F%(}G;|IDik4jJ(9^s)Si%k%9`dtrT)QtLBG-S0PVT%V$(-T0T zKykQyeb5~(TGdVi)vUevioEvaS)oQ4fAe!D)fhoco%77M^Lu-H_8mAcBzrn)0J`}nJob0j?N`6*i0N^|V`d2hgJAdU+4DRq zxseFyu;+c~*fG;_!~=w#mTDLhAbi~4Ly9o8f7To?@!hBD$2*spyE=vsd`0GtGcwMZ z1;NF#M^@Q}Q7Qltn|qgTL;tO{<=fbwH+CW2wd67YZb^$ACIEd4hoZcyeWKKEXpox< z(iOqJpEMS-HRw7Rt;&<+%xIC=h)24g#lbPgDRi~47w{!2a1X&{H^dMysnib5$c3Uu znczc@fd6d8J>W*tfLDMq!XxG~x6&_xZH_z-D;5sf_}SB^m!}WXprSrncQ=?v7_mWF ziKY70Cbd6+8-wY~x8f`VE*gFFYNh$l556ZD7%oA5z_tS`tt#p=XPuS^c&iq_Oyvx$ zS9~3GV5QxyRdav@R5Qs7I}R&e!8Tm8Gse5rV_=6LD}jZ+h?dPnB;*Rq<1MPkx*O6b znxjG{Hr;l<72GW}YdBH9_~oo>%g;-t0~bJzSyI|WCpC*myW&e8c^=BVKu$}sA8VxX zE&J8pU+LpRN-N>!xK@Bxk!EHoKny1EE$InfR%~%SY{NGsH?N}-kmvSV*u8kX+|KT< z|4o!$F4cA*0CKm==PIU~N1#q0b!%;F%lnq+zA~xNP+uQd;WRxin!Ek|txeN%G^}(} zuZj1iqY^8DQ|qO;ULHL71*#^WK@b?jxBQ((X82&*bdRHTGpD&M7MnJ;FNEvi)*4C2 z1XCQL4+PC;;w*Ca?pMg%*6~4z zHbuY!OU+h)!7J?S?Dy=p*JioKTmX}X3+bdw?tDI21l0Fi_lL32A;L?k6KBp?SS_qL zy>ptzzM8^K?=AfNxLhem;l-_+UmKggOD56MNq2a4bty^f^CU#Mdc24}Nk&V$%fN#t zl-Pc|c_7T?Z-d)8O;J_{ixiUJy{CslRAy!-Ry+UqP>px))eBpduZU5BOLbv;rl@Vy zm&W+=V;NYgf#wc6>9}a~xqwT=>B?bZTNRpvjM=y#QY3~+E5`-)t8*1+W>m5}S~YDu|J4cib1c}9*Z=?^#O z-pyCuVuafjR8hJv48KnsVt_Z?nMC65YuQoGvh0s~usTMJxdLpD;1Q*1sR$LP^VeFI zRSp{IQ(uKON?gc}I{-AW_Ev5r`ihumkgHEjaB*)TFtKG!i0wdGf%-(pA5OHLG`wpP9;CPre9C@r2XEP5oY^;2aQDdu_>8bYD~(VEMfeAWI&K4Xa?xp| z!&pV$1?n4tZhp`^{yt^_yt&|S-@eT(q^~SOIJSW8kKFH5VggOuezZQw|J8H>;Gn49 z$ z{_@Gn&?C*#vuuWM61hj3dH?dKC=FjP&R>`~o z$f(xQm00N5I+M%E}e(6Rflv?-3>-VyNrN;!-q7`|qIne^J>O%Z7 zRs^gU5NZG7Tuv2eq1p350O2vJ(B4{~=amm&teO;w`@qcn_$IpZ48Kw=EY_Wi0IqHV zGJu8EzP+svaC3x)ycY7aqN8Im(9zz8@9151i1LBM6k2_~$ZiK54KZ0GiiHv#%&jTo z+^RMGOx=Pzxf2PlzbM)NRI#iem$WXf!-Ok=MvzI_T<)AD$@XWTawB`Y^WG7QBXvzh zZIZZK&SvIJdAgQGzznj&lQgZ|#P7P-SGj0`!NQa!VQRumLebY{d%FcVau zDG1pPMB=>vo49@M835Oa?yo!)5!tnC^ZhM=yU>savt89{I4E?P?du{gbFqEmnwpx= zjj%A6z2vA&*RNubs&p#qpLIIRw@!se2T_gRq{fTSea#x=&s%?-C``*|IB_g0*%R~) z9UbCDAF0#;cC}az6oGeE^Sm$hOj|we(*Fb-XACkBm>cgVcTBuCKOf$kgw48E(A=!D z%+JJBdzc<2&%70CF_8{}mN+U?U0vORMhHDqS){wx@h|AOk`u<2TU#sLf%VhR?? zIk9UOJTDW+Jxs6RT{)Hs0c`UbEkD#SA5y0|dD09N1j17D-@t--7{9ssPq31FJue62 zAll&wxv>Q|EpQnULtpZvZHSB2>0Y6hvWn9W{7+x9H!=Yhu)O?X1RmFR!~*Hg)N^dD z0&P*r&24M{G(5nI4__%`m1SVAqd$Qp!hhTmY6xzl3ITUhg*$gHD7V9Qo}A$V(9=^J z(#DM!pNsFTc81*&zN8RSoOLsct)Zbov&8DA2~GNdXq>ZgbR9Q;fOmAeo(wU=?z>M9 z-d33)95~I5-C1j@#3$sEMyTYAv?i+pJpWhz=0E(~W+|eL_Nn3*cZ_s0lWI0LH^bcU z2;vM5uw=LMucco`y%~Pobp7$r;=pk_nn3i#ayE{sw{TmnQWRrtT5e@cg}mX;r@F!JRktFW-{g$ox( z0mjh$I?E#GzW!)KUKPrj^7n39q*2*Q07zX#$k7YIoiK0p;9khHv$u-=LOA>5JLQkth2L|<^1`pasiBSp4pHT()4t2>jy_xrn-XT<9VPQY+jEA z!8xjhnUyuwht(wl?Vy5gj?rng^i#W#@J;cEIEYGSW@ZXd5;!pQqS1cpHzpoGPK)}8 zudhf!anI+qvvDG#H|J-bj$tLsjKtg@AKD3@Vm`85Oh=16ewGhIz@$ z2PD@A99~60)~V&pr6gEwEd0ba+5MV`_1az4d7=bDtNPorq(Ak_!T+LO`9IN9h}z5Q zZ_Uij8=_WTTrhDF&6(_l#~flpF6eAbpQ7QVbPw3B;D9Iaw;+}i%-Q6ECY>N51TLH8o6PogRjI1E*jJHY;xt*qYcY!B_XO0G>v z?dySUq572jMQG^BOepGE*x8lAZ&9ADid7TiZMhpBw&COz-5zDK#G@T786U6z@z)-R zIcm#tsn?%JczH*#=PFPDyZJla--1SpXD?oaZftCTwKAfAH!`H>43AzR%h2|p&Jd9r zY4cMpWpk*7gdxbT1>pqd1iayI05F%P(dL1?}Wyq4~uSK(;6CLQqU({+Bn zZzFzZZ7)0Qzog}t0}Z$)>RF-_X@1@W0`RpU*xO7pd8nD06p`Kx`7m$ivOME;lEo!^$64eLz|`S-IJuke zX&TQN(Bo-&&{7$^1Ng) zA2Dw|6)?=tNB@g%Adp{|Ovg}3YIOXm%i-<@c#yjpq)kbwvhCc~Ufv23{lB4AQ^0P} z5STwF;K`Cx5Wyo(_<>CSTN|#QKYw;6$sUT5aDSWUPwGn6RU3K@wQ+H0T4RK6NbW3L zGohj2Xn+5hVhs8u$7!rUnkcsHl0$Vr$oOkc9Rz>Ksr8_Uu&}UP*VNKdJzxk;rG{YX zRP)fOzsVrteo`!$3R+e*54k!A+GkAZ{RSitO~{QR&NI-_Og8$?E89;h^{V=`LI(T8 zv&RvJw{ELeQ~f1bAZtEPNJxOr&Cw)>|5?lmoUhp*WH%>M+7^om3^cXUMO$S!M69*x zhIDMmZK0Y$rPg?_Ez4#3L&Jik(oZZy9_@Wh+DFFzA#T}(bMyo1iIf6yu-{m?UDtFl zbll|@ZWpcizzDF<7|qwP^e};3$AdPCL8&Nznp@Yj;xdv8Wd7&btwmnrO5xrD=h=R9 ziM4Nc`?kRZ)&oae52!?Wt3kQF`7Q7e8MwIAgXiJwdqtv-lYtm$N0o0n^e$|4NH+j% zDsyns>Zi@gV&$kzX`!J48(k0eOi{qi4E6Ujb#-;y#n#dfYKpepq4n%s_;=fPsvz7? zeC+(-?7LZQWs`8EPzzXbnqoOvMv`3o!lqg@}VPm4M4K=%Yz|1l^l8^r_pG(&o!>>ef$ zW$Zl$3dS6_PVau>9Y{(_Ds02gmh=|VL9@y!xN92s2Ixp0`q*zT=eZlHtMjCP=AQ;1 zn_xak0xcN-tG{ZfO$S&qgT5bM(RIkAA9jwfAB}T~7CRQsUYGN0(C^35}r1drV zG3mF^snG~ZV9Tu>eGJ-(171{*k8r>rUl<%RXiRucc>z^PKR!8nPYRk9E>de&Z~Quh znzp0>Kw|Ti%+L}a0iNZRcZwP@{L5CvgEihYz_{hn^?=|PT-PtPHL5oE0=YaD8Zw|o zpS#!4;bn8INRj%yMvVIEuG8E=A3b;sZw%dMOUC1&9dl_iIiL~52YaI|b6eCJ0YOa#O6mKb)J(_jD%-cG>c+WLGR9-TC@c9(m09vPq!J0KL zjWx9_(_Q9Gc{m+zU>gLX*o2DD#N(itvhefsw-r`HR~@u&GIyolh=3}(35;D;5zHr9 zS=mz?j88yUg^k@JiP9K6SQ2APO}r7Pv499$b!W^-XP(!ViU4j#{cEKP+Q`(DWTb}J z<6S@a{61@2b`@5Y4cC)>+^%EE%0zN^UgZ9;H1!~+z6Nvwwr)x29UNcw_|xr1j+m1W z9h2tK-(H0~#LzKQ;h~N7%62@G@&5g_k8yr%IjWrt;YFJcXtM+?$c(2a4-->j?RXQ} zq0^Cc964c~bkg9|>YO@b0R}zmR7diYzTSnkt{YMf(5j_jigeSJf-pBG8eWu^3%WSp zgJ0KaMI9}ik?ugpi8z;N^$heU9h=bmp?Zuv7YGOSt@-i^M(R(XjG8^l7TR9HTZt!l zEGcSrz41H^=DM7N3Q2By&oRQg4${L+5mNh`ZEQREVsJy?hXYYR0(hYUgq|ZSs~~}5 zi(=k1h*5iU*C9KJa_{qAhEve2ND-ns%IRV+A=5-FlOCK zF>FjKjrgX$R>lQ7vb^y-m6K^u>H1*J^m%MCHY=mpPX9K3YoY~AB2xLbx8O3Mg7v}5 zwGV}s=R|URy9i(&oQnz+yZhyZ_>hi(Wvpv3d&Gm za3gM5_4D&nMd#dQ_I&G>aJCe~JTQWCL*ZeWkD-E5zZKe{_KVUl@ZMWaz+QfMq=Lix zWzehQOBMq8ok-ewwe1hP$8&NHLSehRS{x-jk*4}y_VvRXj{^kRaa0+C^Lwx=<*gkL z+3=o%WSK_q-bTrNzBW`t-W^-4C(g&$DWGTv-Gx3KPWJ8cxmY2;Sryg$?Wf>?zeWW`*=>j$s+hxg3NCIGS`xOn(>1XM;f2{ z`!_kgI?~g==Ovmi8}rF&xZ&LD{=3M=$b0l-mzlvQ(0Y1aBtQqN_p$?mqX8{S(XC&RqnqUPk3dnUEew6~bem+Jg4$B|Hz z9$yA+L|1ttrwc!$yOw3$c|R*TjCY1P!*jiuH?1eu?QzfGZ_hYguOkG0V2cEMlwrPb z_8O16SLu~UW9<5VVEV)HL8MV?9Q{@%ro*mcDs5)`u3Z=Y82knzwPg(%c+6%B6>?KM jDJm`g=imB0HQM7nt>WjC+xvtcjzbh>)vo8tm_GU+75c;I literal 0 HcmV?d00001 From 4aba12b5f9070eabe39fda0cb03fe343c8ed3244 Mon Sep 17 00:00:00 2001 From: Konrad Bodzioch Date: Wed, 7 Jan 2026 19:51:27 +0100 Subject: [PATCH 4/4] jax notebook first pass working --- .../impl_jax/methods/moments_methods.py | 183 ++++++++++++------ .../impl_jax/methods/physics_methods.py | 15 +- PySDM/backends/impl_jax/storage.py | 33 ++-- .../impl_numba/methods/moments_methods.py | 9 +- PySDM/particulator.py | 1 - ...volume_versus_radius_logarithm_spectrum.py | 2 - jax_test.ipynb | 128 +++--------- readme.png | Bin 14636 -> 13374 bytes 8 files changed, 178 insertions(+), 193 deletions(-) diff --git a/PySDM/backends/impl_jax/methods/moments_methods.py b/PySDM/backends/impl_jax/methods/moments_methods.py index 147191db2e..aabcc1b634 100644 --- a/PySDM/backends/impl_jax/methods/moments_methods.py +++ b/PySDM/backends/impl_jax/methods/moments_methods.py @@ -2,7 +2,7 @@ CPU implementation of moment calculation backend methods """ -from functools import cached_property +from functools import cached_property, partial import time import jax @@ -104,8 +104,10 @@ def moments( # @jax.jit def _spectrum_moments_body(self): # @numba.njit(**self.default_jit_flags) + # @partial(jax.jit, static_argnums=(6,)) + @jax.jit def body( - *, + # *, moment_0, moments, multiplicity, @@ -118,52 +120,54 @@ def body( x_attr, weighting_attribute, weighting_rank, + bin_to_count, + idx_i + # indices + # truth_table, ): - t1 = time.time() - # pylint: disable=too-many-locals - moment_0.at[:, :].set(0) - moments.at[:, :].set(0) - # print(f"initial moments setter: {time.time() - t1}") - # print(length) - for idx_i in range(length): # pylint: disable=not-an-iterable - i = idx[idx_i] - t4 = time.time() - # print(f"x_bins type: {type(x_bins)}") - # print(f"x_attr type: {type(x_attr)}") - for k in range(x_bins.shape[0] - 1): - # print(f"x_attr[i]: {x_attr[i]}") - if x_bins[k] <= x_attr[i] < x_bins[k + 1]: - # if x_bins[k] != x_bins[k + 1] and x_bins[k] == x_bins[k+1]: - t2 = time.time() - atomic_add( - moment_0, - (k, cell_id[i]), - multiplicity[i] * weighting_attribute[i] ** weighting_rank, - ) - atomic_add( - moments, - (k, cell_id[i]), - ( - multiplicity[i] - * weighting_attribute[i] ** weighting_rank - * attr_data[i] ** rank - ), - ) - # print(f"adding... : {time.time() - t2}") - break - # print(f"oop: {time.time() - t4}") - # print(f"After loop 1 checkpoint: {time.time() - t1}") - for c_id in range(moment_0.shape[1]): - for k in range(x_bins.shape[0] - 1): - moments.at[k, c_id].set( - moments[k, c_id] / moment_0[k, c_id] - if moment_0[k, c_id] != 0 - else 0 - ) - # print(f"spectrum moments body: {time.time() - t1}") + # moment_0 - 1 + # moments - 1 + # x_attr - 0 + # multiplicity - 0 + # weighting_attribute - 0 + # attr_data - 0 + + i = idx[idx_i] + +# (k[0] > 0 & ~(x_bins[k[0]-1] <= x_attr < x_bins[k[0]])) | + moment_0 = moment_0.at[bin_to_count, cell_id[i]].add(multiplicity[i] * weighting_attribute[i] ** weighting_rank) + moments = moments.at[bin_to_count, cell_id[i]].add(multiplicity[i] * weighting_attribute[i] ** weighting_rank * attr_data[i] ** rank) + + + # for k in range(x_bins.shape[0] - 1): + # if (x_bins[k] <= x_attr) & (x_attr < x_bins[k+1]): + # k > 0 and not (x_bins[k-1] <= x_attr < x_bins[k]) or k==0 + # moment_0 = moment_0.at[k].add(jax.numpy.multiply(multiplicity, weighting_attribute) ** weighting_rank) + # moments = moments.at[k].add(jax.numpy.multiply(jax.numpy.multiply(multiplicity, weighting_attribute) ** weighting_rank, attr_data ** rank)) + # break + + # Thing 2 (moments = this thing in another func or sth, if we even want to parallelize it) + # np.divide(a, b, out=np.zeros_like(a), where=b!=0) + + return moment_0, moments + # # This part below is only needed to scale it down, so focus on the first part working first + + # for k in range(x_bins.shape[0] - 1): + # moments = moments.at[k,:].divide(moment_0[k,:]) + # # moments = jax.numpy.divide(moments, moment_0) + return body + # @jax.jit + # def generate_calculate_indices(): + # for k in range(x_bins.shape[0] - 1): + # if (x_bins[k] <= x_attr) & (x_attr < x_bins[k+1]): + # moments + # pass + # @jax.jit + + def spectrum_moments( self, *, @@ -182,17 +186,82 @@ def spectrum_moments( ): assert moments.shape[0] == x_bins.shape[0] - 1 assert moment_0.shape == moments.shape - return self._spectrum_moments_body( - moment_0=moment_0.data, - moments=moments.data, - multiplicity=multiplicity.data, - attr_data=attr_data.data, - cell_id=cell_id.data, - idx=idx.data, - length=length, - rank=rank, - x_bins=x_bins.data, - x_attr=x_attr.data, - weighting_attribute=weighting_attribute.data, - weighting_rank=weighting_rank, + # truth_table = [] + + # for k in range(x_bins.shape[0] - 1): + # truth_table.append((x_bins[k] <= x_attr) & (x_attr < x_bins[k+1])) + + + # indices = jax.numpy.argwhere(truth_table) + + # print("done") + @jax.jit + def spectrum_moments_helper(x_bins, x_attr, idx, idx_i): + def cond_fun(k): + return ((x_bins[k] > x_attr[i]) | (x_attr[i] > x_bins[k+1])) & (k < x_bins.shape[0] - 1) + i = idx[idx_i] + bin_to_calculate = jax.lax.while_loop(cond_fun, lambda k: k+1, 0) + return bin_to_calculate + + moment_0.data = moment_0.data.at[:, :].set(0) + moments.data = moments.data.at[:,:].set(0) + idx_idxs = jax.numpy.arange(length-1) + + # maybe vmap just the indexes??? + count_bins_func = jax.vmap(spectrum_moments_helper, (None, None, None, 0)) + bins_to_count = count_bins_func(x_bins.data, x_attr.data, idx.data, idx_idxs) + mapped_spectrum = jax.vmap(self._spectrum_moments_body, (None, None, None, None, None, None, None, None, None, None, None, None, 0, 0)) + + moment_0.data, moments.data = mapped_spectrum( + moment_0.data, + moments.data, + multiplicity.data, + attr_data.data, + cell_id.data, + idx.data, + length, + rank, + x_bins.data, + x_attr.data, + weighting_attribute.data, + weighting_rank, + bins_to_count, + idx_idxs ) + + # This sum takes too much time ??? + moments.data = moments.data.sum(0) + moment_0.data = moment_0.data.sum(0) + # moment_0.data, moments.data = + # self._spectrum_moments_body( + # moment_0.data, + # moments.data, + # multiplicity.data, + # attr_data.data, + # cell_id.data, + # idx.data, + # length, + # rank, + # x_bins.data, + # x_attr.data, + # weighting_attribute.data, + # weighting_rank, + # # indices, + # # truth_table=truth_table, + # ) + + # print(moments.data) + # return self._spectrum_moments_body( + # moment_0=moment_0.data, + # moments=moments.data, + # multiplicity=multiplicity.data, + # attr_data=attr_data.data, + # cell_id=cell_id.data, + # idx=idx.data, + # length=length, + # rank=rank, + # x_bins=x_bins.data, + # x_attr=x_attr.data, + # weighting_attribute=weighting_attribute.data, + # weighting_rank=weighting_rank, + # ) diff --git a/PySDM/backends/impl_jax/methods/physics_methods.py b/PySDM/backends/impl_jax/methods/physics_methods.py index 75b98c9af8..f20f97ab4b 100644 --- a/PySDM/backends/impl_jax/methods/physics_methods.py +++ b/PySDM/backends/impl_jax/methods/physics_methods.py @@ -21,21 +21,16 @@ def _volume_of_mass_body(self): ff = self.formulae_flattened # @numba.njit(**self.default_jit_flags) - # @jax.jit + @jax.jit def body(volume, mass): - t1 = time.time() - for i in range(volume.shape[0]): # pylint: disable=not-an-iterable - # volume.at[i].set(ff.particle_shape_and_density__mass_to_volume.py_func(mass[i])) - volume = volume.at[i].set(ff.particle_shape_and_density__mass_to_volume(mass[i])) - # print(f"Mass[i]: {mass[i]}, Volume: {volume[i]}") - print("Post loop: ", volume) + volume = ff.particle_shape_and_density__mass_to_volume(mass) + # for i in range(volume.shape[0]): # pylint: disable=not-an-iterable + # # volume.at[i].set(ff.particle_shape_and_density__mass_to_volume.py_func(mass[i])) + # volume = volume.at[i].set(ff.particle_shape_and_density__mass_to_volume(mass[i])) return volume - # print(f"volume of water mass run time: {time.time() - t1}") return body def volume_of_water_mass(self, volume, mass): - print("Pre volume: ", volume.data) volume.data = self._volume_of_mass_body(volume.data, mass.data) - print("Post volume: ", volume.data) diff --git a/PySDM/backends/impl_jax/storage.py b/PySDM/backends/impl_jax/storage.py index 3f8b6ca6f8..642f17c2fb 100644 --- a/PySDM/backends/impl_jax/storage.py +++ b/PySDM/backends/impl_jax/storage.py @@ -22,18 +22,25 @@ class Storage(StorageBase): BOOL = jnp.bool_ def row_view(self, i): - print('rowview') return Storage( StorageSignature(self.data[i], (*self.shape[1:],), self.dtype) ) + def ravel(self, other): + print("Begin ravel") + print(self.data) + if isinstance(other, Storage): + self.data = other.data.ravel() + else: + self.data = other.ravel() + print(self.data) + print("End ravel") + def at(self, index): - print('at') assert self.shape == (1,), "Cannot call at() on Storage of shape other than (1,)" return self.data[index] def __imul__(self, other): - print('imul') if hasattr(other, "data"): impl.multiply(self.data, other.data) else: @@ -41,7 +48,6 @@ def __imul__(self, other): return self def __itruediv__(self, other): - print('itruediv') if hasattr(other, "data"): self.data[:] /= other.data[:] else: @@ -49,20 +55,16 @@ def __itruediv__(self, other): return self def download(self, target, reshape=False): - # print('download') - t1 = time.time() if reshape: data = self.data.reshape(target.shape) else: data = self.data - target = jnp.asarray(data) - print(f'download: {time.time() - t1}') - # np.copyto(target, np.asarray(data), casting="safe") + # target = jnp.asarray(data) + # Not sure here? + np.copyto(target, np.asarray(data), casting="safe") @staticmethod def _get_empty_data(shape, dtype): - # print('get_empty_data') - t1 = time.time() if dtype in (float, Storage.FLOAT): data = jnp.full(shape, jnp.nan, dtype=Storage.FLOAT) dtype = Storage.FLOAT @@ -75,17 +77,14 @@ def _get_empty_data(shape, dtype): else: raise NotImplementedError() - print(f'get_empty_data: {time.time() - t1}') return StorageSignature(data, shape, dtype) @staticmethod def empty(shape, dtype): - print('empty') return empty(shape, dtype, Storage) @staticmethod def _get_data_from_ndarray(array): - print('get_data_from_ndarray') return get_data_from_ndarray( array=array, @@ -95,22 +94,17 @@ def _get_data_from_ndarray(array): @staticmethod def from_ndarray(array): - print('from_ndarray') result = Storage(Storage._get_data_from_ndarray(array)) return result def urand(self, generator): - print('urand') generator(self) def upload(self, data): - print('upload') self.fill(data) def fill(self, other): - t1 = time.time() - print('fill') if isinstance(other, Storage): # self.data[:] = other.data # self.data.at[:].set(other.data) @@ -119,4 +113,3 @@ def fill(self, other): # self.data[:] = other # self.data.at[:].set(other) self.data = other - print(f'fill: {time.time() - t1}') diff --git a/PySDM/backends/impl_numba/methods/moments_methods.py b/PySDM/backends/impl_numba/methods/moments_methods.py index b10368c009..20f217e320 100644 --- a/PySDM/backends/impl_numba/methods/moments_methods.py +++ b/PySDM/backends/impl_numba/methods/moments_methods.py @@ -118,16 +118,14 @@ def body( weighting_attribute, weighting_rank, ): + # pylint: disable=too-many-locals moment_0[:, :] = 0 moments[:, :] = 0 - print(length) for idx_i in numba.prange(length): # pylint: disable=not-an-iterable i = idx[idx_i] - t4 = time.time() for k in range(x_bins.shape[0] - 1): if x_bins[k] <= x_attr[i] < x_bins[k + 1]: - # print("adagio") atomic_add( moment_0, (k, cell_id[i]), @@ -143,7 +141,7 @@ def body( ), ) break - print(f"oop: {time.time() - t4}") + return for c_id in range(moment_0.shape[1]): for k in range(x_bins.shape[0] - 1): moments[k, c_id] = ( @@ -172,7 +170,7 @@ def spectrum_moments( ): assert moments.shape[0] == x_bins.shape[0] - 1 assert moment_0.shape == moments.shape - return self._spectrum_moments_body( + self._spectrum_moments_body( moment_0=moment_0.data, moments=moments.data, multiplicity=multiplicity.data, @@ -186,3 +184,4 @@ def spectrum_moments( weighting_attribute=weighting_attribute.data, weighting_rank=weighting_rank, ) + \ No newline at end of file diff --git a/PySDM/particulator.py b/PySDM/particulator.py index a7948912c3..474ecb00fe 100644 --- a/PySDM/particulator.py +++ b/PySDM/particulator.py @@ -381,7 +381,6 @@ def spectrum_moments( weighting_attribute="water mass", weighting_rank=0, ): - print("Attr name: " + attr_name) attr_data = self.attributes[attr] self.backend.spectrum_moments( moment_0=moment_0, diff --git a/PySDM/products/size_spectral/particle_volume_versus_radius_logarithm_spectrum.py b/PySDM/products/size_spectral/particle_volume_versus_radius_logarithm_spectrum.py index c140e43c6a..b9826c4eb1 100644 --- a/PySDM/products/size_spectral/particle_volume_versus_radius_logarithm_spectrum.py +++ b/PySDM/products/size_spectral/particle_volume_versus_radius_logarithm_spectrum.py @@ -35,12 +35,10 @@ def register(self, builder): def _impl(self, **kwargs): vals = np.empty([self.particulator.mesh.n_cell, len(self.attr_bins_edges) - 1]) self._recalculate_spectrum_moment(attr=self.attr, rank=1, filter_attr=self.attr) - for i in range(vals.shape[1]): self._download_spectrum_moment_to_buffer(rank=1, bin_number=i) vals[:, i] = self.buffer.ravel() self._download_spectrum_moment_to_buffer(rank=0, bin_number=i) vals[:, i] *= self.buffer.ravel() - vals *= 1 / np.diff(np.log(self.radius_bins_edges)) / self.particulator.mesh.dv return vals diff --git a/jax_test.ipynb b/jax_test.ipynb index c2a078f3cb..29bb0d529d 100644 --- a/jax_test.ipynb +++ b/jax_test.ipynb @@ -2,7 +2,7 @@ "cells": [ { "cell_type": "code", - "execution_count": 1, + "execution_count": 8, "id": "4c6293f6", "metadata": {}, "outputs": [], @@ -22,7 +22,7 @@ }, { "cell_type": "code", - "execution_count": 2, + "execution_count": null, "id": "fbb45435", "metadata": {}, "outputs": [ @@ -30,74 +30,7 @@ "name": "stdout", "output_type": "stream", "text": [ - "empty\n", - "get_empty_data: 0.051486968994140625\n", - "get_empty_data: 0.01377105712890625\n", - "empty\n", - "get_empty_data: 0.004712820053100586\n", - "get_empty_data: 7.104873657226562e-05\n", - "empty\n", - "get_empty_data: 0.014365196228027344\n", - "empty\n", - "get_empty_data: 7.319450378417969e-05\n", - "empty\n", - "get_empty_data: 0.007919073104858398\n", - "fill\n", - "fill: 5.0067901611328125e-06\n", - "empty\n", - "get_empty_data: 7.295608520507812e-05\n", - "fill\n", - "fill: 5.0067901611328125e-06\n", - "empty\n", - "get_empty_data: 0.011842012405395508\n", - "empty\n", - "get_empty_data: 6.699562072753906e-05\n", - "from_ndarray\n", - "get_data_from_ndarray\n", - "from_ndarray\n", - "get_data_from_ndarray\n", - "from_ndarray\n", - "get_data_from_ndarray\n", - "from_ndarray\n", - "get_data_from_ndarray\n", - "empty\n", - "get_empty_data: 0.012073993682861328\n", - "empty\n", - "get_empty_data: 5.817413330078125e-05\n", - "get_data_from_ndarray\n", - "empty\n", - "get_empty_data: 0.010919809341430664\n", - "empty\n", - "get_empty_data: 0.004091978073120117\n", - "empty\n", - "get_empty_data: 6.985664367675781e-05\n", - "rowview\n", - "upload\n", - "fill\n", - "fill: 5.9604644775390625e-06\n", - "empty\n", - "get_empty_data: 0.011930227279663086\n", - "upload\n", - "fill\n", - "fill: 5.0067901611328125e-06\n", - "empty\n", - "get_empty_data: 8.20159912109375e-05\n", - "upload\n", - "fill\n", - "fill: 3.814697265625e-06\n", - "from_ndarray\n", - "get_data_from_ndarray\n", - "get_data_from_ndarray\n", - "from_ndarray\n", - "get_data_from_ndarray\n", - "empty\n", - "get_empty_data: 4.8160552978515625e-05\n", "dict_keys(['multiplicity', 'signed water mass', 'water mass', 'cell id', 'volume'])\n", - "Pre volume: [nan nan nan ... nan nan nan]\n", - "Post loop: [3.00579754e-18 6.63746385e-18 1.02692410e-17 ... 1.11357449e-12\n", - " 1.16550703e-12 1.25977547e-12]\n", - "Post volume: [3.00579754e-18 6.63746385e-18 1.02692410e-17 ... 1.11357449e-12\n", - " 1.16550703e-12 1.25977547e-12]\n", "[3.00579754e-18 6.63746385e-18 1.02692410e-17 ... 1.11357449e-12\n", " 1.16550703e-12 1.25977547e-12]\n" ] @@ -115,7 +48,7 @@ "radius_bins_edges = np.logspace(np.log10(10 * si.um), np.log10(5e3 * si.um), num=32)\n", "\n", "env = Box(dt=1 * si.s, dv=1e6 * si.m ** 3)\n", - "builder = Builder(n_sd=n_sd, backend=JAX(), environment=env)\n", + "builder = Builder(n_sd=n_sd, backend=CPU(), environment=env)\n", "builder.add_dynamic(Coalescence(collision_kernel=Golovin(b=1.5e3 / si.s), adaptive=False))\n", "products = [ParticleVolumeVersusRadiusLogarithmSpectrum(radius_bins_edges=radius_bins_edges, name='dv/dlnr')]\n", "particulator = builder.build(attributes, products)\n", @@ -126,36 +59,19 @@ }, { "cell_type": "code", - "execution_count": 3, + "execution_count": 10, "id": "3c3564c8", "metadata": {}, "outputs": [ { - "name": "stdout", - "output_type": "stream", - "text": [ - "Attr name: volume\n", - "at\n" - ] - }, - { - "ename": "TypeError", - "evalue": "JAX arrays are immutable and do not support in-place item assignment. Instead of x[idx] = y, use x = x.at[idx].set(y) or another .at[] method: https://docs.jax.dev/en/latest/_autosummary/jax.numpy.ndarray.at.html", - "output_type": "error", - "traceback": [ - "\u001b[31m---------------------------------------------------------------------------\u001b[39m", - "\u001b[31mTypeError\u001b[39m Traceback (most recent call last)", - "\u001b[36mCell\u001b[39m\u001b[36m \u001b[39m\u001b[32mIn[3]\u001b[39m\u001b[32m, line 8\u001b[39m\n\u001b[32m 3\u001b[39m \u001b[38;5;28;01mfor\u001b[39;00m step \u001b[38;5;129;01min\u001b[39;00m [\u001b[32m0\u001b[39m]: \u001b[38;5;66;03m#, 1200, 2400, 3600]:\u001b[39;00m\n\u001b[32m 4\u001b[39m particulator.run(step - particulator.n_steps)\n\u001b[32m 5\u001b[39m pyplot.step(\n\u001b[32m 6\u001b[39m x=radius_bins_edges[:-\u001b[32m1\u001b[39m] / si.um,\n\u001b[32m 7\u001b[39m y=particulator.formulae.particle_shape_and_density.volume_to_mass(\n\u001b[32m----> \u001b[39m\u001b[32m8\u001b[39m \u001b[43mparticulator\u001b[49m\u001b[43m.\u001b[49m\u001b[43mproducts\u001b[49m\u001b[43m[\u001b[49m\u001b[33;43m'\u001b[39;49m\u001b[33;43mdv/dlnr\u001b[39;49m\u001b[33;43m'\u001b[39;49m\u001b[43m]\u001b[49m\u001b[43m.\u001b[49m\u001b[43mget\u001b[49m\u001b[43m(\u001b[49m\u001b[43m)\u001b[49m[\u001b[32m0\u001b[39m]\n\u001b[32m 9\u001b[39m ) / si.g,\n\u001b[32m 10\u001b[39m where=\u001b[33m'\u001b[39m\u001b[33mpost\u001b[39m\u001b[33m'\u001b[39m, label=\u001b[33mf\u001b[39m\u001b[33m\"\u001b[39m\u001b[33mt = \u001b[39m\u001b[38;5;132;01m{\u001b[39;00mstep\u001b[38;5;132;01m}\u001b[39;00m\u001b[33ms\u001b[39m\u001b[33m\"\u001b[39m\n\u001b[32m 11\u001b[39m )\n\u001b[32m 13\u001b[39m pyplot.xscale(\u001b[33m'\u001b[39m\u001b[33mlog\u001b[39m\u001b[33m'\u001b[39m)\n\u001b[32m 14\u001b[39m pyplot.xlabel(\u001b[33m'\u001b[39m\u001b[33mparticle radius [µm]\u001b[39m\u001b[33m'\u001b[39m)\n", - "\u001b[36mFile \u001b[39m\u001b[32m~/Desktop/mgr/PySDM/PySDM/products/impl/product.py:99\u001b[39m, in \u001b[36mProduct.get\u001b[39m\u001b[34m(self, **kwargs)\u001b[39m\n\u001b[32m 98\u001b[39m \u001b[38;5;28;01mdef\u001b[39;00m\u001b[38;5;250m \u001b[39m\u001b[34mget\u001b[39m(\u001b[38;5;28mself\u001b[39m, **kwargs):\n\u001b[32m---> \u001b[39m\u001b[32m99\u001b[39m result = \u001b[38;5;28;43mself\u001b[39;49m\u001b[43m.\u001b[49m\u001b[43m_impl\u001b[49m\u001b[43m(\u001b[49m\u001b[43m*\u001b[49m\u001b[43m*\u001b[49m\u001b[43mkwargs\u001b[49m\u001b[43m)\u001b[49m\n\u001b[32m 100\u001b[39m result /= \u001b[38;5;28mself\u001b[39m.unit_magnitude_in_base_units\n\u001b[32m 101\u001b[39m \u001b[38;5;28;01mreturn\u001b[39;00m result\n", - "\u001b[36mFile \u001b[39m\u001b[32m~/Desktop/mgr/PySDM/PySDM/products/size_spectral/particle_volume_versus_radius_logarithm_spectrum.py:37\u001b[39m, in \u001b[36mParticleVolumeVersusRadiusLogarithmSpectrum._impl\u001b[39m\u001b[34m(self, **kwargs)\u001b[39m\n\u001b[32m 35\u001b[39m \u001b[38;5;28;01mdef\u001b[39;00m\u001b[38;5;250m \u001b[39m\u001b[34m_impl\u001b[39m(\u001b[38;5;28mself\u001b[39m, **kwargs):\n\u001b[32m 36\u001b[39m vals = np.empty([\u001b[38;5;28mself\u001b[39m.particulator.mesh.n_cell, \u001b[38;5;28mlen\u001b[39m(\u001b[38;5;28mself\u001b[39m.attr_bins_edges) - \u001b[32m1\u001b[39m])\n\u001b[32m---> \u001b[39m\u001b[32m37\u001b[39m \u001b[38;5;28;43mself\u001b[39;49m\u001b[43m.\u001b[49m\u001b[43m_recalculate_spectrum_moment\u001b[49m\u001b[43m(\u001b[49m\u001b[43mattr\u001b[49m\u001b[43m=\u001b[49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[43m.\u001b[49m\u001b[43mattr\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mrank\u001b[49m\u001b[43m=\u001b[49m\u001b[32;43m1\u001b[39;49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mfilter_attr\u001b[49m\u001b[43m=\u001b[49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[43m.\u001b[49m\u001b[43mattr\u001b[49m\u001b[43m)\u001b[49m\n\u001b[32m 39\u001b[39m \u001b[38;5;28;01mfor\u001b[39;00m i \u001b[38;5;129;01min\u001b[39;00m \u001b[38;5;28mrange\u001b[39m(vals.shape[\u001b[32m1\u001b[39m]):\n\u001b[32m 40\u001b[39m \u001b[38;5;28mself\u001b[39m._download_spectrum_moment_to_buffer(rank=\u001b[32m1\u001b[39m, bin_number=i)\n", - "\u001b[36mFile \u001b[39m\u001b[32m~/Desktop/mgr/PySDM/PySDM/products/impl/spectrum_moment_product.py:38\u001b[39m, in \u001b[36mSpectrumMomentProduct._recalculate_spectrum_moment\u001b[39m\u001b[34m(self, attr, rank, filter_attr, weighting_attribute, weighting_rank)\u001b[39m\n\u001b[32m 29\u001b[39m \u001b[38;5;28;01mdef\u001b[39;00m\u001b[38;5;250m \u001b[39m\u001b[34m_recalculate_spectrum_moment\u001b[39m(\n\u001b[32m 30\u001b[39m \u001b[38;5;28mself\u001b[39m,\n\u001b[32m 31\u001b[39m *,\n\u001b[32m (...)\u001b[39m\u001b[32m 36\u001b[39m weighting_rank=\u001b[32m0\u001b[39m,\n\u001b[32m 37\u001b[39m ):\n\u001b[32m---> \u001b[39m\u001b[32m38\u001b[39m \u001b[38;5;28;43mself\u001b[39;49m\u001b[43m.\u001b[49m\u001b[43mparticulator\u001b[49m\u001b[43m.\u001b[49m\u001b[43mspectrum_moments\u001b[49m\u001b[43m(\u001b[49m\n\u001b[32m 39\u001b[39m \u001b[43m \u001b[49m\u001b[43mmoment_0\u001b[49m\u001b[43m=\u001b[49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[43m.\u001b[49m\u001b[43mmoment_0\u001b[49m\u001b[43m,\u001b[49m\n\u001b[32m 40\u001b[39m \u001b[43m \u001b[49m\u001b[43mmoments\u001b[49m\u001b[43m=\u001b[49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[43m.\u001b[49m\u001b[43mmoments\u001b[49m\u001b[43m,\u001b[49m\n\u001b[32m 41\u001b[39m \u001b[43m \u001b[49m\u001b[43mattr\u001b[49m\u001b[43m=\u001b[49m\u001b[43mattr\u001b[49m\u001b[43m,\u001b[49m\n\u001b[32m 42\u001b[39m \u001b[43m \u001b[49m\u001b[43mrank\u001b[49m\u001b[43m=\u001b[49m\u001b[43mrank\u001b[49m\u001b[43m,\u001b[49m\n\u001b[32m 43\u001b[39m \u001b[43m \u001b[49m\u001b[43mattr_bins\u001b[49m\u001b[43m=\u001b[49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[43m.\u001b[49m\u001b[43mattr_bins_edges\u001b[49m\u001b[43m,\u001b[49m\n\u001b[32m 44\u001b[39m \u001b[43m \u001b[49m\u001b[43mattr_name\u001b[49m\u001b[43m=\u001b[49m\u001b[43mfilter_attr\u001b[49m\u001b[43m,\u001b[49m\n\u001b[32m 45\u001b[39m \u001b[43m \u001b[49m\u001b[43mweighting_attribute\u001b[49m\u001b[43m=\u001b[49m\u001b[43mweighting_attribute\u001b[49m\u001b[43m,\u001b[49m\n\u001b[32m 46\u001b[39m \u001b[43m \u001b[49m\u001b[43mweighting_rank\u001b[49m\u001b[43m=\u001b[49m\u001b[43mweighting_rank\u001b[49m\u001b[43m,\u001b[49m\n\u001b[32m 47\u001b[39m \u001b[43m \u001b[49m\u001b[43m)\u001b[49m\n", - "\u001b[36mFile \u001b[39m\u001b[32m~/Desktop/mgr/PySDM/PySDM/particulator.py:386\u001b[39m, in \u001b[36mParticulator.spectrum_moments\u001b[39m\u001b[34m(self, moment_0, moments, attr, rank, attr_bins, attr_name, weighting_attribute, weighting_rank)\u001b[39m\n\u001b[32m 384\u001b[39m \u001b[38;5;28mprint\u001b[39m(\u001b[33m\"\u001b[39m\u001b[33mAttr name: \u001b[39m\u001b[33m\"\u001b[39m + attr_name)\n\u001b[32m 385\u001b[39m attr_data = \u001b[38;5;28mself\u001b[39m.attributes[attr]\n\u001b[32m--> \u001b[39m\u001b[32m386\u001b[39m \u001b[38;5;28;43mself\u001b[39;49m\u001b[43m.\u001b[49m\u001b[43mbackend\u001b[49m\u001b[43m.\u001b[49m\u001b[43mspectrum_moments\u001b[49m\u001b[43m(\u001b[49m\n\u001b[32m 387\u001b[39m \u001b[43m \u001b[49m\u001b[43mmoment_0\u001b[49m\u001b[43m=\u001b[49m\u001b[43mmoment_0\u001b[49m\u001b[43m,\u001b[49m\n\u001b[32m 388\u001b[39m \u001b[43m \u001b[49m\u001b[43mmoments\u001b[49m\u001b[43m=\u001b[49m\u001b[43mmoments\u001b[49m\u001b[43m,\u001b[49m\n\u001b[32m 389\u001b[39m \u001b[43m \u001b[49m\u001b[43mmultiplicity\u001b[49m\u001b[43m=\u001b[49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[43m.\u001b[49m\u001b[43mattributes\u001b[49m\u001b[43m[\u001b[49m\u001b[33;43m\"\u001b[39;49m\u001b[33;43mmultiplicity\u001b[39;49m\u001b[33;43m\"\u001b[39;49m\u001b[43m]\u001b[49m\u001b[43m,\u001b[49m\n\u001b[32m 390\u001b[39m \u001b[43m \u001b[49m\u001b[43mattr_data\u001b[49m\u001b[43m=\u001b[49m\u001b[43mattr_data\u001b[49m\u001b[43m,\u001b[49m\n\u001b[32m 391\u001b[39m \u001b[43m \u001b[49m\u001b[43mcell_id\u001b[49m\u001b[43m=\u001b[49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[43m.\u001b[49m\u001b[43mattributes\u001b[49m\u001b[43m[\u001b[49m\u001b[33;43m\"\u001b[39;49m\u001b[33;43mcell id\u001b[39;49m\u001b[33;43m\"\u001b[39;49m\u001b[43m]\u001b[49m\u001b[43m,\u001b[49m\n\u001b[32m 392\u001b[39m \u001b[43m \u001b[49m\u001b[43midx\u001b[49m\u001b[43m=\u001b[49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[43m.\u001b[49m\u001b[43mattributes\u001b[49m\u001b[43m.\u001b[49m\u001b[43m_ParticleAttributes__idx\u001b[49m\u001b[43m,\u001b[49m\n\u001b[32m 393\u001b[39m \u001b[43m \u001b[49m\u001b[43mlength\u001b[49m\u001b[43m=\u001b[49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[43m.\u001b[49m\u001b[43mattributes\u001b[49m\u001b[43m.\u001b[49m\u001b[43msuper_droplet_count\u001b[49m\u001b[43m,\u001b[49m\n\u001b[32m 394\u001b[39m \u001b[43m \u001b[49m\u001b[43mrank\u001b[49m\u001b[43m=\u001b[49m\u001b[43mrank\u001b[49m\u001b[43m,\u001b[49m\n\u001b[32m 395\u001b[39m \u001b[43m \u001b[49m\u001b[43mx_bins\u001b[49m\u001b[43m=\u001b[49m\u001b[43mattr_bins\u001b[49m\u001b[43m,\u001b[49m\n\u001b[32m 396\u001b[39m \u001b[43m \u001b[49m\u001b[43mx_attr\u001b[49m\u001b[43m=\u001b[49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[43m.\u001b[49m\u001b[43mattributes\u001b[49m\u001b[43m[\u001b[49m\u001b[43mattr_name\u001b[49m\u001b[43m]\u001b[49m\u001b[43m,\u001b[49m\n\u001b[32m 397\u001b[39m \u001b[43m \u001b[49m\u001b[43mweighting_attribute\u001b[49m\u001b[43m=\u001b[49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[43m.\u001b[49m\u001b[43mattributes\u001b[49m\u001b[43m[\u001b[49m\u001b[43mweighting_attribute\u001b[49m\u001b[43m]\u001b[49m\u001b[43m,\u001b[49m\n\u001b[32m 398\u001b[39m \u001b[43m \u001b[49m\u001b[43mweighting_rank\u001b[49m\u001b[43m=\u001b[49m\u001b[43mweighting_rank\u001b[49m\u001b[43m,\u001b[49m\n\u001b[32m 399\u001b[39m \u001b[43m\u001b[49m\u001b[43m)\u001b[49m\n", - "\u001b[36mFile \u001b[39m\u001b[32m~/Desktop/mgr/PySDM/PySDM/backends/impl_jax/methods/moments_methods.py:185\u001b[39m, in \u001b[36mMomentsMethods.spectrum_moments\u001b[39m\u001b[34m(self, moment_0, moments, multiplicity, attr_data, cell_id, idx, length, rank, x_bins, x_attr, weighting_attribute, weighting_rank)\u001b[39m\n\u001b[32m 183\u001b[39m \u001b[38;5;28;01massert\u001b[39;00m moments.shape[\u001b[32m0\u001b[39m] == x_bins.shape[\u001b[32m0\u001b[39m] - \u001b[32m1\u001b[39m\n\u001b[32m 184\u001b[39m \u001b[38;5;28;01massert\u001b[39;00m moment_0.shape == moments.shape\n\u001b[32m--> \u001b[39m\u001b[32m185\u001b[39m \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[38;5;28;43mself\u001b[39;49m\u001b[43m.\u001b[49m\u001b[43m_spectrum_moments_body\u001b[49m\u001b[43m(\u001b[49m\n\u001b[32m 186\u001b[39m \u001b[43m \u001b[49m\u001b[43mmoment_0\u001b[49m\u001b[43m=\u001b[49m\u001b[43mmoment_0\u001b[49m\u001b[43m.\u001b[49m\u001b[43mdata\u001b[49m\u001b[43m,\u001b[49m\n\u001b[32m 187\u001b[39m \u001b[43m \u001b[49m\u001b[43mmoments\u001b[49m\u001b[43m=\u001b[49m\u001b[43mmoments\u001b[49m\u001b[43m.\u001b[49m\u001b[43mdata\u001b[49m\u001b[43m,\u001b[49m\n\u001b[32m 188\u001b[39m \u001b[43m \u001b[49m\u001b[43mmultiplicity\u001b[49m\u001b[43m=\u001b[49m\u001b[43mmultiplicity\u001b[49m\u001b[43m.\u001b[49m\u001b[43mdata\u001b[49m\u001b[43m,\u001b[49m\n\u001b[32m 189\u001b[39m \u001b[43m \u001b[49m\u001b[43mattr_data\u001b[49m\u001b[43m=\u001b[49m\u001b[43mattr_data\u001b[49m\u001b[43m.\u001b[49m\u001b[43mdata\u001b[49m\u001b[43m,\u001b[49m\n\u001b[32m 190\u001b[39m \u001b[43m \u001b[49m\u001b[43mcell_id\u001b[49m\u001b[43m=\u001b[49m\u001b[43mcell_id\u001b[49m\u001b[43m.\u001b[49m\u001b[43mdata\u001b[49m\u001b[43m,\u001b[49m\n\u001b[32m 191\u001b[39m \u001b[43m \u001b[49m\u001b[43midx\u001b[49m\u001b[43m=\u001b[49m\u001b[43midx\u001b[49m\u001b[43m.\u001b[49m\u001b[43mdata\u001b[49m\u001b[43m,\u001b[49m\n\u001b[32m 192\u001b[39m \u001b[43m \u001b[49m\u001b[43mlength\u001b[49m\u001b[43m=\u001b[49m\u001b[43mlength\u001b[49m\u001b[43m,\u001b[49m\n\u001b[32m 193\u001b[39m \u001b[43m \u001b[49m\u001b[43mrank\u001b[49m\u001b[43m=\u001b[49m\u001b[43mrank\u001b[49m\u001b[43m,\u001b[49m\n\u001b[32m 194\u001b[39m \u001b[43m \u001b[49m\u001b[43mx_bins\u001b[49m\u001b[43m=\u001b[49m\u001b[43mx_bins\u001b[49m\u001b[43m.\u001b[49m\u001b[43mdata\u001b[49m\u001b[43m,\u001b[49m\n\u001b[32m 195\u001b[39m \u001b[43m \u001b[49m\u001b[43mx_attr\u001b[49m\u001b[43m=\u001b[49m\u001b[43mx_attr\u001b[49m\u001b[43m.\u001b[49m\u001b[43mdata\u001b[49m\u001b[43m,\u001b[49m\n\u001b[32m 196\u001b[39m \u001b[43m \u001b[49m\u001b[43mweighting_attribute\u001b[49m\u001b[43m=\u001b[49m\u001b[43mweighting_attribute\u001b[49m\u001b[43m.\u001b[49m\u001b[43mdata\u001b[49m\u001b[43m,\u001b[49m\n\u001b[32m 197\u001b[39m \u001b[43m \u001b[49m\u001b[43mweighting_rank\u001b[49m\u001b[43m=\u001b[49m\u001b[43mweighting_rank\u001b[49m\u001b[43m,\u001b[49m\n\u001b[32m 198\u001b[39m \u001b[43m\u001b[49m\u001b[43m)\u001b[49m\n", - "\u001b[36mFile \u001b[39m\u001b[32m~/Desktop/mgr/PySDM/PySDM/backends/impl_jax/methods/moments_methods.py:138\u001b[39m, in \u001b[36mMomentsMethods._spectrum_moments_body..body\u001b[39m\u001b[34m(moment_0, moments, multiplicity, attr_data, cell_id, idx, length, rank, x_bins, x_attr, weighting_attribute, weighting_rank)\u001b[39m\n\u001b[32m 135\u001b[39m \u001b[38;5;28;01mif\u001b[39;00m x_bins[k] <= x_attr[i] < x_bins[k + \u001b[32m1\u001b[39m]:\n\u001b[32m 136\u001b[39m \u001b[38;5;66;03m# if x_bins[k] != x_bins[k + 1] and x_bins[k] == x_bins[k+1]:\u001b[39;00m\n\u001b[32m 137\u001b[39m t2 = time.time()\n\u001b[32m--> \u001b[39m\u001b[32m138\u001b[39m \u001b[43matomic_add\u001b[49m\u001b[43m(\u001b[49m\n\u001b[32m 139\u001b[39m \u001b[43m \u001b[49m\u001b[43mmoment_0\u001b[49m\u001b[43m,\u001b[49m\n\u001b[32m 140\u001b[39m \u001b[43m \u001b[49m\u001b[43m(\u001b[49m\u001b[43mk\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mcell_id\u001b[49m\u001b[43m[\u001b[49m\u001b[43mi\u001b[49m\u001b[43m]\u001b[49m\u001b[43m)\u001b[49m\u001b[43m,\u001b[49m\n\u001b[32m 141\u001b[39m \u001b[43m \u001b[49m\u001b[43mmultiplicity\u001b[49m\u001b[43m[\u001b[49m\u001b[43mi\u001b[49m\u001b[43m]\u001b[49m\u001b[43m \u001b[49m\u001b[43m*\u001b[49m\u001b[43m \u001b[49m\u001b[43mweighting_attribute\u001b[49m\u001b[43m[\u001b[49m\u001b[43mi\u001b[49m\u001b[43m]\u001b[49m\u001b[43m \u001b[49m\u001b[43m*\u001b[49m\u001b[43m*\u001b[49m\u001b[43m \u001b[49m\u001b[43mweighting_rank\u001b[49m\u001b[43m,\u001b[49m\n\u001b[32m 142\u001b[39m \u001b[43m \u001b[49m\u001b[43m)\u001b[49m\n\u001b[32m 143\u001b[39m atomic_add(\n\u001b[32m 144\u001b[39m moments,\n\u001b[32m 145\u001b[39m (k, cell_id[i]),\n\u001b[32m (...)\u001b[39m\u001b[32m 150\u001b[39m ),\n\u001b[32m 151\u001b[39m )\n\u001b[32m 152\u001b[39m \u001b[38;5;66;03m# print(f\"adding... : {time.time() - t2}\")\u001b[39;00m\n", - "\u001b[36mFile \u001b[39m\u001b[32m~/Desktop/mgr/PySDM/PySDM/backends/impl_numba/atomic_operations.py:102\u001b[39m, in \u001b[36matomic_add\u001b[39m\u001b[34m(ary, index, value)\u001b[39m\n\u001b[32m 93\u001b[39m \u001b[38;5;250m\u001b[39m\u001b[33;03m\"\"\"\u001b[39;00m\n\u001b[32m 94\u001b[39m \u001b[33;03mAtomically, perform `ary[i] += v` and return the previous value of `ary[i]`.\u001b[39;00m\n\u001b[32m 95\u001b[39m \n\u001b[32m (...)\u001b[39m\u001b[32m 99\u001b[39m \u001b[33;03mThis should be used from numba compiled code.\u001b[39;00m\n\u001b[32m 100\u001b[39m \u001b[33;03m\"\"\"\u001b[39;00m\n\u001b[32m 101\u001b[39m orig = ary[index]\n\u001b[32m--> \u001b[39m\u001b[32m102\u001b[39m \u001b[43mary\u001b[49m\u001b[43m[\u001b[49m\u001b[43mindex\u001b[49m\u001b[43m]\u001b[49m += value\n\u001b[32m 103\u001b[39m \u001b[38;5;28;01mreturn\u001b[39;00m orig\n", - "\u001b[36mFile \u001b[39m\u001b[32m~/Desktop/mgr/PySDM/.venv/lib/python3.11/site-packages/jax/_src/numpy/array_methods.py:617\u001b[39m, in \u001b[36m_unimplemented_setitem\u001b[39m\u001b[34m(self, i, x)\u001b[39m\n\u001b[32m 613\u001b[39m \u001b[38;5;28;01mdef\u001b[39;00m\u001b[38;5;250m \u001b[39m\u001b[34m_unimplemented_setitem\u001b[39m(\u001b[38;5;28mself\u001b[39m, i, x):\n\u001b[32m 614\u001b[39m msg = (\u001b[33m\"\u001b[39m\u001b[33mJAX arrays are immutable and do not support in-place item assignment.\u001b[39m\u001b[33m\"\u001b[39m\n\u001b[32m 615\u001b[39m \u001b[33m\"\u001b[39m\u001b[33m Instead of x[idx] = y, use x = x.at[idx].set(y) or another .at[] method:\u001b[39m\u001b[33m\"\u001b[39m\n\u001b[32m 616\u001b[39m \u001b[33m\"\u001b[39m\u001b[33m https://docs.jax.dev/en/latest/_autosummary/jax.numpy.ndarray.at.html\u001b[39m\u001b[33m\"\u001b[39m)\n\u001b[32m--> \u001b[39m\u001b[32m617\u001b[39m \u001b[38;5;28;01mraise\u001b[39;00m \u001b[38;5;167;01mTypeError\u001b[39;00m(msg.format(\u001b[38;5;28mtype\u001b[39m(\u001b[38;5;28mself\u001b[39m)))\n", - "\u001b[31mTypeError\u001b[39m: JAX arrays are immutable and do not support in-place item assignment. Instead of x[idx] = y, use x = x.at[idx].set(y) or another .at[] method: https://docs.jax.dev/en/latest/_autosummary/jax.numpy.ndarray.at.html" - ] + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAj4AAAHECAYAAAA50yBOAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjcsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvTLEjVAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAMbVJREFUeJzt3Ql0VGWa//EnKySYBJA1sqqI7FsgQsuiBgFtFARE7Zaw2Q3jxtA0I6OEBplhZBMUMKCNEB0QkUWPDaggiI0RIUADMs1iI4QlBGRfE5L8z/N2V/5JSCCpupW6qfv9nHMPubfurftWUVC/vGtATk5OjgAAADhAoK8LAAAAUFoIPgAAwDEIPgAAwDEIPgAAwDEIPgAAwDEIPgAAwDEIPgAAwDEIPgAAwDEIPgAAwDEIPgAAwDEIPkXYuHGj9OzZU6KjoyUgIEBWrlxZouuvXr0qAwcOlGbNmklwcLD06tXrhnOWL18uXbt2lapVq0pkZKS0b99evvjiCwtfBQAAyIvgU4RLly5JixYtZPbs2W5dn5WVJWFhYfLSSy9JXFxckeFKg8+qVaskJSVFHnjgARO2tm/f7mHpAQBAYQJYpPTWtMZnxYoV+Wptrl27Jq+++qosXrxYzp49K02bNpU33nhDunTpcsP1WvOj5xSn1qhJkybSv39/SUhIsPx1AADgdNT4uOmFF16Q5ORk+eijj2Tnzp3Sr18/6d69u+zfv9/t58zOzpYLFy5I5cqVLS0rAAD4J4KPGw4fPizvv/++LF26VDp27Ch33XWXjBo1Su6//35z3F1Tp06VixcvypNPPmlpeQEAwD8F/+tPlMCuXbtMH5577rkn33Ft/rr99tvdes5FixbJ+PHj5dNPP5Vq1apZVFIAAJAXwccNWisTFBRkOiTrn3nddtttJX4+bS4bOnSoqUEqqiM0AADwHMHHDa1atTI1Punp6aapyxPaOXrw4MEm/Dz66KOWlREAANyI4HOTWp0DBw7k7h88eFB27NhhOh5rE9dvfvMbGTBggEybNs0EoZMnT8q6deukefPmuQFmz549kpGRIadPnzadlvV61bJly9zmrfj4eJk5c6bExsZKWlqaOa7D4KOionzyugEA8GcMZy/Chg0bzLw6BWlQWbBggWRmZsrEiRMlKSlJjh49KlWqVJH77rvP9NPRSQtVvXr15NChQzc8h+st16Hv33zzTZH3AAAA1iL4AAAAx2A4OwAAcAyCDwAAcAw6Nxcye/KxY8ckIiLCLFUBAADsT3vu6EAiXVw8MLDoeh2CTwEaemrXru3rYgAAADekpqZKrVq1inyc4FOA1vS43rjIyEhfFwcAABTD+fPnTcWF63u8KASfAlzNWxp6CD4AAJQtt+qmQudmAADgGAQfAADgGDR1AQBgAV3DUWf1h3eEhITcsDC4Owg+AAB4OIxa11o8e/asr4vi9ypWrCg1atTwaLoZgg8AAB5whZ5q1apJeHg4c8B5KVxevnxZ0tPTzX7NmjXdfi6CDwAAHjRvuULP7bff7uvi+LWwsDDzp4Yffb/dbfaiczMAAG5y9enRmh54n+t99qQvFcEHAAAP0bxVdt5ngg8AAHAMgg8AAHAMgg8AAHAMgg/8c9hjxnW3N70eAPxdly5dZMSIEV6/j/6fmpCQYIag68isuLg42b9/v/gKw9nhV/QfWN/EZEk5dMbt54ipW0mWDmtPZ0UAsMDkyZPlrbfekoULF0r9+vVl7Nix0q1bN9mzZ4+UL19eShs1PvArVzKzPAo9auuhM+Z5AMAXtc6lUWM9cOBA+eabb2TmzJnmlzzdfv75Z6+8FzNmzJDXXntNHn/8cWnevLkkJSXJsWPHZOXKleacjIwMeeGFF0yNkAahunXryqRJk8RbqPGB39r6WpyEhxZ/gqvLGVkSM3GtV8sEwP/pL06NE77wyb33TOgm4aG3/mqfOXOm7Nu3T5o2bSoTJkwwx6pWrVroucOGDZMPP/zwps938eLFQo8fPHjQzGytzVsuUVFREhsbK8nJyfLUU0+Z2qDPPvtMPv74Y6lTp46kpqaazVsIPvBbGnqK8x8AADhNVFSUhIaGmgkBde2rm9FgNGrUKLfuo6FHVa9ePd9x3Xc9dvjwYWnQoIHcf//9puZJa3y8iW8FAAAsFBYSZGpefHVvq1WrVs1s3qLNbl27dpWGDRtK9+7d5de//rU8/PDDXrsfwQcAAAtprYU/1TYP86Cpy1WbdOLEiXwLi+p+y5Ytzc+tW7c2TWKrV6+WtWvXypNPPmmaxj755BPxBv/5mwEAAMUWGhpqFlm9FU+aunQUl4afdevW5Qad8+fPy+bNm2X48OG550VGRkr//v3N1rdvX1Pzc/r0aalcubJYjeADAIAD1atXzwQQHc112223mZARGBhoaVOX1n7pXEETJ040/Xhcw9mjo6OlV69e5pzp06eb2qBWrVqZ+y9dutSEpYoVK4o3EHwAAHCgUaNGSXx8vDRu3FiuXLlimps0DFlt9OjRcunSJfnd734nZ8+eNZ2Y16xZkzuHT0REhJnrRyc1DAoKkrZt28qqVasKDWFWIPgAAOBA99xzjxlS7m1a66PNZa5h8wU999xzZistTGAIAAAcg+ADAAAcg+ADAAAcg+ADAAAcg+ADAICHirs4KHz/PhN8AABwU0hIiPnz8uXLvi6KI1z+1/vset/dwXB2AADcpPPO6ER76enpZl8X/dTh27CW1vRo6NH3Wd9vfd/dRfABAMADrvWoXOEH3qOh51aryd8KwQcAAA9oDY8uuaDLOmRmZvq6OH4rJCTEo5oeF4IPAAAW0C9lK76Y4V10bgYAAI5B8AEAAI5B8AEAAI5B8AEAAI5B8AEAAI5B8AEAAI5B8AEAAI5B8AEAAI5B8AEAAI5B8AEAAI5B8AEAAI5B8AEAAI5B8AEAAI5B8AEAAI5B8AEAAI5B8AEAAI5B8AEAAI5B8AEAAI5B8AEAAI5B8AEAAI5B8AEAAI5B8AEAAI5B8AEAAI5B8AEAAI5B8AEAAI5B8AEAAI5h2+AzadIkadu2rUREREi1atWkV69esnfv3ltet3TpUrn33nulfPny0qxZM1m1alWplBcAANifbYPPN998I88//7x8//338tVXX0lmZqY8/PDDcunSpSKv+e677+Tpp5+WIUOGyPbt201Y0m337t2lWnYAAGBPwWJTa9asybe/YMECU/OTkpIinTp1KvSamTNnSvfu3eWPf/yj2X/99ddNaJo1a5YkJiaWSrkBAIB92bbGp6Bz586ZPytXrlzkOcnJyRIXF5fvWLdu3czxoly7dk3Onz+fbwMAAP6pTASf7OxsGTFihPzqV7+Spk2bFnleWlqaVK9ePd8x3dfjN+tLFBUVlbvVrl3b0rIDAAD7KBPBR/v6aD+djz76yPLnHjNmjKlNcm2pqamW3wMAANiDbfv4uLzwwgvy+eefy8aNG6VWrVo3PbdGjRpy4sSJfMd0X48XpVy5cmYDAAD+z7Y1Pjk5OSb0rFixQr7++mupX7/+La9p3769rFu3Lt8x7dysxwEAAILt3Ly1aNEi+fTTT81cPq5+OtoPJywszPw8YMAAueOOO0w/HfXyyy9L586dZdq0afLoo4+aprGtW7fKvHnzfPpaAACAPdi2xuedd94xfW66dOkiNWvWzN2WLFmSe87hw4fl+PHjufsdOnQwYUmDTosWLeSTTz6RlStX3rRDNAAAcI5gOzd13cqGDRtuONavXz+zAQAAlJngA2fT4HslM6vE113OKPk1AADnIPjAlqGnb2KypBw64+uiAAD8jG37+MC5tKbH09ATU7eShIUEWVYmAIB/oMYHtrb1tTgJDy15gNHQExAQ4JUyAQDKLoIPbE1DT3goH1MAgDVo6gIAAI5B8AEAAI5B8AEAAI5B8AEAAI5B8AEAAI5B8AEAAI5B8AEAAI5B8AEAAI5B8AEAAI5B8AEAAI5B8AEAAI5B8AEAAI5B8AEAAI5B8AEAAI5B8AEAAI5B8AEAAI5B8AEAAI5B8AEAAI5B8AEAAI5B8AEAAI5B8AEAAI5B8AEAAI4R7M5Fn332WYmv6dq1q4SFhblzOwAAAN8Fn169epXo/ICAANm/f7/ceeed7twOAADAt01daWlpkp2dXawtPDzcmtICAACUdvCJj48vUbPVb3/7W4mMjHTnVgAAAL5t6nr//fdLdP4777zjzm0AAAAsxaguAADgGG7V+BSUmZlp+vxcvnxZqlatKpUrV7biaQEAAOxR43PhwgXThNW5c2fTf6devXrSqFEjE3zq1q0rzz33nGzZssXa0gIAAJR28Jk+fboJOtrXJy4uTlauXCk7duyQffv2SXJysowbN06uX78uDz/8sHTv3t0MZQcAACiTTV1ak7Nx40Zp0qRJoY+3a9dOBg8eLImJiSYcffvtt9KgQQNPywoAAFD6wWfx4sW5fXu0RkcDTmHBply5cjJs2DDPSggAAGCHUV0hISGyc+dOq8oCAABg7+HsOjnhn//8Z2tKAwAAYOfh7NqJef78+bJ27Vpp06aNVKhQ4YaO0AAAAH4RfHbv3i2tW7c2P+uoroKLkwIAAPhN8Fm/fr01JQEAALBrH5+EhARJSUmxtjQAAAB2DD5HjhyRHj16SK1atWT48OGyevVqycjIsLZ0AAAAdgg+2qFZ1+fSOX0iIiJkxIgRUqVKFenTp48kJSXJ6dOnrSwnAACAb4ezBwYGSseOHWXy5Mmyd+9e2bx5s8TGxsrcuXMlOjpaOnXqJFOnTpWjR496XlIAAABfz+OTly5SOnr0aNm0aZOkpqZKfHy8Wa7CNdMzAABAmQ0+umTFQw89VOgipLpK+5AhQ+TTTz+VUaNGeXIbAAAAS7BkBQAAcAyWrAAAAI7BkhUAAMAxvLpkBQAAgJ2wZAUAAHAMt4LPyJEji3WeLlI6bdo0d24BAABgj+Czffv2fPvbtm0zfX0aNmyY2+QVFBRk+vwAAACU6eCTt3lLOy/rkhULFy6USpUqmWNnzpyRQYMGmVmdAQAA/GY4uzZlTZo0KTf0KP154sSJNHMBAAD/Cj7nz5+XkydP3nBcj124cMHTpwcAALBP8Ondu7dp1lq+fLkcOXLEbMuWLTPLVTzxxBPWlBIAAMAOw9kTExPNWlzPPPOMWbvLPGlwsAk+U6ZMsaKMAAAA9qjxCQ8Plzlz5sgvv/xiRnvpdvr0aXOs4CzOJbVx40bp2bOnREdHm6HxK1euvOn5GzZsMOcV3NLS0jwqBwAA8A8e1/i4aMhp3ry5WOnSpUvSokULGTx4cImazfbu3SuRkZG5+9WqVbO0XAAAwOHBxxt69OhhtpLSoFOxYkWvlAkAADi4qcuOWrZsKTVr1pSuXbvKpk2bbnrutWvXzMi0vBsAAPBPfhV8NOxoZ2sdVaZb7dq1pUuXLmZm6aLoHERRUVG5m14DAAD8k8dNXYcPHzZhQTsR55WTkyOpqalSp04dKS26ZIZr2QzVoUMH+emnn+TNN9+UDz74oNBrxowZk2/tMa3xIfwAAOCfPK7xqV+/fqETGOrILn3M19q1aycHDhwo8vFy5cqZjtB5NwAA4J88Dj5as1OwtkddvHhRypcvL762Y8cO0wQGAADgdlOXq3lIQ8/YsWPNfD4uWVlZsnnzZtPJ2BManvLW1hw8eNAEmcqVK5smNG2mOnr0qCQlJZnHZ8yYYWqZmjRpIlevXpX33ntPvv76a/nyyy89KgcAAHB48NGJCl01Prt27ZLQ0NDcx/RnnX9HZ3T2xNatW+WBBx64IWzFx8fLggUL5Pjx46aPkUtGRob84Q9/MGFIg5jOK7R27dp8zwEAAJzL7eCzfv1686eu0zVz5kyv9I3REVkarIqi4Sev0aNHmw0AAMAro7ref/99T58CAADAvsFHm5xef/11s0xF3qHghZk+fbq7ZQMAAPB98NH+Pa6V2F19fQpT2GgvAACAMhV8XP17Cv4MAABgZ361ZAUAAIDXV2dft26d2dLT0yU7OzvfY/Pnz7fiFgAAAL4PPuPHj5cJEyZITEyMmSGZfj0AAMBvg4+uhq7z6Tz77LPWlAgAAMCufXx0tmRdBR0AAMDvg8/QoUNl0aJF1pQGAADAzk1duhjovHnzzJpYujZWSEhIvseZwBAAAPhN8Nm5c2fuKuy7d+/O9xgdnQEAgF8FHyYwBAAAZQUTGAIAAMfwuMZH5/C5mYSEBE9vAQAAYI/gs2LFinz7unjpwYMHJTg4WO666y6CDwAA8J/gU9jq7OfPn5eBAwdK7969PX16AAAAe/fxiYyMNEtZjB071htPDwAAYK/OzefOnTMbAACA3zR1vfXWW/n2c3Jy5Pjx4/LBBx9Ijx49PH16AAAA+wSfN998M99+YGCgVK1aVeLj42XMmDGePj0AAIB9go+O4AIAAPDbPj6HDx8u0flHjx515zYAAAC+Dz5t27aV3//+97Jly5Yiz9GOze+++640bdpUli1b5kkZAQAAfNfUtWfPHvmv//ov6dq1q5QvX17atGkj0dHR5uczZ86Yx3/88Udp3bq1TJ48WR555BFrSgsAAFDaNT633367TJ8+3YzemjVrljRo0EBOnTol+/fvN4//5je/kZSUFElOTib0AAAA/+jcHBYWJn379jUbAACA305gqGtwaa0OAACA3wefI0eOmAkKa9WqJcOHD5fVq1dLRkaGtaUDAACwQ/CZP3++pKWlyeLFiyUiIkJGjBghVapUkT59+khSUpKcPn3aynICAAD4dq0unaW5Y8eOZuTW3r17ZfPmzRIbGytz5841o7w6deokU6dOZR4fAADgf4uUNmrUSEaPHi2bNm2S1NRUs2zFt99+a2qFAAAAytyoLp2nRxcirVy5spw8edIEm4YNG0qTJk3ynafrdQ0ZMsRsAAAAZa7G57333jOTFcbExMg777wjvXv3lnXr1slTTz1lHgMAAPCbGp+33nrLzMh85coVqVOnjlmgVGt2dHmKzp07y9ChQ71XUgAAgNIMPsHBwWbSQt3uvvtuE3pUVFSUBAQEeFoWAAAA+zR1BQUFydWrV83P33zzTe7xixcvWl8yAAAAXwaftWvXSrly5XJreVwuX74s8+bNs7psAAAAvmvqyht28qpWrZrZAAAA/HaRUhdt/tq5c6ekp6dLdnZ2vscee+wxK24BAADg++CzZs0aGTBggJw6deqGx7TDc1ZWlqe3AAAAsMfMzS+++KL069dPjh8/bmp78m6EHgAA4FfB58SJEzJy5EipXr26NSUCAACwa/Dp27evbNiwwZrSAAAA2LmPz6xZs0xTl67Z1axZMwkJCcn3+EsvveTpLQAAAOwRfHTl9S+//FLKly9van7yzuCsPxN8AACA3wSfV199VcaPHy+vvPKKBAZ63HIGAADgNR4nlYyMDOnfvz+hBwAA2J7HaSU+Pl6WLFliTWkAAADs3NSlc/VMnjxZvvjiC2nevPkNnZunT5/u6S0AAADsEXx27dolrVq1Mj/v3r0732N5OzoDAACU+eCzfv16a0oCAABg1z4+CQkJkpKSYm1pAAAA7Bh8jhw5Ij169JBatWrJ8OHDZfXq1WaEFwAAgN8Fn/nz50taWpqZwDAiIkJGjBghVapUkT59+khSUpKcPn3a2pICAAD4cji7zt3TsWNHM6pr7969snnzZomNjZW5c+dKdHS0dOrUSaZOnSpHjx71tJwAAAAes3TWwUaNGsno0aNl06ZNcvjwYTPHj67hpbVCAAAAZX5UV1GqVasmQ4YMMRsAAIBfBJ+RI0cWelzn8NGFSxs0aCCPPfaYVK5c2dNbAQAA+Db4bN++XbZt22ZmcG7YsKE5tm/fPgkKCpJ7771X5syZY8LRX//6V2ncuLGntwMAAPBdH5/HH39c4uLi5NixY2ZeH910qHvXrl3l6aefNh2btZPzv//7v3t6KwAAAN8GnylTpsjrr78ukZGRuceioqLkT3/6kxntFR4e7vZkhxs3bpSePXuaEWLadLZy5cpbXrNhwwZp3bq1lCtXTu6++25ZsGBBie8LAAD8k8fB59y5c5Kenn7D8ZMnT8r58+fNzxUrVnRrcsNLly5JixYtZPbs2cU6/+DBg/Loo4/KAw88IDt27DBzCw0dOtQsoAoAABBsRVPX4MGDZdq0adK2bVtzbMuWLTJq1Cjp1auX2f/hhx/knnvuKfFz68zQuhVXYmKi1K9f35TFNbxe+xa9+eab0q1btxLfHwAA+BePg49OVqj9d5566im5fv36P580ONjM4aOBQ2kn5/fee0+8LTk52fQ3yksDj9b8FOXatWtmc3HVUgEAAP/jcfC57bbb5N133zUh5x//+Ic5duedd5rjLi1btpTSoEtoVK9ePd8x3dcwc+XKFQkLC7vhmkmTJsn48eNLpXwAAKAM9vHZuXOnZGdn5zumQad58+Zmyxt61I8//phbG2Q3Y8aMMf2UXFtqaqqviwQAAOwUfFq1aiW//PJLsc9v3769WcLC22rUqCEnTpzId0z3dcRZYbU9Skd/6eN5NwAA4J/caurKycmRsWPHmqHqxeHOiC53aMBatWpVvmNfffWVOQ4AAOBW8NEJCXU19uLS4FFUjcvNXLx4UQ4cOJBvuLoOU9flL+rUqWOaqXSCxKSkJPP4sGHDZNasWWahVB1p9vXXX8vHH38sf/nLX0p8bwAA4H/cCj46SWBp2Lp1q5mTp+C6YDpiTCcmPH78eL4mNB3KriFHR5nNnDlTatWqZUaTMZQdAAB4dXV2K3Tp0sU0qxWlsFmZ9RpdPwwAAMDymZsBAADKCoIPAABwDIIPAABwDI+CT2Zmpjz00EOyf/9+60oEAABgx+ATEhJiZnEGAABwRFPXb3/7W/nzn/9sTWkAAADsPJxd1+CaP3++rF27Vtq0aSMVKlTI9/j06dM9vQUAAIA9gs/u3buldevW5ud9+/bleywgIMDTpwcAALBP8Fm/fr01JQEAAPAyhrMDAADHsGTJinXr1pktPT1dsrOz8z2m/X8AAAD8IviMHz9eJkyYIDExMVKzZk369QAAAP8NPomJiWax0GeffdaaEgEAANi1j09GRoZ06NDBmtIAAADYOfgMHTpUFi1aZE1pAAAA7NzUdfXqVZk3b56ZwLB58+ZmGYu8mMAQAAD4TfDRtbpatmyZO5lhXnR0BgAAdsIEhkAhLmdkuXVdWEgQgR8A/H0eH8DfxExc6951dSvJ0mHtCT8A4E/BZ+TIkcU+lz4+KCu0tkaDy9ZDZ9x+Dr32SmaWhIfyOwUA2JFb/ztv3769WOfxWy/KEv28am2NBhd3msbcrSUCANg8+NCvB8WRk5PjdojwZfihtgYA/Bf/w8NroadvYrKkeNBsBACA1ejjA6/Qmh5PQ4/2t9F+NwAA2KqPz7Zt2+T69evSsGFDs79v3z4JCgqSNm3aWFNKlGlbX4uT8NCSBxiGhgMAbNfHR2t0IiIiZOHChVKpUiVz7MyZMzJo0CDp2LGjdSVFmaWhh34zAAC/WKtr2rRpMmnSpNzQo/TniRMnmscAAAD8JvicP39eTp48ecNxPXbhwgVPnx4AAMA+wad3796mWWv58uVy5MgRsy1btkyGDBkiTzzxhDWlBAAAsIDHHS8SExNl1KhR8swzz0hmZuY/nzQ42ASfKVOmWFFGAAAA39b4JCQkSEpKioSHh8ucOXPkl19+MaO9dDt9+rQ5VqFCBWtLCwAA4IsaH23S6tGjh4SGhkrPnj3l8ccflwcffNDsAwAA+FWNz/z58yUtLU0WL15shrO//PLLUqVKFenTp48kJSWZWh8AAAC/6dwcGBho5uqZPHmy7N27VzZv3iyxsbEyd+5ciY6Olk6dOsnUqVPl6NGj1pUYAADAV6O68mrUqJGMHj1aNm3aJKmpqRIfHy/ffvutqRUCAADwNa9Np1u1alUzsks3AAAAvwk+V69elZ07d0p6erpkZ2fne+yxxx6z4hYAAAC+Dz5r1qyRAQMGyKlTp254TBeYzMrK8vQWAAAA9ujj8+KLL0q/fv3k+PHjprYn70boAQAAfhV8Tpw4ISNHjpTq1atbUyIAAAC7Bp++ffvKhg0brCkNAACAnfv4zJo1yzR16bD1Zs2aSUhISL7HX3rpJU9vAQAAYI/go3P0fPnll1K+fHlT86Mdml30Z4IPAADwm+Dz6quvyvjx4+WVV14xMzkDAADYlcdJJSMjQ/r370/oAQAAtudxWtFlKZYsWWJNaQAAAOzc1KVz9egipV988YU0b978hs7N06dP9/QWAAAA9gg+u3btklatWpmfd+/ebUWZAAAA7Bl81q9fb01JAAAA7Bh8dKbm4tDh7NOmTXPnFgAAAPYIPtu3b8+3v23bNrl+/bo0bNjQ7O/bt0+CgoKkTZs21pQSAADAV8Enb/OWdl6OiIiQhQsXSqVKlcyxM2fOyKBBg6Rjx45WlBEAAMAew9m1KWvSpEm5oUfpzxMnTqSZCwAA+FfwOX/+vJw8efKG43rswoULnj49AACAfYJP7969TbPW8uXL5ciRI2ZbtmyZDBkyRJ544glrSgkAAGCH4eyJiYkyatQoeeaZZyQzM/OfTxocbILPlClTrCgjAACAPYJPeHi4zJkzx4Scn376yRy76667pEKFClaUDwAAwD7Bx0WDji5ZAQAAYFcsqQ4AAByD4AMAAByD4AMAABzD9sFn9uzZUq9ePSlfvrzExsbKDz/8UOS5CxYsMOuD5d30OgAAANsHnyVLlpgFUceNG2fWA2vRooV069ZN0tPTi7wmMjJSjh8/nrsdOnSoVMsMAADsy9bBR9cBe+6558wEiY0bNzZzBunw+fnz5xd5jdby1KhRI3erXr16qZYZAADYl22DT0ZGhqSkpEhcXFzuscDAQLOfnJxc5HUXL16UunXrSu3ateXxxx+XH3/88ab3uXbtmll2I+8GAAD8k22Dz6lTpyQrK+uGGhvdT0tLK/Sahg0bmtqgTz/9VD788EPJzs6WDh06mGU0iqILrEZFReVuGpgAAIB/sm3wcUf79u1lwIAB0rJlS+ncubNZP6xq1aoyd+7cIq8ZM2aMnDt3LndLTU0t1TIDAIAyOHOz1apUqSJBQUFy4sSJfMd1X/vuFEdISIi0atVKDhw4UOQ55cqVMxsAAPB/tq3xCQ0NlTZt2si6detyj2nTle5rzU5xaFPZrl27pGbNml4sKQAAKCtsW+OjdCh7fHy8xMTESLt27WTGjBly6dIlM8pLabPWHXfcYfrpqAkTJsh9990nd999t5w9e9YsnKrD2YcOHerjVwIAAOzA1sGnf//+cvLkSUlISDAdmrXvzpo1a3I7PB8+fNiM9HI5c+aMGf6u51aqVMnUGH333XdmKDwAAEBATk5Ojq8LYSc6nF1Hd2lHZ50MEe65nHFdGid8YX7eM6GbhIfaOmN7zGmvFwDK6ve3bfv4AAAAWI3gAwAAHIPgAwAAHIPgAwAAHIPgAwAAHIPgAwAAHIPgAwAAHIPgAwAAHIPgAwAAHIPgAwAAHIPgAwAAHIPgAwAAHIPgAwAAHIPgAwAAHIPgAwAAHIPgAwAAHIPgAwAAHIPgAwAAHIPgAwAAHIPgAwAAHIPgAwAAHIPgAwAAHIPgAwAAHIPgAwAAHIPgAwAAHIPgAwAAHIPgAwAAHIPgAwAAHIPgAwAAHIPgAwAAHIPgAwAAHIPgAwAAHCPY1wWAveXk5MiVzKwSX3c5o+TXAADgbQQf3DT09E1MlpRDZ3xdFAAALEFTF4qkNT2ehp6YupUkLCTIsjIBAOAJanxQLFtfi5Pw0JIHGA09AQEBXikTAAAlRfBBsWjoCQ/l4wIAKNto6gIAAI5B8AEAAI5B8AEAAI5B8AEAAI5B8AEAAI5B8AEAAI5B8AEAAI5B8AEAAI5B8AEAAI5B8AEAAI5B8AEAAI5B8AEAAI5B8AEAAI5B8AEAAI5B8AEAAI5B8AEAAI5B8AEAAI5B8AEAAI5B8AEAAI4R7OsCwLtycnLkSmaWW9deznDvOgAA7Irg4+ehp29isqQcOuProjiKu4ExLCRIAgICLC8PAOD/I/j4Ma3psSL0xNStZL6UUTwxE9e6d13dSrJ0WHvCDwB4EcHHIba+Fifhoe6FF2oiivceaXDZ6kHQ1Gs1rIaH8s8SALzF9v/Dzp49W6ZMmSJpaWnSokULefvtt6Vdu3ZFnr906VIZO3as/Pzzz9KgQQN544035JFHHhGn09DDF6r3aDDU2hp3+lNp05i7tUQAAD8a1bVkyRIZOXKkjBs3TrZt22aCT7du3SQ9Pb3Q87/77jt5+umnZciQIbJ9+3bp1auX2Xbv3i1lva/O5Yzrbmx0Ti7t8KPhsuQbzYgAUFoCcvRb1aZiY2Olbdu2MmvWLLOfnZ0ttWvXlhdffFFeeeWVG87v37+/XLp0ST7//PPcY/fdd5+0bNlSEhMTi3XP8+fPS1RUlJw7d04iIyN9PrJK/3b6JSbLnuPnPSrDngndqPGxKQ2pjRO+8LhJEgDKijAvdKEo7ve3bb8JMzIyJCUlRcaMGZN7LDAwUOLi4iQ5ObnQa/S41hDlpTVEK1euLPI+165dM1veN85qGnpcX2y+QOfksoMmLwBOsMeHv4zbNvicOnVKsrKypHr16vmO6/7f//73Qq/RfkCFna/HizJp0iQZP3682F3jmpH/GvFT8mvpnOz/HaMBAGU8+JQWrVHKW0ukNT7anGb1F5umW0+fg/DinzzpGA0AZVGYD1shbBt8qlSpIkFBQXLixIl8x3W/Ro0ahV6jx0tyvipXrpzZSqPTK1AUPiMA4PBRXaGhodKmTRtZt25d7jHt3Kz77du3L/QaPZ73fPXVV18VeT4AAHAWW/+KqU1Q8fHxEhMTY+bumTFjhhm1NWjQIPP4gAED5I477jD9dNTLL78snTt3lmnTpsmjjz4qH330kWzdulXmzZvn41cCAADswNbBR4ennzx5UhISEkwHZR2WvmbNmtwOzIcPHzYjvVw6dOggixYtktdee03+8z//00xgqCO6mjZt6sNXAQAA7MLW8/j4gjfm8QEAAPb4/rZtHx8AAACrEXwAAIBjEHwAAIBjEHwAAIBjEHwAAIBjEHwAAIBjEHwAAIBjEHwAAIBjEHwAAIBj2HrJCl9wTWStM0ACAICywfW9fasFKQg+BVy4cMH8Wbt2bV8XBQAAuPE9rktXFIW1ugrIzs6WY8eOSUREhAQEBOR7rG3btrJly5abXl/UOZpENUylpqb6zRpgxXk/ytK9rXhOd5+jJNcV91xPPq+Kz6z97+urz2xJr+Ez6z4+s8W/XuOMhp7o6Oh8C5gXRI1PAfpm1apVq9DHgoKCbvmP6Vbn6GP+8g+yOO9HWbq3Fc/p7nOU5LrinmvF51XxmbXvfX31mS3pNXxm3cdntmTX36ymx4XOzSXw/PPPW3KOv/Dla/XGva14TnefoyTXFfdcPq/2eb3euq+vPrMlvYbPrPv4zFpfLpq6SolWwWoSPXfunN/8JgL/xmcWZQ2fWRQHNT6lpFy5cjJu3DjzJ1AW8JlFWcNnFsVBjQ8AAHAManwAAIBjEHwAAIBjEHwAAIBjEHwAAIBjEHwAAIBjEHxsonfv3lKpUiXp27evr4sC3JIuCdClSxdp3LixNG/eXJYuXerrIgFFOnv2rMTExEjLli2ladOm8u677/q6SPAhhrPbxIYNG8waIwsXLpRPPvnE18UBbur48eNy4sQJ80WSlpYmbdq0kX379kmFChV8XTTgBllZWXLt2jUJDw+XS5cumfCzdetWuf32231dNPgANT42ob8968KoQFlQs2ZNE3pUjRo1pEqVKnL69GlfFwsocn0nDT1KA5D+vs/v/M5F8LHAxo0bpWfPnmZFWF3RfeXKlTecM3v2bKlXr56UL19eYmNj5YcffvBJWQGrP7MpKSnmN2pdFRuw6+dVm7tatGhhFqH+4x//aMI6nIngYwGtOtV/UPoPrzBLliyRkSNHmqnUt23bZs7t1q2bpKenl3pZASs/s1rLM2DAAJk3b14plRxOZMXntWLFivK3v/1NDh48KIsWLTJNtXAo7eMD6+hbumLFinzH2rVrl/P888/n7mdlZeVER0fnTJo0Kd9569evz+nTp0+plRXw5DN79erVnI4dO+YkJSWVannhbJ78H+syfPjwnKVLl3q9rLAnany8LCMjwzQFxMXF5R4LDAw0+8nJyT4tG+DuZ1a/fwYOHCgPPvigPPvssz4sLZyuOJ9Xrd3RwSNKV27XprOGDRv6rMzwLYKPl506dcr0f6hevXq+47qvo2Fc9B9pv379ZNWqVaYNmlAEO39mN23aZJoXtK+FdnLWbdeuXT4qMZysOJ/XQ4cOSceOHU0TmP754osvSrNmzXxUYvhasK8LgH9au3atr4sAFNv9998v2dnZvi4GUCzt2rWTHTt2+LoYsAlqfLxMRw7oUMqCHel0X4cBA3bDZxZlCZ9XlBTBx8tCQ0PN5G7r1q3LPaa/Ket++/btfVo2oDB8ZlGW8HlFSdHUZYGLFy/KgQMHcvd1uKRWq1auXFnq1KljhlnGx8ebKdO1ynXGjBlmeOagQYN8Wm44F59ZlCV8XmEpXw8r8wc6DF3fyoJbfHx87jlvv/12Tp06dXJCQ0PN0Mvvv//ep2WGs/GZRVnC5xVWYq0uAADgGPTxAQAAjkHwAQAAjkHwAQAAjkHwAQAAjkHwAQAAjkHwAQAAjkHwAQAAjkHwAQAAjkHwAQAAjkHwAWC5gIAAWblyZbHO/dOf/iQtW7YUu1iwYIFUrFix1Mqn99P3S7cRI0aIt+nrcd1P17QCnIbgA8BtRYWC48ePS48ePcQfjBo1Kt/K394QGRlp3rPXX39dSuP16L1q1arl9XsBdsTq7ABKTJf4y8rKKvLxGjVqiC9lZGRIaGioJc912223mc2btPaltN4z1+sJCgoqlfsBdkOND+DnunTpIi+88ILZoqKipEqVKjJ27FgTXlw++OADiYmJkYiICPMF/Mwzz0h6enru4xs2bDBfzqtXr5Y2bdpIuXLl5MMPP5Tx48fL3/72t9ymE222Kayp68iRI/L0009L5cqVpUKFCuZemzdvLrLM7733njRq1EjKly8v9957r8yZM6dYr1GbivT1devWzRyfPn26NGvWzNyzdu3a8m//9m9y8eLFfNdqmevUqSPh4eHSu3dv+eWXX25aq6X3Ktgk1atXLxk4cGDuvpa3QYMGpvzVq1eXvn37ihXNhdoE53qPf/75Z3POxx9/LB07dpSwsDBp27at7Nu3T7Zs2WLeYw04WvN28uTJEt8f8FfU+AAOsHDhQhkyZIj88MMPsnXrVvnd735nvuyfe+4583hmZqZpZmnYsKEJPCNHjjRf5KtWrcr3PK+88opMnTpV7rzzTvOl/oc//EHWrFkja9euNY9rsCpIg0bnzp3ljjvukM8++8wEq23btkl2dnahZf3f//1fSUhIkFmzZkmrVq1k+/btppwaXuLj42/6GocPHy6bNm3KPRYYGChvvfWW1K9fX/7xj3+Y4DN69OjcIKXhS9+XSZMmmfCir2XcuHHiCX1/X3rpJRMmO3ToIKdPn5Zvv/1WvEXLq3119O9z8ODBJrRqgJ05c6YJc08++aR5P9955x2vlQEoSwg+gANobcebb75pagg03Ozatcvsu4KPfmG6aKjRsKC1Bxpa8jbzTJgwQbp27Zq7r48FBwfftJlm0aJFpsZBayG0xkfdfffdN/0inzZtmjzxxBNmX0PLnj17ZO7cuTcNPlrDMnny5HzH8tbM1KtXTyZOnCjDhg3LDT4aDrp3727CkLrnnnvku+++MwHIXYcPHzYh7de//rUJIHXr1jUBzpt9dlw1XC+//LKpWdM+Sb/61a/MMQ12rloiADR1AY5w3333mdDj0r59e9m/f39uP52UlBTp2bOnqTXQL2utoXF9ieelzScltWPHDvPF7wo9N3Pp0iX56aefzJe1qy+KbhpY9PjNaBNcQVoT9dBDD5naJn1dzz77rGnKunz5snn8//7v/yQ2NjbfNfreeEKDoYYdDZB6P63Bct3PG5o3b577szarKW3ey3ssb7Ml4HQEH8DhNGxojYGOLNIvaa2ZWbFiRW4n4by0JqOktO9Jcbn637z77rsmMLm23bt3y/fff3/TawuWTfvAaK2LBoNly5aZcDd79uxCX1dJaPNZ3v5RrqZCFw1Y2pS3ePFiqVmzpmlmatGihZw9e1Y8VViH8pCQkNyfXeG24LGimhUBJyL4AA5QsCOxhghtGtKRPX//+99NLcj//M//mE6y2pm4uDUEOnLqZqO7lAYPDS/a1+VWtHYiOjra9MfR5rC8mzZ5lYQGHf3C12YzrfHSZqxjx47lO0c7UBf23txM1apVzXBwF339Gszy0ua/uLg40/S2c+dOE8K+/vprKakTJ07k/qzNhQU7ZgMoOYIP4ADaZKUdlvfu3WtqIt5++23TH0Rp85YGGD2mgUM7IBd3PhntN3Pw4EETbE6dOiXXrl274Rztc6J9gLTzsHY81ntoDUxycnKhz6kjxbSzsfYz0hFK2h/p/fffNyO0SkLDktbEuF6XdjZOTEzMd452Qtb+PNphW5v+tEP1rfr3PPjgg/KXv/zFbBoatUN13tqczz//3JRd35NDhw5JUlKSCWDat6qktB+WBjFtktP7KP07LDjyDEDxEXwABxgwYIBcuXJF2rVrJ88//7wJPTqyy1WDoZ1fly5dKo0bNzY1PxoEiqNPnz6mc/ADDzxgnkdDVUEaqr788kupVq2aPPLII6b/id6jqHlkhg4daoaza9jRc7W/kZavpDU+2rykYemNN96Qpk2bmmY8DVR5aU2QNqtpJ2c9X8v52muv3fR5tSO4drLW91TLpn159PXnHXK+fPlyE5C0RknDlr4vTZo0kZLS9/app54yfat0wkH9u9PmuoI1TACKLyCnYGM1AL+i887oPDQsT2BPGup09FnBPkDaN0f7WmlNmTdobZ3etzSWyQDshBofAPCxc+fOmdFr//Ef/+H1e/33f/+3uVfBEXuAUzCPDwD4kDYX3n///ebnvIujeovOY6STGiptngSchqYuAADgGDR1AQAAxyD4AAAAxyD4AAAAxyD4AAAAxyD4AAAAxyD4AAAAxyD4AAAAxyD4AAAAcYr/B2C5m3j1gZ4jAAAAAElFTkSuQmCC", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" } ], "source": [ @@ -180,7 +96,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 11, "id": "d5e22881", "metadata": {}, "outputs": [ @@ -190,7 +106,7 @@ "31" ] }, - "execution_count": 4, + "execution_count": 11, "metadata": {}, "output_type": "execute_result" } @@ -200,6 +116,22 @@ "\n", "31" ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "336e133f", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "7d6cd768", + "metadata": {}, + "outputs": [], + "source": [] } ], "metadata": { diff --git a/readme.png b/readme.png index 1ec42f09932c62c40e69927327a733526661c1d3..d9a997fb31f170c33b3079f8e6ecc7142975351a 100644 GIT binary patch literal 13374 zcmd^lbySvn*XB(Kq9Ps<5hPVWTImLrZc*t}$%pRdu|=hkE~S-{Zk}V&ARUs@-CeV9 z-}9dD{l1xZ)|$0u&6-*B$2reA-2MCQeeLVI_V!m-l_NbteF8xcQU!Tw4Fti*A_(5d zaU%E&f7`$~{D?TqXglArGjnz`aWq9#Oq}hl?VPPG?=!fXIyza}+46E-zs4oV!C>L+ zZ0{t(&297d8C-Ub=G@AkF5iW_kl4%XI3Wm`3HlFShGe=Wg7AJ&kiK!(J#Jx;#QE;M zgxVFqx|_$DIr3#SDMG}h7;e5!BEqvdM<~4zk=I8s_FXzARRnz^Cn!>z(g>Kn6Fo+SsGjV z^YMV+(E+pbBHSxHl^8(;R_a??S~i^-D+9Wb9GwDYHO`?TOwV|imeJdZ#>mDpr>PXU za`pQ~JNY0waa$214E{#c`_4k7G{@`RrZ4G%fABM&VJEiGYy5CN4|K-oa$&Db(w!R5mSkld58e^_9Y9+at>9U~_ph#(p6QABBP ztOv`xni}0MlUgf9aD=>hb1GWI<-^zS-y7<~SW-2!bsby9x1xN#MRp3T+Ee^#L|y08 zw(JjB+F2h|@6LY2&xQR^KlhJBx<+~mL9nLBv5;;Nj0k*(_Ve4g)Fh=(2m5QSXK84x zyN@6(|uUNFyA_vhkpt2?=rBUN$K)>X9Xs(pK^7)mN=C{H4-LjD!cwbfh=x z*;O%JzkWUcLG9^G-J%ccvt3F#hE=mO?>%T`@xJI%4-tL;s`U9&0h4dg@^}7+6B%5} z69^7NN-J^~HjurbIbO79zGLQI44 zRTXHi=IZ)zFYM*Z*yW9iXvsCGVys$1i7H+c)z6b%)JyzL;j~TDP1CNCRLsrC&V=ge)WEOpmdjx z)K3w!;n*>i(-TN9Y0KAfD>6 zHrZa`X1B9o=f_Saz8?HDXR3KSfEdBPu;9(={^~K?Z<@`gg3Iau27nNK8^^UP#=+4f zw)o@S6H7M38s94Gp-LS<7m|*L`}Hh|;WJZtcxbSawV5@&`mvj&f{Tu>8PfU$_fUxy z;AgzX#l^+4oSK>%-G+nkYZ2zUazz@D(7N~$0dY(Wj@P_JzBz_p$=N{f?NP;P1aS`! z-0ou7>`;@e=PVhx|5?*lDf?|2iT0Pg7Jea@*}EOBY6PQCNyx82>OqK`2S$r~6=%{~ z4_8-}9X7mWSFy-4tn$p1OtwLzIZ0jjBP9g|6Q^ErIJ;I(sVVgZe%k?|GUu7l{ZTT% zV5g~;=s)Q$3Rnl@a8hfY_V#(v?CLBpHa9lj!m8DmyUb~Cl@R{pPHuxPkp=EpK9Tt+)rz1pMX60*q55`CSd{W~%1baKNSg{lV8)ey%7LcV^(h{-? z`0O64MTpC--c?NSS`^;d{AoAwn!7p}ih#%TyBFi_sY=HZ)|v%}>+YzDGhLQ_6D;nv z8TBLMCxjmjrzmV$$Ntxo*1dU*t!qq55m7O?r%wsXma6ww#myH7%6{j{82yn&Le1gs z9J5$Utf1Qc`}b*Bq4r=|++|2d0$LI!?GCeb3d|nWeyy!ldid~Rf)C^=&)!MLg?@eo zC8fp!fc?9ChLlJ_e*VLzsOzh-!~3$jCAJ154)yi*8YQ-aVteali`-ecCPh|Vj>XSO zFGg7Q!9{T+E_mvD_*RnDK@z(eO~YRM+be0VZADDX%n>{Mu$u$T?7K+0*wOyl>dsVR zBTA5NJRy2UrvQJ<&B$(VtIi_%Ypekw(6(kfR4( zPUX%smmC}%JT_+-&J2&PhS#_9(aiy&wV$_eHE)mjbM## zkzf0tWb?mG^Z%^}_6K|(wMy(|^-itN_uct$H}hA-Dy&Rz%K}~=(A;k#L$$SXr|Bg) zN``$ael`FtGv<|x!G?u}F|o6Ml-E;N4CH*xVGv=4G zLb&xyjbmIZ+?H&2xA-mFrF9D}Oz@4zz6VSICccdqbrZP_ToRcp6GvJS#D5Y>ES$VR zMa2sH{I<$#tKr?P7gvDyGONT~>B%*Ud_m`l@?M)XECjvOF+1T`8;3hTC^LdLjDN-X5yMlG&pjNs&glkNRc#Z>aej22rY& zm+&ANdvI`II$Z6eY5=tbnJ}Tg5I?@t>UUTW770s$Wl4Y3_9xwIO74i}SV5)QK-yHR zt`CBw7vO3voA64GIbhlRcEOBnCTlP9JC z2O1wWIjOE0V^P8y&1)Wa6(Svpiz7U?&)1AHo)XmMGo7|6W=~CfamFpZx`#ko1^X z%HYJ}pTzXxz}1J>br<@JOuxST15NS@8ZI-23>h^?wz$=`W01PG4W? z%)DcH2yA=#I}#NU!JL_yDY($%L&>FYpASh6+s}u9uo>pt0=B5D8Yg5@eD9+sd@h~e zz5~l^2z_Htx%5Y9;9w3w#jKO<&hgCBqa47hQPsiI4L2xJ;p+i>C-;D>(})hEJr>X zR+a0anafZth!-gbaz~@a*9VoZz-I=3>jit!+)O+UE>a(nt4v~I`Z)#_%u7{U3G35s zx6pu~6>+ZHq7hdu8R+PT@43qtNQ47EtMuhLX|wU`0WcDPS>LASw;I>hab^|iadiFw zNw_)|B-s#aSItWAy?mp3YWhld7>k-~)*HA3MKwt&)YlzV!tBeBHu%-%2FeO(gdJ0G zJ+(pf+*7TIt$loju(a@`Zp2DuIH1UH_Sd)9XD?j1no#VvWZ1!{9}J~BmK5UclG7D5 zOk(}`_f?{<{nIIvsqRn1w-*YfJ(H-IWsN zEi`xLdOzG_FZ0~^Fhj6C5oy$rrV^`o4`L9z+7g0fkV{lpAzMqs+<16+!nrdK;%n5j zjHT(G_{H3_nA%PkBJ0L+r-o zk%4^E*E?E0cnP9zwh8``jjsYny0AKDNFS92Toi1vhO>$t3du(HmTH#9C*c1^{E-Zm zPq{URBLItTQDe{_{nh8DUxee!?9pF&bm<7`5;+(D=eQ9)b)a)DS?g0duy9e?H zf_QC3e65LB)$lxeuUUeOk9N)r{ay7+up78AJp1vwkbLVt9u{sAIIbpuXH?~R4~p@_ z?t!i%H-N#l$`~VJXDukglDzsb&=FxH0&?bJ-=UVy!JL5h*EGc1E6% zrKlxJsGf(^$LT^q3l4kpwyuFkqX0R`#KxxhJ&1nkpZSY){!s&WSzeXEq7mA?f&$|2 zS~enu%y(vn?RQ8Xgv~=3;~qF8brnn1?5KVNq$sQS4VzA0aCEqI#0Bsb4h(aq z_fzuI=g(tm&07=rbHN1KYJYc&7)4^S{pBx4Hy!*!;3e}_Z2i-3@0KW(h2ZUN_s{fjz91=nrbgb`PcK~*&Zim?)ils8novQVg)2#8Jv_x5TFJ59bJ z{|5KOhGFpLKF|-AyYOE2A)^z!3RXujXi%L}`*&0n-orjy^-a;d(Ws$-)59ML_%Fpm z=ky}fSAn9ZaI1`@WUNa<&HTJ|rcOaxN?KYBwx`m=$r4m;!kDx8$RC$wqaa!kibA60 zs;wd3Fo_p0UNnj>RbD>W+lqwf_5eQDj6pv%pBh&(@d+lVcD4xLL=_k*||QG zZ#P)s78!C`u6(W?Jx>~0Ut8l1lbmy%ig$}yPXBQ4BQqF%i{%_}RIH0|er}FkE7y=e zs0Lh{xC2UVL+2hQbbacPQjayWO#L!P)$GJZSR>Q6_cDK97ALwsp_#5qi_&o1$@E9{ zK>{b2ewL_dXlPtvXJ>E4pyc{D*5(C5We4Qsd02Nw7%*BTcLo>S)TGSvySg#SY2kZ} zKBIj2Bk2aP*N0X|j&@Ha;ng5>0N+EWr!6w|N^-&?A#Abf_#?PS1g3M{*{X&REJ(sU zfsoyhSO;Mi2EPRLD$#sA{Pmh!M27|rvkOxrk=#aHu>v-s5Po;aqr7C!_$VV#8kX~s z?b(mH4deA;6K}N)LlocA*!V!+;Eo*lo~5RzH-bDO@q>HNja!_34X8)Kt|{c90Do}L zC()kTjfJ4K9w@!tUu@&m!GMC5=--b-0ThUWU=SINQ1DdI0KrYTY=r55ncmpR5<(Q? zln6$k6fg-3_bC(7($H```2HqDx5%m(Y>oCZXj}Wm%uw^bf`rko|xUdgnQ`jc@+)JKFM70-hJo0lW>%PqRR<-TNR& z;o-xSeXs(}45bF0Uqu1tHGr}-CBjY{2m9?r`Z@glfR?E1w+#IbTu=tkrnv8Tsp{lERQQy85=cgy z5H-@8D9O#r$||^^rlRt;u#j)9)N^B&Ddx&9P^tDD1JP(nKi`l*Nze5RvDLcEjsRJW zLseeQ15Yo@1%^W!S|F7Z5qjf;Zo~X0OG|!76@6}d3nj4! zTu_h#!(}4VK7!^$g=_5lT>Y|3U_qA6!~tiCZvNI~4Tz&2Qd&he0<`RR@=B0;!$R`w zEiR{5gOrKu|8vz0ECjP3uP~V2ym#rZA%oNp@s0O37l1_gvaYTelsjKn`ckF`pD_*u z6GX2Q?YlUmg+yve-%6^wU;Ssv$p;NpVareS5XQyFbCNSFU%OLzPDVi%iF``b4F3K^ z-xl-=q7II=&D?i}z|W`~phH~lPeCZ3-*G23KGub!n1oDlgpw@QQ8~TNRi34vEE}M_ z5UiID^&z$sFI0Zgud%3H0R92*<^ zyHul8D=ZS;NJ=7#2*og-K<;Cdl$7Sh;o;}%8yXtcVm8TO$58DBWUw)u_0AO_y|X>J zYBea*+L2*VpZH<@S&zwR^SH`?6b9mJ;C`38F3gVq^$9|b}tJq$j z=j0*3QDXD6n}aO!BKZw?3wIr)vk3d;^JhoqXDtDED3r7qU$Q`#NJxT=VZ??j-8^Et zlI^g%)0+UPeq_m@s2mX|%8;?H4Zmvdnx#78Q*Q9wzTuU%0(L*d^&WQqOh zGH>lae5TOA|D`*r=}Ge=2skJ@(IYii^H| z{aOcbGXbJA;?0{ke&=wlt%Nu#UNfP8&_!X5e+ug0q$CENJmW|EtBu@d&Uv=xEpgZP zK4@moO8=H_Q}U!q0vnq%9dh{i_-JEr72vs^fK1Fd0rnPm5`&~_~@96_)3eR}XhT3WgezKaPgW^f!OqM@l4 zC+5+fc?q&YyT0u5} z2{!T=?xR+20}5G8OfvcT`Ky(i)5$?1259A<185ht?q#_3{8UIsy4vu1OwL$Lbo9)_ z-*RBwf*wa{=zare;S~i1g|toXE!l0j=kjO7eYTU;KQJ%?JZcdS)VgtJLr&8OsTbQ0 z{xcCIun;G_DmalrCD_%k!^4uTbKKS2@58>qhvbDOI0$HwsR36qX6WX_u+c?;A92mNH4ul_n-ZpXD&H4 zXg>#5CKXT_2rp{R|NS)kqr-#tY9Ft)ZJC{y>}r_huP_Mn;Tv%AGSiGj4yG7xmiAI(1P% z85U@A@jaY{o3S)5^ml{uAUvuzGX^cs_I;?XBLR5VI5|}`?%c8K&)wtFFCF+0T7ymkEn(rDQ654_F*qupdx_o(z{@C&wFP{DPH&p+O zc_Q3jflT5hdAOb2=T=MWCJ(?F&jApK?V=&Ya#3$P|29_KJ?N$2P2K7KwE^gU4~diX ztL6ir4WEF335ZfPm_PVI2Kf*Db^L_3(x%Yu5X~jQg1Z_FG&GveGE~Q|6vS}*M=frM z8`kBx;NUX~ii!<=0NRrl(Qtv-pBVh49Nzu;g4P)SG}r~pNt#YebC9C?Bf`3KK(wpk zvTX;;10$$l>BuWDA!+I9O3)LRmX)2H{|HIrK}Sz-G~xkS77ja(S`HJRYU#q3YL159 zokig_5S)34X(cnWNfmzgT`^%{;bEs3K|5P3VoB&Z?<{{I1FKRApxY`UoxKiHrpfbr z8EUP+`p%ATtl&J;F7x~p)g@pC?NwfG1riRZV>SBu*|7v(n|?m%0_(P<@uSujMF2bs zs3wYsaUw2kt-d}vK6|e|4-Xslm)LR31=6rCN3D$2HiGDMD*pQQt41U(e3hl)9|!?} z$Aj9yvy_yP1!}R*8)e<=6wd@0mxgPMpzaurTz|{1C2;wLlilyf$aUkp6362+g9!cn zV z4apH=mNsdFgLD)w`{I6jzqhYmJil@!Mcf}{%Tjf^Ie${)^w!V z8yyZ2H+cCz>#fyiZFZKjOB{$V%Mv!1hE9Ca1Yt?zE%-_fysSmxqy$3Yl!158?!-@LvF@9+g zc)7V2_>iLi=SOYO54>sw8$F>YYxb$NBF=Ch#*(n+XOIw z0rEGD57VtYKymm*>a)X@Zp8iPcd64Xw~~nF zf7$Wn5+I7LEP5S}DPO`M2-=;k^jK?J%gW-g1fpNL{P`r5s@wtRebB=bpo!3q*wlQ* zRD#zP4i06|F|VAPK{lbIBY>ucFc%~yadwM->eZ=N{*m;d?TYuF)2FS!C(5H~Pa1Ok z^<9I%gFs5=S745erS0y4L5TY|yY2OLXVlSw6rK8wfB=Px5IW$X`et8mQo)?+fr+CL z7cgnfiU$OPMJ=|S4nnsGi5%~q;yXr=hevUCag<_s8fqjHFYldMAY$e<5=TFW4o5m=?S@Pu3db9l$YUo^RFS~=DN=|sno zozX*}lo48Q@J}B<9*k0XeSQ5F30lJUtrdV~OTy)nTPw5BO1--54#;OW%HMUopYI~{ z(oth_u?9S+sN2h~n5V61!V;<3n#5Y32(~-y%74>gUxz6+#edN$9RyOon;NL0n+Y{vdqc9MG zLR3U7i8?EL9?fFd1e_RH6?erf_YXG%C9rVct0MReaP&V>>$FM-G(5eHUNg_1x=xaWWBVB=U43P|~!Y z|NF((#@kpLHnQr>V1$O+y7O(UQL2*726QlVi>+CKD87Xf>sc-@^{6FYlm)`;?D_M| zs7DKCL{End{=j2PW40vtLk#~ZXt1y)dawGYza}`x#>%Py5pL@Hj+uo8ZQf1%cy~*q z$m*UNC26UGnOT~Cg{#Wr$B)q=0A|8P=uO*JZD_*5t9LpPO)1@280VUYc`X^V=|X;Z z7aV~UbXp9i&(tVSUbR2ah6$S|N!=mpYWQ~2h;#3MKY$n;2)4_AVP>L1;S^WryaX5- z|Ln};v+n)zDD;c0sY8twuRYj2XoG&WR`RfIzYX90B7V}BLWLWCdDCP``&rtbZoW@W zZaymmDo>#HFO;O=W>=5c%Ie-Fg)BpR>+1C;$@!&5peVTf4(5y|z%52?|Bh#P%Od^6 zXA<5S9mVitR_CE{>kq%E!m&1GM=n(FM8Cl$`YN_O6s`Bnn?o z=o68=Zo1J4n7q5h}5aZOZpA*u?hlVe_%>9@C?xojFa%`Ul00# z`a1`^`f_UJ(_|lEInXg1WqOLwj16nz6+U~e9rI~EpPa9~wB6fU+EVv7DR@7S3%eTz z7!aLia7UTOpKIg! zV*A&#-kWe51DrMfp8ES=AEJAZzT}V+;F;v25_od(kQP(V($Sm8SW@KcX6!kjT z==tW~n%rXr-Rz&a$8b8nr;W(4rp`t6O-)VsMo?Q4+%-Hs%Lxey!|S~U9UD2#oLI60 zoMuIxeCUjfjO22Nal!x*b#2LVWaz|eGvpY5=Dm**%Fq*hDild3u`LB$nn_5gmoOwf z(N`C^furGYxywfqR`5C^(b@#TIPrx{G8{F9jE;7D;+SP@TAFf6wxQL&+IdG!0o@Mh zh)n82$Kz^V6WSsz1dR@Q>G#mkv6JV|&G1^0fAFn_>9FVnVAl~4LQxn2>t6N4ox~#@ z*xslqeNHE^>8@P8D!;q61pUHQrFTVfp4k^;t@bwNSRm4yJ#r>&y0fBL{UB#d!Q$6A z-#p#Vx3C-Nw>!h))Kp9F#{oye*zSFh2Xl2C04d*TN_CAE^P?lCx3jbJz7f6Dw_BtP zSg{B|GISmZW{p&#Wz6Ka{e@n6+GDfNym=0tyx92!9VO;%8R(x%6+q~VEv*5+se>q+ zfFUQr(Mj+#BSF)ha3g{ zRHJ!hJU8c~*Yu6n*4EeycN}3BDUC=GlUhimKIkzO|BNUG{OfuU;&aZflCjO%--F zt)oh`h!V>6xQ^iP(CN67;#=SIE$M}`pLJjh!Js2|gr0ul(iGci8_kwMzxJ-Di|=pT zi7iFDQ5$nTg3Y<0xx)ZU6)YzNm%oFDH1&yne7H0>Y0?!LdWIZl2)anDZxPkbe211~ z7uK@5^VPwW-(fnHa?%s7SGo^UF58^ENSeimwS3u8ZC@jc`FQizEfqg;8|IcO%RBZC z79$-a$<6h;1Rt6!B#-u;5*|1wDfkYnm{|4X*reTdXlu>I;}^6!AN1x8pGrIU_`D|T zjodZizhFG@F06-ajRSHL8J9%}gI@?t%;EyUpo)pf1X%sKPI;$){^JeREz9NhN7OR$ zL8^OQw9tDbtOQr05$pzA5!t;IOZpqM58n~+tyAk>rUXW}-IuAXDJpoJ;tjV5LI}4{{YBzyNhQv~u)A zfB-c@gqeR%3?y{zvY26b1k^|onoue#Nks48Uw|fgMbv1JhZ63exzGgMuxs(d+at;06L5ZF5 zu(YLA)q$*|g+i*m@6>XT(|v-j>UhSS@g=HLpC0t(t}1-V4%Dpv6qK|%*ijv&>oj{> z*Xv`C%l#c;q$LR`BU>;K%-FV;qut!2*-gwrt>6BOvhIM(m3Vi(n+T306i)?`hoO5U z++4euw#G0!{zEZ#0neQ*1|1LM_;;gu46ZvyU)2o<)N6h)K>T2>-2k8zX2tuW(a9tBY zhzSTnG<<*@KEZShjKYhkqpYr@mL1m7<+i;!qI%o$zO|jBwdEZqXLEZ8OFLVBo-3Dm zgt(b393AgFi1PB<{CxzEo&7ytl>%mMxC-Tc1w97@IeHuYN0cFzZiyhmk&4%^YP-hH z{-E+`Tix4V?6rOyWjOxkisShQpRPM=D`6B<4Z|k$4Oi2(`-AL{P>am^q{#T6P~*(Z z7|P`M6&TLvqm_H%+V}JJVJ~hx_-z31Tm#0W`Zxw|Nj^NJ5v&cn9wAO zH#Ikh4tZ^>zR}LS;N|62&2t1nIDE3Rxkd0xDL#IFO1=l3dwh>0%F-WAbaf;1{D~0c z`R$=<&(6fFBn8Kb@=1%w+m^}R=$E|7^Cv-&gXdXUKWh)WdIfy1#}M6L{r!~>zhtJU zpb+?&;%MWq&;F?frM6_o4s!z{QI~H&s$6ar3l0rsi7;}G`t<2jeeFY%Lrynt6y^Cz zp=&(mpYZudsUbC@rKOd6GwT*xE_Zr))FtV~J z9IGA0sY3U4+f2QpG?AXN2X}eHk0D5+b-u%-|5J~xsr1gbHU?1w zmciY5CMLPnTUN9p_9?L<4qtj@;ZYVY_I4wu!+Lf+ms^BR(9>IQ9b%Ml(`NKs(b(Kw z9~zr29cnFYOSBac3nE6;yEzqF6Qp=)X=$r@p5r#Rx0R)(312;T_MVmysc&yh(Tx0X z849!d`X=hGR=Qf4ag|#`a}>sKv_5Zl1;aO_r0f~04Euk3@W zz`$etW&D>fxA)3&hGbI_wo=>EOm1vwwM1k1y`yG_YD}U9EbTim$?^f=larJFenBLN zdgSvb67Fk}{hVpZXEkk(WY^NwWT<-I-eU z4XhE7cg37%E}I1Ju=yT7c#!Wbg2<-6zJ4w77QMyn$x3;H#f!JHIo#LhFbVIR*JW&T zpM9bm6}0I$6HLuCImgb%hV30(ogb3g-l_K7ZcL3ZzDPoTI55F$tD%{BW2NgxjHTI< z)gBp=D4>3J?6Bm-qw?`Cy&A8oV%s16)=Du}s|9K^*~M+X?1wMrtZx6#${YRkj2bJj zaP6oBi_tshFdY-mulMF*k!vY)Es5v86(n1(x3vwFJMfA+w_tiBkD(db=LxuWSeBVLJJG`G<#%N z|ESbXzMkmKQXOMFt9zokC2zu;eYP7hf~sF$&@!+8a7hpm^xpHxG$_sG)+@Zu=&^8> zKS^hM!+WJe-hOA}!SX=+>uW92H%?L!9qSPdiPuP4Nbi}Kur+393EU4uPV%ly)|{_a z&rI@uz32M1J-sWX#G*rnO(oLw)z@4ju{Hk?#JhMvbeFk4gsmQJrmsm1*`%8Yw-Uc-2-V#fNl@W^1R*lD5_KIFQUU4O)b zlamv6Q2y6}YP%CpkzFSwht@2?W*1ci~iN@egPglRaxIH9_Aj4XszjPjd%d6S5 zy4Bg4AJ(ghgAw&y&UTFqH^_*vNDBEB4}PyDoic8@V+VA&1xJ zulYq%vC5kP1IRy&&Ot37sWgHhO~()j4(@9R>A!tJ0jYGe;=rw82c9$NffU-=x*a>> zo?FFv{&&dQCT2DVg3d0(ZHAG!GF_c)VULaczN3_sSC$7~ui*9scBUq>6H2}2%Kn`F zWU|n_k(&6Y1l^u^f)I}!+^6=u{zy8m@3&j|)YMcp&nui|5iiD@j2XS;Wh!>klN%=< zqeJ`WLH)O*{ja`wfg>LN{WU8U*PDBSj?St_->}e3@pT@+K@KS?se-F(nNAagz53}g zHY6=AO)2YEjv2(AY{eX4M`sO8=V#D-mYf%?aW}6!Izh3C#I8N*uq~h7=mnx!-N`st?!;WKyKQ;$?aG%msV>= z*jdKh=i$08x_Ts$K|!~pJ?jF_@d8LWz0UK0`4UBXbGzHOZyR4?hogjaX^|W8;+)gH zc@aN9K8fj<`T;S_B(=AZjNe)9HJQNwfW*csA|m4cz0@vY#GC?##mN;zfG4_W@YcG) z$$GpwifzMl^*!k-wL;7BWghuF<=~_4S!P<`_YnVA$G9AN{CQ zXMa2-p)lLq+0=Y3LuT2LB%3=u-JQt{8;^zoZgKx?FR8RmXVq$tLD710a&k*sMg~We zm`g6bxU2dJYHCki8FcF2V~Fza*1ZZ23k#cHrXxPI1aTZi1v_sxo1PTl)|4Mj*8h(Z zGRg%fBEr&~qrO~0ZkBe3$QWZec&T(QUbr~CW9^OHxmsa&juV*PAB1m0%8l@|+ zn$8lO4qySMMm>UCPjh(locB(FBOsksZ|)@^Dq(B0-#o^TOW%cL0&LyswJPbf)_E9mMq)Tb8e34tOm!urEpNFH4b1pA;qa zBdr+sK93M&H)dD&&9KwA^h?`Z15U{2y|*)i9fHk?9}ClA zf#en(9NfI96m~Ic&lAqK9a}{pj+<2_3k!@0AYzy0ZW_NPZVp;3jAn=5- z;rEY&BO@bsce3@0Mt3%sm%iN%h-&%hAPKWzfR{j=r+aeh9kONz%C&OKem%cn)LU>@ zzM-+Paeed9wI6S8W<6ioLT~m{;heOM@mKuX!6IKDA7Kp^_z(x@TZLF@L_6#+Mbt!= z-~o%()nW)vFMlMt?ymM$-8v0N=}J-^hMT>lPatQLkod7eiGp!AL~d<~%byc)C`wpH zJKL3yo$jj1doc2)Gh2_J8aa!xV!7sbIOv?$#`|X|HtDua6Wxat%eG)?B5h=4(-6#x zeKgJ>=B(i1QTgJmP-}CY3nbe5IF~`oGy+^(TEv#PC%OtUc~ipPO1gN0U2ksOVR~Wp zU%!5-C>R!DUkX@ue)|3y)-XT3mn1Q_XnJZp9pgc_vUJe(>eZ`!j+5$kG=LekZm(Vv z=js+DSr1nBWqyZCCvvNhB!`IC0E(y$o`Ats$CArYxC6bF59c-EvFR@kY7FBt5ikX; zMR*PppaKY+&VAP^-lI_~vbnh_qIVL$z=_ZyI=O~ZqqFd4I7^an<;s=dz`)^8Iq+Wn zH1Y7>?sh5wI9hnn3R*{cEPkd54i6vylmkN;_D3>Eq#Som;;8dsSMc5WbRhF)R`|+H zpVoL=0xx8Hhom!bo%{B5kq1v*RCscP!K^byvCwndSr4M@`t|GeXYDG_zkdDN9>{)E zD5usaY&*46P>I$my=oPxw@Y9M z=g*(#a~M}a>#dea(;u%G#IUc7$}TEsWEq97w{ zxT#~r3Jm^|1}zM{kwl&(@rhqm_$7nroRs%%%oxl~{Sqryerz*G_on;WOwhaeY6C!v zTYbfrAy7BSkRN5_ZVJD|6}AMYhu+MM!j%4D4oVIyV-Y6t5?f_@McAblp`rRF-WqOL z^R+vQ@X*}ZMBImfb+hvE@%3Jz@^TZGnrMyVXO42Ph3kyHEFrqg$*K6eCnsh)at@av z1`ihBJ4x^NH_~3bgvKcq;n}likZYYf9!R9Mg*l<6-RS7(j1A0FH4khh>ApFV49~eE zjN%+lzq=cp4CRTZ9iZ@h6_yGl6zbyFhINiLKtMDzHB-g&Iyyd7!01Bnz3M@%R1?%X z>|X<4P`rOXKSVYN%~Q`$aoIDV6<6~E7^jd!3i)MG(f*z81tB23{O1y%^larAl&XFB z_)*#MC{@rq=h;R~i@(~uW)W7^@>~9E%W8@H!-VDpsUdMGK#3?x&x#jx7$x!}^c*JM zyCohA?wv}NXTK!=!|G` zQc^Q%TmMRIQNn9M5!^TT-`aZG9o3 z!2(7RK~@XJg)$Ied?l2uhHY1v;G|k#qQGZ~a5ayBKQG)P7F3C?M7n{Ik!R`}A&5=6 zV1yvTILxM}r#Hn*xU(=ZN#|BBDgEJZ--(9-A$Uc@A)wK+)^UFD$3H@K!V`h$vI>8C zR?Im&&mX-MSptQ0SZu5{kZYtYn{wl;j!8>NLRc6p#Zksk@7-yW&QxU@_CJSSg1bVb z#NujXd=F~AbGCW*{CRLkvb@;e!0+^l7cOHsNOAefl>&u{bHev)apv$K1QMc!t?9go z{bldHfSrYpM>pe>Nu;rD@s}>?7Y7DNZq!kM112qXA9yYb? zN_z!`7*Mzgs5aEd7PhZ!@)Z{XqHTq4k$cXZv z6xDLQjfzen5U5K7jfuy3OG`_&0Xq1$`T)UUyg4W(g#`};DCJ)#sY|9ttI5}%9Koi{ zLx7t9sBk)Q{P=mh%0)lgQy%|FP>~d9?RSZ2K?Vww5RGKHhW-+qIbZ}Y zyM}fl8OuBPXqNzbkqanCjar_UsbE}4y!v=F^Zk~w_pZ~_VAZ>9-2xd@ z+*mfWG$_0O`ntS)Q`v(NGpIp2-`wPy8YnNKm-G;BTc+XGWonb$&KB99Nbn9X9dv7| z{4%t+R;tMhdK2N)GE%+uJ1gz&JdlUZ{M*<1itYtuPIZ|N*#aQ~N_jtX@CJJfBcIt* zG%?N8eg52VFToQ_o@mkjuDx`~i%V$GHF;(Bd+!XcrJ+uOu)xy2?-{dgUacl9e8qPW@ZN}(GE|e+bq071NR@7f9X=G1nqEhXSq%9 zC2WVxT4rlPR_48$3Up}!H*#b{P$7LTerevrM!*(jxgw(QZJM{Of-)N1qW+BpNCzW7~6~RfUfcA(W64#uDDu)z+DAhDFFCKvu#$2D0yXm8&3wFr2ME=$b=*`QBtxrdpgRmP=JM zv$#hX zVhBFcf39{38JCao6|x(&mYcw0vEy!YgYNT}dT`|AtS1w~a?uznwdJPEx!L zSIoBbAu#>Cre(-|<6^lktt{Y5>c&9gK!3dm317%_o+pLmTJbrP2bL|dYHT@Ot;M+B z*5<`@j9a|Md$%DNS`60kD-=@v_GJKK*w$Jo?4$IrO3YC{f>pB}wLm;}7VF0Lw&x>& zxBHV@hicCcR<)RQZS_~o-3RcajE&!BVPhMAvv`&t`+_oVODiIwOQVr)Z*MWcJJhNi zuF^oZmH(Zn^U?sBTts`q%%0hN~{&BqsSP>U)nB#uJIe{@t)^MLU*L6hqAV3e7vAl=LTmIQgj=4!y5|3>2FC>*252{HVb6 zI*3sCz^ApWME60sR2zh3hE)G1$!IEKWnmG1zy)tZeUUQneXl7f7X>WZr9lZ~Lp3$p z<7i@wAg>j17+28L430Wb@vmWt?o8Zv9tuk*7ybh*0nRo8C~#|~sRWTw91*_B0$1kR zw3Zu(lCvR0BiRf%F%(}G;|IDik4jJ(9^s)Si%k%9`dtrT)QtLBG-S0PVT%V$(-T0T zKykQyeb5~(TGdVi)vUevioEvaS)oQ4fAe!D)fhoco%77M^Lu-H_8mAcBzrn)0J`}nJob0j?N`6*i0N^|V`d2hgJAdU+4DRq zxseFyu;+c~*fG;_!~=w#mTDLhAbi~4Ly9o8f7To?@!hBD$2*spyE=vsd`0GtGcwMZ z1;NF#M^@Q}Q7Qltn|qgTL;tO{<=fbwH+CW2wd67YZb^$ACIEd4hoZcyeWKKEXpox< z(iOqJpEMS-HRw7Rt;&<+%xIC=h)24g#lbPgDRi~47w{!2a1X&{H^dMysnib5$c3Uu znczc@fd6d8J>W*tfLDMq!XxG~x6&_xZH_z-D;5sf_}SB^m!}WXprSrncQ=?v7_mWF ziKY70Cbd6+8-wY~x8f`VE*gFFYNh$l556ZD7%oA5z_tS`tt#p=XPuS^c&iq_Oyvx$ zS9~3GV5QxyRdav@R5Qs7I}R&e!8Tm8Gse5rV_=6LD}jZ+h?dPnB;*Rq<1MPkx*O6b znxjG{Hr;l<72GW}YdBH9_~oo>%g;-t0~bJzSyI|WCpC*myW&e8c^=BVKu$}sA8VxX zE&J8pU+LpRN-N>!xK@Bxk!EHoKny1EE$InfR%~%SY{NGsH?N}-kmvSV*u8kX+|KT< z|4o!$F4cA*0CKm==PIU~N1#q0b!%;F%lnq+zA~xNP+uQd;WRxin!Ek|txeN%G^}(} zuZj1iqY^8DQ|qO;ULHL71*#^WK@b?jxBQ((X82&*bdRHTGpD&M7MnJ;FNEvi)*4C2 z1XCQL4+PC;;w*Ca?pMg%*6~4z zHbuY!OU+h)!7J?S?Dy=p*JioKTmX}X3+bdw?tDI21l0Fi_lL32A;L?k6KBp?SS_qL zy>ptzzM8^K?=AfNxLhem;l-_+UmKggOD56MNq2a4bty^f^CU#Mdc24}Nk&V$%fN#t zl-Pc|c_7T?Z-d)8O;J_{ixiUJy{CslRAy!-Ry+UqP>px))eBpduZU5BOLbv;rl@Vy zm&W+=V;NYgf#wc6>9}a~xqwT=>B?bZTNRpvjM=y#QY3~+E5`-)t8*1+W>m5}S~YDu|J4cib1c}9*Z=?^#O z-pyCuVuafjR8hJv48KnsVt_Z?nMC65YuQoGvh0s~usTMJxdLpD;1Q*1sR$LP^VeFI zRSp{IQ(uKON?gc}I{-AW_Ev5r`ihumkgHEjaB*)TFtKG!i0wdGf%-(pA5OHLG`wpP9;CPre9C@r2XEP5oY^;2aQDdu_>8bYD~(VEMfeAWI&K4Xa?xp| z!&pV$1?n4tZhp`^{yt^_yt&|S-@eT(q^~SOIJSW8kKFH5VggOuezZQw|J8H>;Gn49 z$ z{_@Gn&?C*#vuuWM61hj3dH?dKC=FjP&R>`~o z$f(xQm00N5I+M%E}e(6Rflv?-3>-VyNrN;!-q7`|qIne^J>O%Z7 zRs^gU5NZG7Tuv2eq1p350O2vJ(B4{~=amm&teO;w`@qcn_$IpZ48Kw=EY_Wi0IqHV zGJu8EzP+svaC3x)ycY7aqN8Im(9zz8@9151i1LBM6k2_~$ZiK54KZ0GiiHv#%&jTo z+^RMGOx=Pzxf2PlzbM)NRI#iem$WXf!-Ok=MvzI_T<)AD$@XWTawB`Y^WG7QBXvzh zZIZZK&SvIJdAgQGzznj&lQgZ|#P7P-SGj0`!NQa!VQRumLebY{d%FcVau zDG1pPMB=>vo49@M835Oa?yo!)5!tnC^ZhM=yU>savt89{I4E?P?du{gbFqEmnwpx= zjj%A6z2vA&*RNubs&p#qpLIIRw@!se2T_gRq{fTSea#x=&s%?-C``*|IB_g0*%R~) z9UbCDAF0#;cC}az6oGeE^Sm$hOj|we(*Fb-XACkBm>cgVcTBuCKOf$kgw48E(A=!D z%+JJBdzc<2&%70CF_8{}mN+U?U0vORMhHDqS){wx@h|AOk`u<2TU#sLf%VhR?? zIk9UOJTDW+Jxs6RT{)Hs0c`UbEkD#SA5y0|dD09N1j17D-@t--7{9ssPq31FJue62 zAll&wxv>Q|EpQnULtpZvZHSB2>0Y6hvWn9W{7+x9H!=Yhu)O?X1RmFR!~*Hg)N^dD z0&P*r&24M{G(5nI4__%`m1SVAqd$Qp!hhTmY6xzl3ITUhg*$gHD7V9Qo}A$V(9=^J z(#DM!pNsFTc81*&zN8RSoOLsct)Zbov&8DA2~GNdXq>ZgbR9Q;fOmAeo(wU=?z>M9 z-d33)95~I5-C1j@#3$sEMyTYAv?i+pJpWhz=0E(~W+|eL_Nn3*cZ_s0lWI0LH^bcU z2;vM5uw=LMucco`y%~Pobp7$r;=pk_nn3i#ayE{sw{TmnQWRrtT5e@cg}mX;r@F!JRktFW-{g$ox( z0mjh$I?E#GzW!)KUKPrj^7n39q*2*Q07zX#$k7YIoiK0p;9khHv$u-=LOA>5JLQkth2L|<^1`pasiBSp4pHT()4t2>jy_xrn-XT<9VPQY+jEA z!8xjhnUyuwht(wl?Vy5gj?rng^i#W#@J;cEIEYGSW@ZXd5;!pQqS1cpHzpoGPK)}8 zudhf!anI+qvvDG#H|J-bj$tLsjKtg@AKD3@Vm`85Oh=16ewGhIz@$ z2PD@A99~60)~V&pr6gEwEd0ba+5MV`_1az4d7=bDtNPorq(Ak_!T+LO`9IN9h}z5Q zZ_Uij8=_WTTrhDF&6(_l#~flpF6eAbpQ7QVbPw3B;D9Iaw;+}i%-Q6ECY>N51TLH8o6PogRjI1E*jJHY;xt*qYcY!B_XO0G>v z?dySUq572jMQG^BOepGE*x8lAZ&9ADid7TiZMhpBw&COz-5zDK#G@T786U6z@z)-R zIcm#tsn?%JczH*#=PFPDyZJla--1SpXD?oaZftCTwKAfAH!`H>43AzR%h2|p&Jd9r zY4cMpWpk*7gdxbT1>pqd1iayI05F%P(dL1?}Wyq4~uSK(;6CLQqU({+Bn zZzFzZZ7)0Qzog}t0}Z$)>RF-_X@1@W0`RpU*xO7pd8nD06p`Kx`7m$ivOME;lEo!^$64eLz|`S-IJuke zX&TQN(Bo-&&{7$^1Ng) zA2Dw|6)?=tNB@g%Adp{|Ovg}3YIOXm%i-<@c#yjpq)kbwvhCc~Ufv23{lB4AQ^0P} z5STwF;K`Cx5Wyo(_<>CSTN|#QKYw;6$sUT5aDSWUPwGn6RU3K@wQ+H0T4RK6NbW3L zGohj2Xn+5hVhs8u$7!rUnkcsHl0$Vr$oOkc9Rz>Ksr8_Uu&}UP*VNKdJzxk;rG{YX zRP)fOzsVrteo`!$3R+e*54k!A+GkAZ{RSitO~{QR&NI-_Og8$?E89;h^{V=`LI(T8 zv&RvJw{ELeQ~f1bAZtEPNJxOr&Cw)>|5?lmoUhp*WH%>M+7^om3^cXUMO$S!M69*x zhIDMmZK0Y$rPg?_Ez4#3L&Jik(oZZy9_@Wh+DFFzA#T}(bMyo1iIf6yu-{m?UDtFl zbll|@ZWpcizzDF<7|qwP^e};3$AdPCL8&Nznp@Yj;xdv8Wd7&btwmnrO5xrD=h=R9 ziM4Nc`?kRZ)&oae52!?Wt3kQF`7Q7e8MwIAgXiJwdqtv-lYtm$N0o0n^e$|4NH+j% zDsyns>Zi@gV&$kzX`!J48(k0eOi{qi4E6Ujb#-;y#n#dfYKpepq4n%s_;=fPsvz7? zeC+(-?7LZQWs`8EPzzXbnqoOvMv`3o!lqg@}VPm4M4K=%Yz|1l^l8^r_pG(&o!>>ef$ zW$Zl$3dS6_PVau>9Y{(_Ds02gmh=|VL9@y!xN92s2Ixp0`q*zT=eZlHtMjCP=AQ;1 zn_xak0xcN-tG{ZfO$S&qgT5bM(RIkAA9jwfAB}T~7CRQsUYGN0(C^35}r1drV zG3mF^snG~ZV9Tu>eGJ-(171{*k8r>rUl<%RXiRucc>z^PKR!8nPYRk9E>de&Z~Quh znzp0>Kw|Ti%+L}a0iNZRcZwP@{L5CvgEihYz_{hn^?=|PT-PtPHL5oE0=YaD8Zw|o zpS#!4;bn8INRj%yMvVIEuG8E=A3b;sZw%dMOUC1&9dl_iIiL~52YaI|b6eCJ0YOa#O6mKb)J(_jD%-cG>c+WLGR9-TC@c9(m09vPq!J0KL zjWx9_(_Q9Gc{m+zU>gLX*o2DD#N(itvhefsw-r`HR~@u&GIyolh=3}(35;D;5zHr9 zS=mz?j88yUg^k@JiP9K6SQ2APO}r7Pv499$b!W^-XP(!ViU4j#{cEKP+Q`(DWTb}J z<6S@a{61@2b`@5Y4cC)>+^%EE%0zN^UgZ9;H1!~+z6Nvwwr)x29UNcw_|xr1j+m1W z9h2tK-(H0~#LzKQ;h~N7%62@G@&5g_k8yr%IjWrt;YFJcXtM+?$c(2a4-->j?RXQ} zq0^Cc964c~bkg9|>YO@b0R}zmR7diYzTSnkt{YMf(5j_jigeSJf-pBG8eWu^3%WSp zgJ0KaMI9}ik?ugpi8z;N^$heU9h=bmp?Zuv7YGOSt@-i^M(R(XjG8^l7TR9HTZt!l zEGcSrz41H^=DM7N3Q2By&oRQg4${L+5mNh`ZEQREVsJy?hXYYR0(hYUgq|ZSs~~}5 zi(=k1h*5iU*C9KJa_{qAhEve2ND-ns%IRV+A=5-FlOCK zF>FjKjrgX$R>lQ7vb^y-m6K^u>H1*J^m%MCHY=mpPX9K3YoY~AB2xLbx8O3Mg7v}5 zwGV}s=R|URy9i(&oQnz+yZhyZ_>hi(Wvpv3d&Gm za3gM5_4D&nMd#dQ_I&G>aJCe~JTQWCL*ZeWkD-E5zZKe{_KVUl@ZMWaz+QfMq=Lix zWzehQOBMq8ok-ewwe1hP$8&NHLSehRS{x-jk*4}y_VvRXj{^kRaa0+C^Lwx=<*gkL z+3=o%WSK_q-bTrNzBW`t-W^-4C(g&$DWGTv-Gx3KPWJ8cxmY2;Sryg$?Wf>?zeWW`*=>j$s+hxg3NCIGS`xOn(>1XM;f2{ z`!_kgI?~g==Ovmi8}rF&xZ&LD{=3M=$b0l-mzlvQ(0Y1aBtQqN_p$?mqX8{S(XC&RqnqUPk3dnUEew6~bem+Jg4$B|Hz z9$yA+L|1ttrwc!$yOw3$c|R*TjCY1P!*jiuH?1eu?QzfGZ_hYguOkG0V2cEMlwrPb z_8O16SLu~UW9<5VVEV)HL8MV?9Q{@%ro*mcDs5)`u3Z=Y82knzwPg(%c+6%B6>?KM jDJm`g=imB0HQM7nt>WjC+xvtcjzbh>)vo8tm_GU+75c;I