diff --git a/.travis.yml b/.travis.yml index 96e5ac93dd..302f08ad8b 100644 --- a/.travis.yml +++ b/.travis.yml @@ -4,9 +4,17 @@ python: "3.6" install: - pip install -U -r requirements.txt - pip install -U -r requirements-dev.txt + - pip install -U -r requirements-examples.txt script: - python -m pytest --cov-report term --cov=SDM + - | + python -m ipykernel install --user + for i in examples/*.ipynb; do + travis_wait 30 jupyter nbconvert --ExecutePreprocessor.timeout=1800 --to markdown --execute --stdout $i > $i.md.travis; + # TODO! + #diff $i.md.repo $i.md.travis + done; after_success: - codecov diff --git a/SDM/backends/__init__.py b/SDM/backends/__init__.py new file mode 100644 index 0000000000..1f666cac2f --- /dev/null +++ b/SDM/backends/__init__.py @@ -0,0 +1,6 @@ +""" +Created at 24.07.2019 + +@author: Piotr Bartman +@author: Sylwester Arabas +""" \ No newline at end of file diff --git a/SDM/backends/numpy.py b/SDM/backends/numpy.py new file mode 100644 index 0000000000..56a6dcd553 --- /dev/null +++ b/SDM/backends/numpy.py @@ -0,0 +1,137 @@ +""" +Created at 24.07.2019 + +@author: Piotr Bartman +@author: Sylwester Arabas +""" + +import numpy as np + + +# TODO backend.array overrides __getitem__ + +class Numpy: + storage = np.ndarray + + @staticmethod + def array(shape, type): + if type is float: + data = np.full(shape, np.nan, dtype=np.float) + elif type is int: + data = np.full(shape, -1, dtype=np.int) + else: + raise NotImplementedError + return data + + @staticmethod + def from_ndarray(array): + if array.ndim > 1: + result = array.copy() + else: + result = np.reshape(array.copy(), (1, -1)) + return result + + @staticmethod + def shuffle(data, length, axis): + idx = np.random.permutation(length) + Numpy.reindex(data, idx, length, axis=axis) + + @staticmethod + def reindex(data, idx, length, axis): + if axis == 1: + data[:, 0:length] = data[:, idx] + else: + raise NotImplementedError + + @staticmethod + def argsort(idx, data, length): + idx[0:length] = data[0:length].argsort() + + @staticmethod + def stable_argsort(idx: np.ndarray, data: np.ndarray, length: int): + idx[0:length] = data[0:length].argsort(kind='stable') + + @staticmethod + def amin(data): + result = np.amin(data) + return result + + @staticmethod + def amax(data): + result = np.amax(data) + return result + + @staticmethod + def transform(data, func, length): + data[:length] = np.fromfunction( + np.vectorize(func, otypes=(data.dtype,)), + (length,), + dtype=np.int + ) + + @staticmethod + def foreach(data, func): + for i in range(len(data)): + func(i) + + @staticmethod + def shape(data): + return data.shape + + @staticmethod + def urand(data, min=0, max=1): + data[:] = np.random.uniform(min, max, data.shape) + + # TODO do not create array + @staticmethod + def remove_zeros(data, idx, length) -> int: + for i in range(length): + if data[0][idx[0][i]] == 0: + idx[0][i] = idx.shape[1] + idx.sort() + return np.count_nonzero(data) + + @staticmethod + def extensive_attr_coalescence(n, idx, length, data, gamma): + # TODO in segments + for i in range(length // 2): + j = 2 * i + k = j + 1 + + j = idx[j] + k = idx[k] + + if n[j] < n[k]: + j, k = k, j + g = min(gamma[i], n[j] // n[k]) + + new_n = n[j] - g * n[k] + if new_n > 0: + data[:, k] += g * data[:, j] + else: # new_n == 0 + data[:, j] = g * data[:, j] + data[:, k] + data[:, k] = data[:, j] + + @staticmethod + def n_coalescence(n, idx, length, gamma): + # TODO in segments + for i in range(length // 2): + j = 2 * i + k = j + 1 + + j = idx[j] + k = idx[k] + + if n[j] < n[k]: + j, k = k, j + g = min(gamma[i], n[j] // n[k]) + + new_n = n[j] - g * n[k] + if new_n > 0: + n[j] = new_n + else: # new_n == 0 + n[j] = n[k] // 2 + n[k] = n[k] - n[j] + + + diff --git a/SDM/colliders.py b/SDM/colliders.py index 2cac5b744c..3d16cb4688 100644 --- a/SDM/colliders.py +++ b/SDM/colliders.py @@ -5,41 +5,52 @@ @author: Sylwester Arabas """ - -import numpy as np -# import numba +from SDM.backends.numpy import Numpy as backend class SDM: - def __init__(self, kernel, dt, dv): - M = 0 # TODO dependency to state[] - N = 1 - self.probability = lambda sd1, sd2, n_sd: \ - max(sd1[N], sd2[N]) * kernel(sd1[M], sd2[M]) * dt / dv * n_sd * (n_sd - 1) / 2 / (n_sd//2) + def __init__(self, kernel, dt, dv, n_sd): + self.probability = lambda sd1n, sd2n, sd1x, sd2x, n_sd: \ + max(sd1n, sd2n) * kernel(sd1x, sd2x) * dt / dv * n_sd * (n_sd - 1) / 2 / (n_sd // 2) + self.rand = backend.array((n_sd // 2,), type=float) + self.prob = backend.array((n_sd // 2,), float) + self.gamma = backend.array((n_sd // 2,), float) - # @numba.jit() #TODO def __call__(self, state): - n_sd = len(state) - - assert np.amin(state.n) > 0 - if n_sd < 2: - return + assert state.is_healthy() # toss pairs state.unsort() + # TODO (segments) + # state.sort_by('z', stable=True) # state.stable_sort_by_segment() + # collide iterating over pairs - rand = np.random.uniform(0, 1, n_sd // 2) + backend.urand(self.rand) + + backend.transform(self.prob, lambda j: self.probability(state._n[state._idx[2 * j]], + state._n[state._idx[2 * j + 1]], + state._x[state._idx[2 * j]], + state._x[state._idx[2 * j + 1]], + state.SD_num), + state.SD_num // 2) + + backend.transform(self.gamma, lambda j: self.prob[j] // 1 + (self.rand[j] < self.prob[j] - self.prob[j] // 1), + state.SD_num // 2) + + # TODO (potential optimisation... some doubts...) + # state.sort_by_pairs('n') - prob_func = np.vectorize(lambda j: self.probability(state[2*int(j)], state[2*int(j)+1], n_sd)) - prob = np.fromfunction(prob_func, rand.shape, dtype=float) + # TODO (when an example with intensive param will be available) + # backend.intesive_attr_coalescence(data=state.get_intensive(), gamma=self.gamma) - # prob = np.empty_like(rand) - # for i in range(1, n_sd, 2): - # prob[i // 2] = self.probability(state[idx[i]], state[idx[i - 1]], n_sd) + for attrs in state.get_extensive_attrs().values(): + backend.extensive_attr_coalescence(n=state._n, + idx=state._idx, + length=state.SD_num, + data=attrs, + gamma=self.gamma) - gamma = np.floor(prob) + np.where(rand < prob - np.floor(prob), 1, 0) + backend.n_coalescence(n=state._n, idx=state._idx, length=state.SD_num, gamma=self.gamma) - # TODO ! no loops - for i in range(1, n_sd, 2): - state.collide(i, i - 1, gamma[i//2]) + state.housekeeping() diff --git a/SDM/runner.py b/SDM/runner.py index f64f693b44..8ec5e15c1c 100644 --- a/SDM/runner.py +++ b/SDM/runner.py @@ -7,6 +7,7 @@ from SDM.stats import Stats + class Runner: def __init__(self, state, dynamics): self.state = state diff --git a/SDM/state.py b/SDM/state.py index c37435348b..6901aa560c 100644 --- a/SDM/state.py +++ b/SDM/state.py @@ -6,88 +6,124 @@ """ import numpy as np +from SDM.backends.numpy import Numpy as backend class State: - def __init__(self, x, n): - assert x.shape == n.shape - assert len(x.shape) == 1 + def __init__(self, n: np.ndarray, intensive: dict, extensive: dict, segment_num: int): + assert n.ndim == 1 + + # https://en.wikipedia.org/wiki/Intensive_and_extensive_properties + for attribute in intensive.values(): + assert backend.shape(attribute) == backend.shape(n) + for attribute in extensive.values(): + assert backend.shape(attribute) == backend.shape(n) + + self.SD_num = len(n) + self.idx = backend.from_ndarray(np.arange(self.SD_num)) + self.n = backend.from_ndarray(n) + self.keys = {} + self.attributes = {'intensive': {}, 'extensive': {}} + # TODO clean + attributes = {'intensive': State.divide_by_type(intensive), 'extensive': State.divide_by_type(extensive)} + # self.attributes['intensive']['int'] = backend.array((len(attributes['intensive']['int']), self.SD_num), int) + # self.attributes['intensive']['float'] = backend.array((len(attributes['intensive']['float64']), self.SD_num), float) + # self.attributes['extensive']['int'] = backend.array((len(attributes['extensive']['int']), self.SD_num), int) + self.attributes['extensive']['float64'] = backend.array((len(attributes['extensive']['float64']), self.SD_num), float) + + for tensive in self.attributes: + for dtype in self.attributes[tensive]: + idx = 0 + for key, array in attributes[tensive][dtype].items(): + self.keys[key] = (tensive, dtype, idx) + self.attributes[tensive][dtype][idx, :] = array[:] + idx += 1 + + self.segment = backend.from_ndarray(np.full(segment_num, 0)) + self.segment_multiplicity = backend.array((segment_num,), int) + self.segment_order = backend.array((self.SD_num,), int) + + @staticmethod + def divide_by_type(attributes: dict) -> dict: + result = {} + for key, ndarray in attributes.items(): + dtype = str(ndarray.dtype) + if dtype not in result: + result[dtype] = {} + result[dtype][key] = ndarray + return result - self.data = np.concatenate((x.copy()[:, None], n.copy()[:, None]), axis=1) + # TODO: in principle, should not be needed at all (GPU-resident state) + def __getitem__(self, item: str) -> backend.storage: + all_valid = self.idx[0, :self.SD_num] + if item == 'n': + result = self.n[0, all_valid] + else: + tensive = self.keys[item][0] + dtype = self.keys[item][1] + attr = self.keys[item][2] + result = self.attributes[tensive][dtype][attr, all_valid] + return result @property - def x(self): - return self.data[:, 0] - - @x.setter - def x(self, value): - self.data[:, 0] = value + def _n(self): + return self.n[0] @property - def n(self): - return self.data[:, 1] + def _idx(self): + return self.idx[0] - @n.setter - def n(self, value): - self.data[:, 1] = value - - # TODO optimize - def __sort(self, key): - idx = np.argsort(key) - self.n = self.n[idx] - self.x = self.x[idx] - - def sort_by_m(self): - self.__sort(self.x) + @property + def _x(self): + return self.attributes['extensive']['float64'][0, :] - def sort_by_n(self): - self.__sort(self.n) + def sort_by(self, item: str, stable=False): + if stable: + backend.stable_argsort(self.idx, self[item], length=self.SD_num) + else: + backend.argsort(self.idx, self[item], length=self.SD_num) - # TODO optimize def unsort(self): - idx = np.random.permutation(np.arange(len(self))) - self.x = self.x[idx] - self.n = self.n[idx] + backend.shuffle(self.idx, length=self.SD_num, axis=1) - def x_min(self): - result = np.amin(self.x) + def min(self, item): + result = backend.amin(self[item]) return result - def x_max(self): - result = np.amax(self.x) + def max(self, item): + result = backend.amax(self[item]) return result - def moment(self, k, m_range=(0, np.inf)): + # TODO update + def moment(self, k, attr='x', attr_range=(0, np.inf)): idx = np.where( np.logical_and( - self.n > 0, # TODO: alternatively depend on undertaker... - np.logical_and(m_range[0] <= self.x, self.x < m_range[1]) + self['n'] > 0, # TODO: alternatively depend on undertaker... + np.logical_and(attr_range[0] <= self[attr], self[attr] < attr_range[1]) ) ) if not idx[0].any(): return 0 if k == 0 else np.nan - avg, sum = np.average(self.x[idx] ** k, weights=self.n[idx], returned=True) + avg, sum = np.average(self[attr][idx] ** k, weights=self['n'][idx], returned=True) return avg * sum - def collide(self, j, k, gamma): - if self.n[j] < self.n[k]: - j, k = k, j + def get_extensive_attrs(self): + result = self.attributes['extensive'] + return result + + def is_healthy(self): + result = backend.amin(self.n[0][self.idx[0, 0:self.SD_num]]) > 0 + return result + + # TODO: optionally recycle n=0 drops + def housekeeping(self): + if self.is_healthy(): + return + else: + self.SD_num = backend.remove_zeros(self.n, self.idx, length=self.SD_num) + + - gamma = min(gamma, self.n[j] // self.n[k]) - if self.n[k] != 0: #TODO: guaranteed by undertaker - n = self.n[j] - gamma * self.n[k] - if n > 0: - self.n[j] = n - self.x[k] += gamma * self.x[j] - else: # n == 0 - self.n[j] = self.n[k] // 2 - self.n[k] = self.n[k] - self.n[j] - self.x[j] = gamma * self.x[j] + self.x[k] - self.x[k] = self.x[j] - def __len__(self): - return self.x.shape[0] - def __getitem__(self, item): - return self.x[item], self.n[item] diff --git a/SDM/undertakers.py b/SDM/undertakers.py deleted file mode 100644 index c5f4a2a445..0000000000 --- a/SDM/undertakers.py +++ /dev/null @@ -1,20 +0,0 @@ -""" -Created at 07.06.2019 - -@author: Piotr Bartman -@author: Sylwester Arabas -""" - - -class Resize: - def __call__(self, state): - # TODO dependency state items - idx_valid = state.n != 0 - state.data = state.data[idx_valid] - - -class Recycle: - def __call__(self, state): - #TODO: state.sort_by_n() - - raise NotImplementedError diff --git a/examples/SDM b/examples/SDM new file mode 120000 index 0000000000..e6f1e03724 --- /dev/null +++ b/examples/SDM @@ -0,0 +1 @@ +../SDM/ \ No newline at end of file diff --git a/examples/Shima_et_al_2009_Fig2/Fig_2.py b/examples/Shima_et_al_2009_Fig2/Fig_2.py deleted file mode 100644 index ac8cbb44a8..0000000000 --- a/examples/Shima_et_al_2009_Fig2/Fig_2.py +++ /dev/null @@ -1,65 +0,0 @@ -""" -Created at 07.06.2019 - -@author: Piotr Bartman -@author: Sylwester Arabas -""" - -import copy -import numpy as np - -from examples.Shima_et_al_2009_Fig2.plotter import Plotter -from examples.Shima_et_al_2009_Fig2.setups import SetupA - -from SDM.runner import Runner -from SDM.state import State -from SDM.colliders import SDM -from SDM.undertakers import Resize -from SDM.discretisations import constant_multiplicity - - -def test_Fig2(): - with np.errstate(all='raise'): - setup = SetupA() - states, _ = run(setup) - - x_min = min([state.x_min() for state in states.values()]) - x_max = max([state.x_max() for state in states.values()]) - - with np.errstate(invalid='ignore'): - plotter = Plotter(setup, (x_min, x_max)) - for step, state in states.items(): - plotter.plot(state, step * setup.dt) - plotter.show() - - -# TODO python -O -def test_timing(): - setup = SetupA() - setup.steps = [100, 3600] - - nsds = [2 ** n for n in range(12, 15)] - times = [] - for sd in nsds: - setup.n_sd = sd - _, stats = run(setup) - times.append(stats.times[-1]) - - from matplotlib import pyplot as plt - plt.plot(nsds, times) - plt.show() - - -def run(setup): - state = State(*constant_multiplicity(setup.n_sd, setup.spectrum, (setup.x_min, setup.x_max))) - collider = SDM(setup.kernel, setup.dt, setup.dv) - undertaker = Resize() - runner = Runner(state, (undertaker, collider)) - - states = {} - for step in setup.steps: - runner.run(step - runner.n_steps) - setup.check(runner.state, runner.n_steps) - states[runner.n_steps] = copy.deepcopy(runner.state) - - return states, runner.stats diff --git a/examples/Shima_et_al_2009_Fig2/input.json b/examples/Shima_et_al_2009_Fig2/input.json deleted file mode 100644 index f42b9f324a..0000000000 --- a/examples/Shima_et_al_2009_Fig2/input.json +++ /dev/null @@ -1,15 +0,0 @@ -{ - "m_min": r2x(10e-6), // not given in the paper - "m_max": r2x(100e-6), // not given in the paper - - "n_sd": 2 ** 13, - "n_part": 2 ** 23, // [m-3] - "X0": 4/3 * np.pi * 30.531e-6**3, - "dt": 1, // [s] - "dv": 1e6, // [m3] - "b": 1.5e3, // [s-1] - "rho": 1000, // [kg m-3] - - "check_LWC": 1e-3, // kg m-3 #TODO - "check_ksi": n_part * dv / n_sd // TODO -} diff --git a/examples/Shima_et_al_2009_Fig2/plotter.py b/examples/Shima_et_al_2009_Fig2/plotter.py deleted file mode 100644 index cf358bc357..0000000000 --- a/examples/Shima_et_al_2009_Fig2/plotter.py +++ /dev/null @@ -1,74 +0,0 @@ -""" -Created at 09.07.2019 - -@author: Piotr Bartman -@author: Sylwester Arabas -""" - -from matplotlib import pyplot -from examples.Shima_et_al_2009_Fig2.utils import * - -import numpy as np - - -class Plotter: - def __init__(self, setup, xrange): - self.setup = setup - - self.x_bins = np.logspace( - (np.log10(xrange[0])), - (np.log10(xrange[1])), - num=64, - endpoint=True - ) - self.r_bins = x2r(self.x_bins) - - def show(self): - pyplot.show() - - def save(self, file): - pyplot.savefig(file) - - def plot(self, state, t): - s = self.setup - - if t == 0: - analytic_solution = s.spectrum.size_distribution - else: - analytic_solution = lambda x: s.norm_factor * s.kernel.analytic_solution( - x=x, t=t, x_0=s.X0, N_0=s.n_part - ) - - dm = np.diff(self.x_bins) - dr = np.diff(self.r_bins) - - pdf_m_x = self.x_bins[:-1] + dm/2 - pdf_m_y = analytic_solution(pdf_m_x) - - pdf_r_x = self.r_bins[:-1] + dr / 2 - pdf_r_y = pdf_m_y * dm / dr * pdf_r_x - - pyplot.plot( - m2um * pdf_r_x, - kg2g * pdf_r_y * r2x(pdf_r_x) * s.rho / s.dv, - color='black' - ) - - vals = np.empty(len(self.r_bins)-1) - for i in range(len(vals)): - vals[i] = state.moment(1, (self.x_bins[i], self.x_bins[i+1])) - vals[i] *= s.rho / s.dv - vals[i] /= (np.log(self.r_bins[i+1]) - np.log(self.r_bins[i])) - - pyplot.step( - m2um * self.r_bins[:-1], - kg2g * vals, - where='post', - label=f"t = {t}s" - ) - pyplot.grid() - pyplot.xscale('log') - pyplot.xlabel('particle radius [µm]') - pyplot.ylabel('dm/dlnr [g/m^3/(unit dr/r)]') - pyplot.legend() - diff --git a/examples/Shima_et_al_2009_Fig2/setups.py b/examples/Shima_et_al_2009_Fig2/setups.py deleted file mode 100644 index 5343717959..0000000000 --- a/examples/Shima_et_al_2009_Fig2/setups.py +++ /dev/null @@ -1,43 +0,0 @@ -""" -Created at 09.07.2019 - -@author: Piotr Bartman -@author: Sylwester Arabas -""" - -from SDM.spectra import Exponential -from SDM.kernels import Golovin -from examples.Shima_et_al_2009_Fig2.utils import * -import numpy as np - - -class SetupA: - x_min = r2x(10e-6) # not given in the paper - x_max = r2x(100e-6) # not given in the paper - - n_sd = 2 ** 13 - n_part = 2 ** 23 # [m-3] - X0 = 4/3 * np.pi * 30.531e-6**3 - dv = 1e6 # [m3] - norm_factor = n_part * dv - rho = 1000 # [kg m-3] - - dt = 1 # [s] - steps = [0, 1200, 2400, 3600] - - kernel = Golovin(b=1.5e3) # [s-1] - spectrum = Exponential(norm_factor=norm_factor, scale=X0) - - # TODO: rename? - def check(self, state, step): - check_LWC = 1e-3 # kg m-3 - check_ksi = self.n_part * self.dv / self.n_sd - - # multiplicities - if step == 0: - np.testing.assert_approx_equal(np.amin(state.n), np.amax(state.n), 1) - np.testing.assert_approx_equal(state.n[0], check_ksi, 1) - - # liquid water content - LWC = self.rho * np.dot(state.n, state.x) / self.dv - np.testing.assert_approx_equal(LWC, check_LWC, 3) diff --git a/examples/Shima_et_al_2009_Fig2/utils.py b/examples/Shima_et_al_2009_Fig2/utils.py deleted file mode 100644 index 6e7f0f2830..0000000000 --- a/examples/Shima_et_al_2009_Fig2/utils.py +++ /dev/null @@ -1,20 +0,0 @@ -""" -Created at 09.07.2019 - -@author: Piotr Bartman -@author: Sylwester Arabas -""" - -from numpy import pi - -kg2g = 1e3 -m2um = 1e6 - - -def x2r(x): - return (x * 3/4 / pi)**(1/3) - - -def r2x(r): - return 4/3 * pi * r**3 - diff --git a/examples/Shima_et_al_2009_Fig_2.ipynb b/examples/Shima_et_al_2009_Fig_2.ipynb new file mode 100644 index 0000000000..e5ea93990c --- /dev/null +++ b/examples/Shima_et_al_2009_Fig_2.ipynb @@ -0,0 +1,293 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 14, + "outputs": [], + "source": [ + "import copy\n", + "import numpy as np\n", + "from matplotlib import pyplot\n", + "\n", + "from SDM.runner import Runner\n", + "from SDM.state import State\n", + "from SDM.colliders import SDM\n", + "from SDM.discretisations import constant_multiplicity\n", + "from SDM.spectra import Exponential\n", + "from SDM.kernels import Golovin\n" + ], + "metadata": { + "collapsed": false, + "pycharm": { + "name": "#%%\n", + "is_executing": false + } + } + }, + { + "cell_type": "code", + "execution_count": 15, + "outputs": [], + "source": [ + "def x2r(x):\n", + " return (x * 3/4 / np.pi)**(1/3)\n", + "\n", + "def r2x(r):\n", + " return 4/3 * np.pi * r**3\n", + "\n", + "kg2g = 1e3\n", + "m2um = 1e6" + ], + "metadata": { + "collapsed": false, + "pycharm": { + "name": "#%%\n", + "is_executing": false + } + } + }, + { + "cell_type": "code", + "execution_count": 16, + "outputs": [], + "source": [ + "class Plotter:\n", + " def __init__(self, setup, xrange):\n", + " self.setup = setup\n", + "\n", + " self.x_bins = np.logspace(\n", + " (np.log10(xrange[0])),\n", + " (np.log10(xrange[1])),\n", + " num=64,\n", + " endpoint=True\n", + " )\n", + " self.r_bins = x2r(self.x_bins)\n", + "\n", + " def show(self):\n", + " pyplot.show()\n", + "\n", + " def save(self, file):\n", + " pyplot.savefig(file)\n", + "\n", + " def plot(self, state, t):\n", + " s = self.setup\n", + "\n", + " if t == 0:\n", + " analytic_solution = s.spectrum.size_distribution\n", + " else:\n", + " analytic_solution = lambda x: s.norm_factor * s.kernel.analytic_solution(\n", + " x=x, t=t, x_0=s.X0, N_0=s.n_part\n", + " )\n", + "\n", + " dm = np.diff(self.x_bins)\n", + " dr = np.diff(self.r_bins)\n", + "\n", + " pdf_m_x = self.x_bins[:-1] + dm / 2\n", + " pdf_m_y = analytic_solution(pdf_m_x)\n", + "\n", + " pdf_r_x = self.r_bins[:-1] + dr / 2\n", + " pdf_r_y = pdf_m_y * dm / dr * pdf_r_x\n", + "\n", + " pyplot.plot(\n", + " m2um * pdf_r_x,\n", + " kg2g * pdf_r_y * r2x(pdf_r_x) * s.rho / s.dv,\n", + " color='black'\n", + " )\n", + "\n", + " vals = np.empty(len(self.r_bins) - 1)\n", + " for i in range(len(vals)):\n", + " vals[i] = state.moment(1, attr='x', attr_range=(self.x_bins[i], self.x_bins[i + 1]))\n", + " vals[i] *= s.rho / s.dv\n", + " vals[i] /= (np.log(self.r_bins[i + 1]) - np.log(self.r_bins[i]))\n", + "\n", + " pyplot.step(\n", + " m2um * self.r_bins[:-1],\n", + " kg2g * vals,\n", + " where='post',\n", + " label=f\"t = {t}s\"\n", + " )\n", + " pyplot.grid()\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", + "\n", + "\n" + ], + "metadata": { + "collapsed": false, + "pycharm": { + "name": "#%%\n", + "is_executing": false + } + } + }, + { + "cell_type": "code", + "execution_count": 17, + "outputs": [], + "source": [ + "class SetupA:\n", + " x_min = r2x(10e-6) # not given in the paper\n", + " x_max = r2x(100e-6) # not given in the paper\n", + "\n", + " n_sd = 2 ** 13\n", + " n_part = 2 ** 23 # [m-3]\n", + " X0 = 4/3 * np.pi * 30.531e-6**3\n", + " dv = 1e6 # [m3]\n", + " norm_factor = n_part * dv\n", + " rho = 1000 # [kg m-3]\n", + "\n", + " dt = 1 # [s]\n", + " steps = [0, 1200, 2400]\n", + "\n", + " kernel = Golovin(b=1.5e3) # [s-1]\n", + " spectrum = Exponential(norm_factor=norm_factor, scale=X0)\n", + "\n", + " # TODO: rename?\n", + " def check(self, state, step):\n", + " check_LWC = 1e-3 # kg m-3\n", + " check_ksi = self.n_part * self.dv / self.n_sd\n", + "\n", + " # multiplicities\n", + " if step == 0:\n", + " np.testing.assert_approx_equal(np.amin(state['n']), np.amax(state['n']), 1)\n", + " np.testing.assert_approx_equal(state['n'][0], check_ksi, 1)\n", + "\n", + " # liquid water content\n", + " LWC = self.rho * np.dot(state['n'], state['x']) / self.dv\n", + " np.testing.assert_approx_equal(LWC, check_LWC, 3)\n", + "\n" + ], + "metadata": { + "collapsed": false, + "pycharm": { + "name": "#%%\n", + "is_executing": false + } + } + }, + { + "cell_type": "code", + "execution_count": 18, + "outputs": [], + "source": [ + "def run(setup):\n", + " x, n = constant_multiplicity(setup.n_sd, setup.spectrum, (setup.x_min, setup.x_max))\n", + " state = State(n=n, extensive={'x': x}, intensive={}, segment_num=1)\n", + " collider = SDM(setup.kernel, setup.dt, setup.dv, n_sd=setup.n_sd)\n", + " runner = Runner(state, (collider,))\n", + "\n", + " states = {}\n", + " for step in setup.steps:\n", + " runner.run(step - runner.n_steps)\n", + " setup.check(runner.state, runner.n_steps)\n", + " states[runner.n_steps] = copy.deepcopy(runner.state)\n", + "\n", + " return states, runner.stats\n" + ], + "metadata": { + "collapsed": false, + "pycharm": { + "name": "#%%\n", + "is_executing": false + } + } + }, + { + "cell_type": "code", + "execution_count": 20, + "outputs": [ + { + "data": { + "text/plain": "
", + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEOCAYAAABmVAtTAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzs3Xlc1NX6wPHPAYZNQHYXUHFBEZcsyeWallmJaNo1t0xTq5+23rSytO1a11t2M9tNy3KpNJfbTVNDM820tMQyXIIQFQVRdpB9O78/hhkBYRhgNuC8X6955cx85/t9qOSZ8z3nPI+QUqIoiqIoAHbWDkBRFEWxHSopKIqiKHoqKSiKoih6KikoiqIoeiopKIqiKHoqKSiKoih6KikoiqIoeiopKIqiKHoqKSiKoih6KikoiqIoeg7WDqC+fH19ZVBQkLXDUBRFaVKOHj2aJqX0q+u4JpcUgoKCiIqKsnYYiqIoTYoQIsGY49TtI0VRFEVPJQVFURRFTyUFRVEURa/JzSkoLU9JSQmJiYkUFhZaO5RmwdnZmcDAQDQajbVDUWyQSgqKzUtMTMTd3Z2goCCEENYOp0mTUpKenk5iYiKdO3e2djiKDTLb7SMhxKdCiBQhxAkDx9wihDgmhDgphNhvrliUpq2wsBAfHx+VEExACIGPj48adSm1MudIYQ3wPrCupjeFEJ7AciBcSnleCOFvxlhavPW/nGfrsaRrXh/XL4CpAztaIaL6UQnBdNS/S8UQs40UpJQ/AhkGDpkKfCWlPF9xfIq5YlFg67EkTiXnVHntVHJOjYlCqSorK4vly5eb5FxFRUVMnjyZbt26MXDgQM6dO2eS8yqKqVhz9VF3wEsI8YMQ4qgQ4j4rxtIihLbz4FaiSd34PF/OHkRoOw9rh9QkmDIpfPLJJ3h5eXH69GnmzZvHs88+a5LzKoqpWDMpOAD9gdHASOBFIUT3mg4UQswWQkQJIaJSU1MtGWOzIqXk5ZdfZt++fSQmJlo7nCZjwYIFxMfH069fP+bPn9+oc23dupUZM2YAMGHCBL7//nuklJw8eZIBAwbQr18/+vbtS1xcnClCV5R6s+bqo0QgTUqZB+QJIX4ErgP+qn6glPIj4COAsLAwadEom5HU1FSSk5MBOHr0KNDGugE1wMvfnOTUxZy6D6yH0PYe/PPOXrW+v2TJEk6cOMGxY8dqfH/o0KFcuXLlmteXLl3KbbfdVuW1pKQkOnToAICDgwOtW7cmPT2dFStW8MQTT3DvvfdSXFxMWVlZI34iRWk4ayaFrcD7QggHwBEYCLxlxXiavcTERLp27cq5c+e0SSEwwtohNQsHDhww+lgpr/1OI4Rg8ODB/Pvf/yYxMZHx48cTHBxsyhAVxWhmSwpCiA3ALYCvECIR+CegAZBSrpBS/imEiASigXJglZSy1uWrSuPk5ORw5coVXpw3jxUrVnD06FHcm2BSMPSN3lrqM1IIDAzkwoULBAYGUlpaSnZ2Nt7e3kydOpWBAweyY8cORo4cyapVq7j11lst9SMoip7ZkoKU8h4jjnkDeMNcMShXJSUl4eBgz4wZMzhy5Ajffvstt4yzdlRNg7u7e42/9HXqM1IYO3Ysa9euZfDgwWzZsoVbb70VIQRnzpyhS5cu/OMf/+DMmTNER0erpKBYhap91AIkJSWRmppK27btcHNzo3///qSkpFBUVGTt0JoEHx8fhgwZQu/evRs90fzAAw+Qnp5Ot27dWLZsGUuWLAFg48aN9O7dm379+hETE8N996nFeIp1qDIXLcDy5cuR0p2AgAAA+vfvD0Bubi5OTk7WDK3JWL9+vUnO4+zszObNm695feHChSxcuNAk11CUxlAjhWauoKCAlStX4uvri7OzMwD9+vXDzs7O4C0RRVFaJpUUmrkvvviC9PR0/SgBwNXVlZ49e6qkoCjKNdTto2aocp2jqKgCuvzfu1wssMPT8+ox/fv3Z1+uSgqKolSlRgrNkK7OUVZWFnl5eQQGBBLazoNx/a6OFvr3709xcQnFxWqyWVGUq9RIoZkKbedBr8vfsXvDP4nPzaVVq1ZV3u/fvz8cPsCVK7lWilBRFFukRgrNWFxcHB06dLgmIYB2shlQ8wqKolShkkIzFhcXV2u5hFatWuHq6soVNa9QJ1NWSf3xxx+54YYbcHBwYMuWLfrXjx07xuDBg+nVqxd9+/Zl48aN+vfOnj3LwIEDCQ4OZvLkyRQXFwOqDLdiHiopNGOGkgJod+rmqttHdTJlUujYsSNr1qxh6tSpVV53dXVl3bp1nDx5ksjISObOnUtWVhYAzz77LPPmzSMuLg4vLy8++eQTQJXhVsxDJYVmqrS0lIyMjDqTQnFxMRcvXrRgZE2PKUtnBwUF0bdvX+zsqv7V6969u/6/Vfv27fH39yc1NRUpJXv37mXChAkAzJgxg6+//hpQZbgV81ATzc1UQUE+gMGk4ObmBmjLaLdv394icTXatwvg0nHTnrNtHxi1pNa3TVk62xi//vorxcXFdO3alfT0dDw9PXFw0P5VDQwMJClJu9xYleFWzEElhWYqP78AMD4p3HnnnRaJqzmqT0G8uiQnJzN9+nTWrl2LnZ1draW2QZXhVsxDJYVmqqCgADs7O7p06VLrMfb29ri6unL06E8WjKyRDHyjtxZTjRRycnIYPXo0ixcvZtCgQQD4+vqSlZVFaWkpDg4OJCYm6kd1qgy3Yg4qKTRTBQUFdOzYsc6Cd+7u7hzdedRCUTVNpiydXZvi4mL+/ve/c9999zFx4kT960IIhg8fzpYtW5gyZQpr165l3DhtzXNVhlsxB7NNNAshPhVCpAghDDbOEULcKIQoE0JMMFcsLVFBQYFRtw3c3NxITk7Wt+lUrmXK0tlHjhwhMDCQzZs3M2fOHHr10jYN2rRpEz/++CNr1qyhX79+9OvXTz+H8frrr7Ns2TK6detGeno6DzzwAKDKcCvmIWq6L2mSEwsxDMgF1kkpe9dyjD3wHVAIfCql3FLTcZWFhYXJqKgok8ba3ExeeYiffjrIXR4JvP/++waPO34hg4wzx+nduzc+Pj6M6xfA1IEdLRht3f7880969uxp7TCaFfXvtOURQhyVUobVdZzZRgpSyh+BjDoOexz4L5BirjhaopKSEkpLy+ocKYzrF0Bo+9aAtrfCqeQcfSE9RVFaJqvtUxBCBAB/B1ZYK4bmqqCg7pVHAFMHdmTLIzfh9fta/E+sJ7SdhyXCUxTFhllz89rbwLNSyjoXUQshZgshooQQUampqRYIrWkzNino9OnTh5iYGHOGpChKE2HN1UdhwJcVa659gQghRKmU8uvqB0opPwI+Au2cgkWjbILy8/MRQrt71hht27YlJUXdwVMUxYpJQUrZWfdnIcQaYHtNCUGpv4KCApydXdBoNEYd36ZNG7KzsykvL7+m/IKiKC2L2ZKCEGIDcAvgK4RIBP4JaACklGoewYwKCgpwcXEx+nh/f39AO0Fd174GRVGaN3OuPrpHStlOSqmRUgZKKT+RUq6oKSFIKWcasxxVqZuUkoKC/AYmhWJzhdWkmbJK6rJlywgNDaVv376MGDGChISEKu/n5OQQEBDAY489pn/t6NGj9OnTh27duvGPf/xDX94iIyOD22+/neDgYG6//XYyMzNNEqPSsql7Bc3MpUuXKCsrr1dSaNOmDQDFxSXmCqtJM2VSuP7664mKiiI6OpoJEybwzDPPVHn/xRdf5Oabb67y2sMPP8xHH31EXFwccXFxREZGAtpCfSNGjCAuLo4RI0boN68pSmPUmhSEEN5GPDxr+7xiHboyya6u9R8p6Jq3KFWZsnT28OHDcXV1BWDQoEEkJibq3zt69CiXL1/mjjvu0L+WnJxMTk4OgwcPRgjBfffdV2Pp7Moltffv36/fFX399der7npKvRiaU7hY8RAGjrEHbGv7awunSwouLq5Gf6bynIKte/3X14nJMO3y2RDvEJ4dUHuDGnOVzv7kk08YNWoUAOXl5Tz11FN89tlnfP/99/pjkpKSCAwM1D+vXDr78uXLtGvXDoB27drpV5AtXbqUDz74gCFDhpCbm4uzs3OtMShKdYaSwp9SyusNfVgI8buJ41EaKS4uDiG86jVh7ObmhqurqxopNFBDCuJ9/vnnREVFsX//fgCWL19ORESEvj+CjqHS2bUZMmQITz75JPfeey/jx4+vklQUpS6GksJgIz5vzDGKBcXFxeHsNazOXxzV+fv7N4mRgqFv9NZS35HCnj17+Pe//83+/fv1yfvQoUMcOHCA5cuXk5ubS3FxMW5ubjzxxBNVbjFVLp3dpk0bkpOTadeuHcnJyfoR34IFCxg9ejQ7d+5k0KBB7Nmzh5CQEHP86EozVGtSkFIW6v4shPAC2gMFwDkpZXn1YxTbEBcXh+vN4fX+nL+/PzlqpFAjU5bO/v3335kzZw6RkZH6X+IAX3zxhf7Pa9asISoqSj9x7O7uzuHDhxk4cCDr1q3j8ccfB66Wzl6wYEGVktrx8fH06dOHPn36cOjQIWJiYlRSUIxmaKK5tRDiOSHEceAwsBLYBCQIITYLIYZbKkjFOOXl5Zw+fbpeK4902rRpQ7FaklojU5bOnj9/Prm5uUycOJF+/foxduzYOj/z4Ycf8uCDD9KtWze6du2qn4dYsGAB3333HcHBwXz33XcsWLAAgLfffpvevXtz3XXX4eLioj9eUYxh6PbRFmAdMFRKmVX5DSFEf2C6EKKLlPITcwaoGO/ixYv13rim4+/vzzG1JLVW69evN8l59uzZU+cxM2fOZObMmfrnYWFhnDhxbVsSHx+fKpPSOu+9916jYlRaNkO3j24X2hvTgUBWtfeOAqpdl425uvKoYSOFknNqpKAoLZ3BzWtSu/RB1SNqIhqyR0HH398fKaG0VI0WFKUlM6b20WEhxI1SyiNmj0ZplLi4OJycnHByqv+6dO2kZ4La1aw0CZkbN5Gzffs1r3uMGYPX5ElWiKj5MKbMxXDgkBAiXggRLYQ4LoSINndgSv3FxcXRtWvXBn32aqkLdQtJsX0527dTWK0HSGFMTI2JQqkfY0YKaulCExEXF2d0Y53qmtKuZkUBcA4JodNn6/TPE6bfZ8Vomg9DS1KjhBDvAD2By1LKhMoPy4WoGKO8vJz4+PhGJwU1UlCUls3Q7aNBwP/Q9kTYL4TYKYR4QgjR3SKRKUZZ/8t5Jq88xPj39+M5fhG/tBrMqeScep/Hx8cHUCOFmpiySuqKFSvo06cP/fr146abbuLUqVP696Kjoxk8eDC9evWiT58+FBZq94aq0tmKJdWaFKSUpVLKH6SUC6SUA4EHgCvAYiHEb0II0/wtURpl67EkTiXnkJ+v7cvs4uJCaDsPxvULqNd57O3t0Wg0aqRQA1MmhalTp3L8+HGOHTvGM888w5NPPglAaWkp06ZNY8WKFZw8eZIffvhB3zlPlc5WLMngRLMQwl4I8QaAlDJZSvmplHIS2v7KX9Tx2U+FEClCiGt33Wjfv7di4jpaCPGzEOK6hv4QLV1oOw8m+V7k8oaFrJzck41zBjN1YP2L1zo6OqpdzTUwZelsDw8P/Z/z8vL0Nap2795N3759ue467V8DHx8f7O3tVelsxeIMTjRLKcuEEP2FEEJWKtdYUfvopzrOvQZ4H+2u6JqcBW6WUmYKIUYBHwEDjY5cqSIhIQE7OzsCAuo3QqhM46ihxMaXpF569VWK/jRt6WynniG0fe65Wt83densDz74gGXLllFcXMzevXsB+OuvvxBCMHLkSFJTU5kyZQrPPPOMKp2tWJwxq49+B7YKITYDeboXpZRfGfqQlPJHIUSQgfd/rvT0MNqd00oDJSQkEBAQoL/l0BCOGkdyCus/H9HS1bd09qOPPsqjjz7K+vXrWbx4MWvXrqW0tJSDBw9y5MgRXF1dGTFiBP37968ystBRpbMVczImKXgD6cCtlV6TgMGkUE8PAN+a8Hwtzrlz5wgKCmrUORwdHW1+pGDoG721NLTJzpQpU3j44YcB7Qjg5ptvxtfXF4CIiAh+++03pk2bpkpnKxZVZ1KQUs4yZwAV1VYfAG4ycMxsYDZAx46q0VtNEhISGDZsWKPOodFoKCsrIz8/X98yUjFt6ezKe0l27Nih//PIkSP5z3/+Q35+Po6Ojuzfv5958+bRrl07VTpbsahak4IQ4j20I4IaSSn/0diLCyH6AquAUVLKdAPX+gjtnANhYWG1xtRSSSlJSkqiU6dOjTqPo6MjAKmpqY0+V3NSuXT2qFGjeOONNxp8rvfff589e/ag0Wjw8vJi7dq1AHh5efHkk09y4403IoQgIiKC0aNHA9rS2TNnzqSgoIBRo0ZVKZ09adIkPvnkEzp27MjmzZsBbensffv2YW9vT2hoqCqdrdSLoZFCVMU/hwChwMaK5xMxQYVUIURHtLegpksp/2rs+Vqy4uIiysrKGv2LXDcfcfnyZZUUqjFV6ex33nmn1vemTZvGtGnTrnldlc5WLMlQ6ey1AEKImcBwKWVJxfMVwO66TiyE2IB245uvECIR+CegqTj3CuAlwAdYXjFxViqlDGvEz9Ji6TY5mWJOAdCvYlEUpeUxZqK5PeAOZFQ8d6t4zSAp5T11vP8g8KAR11fqUFhYBGCykYJKCorSchmTFJYAvwsh9lU8vxlYZLaIlHrTjRQaOwmvGylcvny50TEpitI0GbP6aLUQ4luubixbIKW8ZN6wlPooLCqkTZs2jd6kZGdnh729PSkpF00UmelIKetcn68Yp9I+VEW5hjEjBSqSwFYzx6I0UGFhYaPnE3Q0Go3NjRScnZ1JT0/Hx8dHJYZGklKSnp6udjkrtTIqKSi2raiwyGSrhRwdHUk5b1tzCoGBgSQmJpKammrtUJoFZ2dntctZqZVKCs1AYWGhyZKCxlFjcxPNGo2Gzp07WzsMRWkR6mzHKYT4zJjXFOsoLi5GSmm6kYLG0eZuHymKYjnG9GjuVfmJEMIe6G+ecJT6MtUeBR1HR0fS0tIoKyszyfkURWlaDLXjXCiEuAL0FULkVDyuACmoSWeboUsKJrt9pNFQXl5ORkZG3QcritLsGOq89pqU0h14Q0rpUfFwl1L6SCkXWjBGxYCiItNsXNNRexUUpWUzVBAvREoZA2wWQtxQ/X0p5W9mjUwxSmFhIQ4ODri7u5vkfGpXs6K0bIZWHz2Jtlz1mzW8J6naX0GxksLCQpOuOVf1jxSlZTNUEG92xT+HWy4cpb4KCwtN2vvA0fFqpVRFUVoeo/YpCCH+BgRVPl5KWVvvZcVCpJQUFRXi5e1lsnM6OGgqSl2okYKitER1JoWKPQldgWOAbp2iBFRSsLKMjAzKyspxdjJtyQI/Pz+VFBSlhTJmpBAGhEpVRcvmnDt3DsDkdWzatGmjbh8pSgtlTFI4AbQFkutzYiHEp8AYIEVK2buG9wXwDhAB5AMz1Yqm+klISABMnxT8/f3VSEGxisyNm8jZvv2a1z3GjMFr8iQrRNTyGLOj2Rc4JYTYJYTYpnsY8bk1QLiB90cBwRWP2cCHRpxTqeRqUnAy6Xn9/f3VSEGxipzt2ymMianyWmFMTI2JQjEPY0YKixpyYinlj0KIIAOHjAPWVdyWOiyE8BRCtJNS1mtE0pIlJCRgb98eBweNSc/bpk0bNVJQrMY5JIROn12dskyYfp8Vo2l5jGmys99M1w4ALlR6nljxmkoKRjp37hxOfl1Mfl5/f3/y8/PJy8ujVatWJj+/oii2y5jVR1fQrjYCcAQ0QJ6U0qOR166pW0qNk9lCiNlobzE1uuVkc5KQkIBzoOmbpbRp0wbQ7lXo0sX0SUdR6qswJqbKiKEwJgbnkBArRtR81TmnUFHvSFf7yBm4G3jfBNdOBDpUeh4I1NgHUkr5kZQyTEoZ5ufnZ4JLNw8JCQlm6aDl7+8PqF3Nim3wGDPmmgTgHBKCx5gxVoqoeat3kx0p5ddCiAUmuPY24DEhxJdo+z9nq/kE4+Xk5JCZmUmIiZPCqeQc3r/iQJt7XuOFvWn4/HGIcf0CmDpQjdCsafNfm9l5Zuc1r0d0iWBi94lWiMhyvCZPUiuPLMiY20fjKz21Q7tvoc49C0KIDcAtgK8QIhH4J9pbT0gpVwA70S5HPY12SeqsesbeopljOeq4fgHA1cqrxcXFnErOAVBJwcp2ntlJbEYsPbx76F+LzYgFaPZJQbEsY0YKd1b6cylwDu3KIYOklPfU8b4EHjXi+koNdEnBycl0y1GnDuzI1IEdKSwsxOWJ4Ty+eDHRvqruoa3o4d2D1eGr9c9nRarvUYrpGbP6SP2fZ4PMtXFNd04PDw/tnIKvyU+vKIoNM9R57QUhhLeB928VQqiZHitJSEjA0dFRX+ra1FSpC0VpmQyNFI4D3wghCoHfgFTAGe0O5H7AHuBVs0eo1CghIcFk3dZq4ufnR2pqqhooKEoLY6ifwlZgqxAiGBgCtANygM+B2VLKAsuEqNTk3LlzZk8KZ86cUUlBUVoYQ+04FwKRUsrfgTjLhaQYIyEhgTFjxnDFTOf39fXl119/NdPZFUWxVYZuH50FnhBCXAf8AXwL7JZSZlokMuUa6385z9ZjSZSXl8OtcznhH0Recg6h7Rq7ufxafn5+pKWlmfy8iqLYNkO3j74EvgQQQlyPtuLpV0IIe7TzCZFSSvVV0oK2HkviVHIOnT21/9mcnZ3p1MZDv7/AlHx9fSkpKaGsrBR7+3rvcVQUpYky6m97xS2k34HXhBAewO3Ag4BKChYW2s6DGR0y+WbeQv578CBDhgw2y3V05USKi0twcVFJQVFaCoN/24UQA9DuMzsihAhFO1qIkVL+F/ivJQJUrhUfHw9A165dzXYNX1/tFHNJSQkuLi5mu46iKLbF0ETzP9E2wnEQQnyHtj7RD8ACIcT1Usp/WyZEpbr4+HhatWqlr2ZqDrqRQklJidmuoSiK7TE0UpiAdj+CE3AJCJRS5ggh3gB+AVRSsJL4+Hi6du2KtqOpeVQeKSiK0nIYKp1dKqUsk1LmA/FSyhyAiv0J5RaJTqmRLimYk0oKitIyGUoKxUII14o/99e9KIRojUoKVnXmzBmzJwU3NzecnJxUUlCUFsbQ7aNhUsoiACll5SSgAWaYNSqlVkVFRRQVFZk9KQgh9MtSFUVpOQztUyiq5fU0QO1qspKCAm11EXMnBdBONherpKAoLYoxTXYWSSkXWSAWxQiFhYWAZZKCr68vCbacFKJWw/Et177eZwKEqYrvitIQhpak2gEfAw1u1CuECAfeAeyBVVLKJdXe7wisBTwrjlkgpby256CiV1BQgIODAx07mr8Tmp+fH3ElxWa/ToMd3wKXjkPbPldfSziofVRPFipRKIpRDI0UvgFOSSkXNuTEFeUwPkC7+zkROCKE2CalPFXpsBeATVLKDys2x+0EghpyvZaioKCATp064eBg/l3Gvr6+lGbZ8EgBtAlh1o6rz2saPVw6rv2nSgqKUidDv1nCaNxehAHAaSnlGQAhxJdo23hWTgoS0FVzaw1cbMT1WoSCwgK6detmkWv5+flRmlaGtnNqExE269pf/qtHWycWRWmCDCWF4cBGIcSDUspfGnDuAOBCpeeJaHdFV7YI2C2EeBxoBdzWgOu0KIUFhRaZTwDdXoVstQLJhsVmxFbp1RzRJYKJ3SdaMSKlqat1n0LFbZ6RwBsNPHdN222rf+W8B1gjpQwEIoDPKuYyqp5IiNlCiCghRFRqamoDw2n6SktLKS0ttVhSUKUubFtElwh6ePfQP4/NiGXnGTUlpzSOwRvTUsqLQoiGjr0TgQ6Vngdy7e2hB9AW2UNKeUgI4Yy2VXyVyW0p5UfARwBhYWFN6F6GaVlyOSroRgqnm0RSyMnJISMjo8prjo6OtG/f3koRmd/E7hOrjAoqjxgUpaHqnK2UUja0udcRIFgI0RlIAqYAU6sdcx4YAawRQvRE2wO65Q4F6mDppGDVkUJNE8a1rCD64osvmD17Nvn5+de8N3fuXJb1lYgaB66KolRnzD6FMOB5oFPF8QJtOe2+hj4npSwVQjwG7EK73PRTKeVJIcQrQJSUchvwFPCxEGIe2ltLM2WTmtW0LF1S6NKli0WuZ9X6R9WXm9awgqi8vJzT8aeZ9to0hg4dyqxZs6oUCTx48CBvv/02j73Qgy5duqi0YCaZGzeRs337Na97jBmD1+RJVohIaQxj1jV+AcwHjlPPmkcVew52VnvtpUp/PgUMqc85W7LCwkIcHR1xdXWt+2AT8PHxAaC42Ep7FSovN109WpsYKlYSFRYWUpL4GxcvFjN//nxeffXVa5bpzpgxAycnJy5cWIOdnR2dLR1/C5GzfTuFMTE4h4ToXyuMiQFQSaEJMiYppFZ8q1esrKCgwKINbxwcHHBwcKCk1AbmFPpM0P8xIyODU39qVzZ7Dv0//vPAf2r8iBCC9957j7jnd5GQkMDn//oXL774okXCbWmcQ0Lo9Nk6/fOE6fdZMRqlMYxJCv8UQqwCvgf09ZCklF+ZLSqlRgUFBXh7e1v0mhqNhpJiG0gKFfsPfv31V4YOHUpoaChbtmypc37Fzs6O7t27I6XkpZdewtHRkWeffdZCQVvf5r8217giSS1dVWpjTFKYBYSgrY6qu30kAZUULKigoIDi4mKcnZ0tel2NRmMzq49SUlK4++67adeuHXv27NHf3qqLQNCjdTEn5weR8ucrZC7bipeXV4sofbHzzE5iM2KvWboKqKSg1MiYpHCdlLJP3Ycp5nTmzBkAi/dL1jhqKCwotOg1a1JaWsqUKVNIS0vjp59+MjohANBnAgIIaVNOVlYW5y+cx6soUfteE0wK5eXlXLx4kbNnz+ofwcHBUMsgsod3D1aHr9Y/V0tXFUOMSQqHhRCh1WoWKRYWHx8PWCEpaDRcyWnoqmTTWbhwIfv27WPt2rXccMMN9ftwxa0nO+BA6ussWLCAK+8Ow80skZpPeXk5jz76KKtWrapx8n/gsoF0797dCpGZR3l5OWlpaSQlJXHx4kWSkpIoLCxk5syZeHh41H0CpUGMSQo3ATOEEGfRzikYtSRVMa2rScGyt48cNY6UlJQgpTRrT2hDNm3axNKlS3nsscfaLx8aAAAgAElEQVS4777GTWDOmTOHxYsXcyHxAj1DepooQvMrKiri5MmTRC+PZubMmQwYMIDOnTvTpUsXOnTowJtvvsmnKZ+Sl5dHXPc47cihCduyZQuzZ88mMzPzmveWL1/O119/TUil1U6K6RiTFMLNHoVSp/j4eBwcOuHgoLHodTUaDVJKcnJyaN26tUWvDXDixAnuv/9+hgwZwptvvtno83l6evLggw+SkvIxnTt3xrIptmF+/PFHjh49Snl5OVu2bOHuu+++5pgXXniBI+uP8OeffxIWFsa6desYN26cFaJtnHIpORMfz8SJExk0aBD33nsv7du3JyAggPbt2/PXX38xZcoUBgwYwOeff87YsWOtHXKzU2vto4paQ+8APYHLUsqEyg/LhagAnD59Gmdny946Am1SAEhLs3yzvXJZzuTJk3F3d2fz5s04Ojqa5LxPPPEEAElJSSY5n7lIKXn33XcZMWIEDg4O3HDDDTUmBB1vb2/69+9P9+7dueuuu3jnnXcsGG3jnT9/nmPHficxKYm5c+eyf/9+HnvsMcaPH8/AgQPp0KEDI0aM4OjRowQHBzNu3DhefvllystVy3hTqjUpAIOA/wG3APuFEDuFEE8IIZrPTcsmJD4+3uLzCXA1KVijEGFCQgKnTp3i008/pV27diY7b1BQEH5+fly8eJHs7GyTndfUXnrpJZ544gkiIiK44YYbjNq06OzszIEDBxg9ejQLFy603sbDetq1axc33HADeXn59OoVyltvvVXrl4COHTty8OBBpk+fzqJFixg/fjylpaUWjrj5MtSjuRT4oeKBEKIdMApYLIQIBg5JKR+xQIwtXmlpKefOneN6KyYFs44UaqhzVHbxGOfPZzNt2jRGjRpl8kt26NCBlJQUVq1axVNPPVVrHIBVlq5u2LCBxYsXc//99/Pxxx/zwO4HjP6ss7Mz7777LiEhIZw7d87mJ5/37NnD6NGj6dWrF/2Du+NqxP/nLi4urF27lj59+vDMM8+wevVq7rBArC2BoZFCFVLKZCnlp1LKSUB/tOUvFAu4cOECpaWlzXekoKtzVEFKybFLZWw768Rbb71llku6u7nj6enJ22+/fXUfRrU4AO3zmhKFGR05coT777+foUOH8uGHH2JnZ/RfU70uXbrw8MMPk5ycXGOhQFsRExPDhAkT6NmzJwcOHDAqIegIIXj66acZPHgwixYtokzdRjKJOv9vE0J8I4TYVvmBtq9yWEWpa8XMdCuPLL1xDSw4p6CrczRrB8syRxD2/mUGPLxcX5TPHDp06EBiYiIbN26sMQ5m7aja/9kCkpKSGDduHG3atOG///1vo+ZRXnjhBezt7Tl79qwJIzSdtLQ0xowZg5OTE9u3b2/QMlMhBEuWLNEuWU1MNEOULY8xX0HOALnAxxWPHOAy0L3iuWJm1tqjAGBvb4+dnbDYnMLp06d58cUXGTt2LJMmmbeYmre3N6GhoSxdutQmWo4WFBRw1113ceXKFb755ht96fKG8vPzo0OHDqSlpXH48GETRWkaRUVFjB8/nsTERLZu3UqnTp0afK5hw4YxevRozl84T4maW2g0Y5LC9VLKqVLKbyoe04ABUspHgXruIlIaIj4+HicnJ5ycnKxyfY1GY5HVR1JK/u///g+NRsPy5cvNvi9CIJg7dy5//PEHUVFRZr1WXaSU3H///Rw9epQvvviCPn1MM0IJDAxEo9HwzDPP2ETiA+3POmfOHA4cOMCaNWsYNGhQo8/52muvUVpaxvnz500QYctmzD4FPyFERynleQAhREe03dEAmsbShiYuPj6ezp2tV/hZo3EkNdn8I4VVq1bxww8/8NFHHxEQEGD26wH8/e9/56GHHuKbb77hRiv8K9YVrEtKSuJ00GmGrxjO/xz/R9FfRSapTWRvb09QUBC7D+xmx44djBkzxgRR11/lngvnL5xn+JmzzBo6lKAdO0nYcbVgX/US3Mbq06cPF9u0ISkpkQsXLtChQ4drrquj+jwYZsxI4SngoBBinxDiB+AAMF8I0Qrt3IJiZvHx8RbrtlYTS4wUikuKefbZZxk2bBgPPvigWa9Vma+vL0OGDGHbNutUh995Zien0k5x5swZvL296dixo8l7Lbdr147g4GAWLFhAWVmZyc5bH7qeC9k52Zw9cxZ/Pz86dQq65jjnkBA8Gpi4gjoHgYRFixZdc12dwpiYGhsCKVcZ045zZ8US1BC0JS5ipJS6CmlvG/qsECIceAdt57VVUsolNRwzCViEtvLqH1LK6i07W6T1v5xn6zHt5qrUPvdi37Ytxck5hLazfM0XjUZj9jmFM2fOcOXKFT788EOLl9MYO3Ys8+fPp7DIB2cny07mSynJT8gn+6NsDh4/SNu2bU1esE4IwauvvsrEiRP5/PPPoY1JT280h+BuhB86hJ3Ggd/37TV5/SJnJ2faBwSwZs0annrqKUJDQ7WvV+r1oPo81M3Qjmb9fIGUskhK+YeU8lilhFDlmBo+bw98gHZvQyhwjxAitNoxwcBCYIiUshcwt8E/STOz9VgSp5JzKCoqoqysDFdXV0LbeTCun2Vuq1Rm7pFCdk42ly5dYt68efq/yJakK5WQnpZu8WufPXeWvNw8Vq1aRdu2bc12nbvvvpvQ0FBWrVpltmvUJS4ujgsXLvD555+braBdp44dcXNz47nnnjPL+VsCQ7ePVgshvIQQ3rU9gE8MfH4AcFpKeUZKWQx8CVQvxvJ/wAdSykwAKWVKY36Y5ia0nQezu+VzecNCltzRlo1zBjN1YEeLx6HRaMjOzjbL7tjS0lLi4uJwcnLipZdeqvsDplTR3rP7T/P45WFfWhdadpLywIEDXDh/gbZt25q9TpEQgqlTp3Lw4EGKiorq/oCJXU5J4fLlFF588UUGDx5stutoNBrmzp3L1q1b1aRzAxlKCq2Bo3U8DHVfCQAuVHqeWPFaZd2B7kKIn4QQhytuN11DCDG7ohZTlDXKLVhTdHQ0AL1797ZaDLq9Cunppv8mvXz5cnJzc+nWtRtubhYsZt1nQpU9CD4+PhxNKqYg2DITsTk5OUyfPh1nZ2e6detmkWtOnjwZ0DYrsqSEhATi/voLDw8Pnn/+ebNfT1dJ98svvzT7tZojQ2Uughp57ppuDFdfE+cABKOtrxQIHBBC9JZSZlWL5SPgI4CwsDDbWFdnIdHR0QQFBVmlQqlO5V3NpqxBlJyczAsvvMCBB7zw9TPfJrUaVfRY0Ll44AC3DBvGpghfzNWPrHJrzJiYGOzuscOtixv29vZmumJV3bp148YbbyQlJUW/OsfcysrKmD59OrOBnj174uBgzILHxunatSsDBgzgyy+/ZHIv632ZaqrM+V8oEaj8f14gcLGGYw5LKUuAs0KIWLRJ4ogZ42pSoqOjTbZmvaHMtav56aefpqioiODgYESN3yEsZ/Dgwfj4+LBt2zYmTjRPWtC1xvQp9+Hy5ct07NSRzn6diegSYbJrxGbEVpmort6Kc8qUKbyf+j4FBQV1nqum5ZxQvyWdr7/+OgcOHGDp7XfgYsEd+ffccw/z5s0jv0sXXF3qLiSoXFX/oirGOwIECyE6CyEcgSlA9XV/XwPDAYQQvmhvJ50xY0xNipTlxMTE0LevdfsZOTqaPin88MMPrF+/nmeeecYm/tI6ODgwevRoduzYYdaKm13cu3Dk6SO0jmzNdw98x+rw1SbrlRzRJaJKAgBtK87KSac+t5CqL+eE+i3pzM3LZdGiRUyaNAn/NpZd8jRp0iSEEKRcVtOU9WVwpCC0awMDpZQXDB1XEyllqRDiMWAX2iWpn0opTwohXgGipJTbKt67QwhxCigD5kspLb8ExEbl5eVTVlZm9aRg6qJ4xcXFPPLIIwQFBbFw4ULYaBsN5MeOHcu6dev46aefuPnmm81yjdOnT5Oenk5kZKT+32tN6vrGX5OJ3SfWmWACAgJo3bo1KSkpRnXTq7ycE4xf0lkuJTExMXh7e2vnjebOM+pzptK+fXtuvvlmUpKT6RQUZNFrN3UGRwpSuy/+64aeXEq5U0rZXUrZVUr574rXXqpICEitJ6WUoVLKPlJKNTNUSV5eHoDVk4Ku25upRgpvv/02f/75J++++65RPQIs5Y477sDR0dFsG9nS0tJISUnh+eefp1+/frUeZ8w3/sbw9/cnPz+f48eP131wA51PSCA3N4+VK1fi4+NjtusYcs8995CfX0Bebq5Vrt9UGTOncFgIcaOUUt3nt7DcvFyLrk6pjRACLy8vk4wUzp8/z8svv8zYsWO58847TRCd6bi7u3Prrbeybds2li5datJNdBkZGcTFxdGqVSuee8bwGnpjvvE3hp+fH3FxcWzYsMEsXzh+++03Es4n0KZNG8Kt2BL07rvvZse/X+VySgq9rBZF02PMnMJw4JAQIl4IES2EOC6EiDZ3YArk5ebRq1cvi6zYqIufn59JRgrz5s1DSmmzrSLHjh3L6dOniY2NNel5586dS0lJCSEhISZrK9pQGo0GLy8vvvzyS5MXySsqKmLGjBk4ahyt/mXGx8cHby8vUlNSVMvOejAmKYwCugK3AncCYyr+qZhZXl6u1W8d6fj6+jZ6pBAZGclXX33FCy+8QJCN3ufVFYwz5S2kb775hs8++4yOFbttbYF/G3/OnTvHL7/8YtLzvvLKK5w4cYLuPXqgsYEvM/7+/hQWFXHo0CFrh9Jk1JkUpJQJaJeOlqDdZ6B7KGZUUlJMcXGJzSSFxo4UCgsLeeyxx+jevfvV9pc2qEOHDlx//fUmSwoZGRnMmTOHPn36NKpngKn5+vji5ORk0g1eR44cYcmSJdx///34eHub7LyN4ePri52dHRs2bLB2KE2GMZ3XHkfbVOc7YEfFQ5UZNLPcXNuYZNZp7Ejh9ddfJz4+ng8++MBqfSGMdeedd/Lzzz+TmZnZ6HP94x//IDU1lTVr1li80J8hDg4OREREsHHjRpNUTi0rL+e+++6jffv2LFu2zAQRmoaDvT0+Pj5s3rzZrEuNmxNjbh89AfSQUvaqWCHUR0ppG7+pmjHdyiNrb1zT0Y0UGnIPOj4+ntdee43Jkydzm+cFWD266qN6X2QrGzFiBFJKfvzxxwafY/Nfm7nz8zv5ucPP3PTeTbyX8h6xGaadp2isKVOmcOnSJfbv39/oc509e5aYmBhWr15t1d33NfH39yclJYV9+/ZZO5QmwZikcAHINncgSlW5ebk4Ojo2uiWjqfj6+lJSUkJOTk69Piel5KGHHsLR0ZE333wTjm+5Ngm07aOtRWQjBg4ciLOzc6N+iWyN3cqZ3DO0cmulv21kymWlpjBmzBicnZ355ptvGnWerOwsEhMTefTRR7nttttMFJ3p+Hh74+HhoW4hGcmYmaAzwA9CiB2AvryilNJ2xojNUF5uHq3cWlk7DD1dckpLSzP+m2DUai7veZ/nA2IIfjaYgN2ztQmhbR+YtcOM0TaOk5MTQ4YM0SaF6wIbdI64uDgK0wr5+p6vbWa0V52rqyvDhg0jMjKSt956q0HnuHLlCjExsbi4uPD666+bOELTsLOzY9y4cWzdupWXIiKsXlLF1hkzUjiPdj7BEXCv9FDMpLS0lLy8PNxa2cZKlVPJOXyW5EObe17jsa9OM3nlIdb/UndZ4uLf1uOSfRoPDw/at2+vfdHGRgW1GT58ONHR0RSX1L9c+KZNm0hNTSUoKMgmE8J1P19mynsnSZh+Hy8Vl7AgL5/YSZPI3Lip3ud66qmnKCosJCQkhFatbOdLTHURERFkZGRw5coVa4di84zpvPayJQJRrvrrr7+QUtrEXzJdUx/dX6aSkhJOJWtvIdXV2yE+Pp6UlDJ8n96D6NW0tg/deuutAGRnZdfrFt7ly5d55JFHCHgywGKVSOur59E0/JPywRu8vL0hPp7i2L/I2b69SqG71IJU0gvSWVSp3MaUjBh8XHzoBHz77bd8/PHH3PO3IbQ2U9McU7n99tuxs7MjIyMDD3fbjtXaak0KQohvMLD0VEo51iwRKfoeCrawpn3qwI5MHdiRs2fP0mX+SF779FMii+turL57924cL1+mU6dOdG5iCQEgLCyMVq1akZmVaXRSkFIyZ84ccnNzCQkJsanVRtWlBLhy/WfrkFIS3qkTK+ztr9n1m16QTn5p1Wqq+aUFUJBORkYGDzzwAL169dL2RrZxPj4+DBgwgIyMDIJq6A2tXGVopLDUYlEoVURHRyNEa5uqC+Trq+13kJaWBnV80crLy+Ohhx7iy3BXOnW0nbX59aHRaBg6dChZWVFGf2bVqlVs3bqVpUuXcsL1hBmjMx0hBCNHjiRzz/eU17CyzNXBhdXhq/XPd703AICHHnqI1NRUtm/fjt1bBlu124zw8HCufPIpJaWGeoMptc4pSCn3G3pYMsiWJjo6GldXV5v6punm5oaTk1PtexWiVuuXmF5a0p9Ph13mhvYa7OzMWZ3dvIYPH05+fr5RbUhjY2OZO3cuI0aMYN48y1YEbazw8HBKy8q4YuTKsoLMAjZv3szixYu54YZa27TbnPDwcCSYZP9Jc1br31hdjaPaHpYMsqWJjo62qZVHoP1GGRAQQGJiYs0HVCw1vXLlComJibRr1w6HgOubxKRybYYPHw5AZpbhXyLFxcXce++9ODs7s3bt2iaXCEeMGIEQ2t3XdSktKiPn4hVuvfVW5s+fb4HoTCcsLAyNgwMZ6XX/nC2ZodtHuma1j1b887OKf94L5JstohYuMzOTCxcuEGgjK48qCwoK4ty5c3SopdVAmX8oA988T06OJ8ePHwEvL8sGaGLXX389hx0cyMrKwlCLmJdeeomjR4/y1VdfERBQvQ257fP09MTD3YOMTMO/LIuKisi+kIUQsG7duiaX/Ozt7fHy9iIjM4Py8vImF7+lGLp9lFBR92iIlPIZKeXxiscCYKQxJxdChAshYoUQp4UQCwwcN0EIIYUQYfX/EZoXXY17W1h5VJ0uKdTm3Llz/Pnnn3zyySd4NfGEANpSEK1btyYrK6vWY/bt28fHv3zM3975G9tctjErchazImfZ3O7lunh7e3PlSq7BUiYvvPACJQWltA5s3SSTH2h/zuLiEv1iDuVaxqTKVkKIm3RPhBB/A+r8jSWEsAc+QFtlNRS4RwgRWsNx7sA/ANOWa2yibGnlUXVBQUEkJyfXWIY4KzuLCxcu8NBDDzFypFHfGZoET09PCgoK9LfNNv+1Wf+Lf/r26Tz8w8O0n9menNZV78fb2u7lunhXFLD77rvvanx/9+7dLF26FFdvF5w8bLt2lSHeXtqfMzIy0sqR2C5jdjQ/AHwqhGiNdolqNnC/EZ8bAJyWUp4BEEJ8CYwDTlU77l/Af4CnjQ26OYuOjsbHx8fqNfdroit3XVRUhIuLi/713NxcYmJicHZ25o033rBSdObh5akd8ezbt4/p06ez88xOfWvM2NhYSkpK6O3am8nXTTZrYxxzc3N3Q6NxIDIykqlTp1Z5LyUlhRkzZhAaGop7O9tZEdcQjo6OuLm5ERkZyYIFtd68aNEMTTQPFkIIKeVRKeV1QF+gn5Syn5TyNyPOHYC2bpJOYsVrla9xPdBBStmiq66u/+U8k1ceYvLKQ+yzuw7fSYv1G8RsiS4pFBYWVnn96aefprCwkJ4hPW1yhNMYrdxa4eDgwN69e/Wv9fDuQe+TvTny9BEe8XmELRO3NOmEACAQeHl5s3v37mtGgvfccw+ZmZls2LABYWc7K+Iaytvbm59++qnedbxaCkO3j2YAR4UQXwohZgKuUsr6FMar6f8e/UJoIYQd8BZQZ3F9IcRsIUSUECLKVM3jbcnWY0mcSs6hvLycvNxc3NzdCG3nod9NbCtqSgqRkZGsXLmSDh062Fx1TFMQCDw9PasUx8vOzubZZ59l/PjxPPnkk1aMzrS8vb25fPkyf/zxh/61K5eusHfvXlasWGEzZdwby9vbm9LS0iqJXrnK0ETzQ1LKG4BFgBewRghxSAjxqhBiWMWcgSGJQOV9/oHAxUrP3YHeaIvtnQMGAdtqmmyWUn4kpQyTUobZStVQUwtt58FD3QtJ/mIBr9zszcY5g+ssI2Fp7du3x8HBQZ8U0tPT9btaOwd1tnJ05uPl6UVCQgJnz56luLiYU6dO0bVrV1avXm1Te0kay7ticcCuXbsAKMwpIi81nzlz5jBz5kwrRmZaHh4euLu7q3mFWhjTeS1GSvmWlDIcbUvOg8BE6p4YPgIECyE6CyEcgSmAvp2VlDJbSukrpQySUgYBh4GxUkrjt5A2M3v37sXBwYGhQ4daO5Qa2dvb07FjRwoLC5Hl5cyYMYO0tDQ+++yzZr28z9PLE9BOwp48dZLS0lL++9//4mHj9X7qy9HRkX79+hEZGUlsbCzZF7LRuGhstp92Q9kJwW233ca3335r8h7VzYFRf5OFEF5CiL5AT+ASsFpKaXD5qJSyFHgM2AX8CWySUp4UQrwihFB1k2rw/fffM2DAANzdbbcIbefOnRknfmR+/Cye9t3HmZd6cv2xF2yuUY4pubq64u/vz9tvv01Odg7du3end+/e1g7LLEaOHMnBgwe56667EAI8O7W2+U55DREeHs758+eJiYmxdig2p87VR0KIfwEz0fZV0M1ASbSjBoOklDuBndVee6mWY2+p63zNWVlZKUeOHOG5556zdigGBQUFES5+JNjhEkl+fk2uJHZDCATdunXj559/ZsjDQ2jTxtBWtqZJVxE1JiIH52BnYmJjcPlbCPaauu4SN03h4eGAdk6sZ8+eVo7GthizJHUS0FVKWf/C8orRsrKyKS8vZ8SIEdYOpaqo1doSFhWea38G7+LLHE+3p+/SPxDNcHK5uiNnszhy5AjA1STYzOgqoiZ8m0BeTB6d7uyEe5k7Pi4+1g7NLDp27EhoaCjffvttk6tVZW7G3D46AXiaO5CWLisrC2dnZwYNGmTtUKqq1D5TIsnIyODYpTIOe49rlquNqktIK+DOd6L0owNDu5ubOpErOPbZMQIDA7E7YUeIdwh+Ls1zYQdoRws//vgj+fmqak9lxiSF14DfhRC7hBDbdA9zB9bSZGZmctNNN+Hs7GztUK5V0T7znew7uPGDFMYeH853Xk17Xb4xsrOzGfPOEQpLytm1axdBQUF1FsdrqorzislJzOaWW27h6aef5uzZsxQUFNR4bH5pgX5X96zIWcRkxJBa0PSWioeHh1NUVMT+/aroc2XGJIW1wOvAEuDNSg/FREpKSsjLy9N3+7JFP/zwA/Pnz+eOO+7ApfvfKCwsqvtDTVhJSQmTJk0i5lIe/330BkJDQ7ntttvIysxqditWYmJiyErIwt7Rnq+++orRo0cDNVdN9XHxwdXBpcpr+aUFpBekWyRWUxo6dCguLi58++231g7Fphgzp5AmpXzX7JG0YLpbEjY3n1ChoLCACRMmEBwczIYNG+j1+MfX7GpuTqSUPP744+zevZtPZvVhRKi2wdCIESP4/o/vyc3NtXKExsvcuImc7VULBvgn5ZMSoC1XkZKSQkREBC8LgVeQF15e2keXLl3IyMy4pvCdn4sffi5+NTbeaWqcnZ0ZPny42q9QjTEjhaNCiNcqyl7coHuYPbIWJDMzE3t7e5tsWFJaVsqJEycoLy9n27ZteHt74+Ts1KyTwmuvvcbKlStZuHAh9w+9uv9SN5JrSk1acrZvp7DassuUAFf+7O9Lbm4uY8eO5dKlS3gFeWLveHWl0ciRI8nKyqqxG1tzEh4eTlxcHPHx8dYOxWYYM1K4vuKflWdAjVqSqhgnMysTT09PHByM+c9hRtVWGkkkheeOkJ9fzObNW+nWrRug/YbVXJPC22+/zfPPP8+9997L4sWLYe2d+vf8/f21fZubUFIAcA4JodNn6/TPF0XOorSwlK133klUVBRbtmxB8/GrVT4zcuRIsrbvICe7PpVtqiqMiSFh+n36PzuH1N3b29JGjRoFaHdxP/LII1aOxjYYs6N5eA0PlRBMJCEhgcKCQtvoP1BppRHA2bNniUosorzX+Cq3tpprUli5ciXz5s3j7rvvZs2aNTXu0vby8iI7O7vWSdimoLS4lL3/2sv+/ftZt24dd9111zXHDB8+HCGEUd3YauIxZkyVJOAcEoLHmDEGPmEd3bp1o2vXruoWUiW1fjUVQhis9CWlXGb6cFoeXVEuT08bWfVbsdJo/fr13Puve5k9ezYrnlpR5RBnJ2eKi4spKCioUkK7KVu7di0PPfQQo0ePZv369fpR22Zy2SnyIHIWAAXuBUgp+fnnn212DsiQ4uJifvj3D1z87SKffvrpNWWydTw8PGjtUXc3ttp4TZ6E1+RJjQnVYsLDw1mzZg1FRUXNcvd2fRkaKbhXPMKAh9GWvQ4AHkLbNEcxgb1796LRaGyq09revXuZNWsWw4YN47333rum6Jtu2ez58+etEZ7Jbdy4kfvvv5+/PfQ3fB7xYc7eOfrllq/YZRAlrq606unTk5xfc9izZ48VI26YkpISpkyZQuKviQx+fDCzZs0yeLy3tze5uXlcunTJQhFaR3h4OHl5eRw8eNDaodgEQ1VSX5ZSvgz4AjdIKZ+SUj4F9Edb8VRpJCkl33//vb7gmi24knuFcePG0b17d77++usam/3okoKh1pxNxebNm5k2bRpDhgyh65iuxGXFVXk/TDrxUrk3q8NXszp8NetGr6NncU++//57K0XcMCUlJUybNo3//e9/DHhoAD0ietT5Ga+Kbmy7d+82d3hWdcstt+Do6KhuIVUwZvVRR6ByiYtiIMgs0bQwsbGxJCcn67t7WVt+Qb6+81tkZGSt8xzNJSl8+OGHTJ48mUGDBrFjxw7s7e3p4d1DnwBWh69mtWzDRCo1DopazWcjMvhP6ClKPr4DVo/WPqJW134hKysrK+POO+9k06ZNvPHGG4SOM26g7+bmhkaj0ZfSbq7c3NwYOsZoXpwAACAASURBVHSoSgoVjEkKnwG/CiEWCSH+ibZk9lrzhtUy6L5tellrpBC1Wv9LrWjlbZQm/g5oV2IYaszu5OSIEKLJJgUpJa+88gqPPPIIo0ePZteuXcZXpj2+hUCHTEBeLXlx6XiVVVu2pKSkhD/++IPvvvuOVatW8fTTxne9FWh7LNTUja25GTVqFCdOnODChQt1H9zMGbP66N/ALCATyAJmSSlfM3dgzZmu/eaHsU50mrWMM5kl1gmkYrVRaWkp0dHR/HG5nIDwJ+jRo65bCwInZ6cmkxQ2/7W5SlmGoR8MZU3BGiIWRPDVV1/h6lq/vsN27a9jzBbBywkDYNYO7eS8DUpISOD3Y7+Tm5fH//73Px544IF6n8Pb25u0tDR+//13M0RoO3RVU5v7qMgYRi2Mr+jJbExfZsUIW48lcfJiNunp6bRv355uVmy9WerXk5tX53LkSBY7duwg6Pbbjfqcs7Mz506dM29wJrLzzE5iM2IJ9gwmNjaWlJQUWge3xr+9PxqNpt7nsxN2DBs2zKbnFU6cOMHIkSP5j5Mz1/XtS9+xdbcw8U/K1+8rAO3eAq+uXQDtL8v+/fubLV5rCw0NJTAwkMjISB588EFrh2NVVt4t1XJ5lOVwcv0C9pw4Qa9evawSQ2mZdoTw669pfPnll9xuZEIAcHF2aTIjBYDObp2JWxzHL4d+4fXXX+fP9n8a/+FLx7W32XR/btuH2267jaeeeoqkpCSs1Um7phIWHmPG8KOzE9OnT8fd3Z3r+/UzamXbn/19gTS8K72m21vQLzaWXbt22Xyvj8YQQhAeHs6mTZsoKSlp0JeF5sKsSUEIEQ68A9gDq6SUS6q9/yTwIFAKpAL3SykTzBmTrUhOTmbIkCFWSwg5OTmcjo7mypUrbNy4kfHjx9fr887Ozly6dMnm9ips/mszO89U6evEydST5JzOITE6kc2bNzNhwgRmRRpejqlXvXFQRTOhERptSZLvv/+e+2r4mCXoSljoNokVxsRw9tw57jp4gLCwML766ivKn3veqHP98bc2/PG3NlVqGumEH/udpUuXkpOTY9L4bU14eDirVq3il19+4aabbrJ2OFZjtqQghLAHPgBuBxKBI0KIbVLKU5UO+x0Ik1LmCyEeBv4DTDZXTLYiKyuLgoIC5syZY5XrZ2dnEx4ezpKQK4SGhuJXz4QAVfcq1D0HYR41JYCoy9oW32FttN1iky8lk/5XOvZx9hw6dIg+fa7e/4/NiK2SHGIzYunhXe1nCZulfVTTp7wcX19f9uzZw33DTfUT1Z+uhEVWVhbHbxmOX1YWkTcOoHtwMOXPPW+S8hJjxoxhyZIlbNu2DdvsHm4aI0aMwN7enh07drTopGDObusDgNNSyjMVXdu+BMZVPkBKuU9KqetwcZgWsv/hYvJFHBwcmDDB8u0rdQkhKipKmxB8G9ZExdlZu/PTmreQdHMFlYW1CeOlwS/x4S0f4rLNhV2zdtE1qiuHVx6ukhAiukRckwB6ePcgokuEUde2s7MjIiKCb775hnJp3ZU5J06c4MYbb+Tz8+cp79iRHj166Et0mKK8xODBg+nQoQMbNmwwRbg2y9PTk1tvvZVNmzY1u/Lo9WHO20cBQOX1XYnAQAPHPwA0+8LmKSkppKWm0T6gvcVvu6SkpDBq1Ciio6PZvHkzfpkfN/hctrJXQbevoLLjx48zYMAAjh8/ztNPP81rr712TbHBid0nMrF74xoFTZkyhXXr1pGRkYGvj2+jztUQEkhKTGREWBienp48tm0rg83wDdfOzo4pU6bw1ltvUTJhAhqH5nu/ferUqcyaNYvDhw8zePBga4djFeYcKYgaXqsx/QohpqEtp/FGLe/PFkJECSGiUlObXoenytauXYuUkvbtLNvrN333m5z9ZyhvXfcXl5b0567Mj6sUv6svR0cnNBqN1ZNCZeXl5SxbtoywsDBSUlLYsWMHb7zxhtmqz9522214e3uTkvL/7Z15fFTV+fC/z0yWyUpYEiEhIaCByI6AgG3FtkBRNqXs+amlpBQrLhWt9ef7UYvVyivv2w9U3tYUV5ZfMCACivVjKVSlyCImrAESliRsgQlkIzPJzJz3j5mMIQaYJDOZSTjffG7m3nPPPee5c+be557nnvM8xT4p/3oUFhayLyeHvPx8Ro8eTU5Ojk9NHjNnzsRmsxHo11+tZ9bapb7b8BsxefJkTCYTq1at8pGEgY8vlUIRkFhnuytwpn4mERkFPA9MVEo1GM5LKZWhlBqilBoSG9t6Y8Y6HA4yMjJo165do8fGN4eDBw9yZO0fSY2pYcCAAXTs4ArG7npp2lS6desWMEqhsLCQ0aNHs2DBAu69917279/Pffd5ZgpqKsHBwUyZMgWz2YzdYfdpXbUopVi1ahX9+vWjrLyMXj17snHjRncMaV8xcOBAevXqRfH5lleAnlLfMys03nwWHR3NhAkTWLNmDTU1fpo/5Gd8aT7aDaSISHfgNDADuMolo4gMAt4ExiqlAvfX5iW2bt1KXl4eIx/s0mJ17tixg3HjxvHxFDAmDCRyvvfi0SYnJ3PixAmvldcUlFIsW7aM559/HpvNxvLly/nlL3/5PSd+vmLmzJnY312J2Wwmzsd1nTx5kieffJINGzZw1113MSQ+njBTWIucq4gwc+ZMLr/1NtbqwAzF6i3PrGlpaWRlZfH555/7/MEiEPFZT0EpZQPmA58Bh4EPlFIHRWShiNTOpHkdiASyRCRbRDb6Sp5AICMjg/bt2+Pz3o7LfYX5/wyn5u8/49PpBoZ1CycyIvLGxzaC5ORkv/YUysvL2bt3L/Pnz2fo0KHk5OQwZ86cFlMI4IzzGxIS4lMTktVq5ZVXXqF37958/vnnLFq0iC+++IIwU8u+k5o5cyYAxcWBbUJqLvfeey/t27e/aU1IPp2noJTaDGyul/ZCnfVRvqw/kCguLmb9+vU8+uijnGkgeItX2b8Wa8Ee9p+oJDIygv79+2MMDmmWqaghkpOTOX/+fIvMVag7/NRms3H8xHEuBV3CUe1gzZo1TJ06tUWVQS1Go5G4uDjOnDlDaWkp7dq1a3aZdSellVy6xLFjx+haVcULw4eT9t57JCYm3qAE39CzZ08OREZSXHzeL/W3FCEhIUydOpWVK1dSUVFBZKR3H6YCHT2j2ces3lnAhuzTFBQU0GHKQvKShnLybBm9u0T7pD673c6J/DyKiir4vyU/ZfUbqwnxwY/60Nkyio19uWXmn5j2t/8QHh7OpIEJzBqW5PW6wDn8NLckl4jKCAoKCrDZbCR0TWDu5LlMG+DfYC5xsXEUFRWxYcMGHnqo+VPZyj7+mMqDBznucGAuKSEsLIwhHTsQmdDVbwqhlri4OPKPH+fYsWOkpKT4VRZfkpaWRkZGBhs2bCAtLc3f4rQoPn5k1WzIPs3B05edSqFDB8LDw+ntI19HFRUVPPDAAxQVFdG1a1fWr1/vk6ecSQMT6N0l2j1XwWKxcOhsGRuyT3u9LnD2DM6ePculI5fYMncLPXb3YP2M9Xz56Jc8OOBBn9TZGKKiozCZTF4Zx3/kyBEOHTrE7osXeajgFOce/Q2T9uUQ2TcwnO7FxjnfnGRmZvpZEt/ywx/+kMTExJvShKR7Ci2Ao6SQyx++xFcHDtC9e3ef1FFUVMTEiRPJycnhrwv7kRCfAEajT+qaNSyJWcOSOH36NF1/dy//PXwJ2zsM9WodWUez+CT/E4qLizl16hQqVmEymdi2bRsjR470al3NRRDi4uL4fPnnXLx4kU6dPJ+zUGsqulJVRUFBAefPnSPVZCKsWzeOb9oYOGFaXZhCQ+kbEUH4qtWcPHoMwelIrzih5UbTtQQGg4FZs2axePFiiouLiYvz9TCCwEErBR9TXFxMSUkJr776qm8Uwp53uPzlck4dOsiSQQ56pw2lY/UZaAE3bfHx8aSkpLBhwwY6TWu6UqjvrsJut/PtRaer5srcSiIiIuge3Z0Hf/IgI3sFlkKoJS4uDrs9l3Xr1nnsvkQpReGqVdjy8thfUYHBICR07Ur7pCRSHnjApwrBIxcfDRA9fjxnzpzhyrFjVFZWEBkRSXFCOIcHd+JnPpPWP6SlpbFo0SI++OAD5s+f729xWgxtPvIhxcXF5OUdIzo62ic/KqUUhZ/+Gc7vJzg4mMF3DHbOQWjm/ANPERFmzZrF1q1bqa6uvvEB16DWXUV1dTUnTpzg66+/pjK3kvCvwln6g6Xk/HcOG9M2MrVX82Yg+5KIiAhSU1M9MiHV1NSQlZXFiBEj+DY7m1yLhfxfPMyYnTv56X+2k5L5Pz4Net8cFx/tp08jZU0m6WfPsKJnT7qteJ/Mx/qQc5dv50n4g379+tGvX7+bzoSkewo+5IknnsBmcE76MXrDlLPnHXeEL7vdzpGjR4g3mCmytef2RQcI9zR6mBeZNWsWf/jDHyguLqZr18a7rrLb7ZSUlFBztoY1z67B4XDwwAMP8MwzzzB8+HAfSOwbBOc4/pdeesnpTruByHW5S5dyZs0HnD93juqaGp4KC+P2mBgi+vVj1sKFN6yjdrZu3e2mOLtrrouPTp06MXr0aDIzM3nllVeaXE5rIC0tjd///vfk5+dz6623+lucFkErBR+xceNGMjMzGfbcJO/NXnZFSquI6s7hw4eprLxCZI/u9LnvKcQPCgGcwxSHDBnC+eLzjVIKJ06cYMWKFSxfvpygtCCCg4N5+umnSU9P57bbbvOhxL5j+vTpvPjii2RmZrJgwQLA6YDwo48+4q233mLOyVOkhoYS3S6aLp270KFjBwTxaMZtQ3m84eyuqfzqV79i8uTJvPXWW84o7q2I+sr1WkSPH8/MmTN57rnnWLp0KUuWLGkB6fyPVgo+oLS0lEceeYR+/fqRlOS9K0ahKKqJ4bYFXxMTE8PKletJakRgHF8xa9YsFu2uoKrqynXznT59mqysLDIzM9m5cycAY8aMIbh3MJ06deK1e1+77vGBTq9evfjRj37En/70J6Kiovjkk0/4xz/+QXV1NSkpKfTo0Z0Ot3Rm2JrGj9zx1mxdb3H//fdz99138/zzzzP2nbE+8y/lbTxVorU+k7pNn8bcuXNZtmwZ6enpV3nabau0jpZsJazeWcCH3xRw4OAB7CMfo+sdd3D4XLlX5iScOHGCy9nZlJaWMm7cON58803fz4z2kOnTp7No9/ucr+cXJ+toFusOrsNsNnPx4kVKS0sBiJwQyU/n/JTY2FhMJhNHSo4QK4FxLk2lqKSKT//+d4xGI2azmV//+tckJCTw6KOPMnXqVIYPH07BQw/7W0yvISIsWbKEO+64g1OnTrUa04qnyrVuT+KVV14hKyuLxx57jK1bt/plkmRLopWCF1m7+yTfnriA5fJlUlNTiYqKoncUzZqT4HA4ePvtt/ntb3/L5mlCamoq655dF1A/zPj4eGJiYiguLsZisbBjxw4++eQTPo3+FBWrsBRYCA8PJzk5mdjY2O+Z0xoTx8Df1A4hdTgclJWVUXKphJLzZ1hzoYSs0q0kJSUxYMAAcnJyWL9+PUOHeneobiAxcOBA0tPT2XJ6C126tJw/r5amY8eOvPrqq8ybN481a9YwY8YMf4vkU7RS8BKFhYVkZ2djtVjImNGHCRMmNK/APe9w+avl5Ofnc1t5OV/OiaF/nGC4pTMEkEKw2Wzs3bsXu9qAoWsBfRYuRSmFhAgRXSLobOjMm3Pf9Nn8jJaiqqqKXbt2wbI3iLxwkUNVV3A4FCJCn7BQHk/uwEtT+nN7fCR2h51du2L456L/YvAHh90Bb9oif/zjH9m6dCv5+fn+FsWnpKenk5GRwYIFCxg/fnybdn2hlYIXyM3NZcyYMdhHPkb//v2ZMGFss8rLy8vDkvkSXYMuUV0dxO2ptxN3SxyCtMhQ07rUn0Ngs9koKyujrKyM0rJSysvKsdvtRAyPAAxQFE5yt2Tat2+P0Wjkvh73tTqF4HA4OHbsGDt37qR83Yd0OXaMispKlFKkhoZyKjiY7J//nHvuuYe7776byw9PxXLiNKZ/l1NAOQBdrB25r6qanWN+RrzrKdoboTEDjbi4OLold+N4/nE2b97cZr2KGo1Gli1bxogRI3j55ZdZtGiRv0XyGVopNJFan0YXL17k6NGj8OMnaJd0e9McormGmlZVVVF0uogzZ84w4BYDpaYkBvx5V4vFXmhIAWSbswGIKImgoqKCK1e+e5kcERnBLbfcQkxMDOctBqrOdsf81w/ILsr2zhDcFqC8vJxDhw6Rk5NDdnY2OTk57Nu3j4qKCgBWJHenm8lEWWIi7aKjadeuHT3vv5/JdezSjpnp4HJgV4uc209QTQUnjh8ntlMngoOD/TpayJd0TejK2TNneeqppxg1ahQhISH+FsknDB8+nNmzZ/PnP/+Z2bNnk9rGFHwtWik0kff/fYgjxZVUnTlGeHg4ffv2JSwsrNHvDxwOB5f/nUHIpSPsKbIiInTu3Jmw5GRiBqdBCygEi8XC0aNHeTf7Xc7YzhBUEkRlZSUWiwWAyzsuE5EfwaBBgxg2bBgjRozgzjvvvKoLPf3NHVwwXODguaVs3bqVUaMCxwGuzWajoKCA/Px88vLyyM3N5fDhwxw+fJiioiJ3vgdv6cwzsbFE9u1LZGQk0dHRGAsKMaWmMnjF+9csv8GXl++Mo/JKJQOf3M7DP/kxy5cv99Xp+R0R4dZbb2XTkU288cYbPPXUU/4WyWe89tprfPjhhzz++ON89tlnAfVuz1topdBI9u/fz3PPPUd29A8JDQ3hlZ/G8vDDDzdqSJ7a8w5Xdr6H2Wzm7Nmz9GpXzb4SI9uSn2Hu3LnEx3s/VKfNZqOwsJD8/HyOHTvGkSNH+Mb2DeZYs/vmb0oyYS20ErYxjP59+9O/f38GDx7MoN8O8sj3y4WaUOL/63/z9OYieuXvAKB3fDQvTujj9fOB71762h0OqqursVotWCxWrBYLFqvzs8pShcVidQdiDweGGo2MDA8jvHMXwnvcSkR4OBGRkaj9+6GmhvDb6nj/bMbTfUR4BE8++SSLFy9m0qRJzX/PFMB07NiRsWPH8uyzz2K321mwYEGbfJcSFxfHyy+/zOOPP84TTzzBokWLWjzWuq/xqVIQkbHAEsAILFdKvVZvfyjwPjAYMAPTlVInfSlTY1FKsXj916zbcwqz2Ux5eTlBHe4hOqEnA5I6MmfODzwqx2KxsH37djZt2sQMywpSY2o4fs5OdHQ01vbJ3PnAPO4alt4sOc1mMydPnmTtsbXsqdiDpfamWGXBYrG4b4wAhggDYT3DMGKkW1U3wsPDiYiIYPLIycx4tfGjK2p7SEfKY7lw4QIpKSlNvikopbhy5Yp7KOuFCxdw/OMzorK/pbq6hprqaqqrq7nN6owAtuvK9+dHhAQHE2oKJSoyitjYOMLCwggLMxFmCiMkNLTBAOIMHUr0+PFenQ/w4osvsnnzZiZOnMjs2bNZvHgxHTp08Fr5gcTq1atJT0/nd7/7HVu2bOG9997zeZhQf/Cb3/yGvLw8li5dypYtW1i5ciWDBg3yt1hew2dKQUSMwDJgNM54zbtFZKNS6lCdbHOAS0qp20RkBrAImO4rmTzh3S/zyNp1nIqKSioqKjCbzRCXAkQRLGa6d+9OfHw8QUFB1zQVlZaWkpuby969e4nIzaK3/SCVlZUYlWKywUDf+GDKwm7ltvmbPJoFrJTi4sWLrNy3km3nt2G1WrFWW7FarFitViwWC1arFYfDAUBEagSEgLXAislkIjIykk6xnQgzhblujmGEhjrdXt/X475muTyopdZz6j//Wcno0bPJze7LI488wl2Jd7F9+3bKy8spKysj7MuvaH9gP3abHZvN5l5qbDZsthrnek0NDoe6qvw7XWa0PZYqgoNDCA0J4Xh4GKeSkrg0dCjx8fEkJiaSlJREYmJiwDy9RUZGsmfPHhYuXMjrr7/O5s2b+ctf/sKUKVPanOmhffv2rF27loyMDJ588kkGDBjA+++/z5gxY/wtmlcxGo0sWbKEcePGMXv2bIYNG8bChQt55plnWs27tOshdZ8evVqwyAjgJaXUz1zbzwEopf5UJ89nrjw7RCQIOAfEqusINWTIELVnz55mybZv3z4KCwv5OPcy35qNLtODlStXrqBinS4WLAX7MRgMxMTE0KF9Ox4Y2JVpdyZRXl6O2Wwm8shabin+AluNjeqaaiwWi3ux2Wzuuu5Jdurdo7YuhJgiCQqLYkNQFdtiwrCb2te7MdZgq3HeFGtqaqiuqaam2rmulHLe7HF6DgVnhKjQ0FBCQkMwBZsIDgomJDiEkKAQRiWO4v6e939Xluvpuu4S+sUXRHzzDcqhcDgcOJQDh2tdORzOtAYWu8OOw+7A7nDgcNix2x047Hbntt3O9X5RtTf3XVeuYDQaCQ4KIqh2CQ4iKCiY4OBggoOdri+Ca7dDQoiZMJ7ODz3UOm6m74xzfs7+xJ307bffkp6ezt69e5k4cSLTpk2jV69e9OrViyg/uSnxBrXeVt8Z+4477cCBA8yYMYODBw8ydepURowYwaBBgxg4cGDAuQNviFMPPtTgaLH6PcmSkhLmzZtHVlYWw4cPZ9KkSfTv7zS9JiQkBNRvVUS+UUoNuWE+HyqFKcBYpVS6a/tBYJhSan6dPAdceYpc2/muPBevVW5TlcIHs4YSWey8mZaet+KweXre6rqb192lGlhVIEYBpVCO2h2qXiZp4EDXsa4/lEIphcP12RTq3qAbQgCD0YDBUH8xYjQYXPuMGI0GjAYjBuPV61cqKzGXmLFaqwkJcSoso9FIdnQ0B+K7YDKZGrxoAulCairPdfmKpJBSCqqvHo12yhrN818aycnJwW63u9PDw8OJiorCaDRiMBgQEfdnoH9HoQ+FYuhswHHO8b19NTU17vOsKqji3Opz7oeA2nOoPcdAOqf7wyP4Wb1BHoNdvetvXGbL73Bdi47vX4fi/ge5Viuvm83NkqtPnz5kZ2c36VhPlYIv3yk01ML1vzVP8iAic4G5gFd8CdmrHThqPLuRNvgzvV6i1Etx/atNFgHsYJQgggxBiEEwiPNmKwbXxYFcfaGIYJCr09yLQRAxfLe/7vZV6/XSDIJdDKi7f8TwcePcvY7Q0FBMJhOhoaFe8WfjcDhYt24dq1atoqKigqqqKqqqqrAcOoT1excXTVZwgcaKFAsTuwOUXpV+tqSSS5ciSExMdPfiapfLly8Dzu+g9ntoDd9Huy/a0W7YtYdii4hzQqNL0dntdux2e0Cf28qKClbWS5varh3johs35PyqZz6HusqK0BROn/ZNdMO63JTmI41Go7nZ8LSn4MsxY7uBFBHpLiIhwAxgY708G4FaL2FTgH9dTyFoNBqNxrf4zHyklLKJyHzgM5xDUt9WSh0UkYXAHqXURuAtYIWI5AElOBWHRqPRaPyET+cpKKU2A5vrpb1QZ90CBG6MRY1Go7nJaHtTDjUajUbTZLRS0Gg0Go0brRQ0Go1G40YrBY1Go9G40UpBo9FoNG58NnnNV4jIBeBUveR21J862nBaJ+CaLjRakIZk80d5nh7nSb4b5bnW/sakB0L7tba28zTv9fI0dl+gth0ERvu1ZNvV3d9NKRV7wxprp9S35gXI8DBtj79lvZZs/ijP0+M8yXejPNfa35j0QGi/1tZ23mi/xu4L1LYLlPZrybZrioxtxXy0ycO0QMHbsjW1PE+P8yTfjfJca39j0/1Na2s7T/NeL09j9wVq20FgtF9Ltl1j62t95qPmICJ7lAe+PzSBiW6/1otuu9ZDW+kpeEqGvwXQNAvdfq0X3XathJuqp6DRaDSa63Oz9RQ0Go1Gcx20UtBoNBqNG60UNBqNRuPmplYKItJDRN4SkbX+lkXTeETkfhH5u4hsEJEx/pZH4zkicruI/E1E1orII/6WR/MdbU4piMjbIlIsIgfqpY8VkSMikicivwdQSh1XSs3xj6Sahmhk+32klPoV8Atguh/E1dShkW13WCk1D5gG6KGqAUSbUwrAu8DYugkiYgSWAfcCvYGZItK75UXTeMC7NL79/pdrv8a/vEsj2k5EJgJfAVtaVkzN9WhzSkEp9QXO0J51uRPIc/UMqoFMYFKLC6e5IY1pP3GyCPhUKbW3pWXVXE1jrz2l1Eal1F1AWstKqrkebU4pXIMEoLDOdhGQICIdReRvwCARec4/omk8oMH2Ax4DRgFTRGSePwTT3JBrXXv3iMhSEXmTeiF7Nf7FpzGaAwhpIE0ppcyAvpkEPtdqv6XA0pYWRtMortV224BtLSuKxhNulp5CEZBYZ7srcMZPsmgaj26/1otuu1bGzaIUdgMpItJdREKAGcBGP8uk8Rzdfq0X3XatjDanFETkf4AdQC8RKRKROUopGzAf+Aw4DHyglDroTzk1DaPbr/Wi265toB3iaTQajcZNm+spaDQajabpaKWg0Wg0GjdaKWg0Go3GjVYKGo1Go3GjlYJGo9Fo3GiloNFoNBo3Wilo2gyu+Aq962wvFJFR18l/j4h83AJy/UJE3nCtzxORh7xQpl1EskUkvvkSgohsFZEKEdFurG9ybhbfR5o2jogEAfcDHwOHAJRSL/iwPqNSyt7Y45RSf/OSCFVKqYFeKgul1I9FZJu3ytO0XnRPQRMQiEiyiOSKyHsiss8VkSvcte8FEdktIgdEJENExJW+TUReFZF/A88CE4HXXU/Qt4rIuyIyxZV3qIj8R0RyRGSXiETVqz/CFSRmt4h8KyLfc63u6llsFZHVwH5X2kci8o2IHBSRuXXyzhaRoy7ZflAn/SURebqO/ENc651E5KRrvY9LxmzXd5HiwfdXUWd9ioi861p/V0T+6pL7uIiMdJ3n4do8Gk1ddE9BE0j0AuYopbaLyNvAb4DFwBtKqYUAIrICGA9sch0To5Qa6dqXAnyslFrr2sb1GQKsAaYrpXaL8rFXLAAAAqlJREFUSDRQVa/u54F/KaV+KSIxwC4R+adSqrJevjuBvkqpE67tXyqlSkQkDNgtIuuAEOAPwGCgFNgKfNuI72EesEQptcolu7ERxzZEe+AnOJXmJpxKKt0l70ClVHYzy9e0IXRPQRNIFCqltrvWVwI/dK3/WER2ish+nDe3PnWOWeNBub2As0qp3QBKqTKXT566jAF+LyLZOF06m4CkBsraVUchADwuIjnA1zi9gaYAw4BtSqkLrsAynshYlx3Af4vIs0A3pVR9BdZYNimnP5v9wHml1H6llAM4CCQ3s2xNG0P3FDSBRH1HXEpETMD/A4YopQpF5CWcN+xa6j/JN4Q0UHZDeX6ulDpyg3zu+kTkHpxBfkYopa64bPK1snniVMzGdw9m7nNSSq0WkZ3AOOAzEUlXSv3LA/lrCa63z+r6dNRZr93W9wDNVeiegiaQSBKREa71mTjj99beLC+KSCQw5TrHlwNRDaTnAvEiMhRARKJcL6br8hnwWJ33FYM8kLcdcMmlEFKB4a70ncA94ozsFwxMvcbxJ3GamKDOeYlID+C4K4jQRqC/B7KE1xl5dQ/NNzlpblK0UtAEEoeBh0VkH9AB+KtS6jLwd5ymj49w+ue/FpnAM64XxbfWJrpMONOBv7hMPZ9zdW8D4GWcT9j7ROSAa/tG/AMIcsn7Mk4TEkqps8BLOM1A/wSuFT96MfCIiPwH6FQnfTpwwGXKSgXe90CWKuAFEfkG53uMH4vIXR4cp9FchXadrQkIRCQZ50vivn4WpVUgIhVKqchrbTexzG3A00qpPc2VT9N60T0FjaZ1UubtyWtAD6DGG+VpWi+6p6DRaDQaN7qnoNFoNBo3WiloNBqNxo1WChqNRqNxo5WCRqPRaNxopaDRaDQaN1opaDQajcbN/weWJv0NBpz8bAAAAABJRU5ErkJggg==\n" + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "with np.errstate(all='raise'):\n", + " setup = SetupA()\n", + " states, _ = run(setup)\n", + "\n", + " x_min = min([state.min('x') for state in states.values()])\n", + " x_max = max([state.max('x') for state in states.values()])\n", + "\n", + "with np.errstate(invalid='ignore'):\n", + " plotter = Plotter(setup, (x_min, x_max))\n", + " for step, state in states.items():\n", + " plotter.plot(state, step * setup.dt)\n", + " plotter.show()\n" + ], + "metadata": { + "collapsed": false, + "pycharm": { + "name": "#%%\n", + "is_executing": false + } + } + }, + { + "cell_type": "code", + "execution_count": null, + "outputs": [], + "source": [ + "# TODO python -O\n", + "def test_timing():\n", + " setup = SetupA()\n", + " setup.steps = [100, 3600]\n", + "\n", + " nsds = [2 ** n for n in range(12, 15)]\n", + " times = []\n", + " for sd in nsds:\n", + " setup.n_sd = sd\n", + " _, stats = run(setup)\n", + " times.append(stats.times[-1])\n", + "\n", + " from matplotlib import pyplot as plt\n", + " plt.plot(nsds, times)\n", + " plt.show()\n" + ], + "metadata": { + "collapsed": false, + "pycharm": { + "name": "#%%\n" + } + } + } + ], + "metadata": { + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 2 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython2", + "version": "2.7.6" + }, + "kernelspec": { + "name": "python3", + "language": "python", + "display_name": "Python 3" + }, + "pycharm": { + "stem_cell": { + "cell_type": "raw", + "source": [], + "metadata": { + "collapsed": false + } + } + } + }, + "nbformat": 4, + "nbformat_minor": 0 +} diff --git a/requirements-dev.txt b/requirements-dev.txt index 9cf8484615..fc522cbbe3 100644 --- a/requirements-dev.txt +++ b/requirements-dev.txt @@ -1,3 +1,4 @@ pytest pytest-cov codecov +matplotlib diff --git a/requirements-examples.txt b/requirements-examples.txt new file mode 100644 index 0000000000..558c1d0b66 --- /dev/null +++ b/requirements-examples.txt @@ -0,0 +1,3 @@ +nbconvert +jupyter_client +ipykernel \ No newline at end of file diff --git a/requirements.txt b/requirements.txt index f9cec45772..4559b7ec25 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,4 +1,3 @@ numpy scipy -matplotlib mpmath diff --git a/tests/debug.py b/tests/debug.py index 61c020699a..c23ca25698 100644 --- a/tests/debug.py +++ b/tests/debug.py @@ -3,13 +3,13 @@ def plot(state): - state.sort_by_m() + state.sort_by('x') - pyplot.plot(state.x, state.n) + pyplot.plot(state['x'], state['n']) pyplot.grid() pyplot.show() -def plot_mn(m, n): - state = State(m, n) +def plot_mn(x, n): + state = State({'x': x, 'n': n}) plot(state) diff --git a/tests/test_colliders.py b/tests/test_colliders.py index 1a191ddfb5..f3f9e7190d 100644 --- a/tests/test_colliders.py +++ b/tests/test_colliders.py @@ -22,25 +22,25 @@ def __call__(self, m1, m2): class TestSDM: @pytest.mark.parametrize("x, n", [ - pytest.param(np.array([1, 1]), np.array([1, 1])), - pytest.param(np.array([1, 1]), np.array([5, 1])), - pytest.param(np.array([1, 1]), np.array([5, 3])), - pytest.param(np.array([4, 2]), np.array([1, 1])), + pytest.param(np.array([1., 1.]), np.array([1, 1])), + pytest.param(np.array([1., 1.]), np.array([5, 1])), + pytest.param(np.array([1., 1.]), np.array([5, 3])), + pytest.param(np.array([4., 2.]), np.array([1, 1])), ]) def test_single_collision(self, x, n): # Arrange - sut = SDM(StubKernel(), dt=0, dv=0) - sut.probability = lambda m1, m2, n_sd: 1 - state = State(x, n) + sut = SDM(StubKernel(), dt=0, dv=0, n_sd=len(n)) + sut.probability = lambda sd1n, sd2n, sd1x, sd2x, n_sd: 1 + state = State(n=n, extensive={'x': x}, intensive={}, segment_num=1) # Act sut(state) # Assert - assert np.sum(state.n * state.x) == np.sum(n * x) - assert np.sum(state.n) == np.sum(n) - np.amin(n) - if np.amin(n) > 0: assert np.amax(state.x) == np.sum(x) - assert np.amax(state.n) == max(np.amax(n) - np.amin(n), np.amin(n)) + assert np.sum(state['n'] * state['x']) == np.sum(n * x) + assert np.sum(state['n']) == np.sum(n) - np.amin(n) + if np.amin(n) > 0: assert np.amax(state['x']) == np.sum(x) + assert np.amax(state['n']) == max(np.amax(n) - np.amin(n), np.amin(n)) @pytest.mark.parametrize("n_in, n_out", [ pytest.param(1, np.array([1, 0])), @@ -49,56 +49,56 @@ def test_single_collision(self, x, n): ]) def test_single_collision_same_n(self, n_in, n_out): # Arrange - sut = SDM(StubKernel(), dt=0, dv=0) - sut.probability = lambda m1, m2, n_sd: 1 - state = State(x=np.full(2, 1), n=np.full(2, n_in)) + sut = SDM(StubKernel(), dt=0, dv=0, n_sd=2) + sut.probability = lambda sd1n, sd2n, sd1x, sd2x, n_sd: 1 + state = State(n=np.full(2, n_in), extensive={'x': np.full(2, 1.)}, intensive={}, segment_num=1) # Act sut(state) # Assert - np.testing.assert_array_equal(state.n, n_out) + np.testing.assert_array_equal(sorted(state._n), sorted(n_out)) @pytest.mark.parametrize("x, n, p", [ - pytest.param(np.array([1, 1]), np.array([1, 1]), 2), - pytest.param(np.array([1, 1]), np.array([5, 1]), 4), - pytest.param(np.array([1, 1]), np.array([5, 3]), 5), - pytest.param(np.array([4, 2]), np.array([1, 1]), 7), + pytest.param(np.array([1., 1]), np.array([1, 1]), 2), + pytest.param(np.array([1., 1]), np.array([5, 1]), 4), + pytest.param(np.array([1., 1]), np.array([5, 3]), 5), + pytest.param(np.array([4., 2]), np.array([1, 1]), 7), ]) def test_multi_collision(self, x, n, p): # Arrange - sut = SDM(StubKernel(), dt=0, dv=0) - sut.probability = lambda m1, m2, n_sd: p - state = State(x, n) + sut = SDM(StubKernel(), dt=0, dv=0, n_sd=len(n)) + sut.probability = lambda sd1n, sd2n, sd1x, sd2x, n_sd: p + state = State(n=n, extensive={'x': x}, intensive={}, segment_num=1) # Act sut(state) # Assert gamma = min(p, max(n[0] // n[1], n[1] // n[1])) - assert np.amin(state.n) >= 0 - assert np.sum(state.n * state.x) == np.sum(n * x) - assert np.sum(state.n) == np.sum(n) - gamma * np.amin(n) - assert np.amax(state.x) == gamma * x[np.argmax(n)] + x[np.argmax(n) - 1] - assert np.amax(state.n) == max(np.amax(n) - gamma * np.amin(n), np.amin(n)) + assert np.amin(state['n']) >= 0 + assert np.sum(state['n'] * state['x']) == np.sum(n * x) + assert np.sum(state['n']) == np.sum(n) - gamma * np.amin(n) + assert np.amax(state['x']) == gamma * x[np.argmax(n)] + x[np.argmax(n) - 1] + assert np.amax(state['n']) == max(np.amax(n) - gamma * np.amin(n), np.amin(n)) @pytest.mark.parametrize("x, n, p", [ - pytest.param(np.array([1, 1, 1]), np.array([1, 1, 1]), 2), - pytest.param(np.array([1, 1, 1, 1, 1]), np.array([5, 1, 2, 1, 1]), 1), - pytest.param(np.array([1, 1, 1, 1, 1]), np.array([5, 1, 2, 1, 1]), 6), + pytest.param(np.array([1., 1, 1]), np.array([1, 1, 1]), 2), + pytest.param(np.array([1., 1, 1, 1, 1]), np.array([5, 1, 2, 1, 1]), 1), + pytest.param(np.array([1., 1, 1, 1, 1]), np.array([5, 1, 2, 1, 1]), 6), ]) def test_multi_droplet(self, x, n, p): # Arrange - sut = SDM(StubKernel(), dt=0, dv=0) - sut.probability = lambda x1, x2, n_sd: p - state = State(x, n) + sut = SDM(StubKernel(), dt=0, dv=0, n_sd=len(n)) + sut.probability = lambda sd1n, sd2n, sd1x, sd2x, n_sd: p + state = State(n=n, extensive={'x': x}, intensive={}, segment_num=1) # Act sut(state) # Assert - assert np.amin(state.n) >= 0 - assert np.sum(state.n * state.x) == np.sum(n * x) + assert np.amin(state['n']) >= 0 + assert np.sum(state['n'] * state['x']) == np.sum(n * x) # TODO integration test? def test_multi_step(self): @@ -106,21 +106,19 @@ def test_multi_step(self): n = np.random.randint(1, 64, size=256) x = np.random.uniform(size=256) - sut = SDM(StubKernel(), dt=0, dv=0) - sut.probability = lambda x1, x2, n_sd: 0.5 - state = State(x, n) + sut = SDM(StubKernel(), dt=0, dv=0, n_sd=len(n)) - from SDM.undertakers import Resize - undertaker = Resize() + sut.probability = lambda sd1n, sd2n, sd1x, sd2x, n_sd: 0.5 + state = State(n=n, extensive={'x': x}, intensive={}, segment_num=1) # Act for _ in range(32): sut(state) - undertaker(state) + # undertaker(state) # Assert - assert np.amin(state.n) >= 0 - actual = np.sum(state.n * state.x) + assert np.amin(state['n']) >= 0 + actual = np.sum(state['n'] * state['x']) desired = np.sum(n * x) np.testing.assert_almost_equal(actual=actual, desired=desired) @@ -130,13 +128,12 @@ def test_probability(self): dt = 666 dv = 9 n_sd = 64 - sut = SDM(StubKernel(kernel_value), dt, dv) + sut = SDM(StubKernel(kernel_value), dt, dv, n_sd) # Act - actual = sut.probability((0, 1), (0, 1), n_sd) # TODO dependency state [] + actual = sut.probability(1, 1, 0, 0, n_sd) # TODO dependency state [] # Assert desired = dt/dv * kernel_value * n_sd * (n_sd - 1) / 2 / (n_sd//2) assert actual == desired - #TODO diff --git a/tests/test_kernels.py b/tests/test_kernels.py index 84b7f13d43..9f3b817c5e 100644 --- a/tests/test_kernels.py +++ b/tests/test_kernels.py @@ -11,18 +11,18 @@ class TestGolovin: - def test_analytic_solution(self): - b = 1.5e3 - x_0 = 4/3 * np.pi * 30.531e-6**3 - N_0 = 2**23 - - sut = Golovin(b) - - from matplotlib import pyplot - x = np.linspace(8e-25, 500e-15, 100) - print(sut.analytic_solution(x=1e-10, t=0.01, x_0=x_0, N_0=N_0)) - pyplot.plot(x, sut.analytic_solution(x=x, t=0.01, x_0=x_0, N_0=N_0)) - pyplot.show() + # TODO optional + # def test_analytic_solution(self): + # b = 1.5e3 + # x_0 = 4/3 * np.pi * 30.531e-6**3 + # N_0 = 2**23 + # + # sut = Golovin(b) + # + # from matplotlib import pyplot + # x = np.linspace(8e-25, 500e-15, 100) + # pyplot.plot(x, sut.analytic_solution(x=x, t=0.01, x_0=x_0, N_0=N_0)) + # pyplot.show() @pytest.mark.parametrize("x", [ pytest.param(5e-10), pytest.param(np.full(10, 5e-10)) diff --git a/tests/test_spectra.py b/tests/test_spectra.py index 968b81ca2f..c68ac1ab6a 100644 --- a/tests/test_spectra.py +++ b/tests/test_spectra.py @@ -58,20 +58,21 @@ def test_size_distribution_n_part(self, scale): # Assert assert_approx_equal(np.sum(sd) * dm, n_part, 2) - - def test_plot(self): - from matplotlib import pyplot as plt - - norm_factor = 1e10 - scale = 1e-13 - sut = Exponential(norm_factor, scale) - - x = np.logspace(-25, -11, 100) - y = sut.size_distribution(x) - - plt.loglog(x, y) - plt.show() - + # TODO optional + # def test_plot(self): + # from matplotlib import pyplot as plt + # + # norm_factor = 1e10 + # scale = 1e-13 + # sut = Exponential(norm_factor, scale) + # + # x = np.logspace(-25, -11, 100) + # y = sut.size_distribution(x) + # + # plt.loglog(x, y) + # plt.show() + + # TODO @pytest.mark.xfail def test_underflow(self): np.seterr(all='raise') # TODO: use with construct diff --git a/tests/test_state.py b/tests/test_state.py index 0a18cdbb8f..c92a19b638 100644 --- a/tests/test_state.py +++ b/tests/test_state.py @@ -2,10 +2,71 @@ from SDM.discretisations import linear from SDM.spectra import Lognormal +import numpy as np +import pytest + class TestState: - # TODO: test copy() calls in ctor + @staticmethod + def check_contiguity(state, attr='n', i_SD=0): + item = state[attr] + assert item.flags['C_CONTIGUOUS'] + assert item.flags['F_CONTIGUOUS'] + + sd = state.get_SD(i_SD) + assert not sd.flags['C_CONTIGUOUS'] + assert not sd.flags['F_CONTIGUOUS'] + + @pytest.mark.xfail + def test_get_item_does_not_copy(self): + # Arrange + arr = np.ones(10) + sut = State({'n': arr, 'a': arr, 'b': arr}) + + # Act + item = sut['a'] + + # Assert + assert item.base.__array_interface__['data'] == sut.data.__array_interface__['data'] + + @pytest.mark.xfail + def test_get_sd_does_not_copy(self): + # Arrange + arr = np.ones(10) + sut = State({'n': arr, 'a': arr, 'b': arr}) + + # Act + item = sut.get_SD(5) + + # Assert + assert item.base.__array_interface__['data'] == sut.data.__array_interface__['data'] + + @pytest.mark.xfail + def test_contiguity(self): + # Arrange + arr = np.ones(10) + sut = State({'n': arr, 'a': arr, 'b': arr}) + + # Act & Assert + self.check_contiguity(sut, attr='a', i_SD=5) + + def test_reindex_works(self): + pass + + @pytest.mark.xfail + def test_reindex_maintains_contiguity(self): + # Arrange + arr = np.linspace(0, 10) + sut = State({'n': arr, 'a': arr, 'b': arr}) + idx = range(len(arr) - 1, -1, -1) + assert len(idx) == sut.data.shape[1] + + # Act + sut._reindex(idx) + + # Assert + self.check_contiguity(sut) def test_moment(self): # Arrange (parameters from Clark 1976) @@ -18,7 +79,8 @@ def test_moment(self): n_sd = 32 spectrum = Lognormal(n_part, mmean, d) - sut = State(*linear(n_sd, spectrum, (mmin, mmax))) + x, n = linear(n_sd, spectrum, (mmin, mmax)) + sut = State(n=n, extensive={'x': x}, intensive={}, segment_num=1) #debug.plot(sut) true_mean, true_var = spectrum.stats(moments='mv') @@ -35,3 +97,21 @@ def test_moment(self): true_mrsq = true_var + true_mean**2 assert abs(discr_mrsq - true_mrsq) / true_mrsq < .05e-1 + + @pytest.mark.parametrize("x, n", [ + pytest.param(np.array([1., 1, 1, 1]), np.array([1, 1, 1, 1])), + pytest.param(np.array([1., 2, 1, 1]), np.array([2, 0, 2, 0])), + pytest.param(np.array([1., 1, 4]), np.array([5, 0, 0])) + ]) + def test_housekeeping(self, x, n): + # Arrange + sut = State(n=n, extensive={'x': x}, intensive={}, segment_num=1) + + # Act + sut.housekeeping() + + # Assert + assert sut['x'].shape == sut['n'].shape + assert sut.SD_num == (n != 0).sum() + assert sut['n'].sum() == n.sum() + assert (sut['x'] * sut['n']).sum() == (x * n).sum() diff --git a/tests/test_undertakers.py b/tests/test_undertakers.py deleted file mode 100644 index b814a738b7..0000000000 --- a/tests/test_undertakers.py +++ /dev/null @@ -1,30 +0,0 @@ -""" -Created at 07.06.2019 - -@author: Piotr Bartman -@author: Sylwester Arabas -""" - -from SDM.undertakers import Resize -from SDM.state import State -import pytest -import numpy as np - - -class TestResize: - - @pytest.mark.parametrize("x, n", [ - pytest.param(np.array([1, 1, 1, 1]), np.array([1, 1, 1, 1])), - pytest.param(np.array([1, 2, 1, 1]), np.array([2, 0, 2, 0])), - pytest.param(np.array([1, 1, 4]), np.array([5, 0, 0])) - ]) - def test___call__(self, x, n): - sut = Resize() - state = State(x, n) - - sut(state) - - assert state.x.shape == state.n.shape - assert state.n.shape[0] == (n != 0).sum() - assert state.n.sum() == n.sum() - assert (state.x * state.n).sum() == (x * n).sum()