diff --git a/.gitignore b/.gitignore new file mode 100644 index 00000000..4757a9d5 --- /dev/null +++ b/.gitignore @@ -0,0 +1,8 @@ +.coverage +__pycache__ +*.egg-info +*.pyc +.ipynb_checkpoints +.mypy_cache +.envrc +docs/.build diff --git a/CHANGELOG.rst b/CHANGELOG.rst index bc3bf756..01c95423 100644 --- a/CHANGELOG.rst +++ b/CHANGELOG.rst @@ -1,12 +1,19 @@ Master Branch ============= +Version 4.0.1 (2020-10-26) +========================== + +Release! + ADDED: * Gaussian Mixture Model: `GaussianMixture`. * Tutorial for how to use `scikit-learn` mixture models to fit a model, and `chaospy` to generate quasi-random samples and orthogonal polynomials. CHANGED: * `chaospy.Trunc` updated to take both `lower` and `upper` at the same time. +REMOVED: + * `chaospy.SampleDist` removed in favor of `chaospy.GaussianKDE`. Version 4.0-beta3 (2020-10-22) ============================== diff --git a/README.rst b/README.rst index acb30b85..131b3fdd 100644 --- a/README.rst +++ b/README.rst @@ -141,8 +141,9 @@ Also a few shout-outs: | `orthopy`_ | Thanks to `Nico Schlömer`_ for providing the implementation | | `quadpy`_ | for several of the quadrature integration methods. | +--------------+--------------------------------------------------------------+ -| ``UQRF`` | Thanks to `Florian Künzner`_ for providing the | -| | implementation for `sample distribution`_. | +| ``UQRF`` | Thanks to `Florian Künzner`_ for providing the initial | +| | implementation of kernel density estimation and | +| | quantity-of-interest distribution. | +--------------+--------------------------------------------------------------+ .. _OpenTURNS: http://openturns.github.io/openturns/latest diff --git a/chaospy/descriptives/quantity_of_interest.py b/chaospy/descriptives/quantity_of_interest.py index d92d875d..e75098f0 100644 --- a/chaospy/descriptives/quantity_of_interest.py +++ b/chaospy/descriptives/quantity_of_interest.py @@ -3,9 +3,7 @@ from functools import reduce from operator import mul import numpy - -from .. import distributions -from ..external import SampleDist +import chaospy def QoI_Dist(poly, dist, sample=10000, **kws): @@ -28,7 +26,7 @@ def QoI_Dist(poly, dist, sample=10000, **kws): Number of samples used in estimation to construct the KDE. Returns: - (numpy.ndarray): + (Distribution): The constructed quantity of interest (QoI) distributions, where ``qoi_dists.shape==poly.shape``. @@ -37,44 +35,16 @@ def QoI_Dist(poly, dist, sample=10000, **kws): >>> x = chaospy.variable(1) >>> poly = chaospy.polynomial([x]) >>> qoi_dist = chaospy.QoI_Dist(poly, dist) - >>> values = qoi_dist[0].pdf([-0.75, 0., 0.75]) + >>> values = qoi_dist.pdf([-0.75, 0., 0.75]) >>> values.round(8) - array([0.29143037, 0.39931708, 0.29536329]) + array([0.29143989, 0.39939823, 0.29531414]) + """ shape = poly.shape poly = poly.flatten() dim = len(dist) - #sample from the inumpyut dist - samples = dist.sample(sample, **kws) - - qoi_dists = [] - for i in range(0, len(poly)): - #sample the polynomial solution - if dim == 1: - dataset = poly[i](samples) - else: - dataset = poly[i](*samples) - - lo = dataset.min() - up = dataset.max() - - #creates qoi_dist - qoi_dist = SampleDist(dataset, lo, up) - qoi_dists.append(qoi_dist) - - # reshape the qoi_dists to match the shape of the input poly - if shape: - def reshape(lst, shape): - if len(shape) == 1: - return lst - n = reduce(mul, shape[1:]) - return [reshape(lst[i*n:(i+1)*n], shape[1:]) for i in range(len(lst)//n)] - qoi_dists = reshape(qoi_dists, shape) - else: - qoi_dists = qoi_dists[0] - - if not shape: - qoi_dists = qoi_dists.item() - - return qoi_dists + #sample from the input dist + samples = numpy.atleast_2d(dist.sample(sample, **kws)) + qoi_dist = chaospy.GaussianKDE(poly(*samples)) + return qoi_dist diff --git a/chaospy/distributions/__init__.py b/chaospy/distributions/__init__.py index b87e2798..78b8c3a1 100644 --- a/chaospy/distributions/__init__.py +++ b/chaospy/distributions/__init__.py @@ -73,7 +73,7 @@ from .operators import * from .constructor import construct from .approximation import * -from .kernel import GaussianKDE, GaussianMixture +from .kernel import * from . import ( baseclass, sampler, approximation, diff --git a/chaospy/distributions/kernel/__init__.py b/chaospy/distributions/kernel/__init__.py index 41061f16..983a633c 100644 --- a/chaospy/distributions/kernel/__init__.py +++ b/chaospy/distributions/kernel/__init__.py @@ -1,3 +1,34 @@ -"""Kernel density estimation.""" +""" +In some cases a constructed distribution that are first and foremost data +driven. In such scenarios it make sense to make use of +`kernel density estimation`_ (KDE). In ``chaospy`` KDE can be accessed through +the :func:`GaussianKDE` constructor. + +Basic usage of the :func:`GaussianKDE` constructor involves just passing the +data as input argument:: + + >>> data = [3, 4, 5, 5] + >>> distribution = chaospy.GaussianKDE(data) + +This distribution can be used as any other distributions:: + + >>> distribution.cdf([3, 3.5, 4, 4.5, 5]).round(4) + array([0.1393, 0.2542, 0.3889, 0.5512, 0.7359]) + >>> distribution.mom(1).round(4) + 4.25 + >>> distribution.sample(4).round(4) + array([4.7784, 2.8769, 5.8109, 4.2995]) + +In addition multivariate distributions supported:: + + >>> data = [[1, 2, 2, 3], [5, 5, 4, 3]] + >>> distribution = chaospy.GaussianKDE(data) + >>> distribution.sample(4).round(4) + array([[2.081 , 3.0304, 3.0882, 0.4872], + [3.2878, 2.5473, 2.2699, 5.3412]]) + +.. _kernel density estimation: \ +https://en.wikipedia.org/wiki/Kernel_density_estimation +""" from .gaussian import GaussianKDE from .mixture import GaussianMixture diff --git a/chaospy/distributions/kernel/gaussian.py b/chaospy/distributions/kernel/gaussian.py index 979f9c27..ec829679 100644 --- a/chaospy/distributions/kernel/gaussian.py +++ b/chaospy/distributions/kernel/gaussian.py @@ -83,8 +83,12 @@ def _mom(self, k_loc, cache): def _lower(self, idx, dim, cache): """Lower bounds.""" + del dim + del cache return (self.samples[idx]-10*numpy.sqrt(self.h_mat[:, idx, idx]).T).min(-1) def _upper(self, idx, dim, cache): """Upper bounds.""" + del dim + del cache return (self.samples[idx]+10*numpy.sqrt(self.h_mat[:, idx, idx]).T).max(-1) diff --git a/chaospy/external/__init__.py b/chaospy/external/__init__.py index b5917e4b..5cfc09d4 100644 --- a/chaospy/external/__init__.py +++ b/chaospy/external/__init__.py @@ -7,4 +7,3 @@ """ from .openturns_ import openturns_dist, OpenTURNSDist from .scipy_stats import ScipyStatsDist -from .samples import sample_dist, SampleDist diff --git a/chaospy/external/samples.py b/chaospy/external/samples.py deleted file mode 100644 index 7375cc3a..00000000 --- a/chaospy/external/samples.py +++ /dev/null @@ -1,157 +0,0 @@ -""" -In some cases a constructed distribution that are first and foremost data -driven. In such scenarios it make sense to make use of -`kernel density estimation`_ (KDE). In ``chaospy`` KDE can be accessed through -the :func:`SampleDist` constructor. - -Basic usage of the :func:`SampleDist` constructor involves just passing the -data as input argument:: - - >>> data = [3, 4, 5, 5] - >>> distribution = chaospy.SampleDist(data) - -This distribution can be used as any other distributions:: - - >>> distribution.cdf([3, 3.5, 4, 4.5, 5]).round(4) - array([0. , 0.1932, 0.4279, 0.7043, 1. ]) - >>> distribution.mom(1).round(4) - 4.25 - >>> distribution.sample(4).round(4) - array([4.4131, 3.3111, 4.9139, 4.1042]) - -It also supports lower and upper bounds defining where the range is expected to -appear, which gives a slightly different distribution:: - - >>> distribution = chaospy.SampleDist(data, lo=2, up=6) - >>> distribution.cdf([3, 3.5, 4, 4.5, 5]).round(4) - array([0.1344, 0.2543, 0.4001, 0.5716, 0.7552]) - -In addition multivariate distributions supported:: - - >>> data = [[1, 2, 2, 3], [5, 5, 4, 3]] - >>> distribution = chaospy.SampleDist(data) - >>> distribution.sample(4).round(4) - array([[1.5286, 2.0468, 2.1125, 1.8947], - [4.402 , 4.1522, 4.4384, 4.5737]]) - -.. _kernel density estimation: \ -https://en.wikipedia.org/wiki/Kernel_density_estimation -""" -import numpy -from scipy.stats import gaussian_kde -import chaospy - -from chaospy.distributions import SimpleDistribution - - -class sample_dist(SimpleDistribution): - """A distribution that is based on a kernel density estimator (KDE).""" - - def __init__(self, samples, lo, up): - samples = numpy.asarray(samples) - self.samples = samples - self.kernel = gaussian_kde(samples, bw_method="scott") - self.flo = self.kernel.integrate_box_1d(0, lo) - self.fup = self.kernel.integrate_box_1d(0, up) - self.unbound = numpy.all(lo == samples.min()) - self.unbound &= numpy.all(up == samples.max()) - super(sample_dist, self).__init__( - parameters=dict(lo=lo, up=up), - repr_args=[repr(samples), lo, up], - ) - - def _cdf(self, xloc, lo, up): - cdf_vals = numpy.array([self.kernel.integrate_box_1d(0, x) - for x in xloc]) - return (cdf_vals-self.flo)/(self.fup-self.flo) - - def _pdf(self, x, lo, up): - return self.kernel(x) - - def _lower(self, lo, up): - return lo - - def _upper(self, lo, up): - return up - - def _mom(self, k, lo, up): - if self.unbound: - return numpy.prod(numpy.mean(self.samples.T**k, -1)) - raise chaospy.StochasticallyDependentError("component lack support") - - -def SampleDist(samples, lo=None, up=None, threshold=1e-5): - """ - Distribution based on samples. - - Estimates a distribution from the given samples by constructing a kernel - density estimator (KDE). - - Args: - samples (numpy.ndarray): - Sample values to construction of the KDE. Either shape - ``(N,)`` or ``(D, N)``, where ``N`` are the number of - samples, and ``D`` is the number of dimension in the - distribution. - lo (float): - Location of lower bound. - up (float): - Location of upper bound. - threshold (float): - Threshold for how low the correlation between two - columns should be before defining them as - stochastically independent. - - Example: - >>> distribution = chaospy.SampleDist([0, 1, 1, 1, 2]) - >>> distribution - sample_dist(array([0, 1, 1, 1, 2]), 0, 2) - >>> q = numpy.linspace(0, 1, 5) - >>> distribution.inv(q).round(4) - array([0. , 0.6016, 1. , 1.3984, 2. ]) - >>> distribution.fwd(distribution.inv(q)).round(4) - array([0. , 0.25, 0.5 , 0.75, 1. ]) - >>> distribution.pdf(distribution.inv(q)).round(4) - array([0.2254, 0.4272, 0.5135, 0.4272, 0.2254]) - >>> distribution.sample(4).round(4) - array([0.3662, 0.6073, 0.9156, 1.0883]) - >>> distribution.mom(1).round(4) - 1.0 - - """ - samples = numpy.atleast_2d(samples) - assert samples.ndim == 2, "samples have too many dimensions provided" - - if lo is None: - lo = samples.min(axis=-1) - else: - lo = numpy.broadcast_to(lo, len(samples)) - if up is None: - up = samples.max(axis=-1) - else: - up = numpy.broadcast_to(up, len(samples)) - - # construct vector of marginals - distributions = [] - for samples_, lo_, up_ in zip(samples, lo, up): - #construct the kernel density estimator - try: - dist = sample_dist(samples_, lo_, up_) - #raised by gaussian_kde if dataset is singular matrix - except numpy.linalg.LinAlgError: - dist = chaospy.Uniform(lower=-numpy.inf, upper=numpy.inf) - distributions.append(dist) - - if len(samples) == 1: - distributions = distributions[0] - - else: - distributions = chaospy.J(*distributions) - - # Attach dependencies to data. - correlation = numpy.corrcoef(samples) - correlation[numpy.abs(correlation) <= threshold] = 0 - if numpy.any(correlation != numpy.diag(numpy.diag(correlation))): - distributions = chaospy.Nataf(distributions, correlation) - - return distributions diff --git a/docs/tutorials/advanced/gaussian_mixture_model.ipynb b/docs/tutorials/advanced/gaussian_mixture_model.ipynb index d4e674f9..7c785019 100644 --- a/docs/tutorials/advanced/gaussian_mixture_model.ipynb +++ b/docs/tutorials/advanced/gaussian_mixture_model.ipynb @@ -77,7 +77,7 @@ "`chaospy` supports Gaussian mixture model representation, but does not provide an automatic method for constructing them from data.\n", "However, this is something for example `scikit-learn` supports.\n", "It is possible to use `scikit-learn` to fit a model, and use the generated parameters in the `chaospy` implementation.\n", - "For example, let us consider the (Iris example from scikit-learn's documentation)[https://scikit-learn.org/stable/auto_examples/mixture/plot_gmm_covariances.html] and its \"full\" implementation and the 2-dimensional representation:" + "For example, let us consider the [Iris example from scikit-learn's documentation](https://scikit-learn.org/stable/auto_examples/mixture/plot_gmm_covariances.html) (\"full\" implementation in 2-dimensional representation):" ] }, { diff --git a/docs/tutorials/advanced/kernel_density_estimation.ipynb b/docs/tutorials/advanced/kernel_density_estimation.ipynb index 5f21faa8..965be551 100644 --- a/docs/tutorials/advanced/kernel_density_estimation.ipynb +++ b/docs/tutorials/advanced/kernel_density_estimation.ipynb @@ -59,7 +59,7 @@ ], "source": [ "import chaospy\n", - "distribution = chaospy.GaussianKDE(samples)\n", + "distribution = chaospy.GaussianKDE(samples, h_mat=0.05)\n", "distribution" ] }, @@ -70,7 +70,7 @@ "outputs": [ { "data": { - "image/png": "iVBORw0KGgoAAAANSUhEUgAAAgAAAAF4CAYAAAAi4UHLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALiAAAC4gB5Y4pSQAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3XuM3eV95/HPd8Yz+JbgiJrNGtvQbrQES9iEQrsRuXkVqshqqmykCCrCioWQbuVo26jtokZV/4hWq6ZahSqpEyWSCQpQmQiRXa3qahM2qxhUpJWVxoBNyCKVYMcWnpCYBHybGX/3jznHc+bMufzOOb/Lc3m/JEuemePxw+H4vL/Pcy5j7i4AAJCXqaYXAAAA6scAAABAhhgAAADIEAMAAAAZYgAAACBDDAAAAGSIAQAAgAwVGgDM7MtmdsLMFgZcZreZHTOzl83sITNbU94yAQBAmYqeADwu6Tf7fdHMpiXtl/QJd3+XpI2S7p58eQAAoAqFBgB3f8bdXxtwkVslnXD3o62P90v6+KSLAwAA1SjrmH6rpOMdH7/a+twKZrZX0t6OT91gNlPSEgCkbq02NL2EqKybtqaXgAr9fOEXv3L3t4/758t8nL7zhwr0vNW5+z5J+9ofT03N+vrZa0tcAoAU3Ty1u+klROOWTeuaXgJq8uCpL52Y5M+XNQAcl7S94+OtkiZaGABIxL8Ioo9xlDUAHJZ0jZntcPdjku6V9GRJ3xtAhgj/cIQfkyj6MsCvmdkJSdOtlwN+zcxuMbODkuTui5Lul/SEmb0s6aykRypbNYCkEf/+btm07vIvYBKFTgDc/Q/6fGlPx2W+J2lHGYsCkCfC3x/BR9l4sx4AQSD+qxF9VIkBAEDjiP9KhB91YAAA0BjCvxLhR50YAAA0gvgvI/zVufmqXzW9hJ5+8Prbml4CAwCAehH+ZYR/MqHGvYgia696SGAAAFAb4r+E8BcXc+Qn1e+/vazBgAEAQC2IP+EfJOfQj+rydXVqsu/DAACgcrnHn/CvRvCbxwAAoDKEn/BLxD5UDAAAKpFz/Ak/0Y8BAwCA0hH//BD8+DAAAChVrvHPMfxEP24MAABKk2P8cws/0U8HAwCAieUYfimf+BP9NDEAAJhIjvHPIfxEP30MAADGRvzTQvTzwgAAYCy5xZ/wIzUMAABGRvzTQPjzxgAAYCQ5xT/F8BN9tDEAACiM+MeL8KMbAwCAQoh/nAg/+mEAADBULvEn/MgJAwCAgYh/XAg/imIAANAX8Y8H4ceopppeAIAwEf94EH+MgxMAAKsQ/zgQ/mJ2bTne2N995OS2xv7uYRgAAKyQQ/wJf3qajPwgg9bV9HDAAADgMuIfvtzjH2rox9Hvv6WuwYABAIAk4h+6XMOfUvCL6v5vrmogYAAAQPwDl1P8cwz+MJ3XSZnDAAMAkDniH65cwk/0i1txXb0w2fdiAAAyRvzDlXr8iX7zGAAAJIv4h4Xoh4UBAMhU6rv/GOOfYviJfrgYAIAMEf/wpBZ/wh8+BgAgM8Q/PKnEn+jHhQEAyAjxD08K8Sf8cWIAADJB/MNC+NE0fhogkAHiHxbijxBwAgAgasS/XoQ/HQwAQOJS3/3HJOb4E/708BAAkLDU4x/T7p/4IzScAACJIv7hiDX+hD9tnAAACSL+4SD+CBUnAEBiiH84Yow/4S/PeVujMzPrtWn+rNb6QtPLWYUBAAAqQPzztSjTw1vep4Obd8klmaQ9c0d0z8lnNC1venmXMQAACWH3j3ER//IsxX+nLk4tJ/bg5p2SpPtOPt3UslbhOQBAIoh/OGLb/RP/8py3NTq4eZcuTs2s+PzFqRkd3LxT5y2cfXc4KwEwNuIfjpjiT/jLd2Zmfd9DfpfpzMx6vfPiL2tdUz+cAACRSz3+MSH+2DR/VtbnaybXpvmzta5nEAYAAEGLZfdP/CFJa31Be+aOaPbS/IrPz16a156554J6NQAPAQARS333H0v8Y0L8q3fPyWckLT3xz2UyufbMPXf586FgAAAilXr8YxLL7p/412NarvtOPq27Tj3L+wAAKFcO8Y9l9x9L/FG/tb4QzBP+euE5AEBkiD/Gwe4f3RgAgIgQ/7DEsvsn/uiFAQCIRA7xjwnxR+wYAIAI5BL/mHb/MSD+GKTQAGBmu83smJm9bGYPma1+L0Mz+4iZ/bD16x/N7IbylwsAzYtl9w8MMnQAMLNpSfslfcLd3yVpo6S7e1z065J+391vkvSwpM+XuE4gW+z+wxJL/Nn9Y5giJwC3Sjrh7kdbH++X9PEel3NJb2/9/kpJpyZfHpA34o9xEH8UUeR9ALZK6rw1vdr6XLe7Jf29mZ2T9Jak2yZfHpCvXOIfk1h2/0ARRZ8E2PnDjVb9nIPWwwR/Lmm3u2+T9NeSvtnjcntbzyU4ZmbH3C+Ns2YgeTnFn91/udj9o6giA8BxSds7Pt4q6UTXZd4j6e3u/nzr40clfaj7G7n7Pnff0f5lxosQgG45xT8mMez+if9orr3xJV1740tNL6MxRR4COCzpGjPb4e7HJN0r6cmuy/xU0vVmdo27/1TSRyS9WO5SgfTlFn92/2hKZ/g7f/+T569vYjmNGDoAuPuimd0v6Qkzm5V0SNIjZnaLpM+7+x53P2VmfybpO2a2IOlNSfdVunIgMcQ/XDHs/lHcoF1/99dSHgjM3YdfqiJTU7O+fvbaxv5+IBS5xV9iACgbx//DTXLcH+IgsPM7T73o7jvG/fM8CA80jPiHLYb4Y7hJH+tP8bkC/DhgoCE5hl+KK/6xYPc/WIrxLgMnAEADco1/bNj9x63sZ/mnNkgwAAA1yzn+7P5Rl9RiXQUGAKBGOccf1eD4v7cQn7QXGp4DANSA8Me3++f4P37tIaCM04AUBwoGAKBixD+++CMtkwwCKYa/jQEAqBDxjxO7/zQVHQRSjn4nBgCgAoR/Gbt/hKbXIJBL9DsxAAAlI/6oC08AnEyO0e/EAACUhPCvFuPun+N/5IIBAJgQ4e8txvgjbudtjc7MrNem+bNa6wtNLyd4DADAmAg/EIZFmR7e8j4d3LxLLskk7Zk7ontOPqNpNfcD70LHAACMgfgPFuvun+P/OC3Ff6cuTi0n7eDmnZKk+04+3dSygsc7AQIjuHlqN/EfItb4I07nbY0Obt6li1MzKz5/cWpGBzfv1Hljn9sP1wxQANEvhvijbmdm1vc95HeZzsys1zsv/rLWNcWCAQDog+jnheP/OG2aPyvr8zWTa9P82VrXExMeAgC6cMw/Hnb/aMJaX9CeuSOavTS/4vOzl+a1Z+45Xg0wACcAgNjtTyr2+LP7j9s9J5+RtPTEP5fJ5Noz99zlz6M3BgBkjfBPLvb4I37Tct138mnddepZ3gdgBAwAyA7RLw/xR0jW+gJP+BsBAwCyQPTRD8f/yBUDAJJF9KvF7h+IGwMAkkHw60P8gfgxACBaBL8ZKcWf43/kjAEA0SD4zUsp/ik4cnKbdm053vQyECkGAASJ2IeH+ANpYQBA44h9+FKMP8f/yB0DAGpF7OOTYvwBMACgIoQ+DcQ/fDwPAONiAMBECH2aUg8/x/8AAwAKIvT5SD3+AJYwAOAyIo8c4p/i7p+HATAOBoCMEHj0k0P4AazEAJAYIo9REX8gTwwAESLyKEOO4U/x+L+NhwEwKgaAgBF6VCHH8ANYjQGgYUQedco5/inv/ts4BcAoGABqQujRpJzDL+URf2BUDAAlIvIITe7hzxGnAOk7cnJbKd+HAWAMhB4hI/orsftHKsoKfxsDwBDEHrEg/GjjFCAtZYe/jQGghdAjVoS/v5x3/wwBcasq+p2yHACIPWJH9IfLOf5tDAHxqSP8bVkMAAQfKSD6QLrqDH9bcgMAsUdKiP542P0v4xQgXE1Ev1PUAwCxR2oI/uSI/2oMAWFpOvxtUQ0ABB+pIfjlIv79MQQ0K5Todwp2ACD2SBHBrw7xH44hoF4hRr9TUAMA0UdqCH49iH9xDAHVCj36nRodANZqA9FHMoh9M4j/6BgCyhNT8LsFdQIAxIDQh4P4j68dLgaB0cUc/U4MAEAfhD5chL88nAYMl0rwuzEAIGtEPj7Ev3ycBqyUavC7MQAgWcQ9LYS/erkOArkEvxsDAKJE3PNB+OuX8iCQa+x7YQBAcIg72oh/s7pjGdtAQOwHYwBAbQg7iiL8YeoMaijDAJEfHwMAJkbYURbCH49B4S1rOCDu1So0AJjZbkn7JM1KOiTp0+6+0HWZjZK+IunfSLok6UF3/1q5y0WdCDvqQvjTQrjjMHQAMLNpSfslfdTdj5rZtyTdLekbXRf9oqQX3P3fm5lJ2lz6alEKwo4QEH2gWUVOAG6VdMLdj7Y+3i/pM+oYAMzsbZJ+V9J2SXJ3l3S63KWiCOKOkBF9IBxFBoCtkjof0Hm19blOvyHpNUl/a2a/1br8H7n7K2UsEssIPGJD9IEwFX0SoHf83np8fUbSTZI+5+7/0cw+JekhSf+280JmtlfS3uU/RMx6IfJIAeEHwlZkADiu1tF+y1ZJJ3pc5ufu/g+tj/9O0n/r/kbuvk9LTyaUJK2feod3XyYXRB6pIfhAXIoMAIclXWNmO9z9mKR7JT3ZeQF3f83MXjCzW9z9sKTbJR3t8b2yQ+iRMqIPxGvoAODui2Z2v6QnzKz9MsBHzOwWSZ939z2ti/6hpP1mtkHSGUmfqmrRISL0SB2xB9JS6DkA7v49STu6Pn1Y0p6OyxyT9N7ylhYuYo/UEXsgfbwT4ACEHjkg9kCeGABaiD1SR+iBNPzg9beV8n2yHQAIPlJD4IF0lRX9TtkMAAQfsSPwQH6qCH9bsgMAwUdMiDuATlWGvy2JAYDYIyTEHMC46gh/W5QDAMFH1Yg4gDrVGf62aAYAoo8iCDeAmDQR/rZgBwCCnx/iDSAXTYa/LagBgOjHj4gDQH8hhL+t0QFg3bQR/cARdACYXEjhbwvqBAD1IOoAUI8Qw9/GAJAY4g4AzQs5/G0MAJEh8AAQrhjC38YAECAiDwBxiSn8bQwADSHyABC/GMPfxgBQMUIPAOmJOfxtDAAlIvYAkLYUwt/GADAmYg8AeUkp/hIDQCHEHgDylVr42xgAuhB7AICUbvjbsh8ACD4AoFvq8ZcyHAAIPgCgnxzC35b8AEDwAQBF5BR/KdEBgOij264tx5teAiJy5OS2ppeAGuUW/rYkBgCCHz4CjJgUub0yJKQh1/hLEQ8ARL8+xBtYrde/C4aCuOQcfymyAYDol4eoA+VjKIhD7uFvC3oAIPijI+xAWNr/JhkEwkD8lwU3ABD9wQg8EKfOf7sMA80g/isFMQAQ/dUIPZAuTgXqR/xXa3QAWDt9Kev4E3kgb5wK1IP49xbECUAOiD2AQXZtOc4QUAHi3x8DQAWIPYBxMASUi/gPxgAwIWIPoEw8P6AcxH84BoAREHsAdeE0YHzEvxgGgAEIPoAmMQSMjvgXxwDQgeADCA1DAKqS/QBA9AGEjucFFJPT7v/wmXMTf4/sBgCCDyBWnAb0R/xHl8UAQPQBIF3EfzzJDgBEv3zX3vhS00sI0k+ev77pJSAjnALkq8z4SwkNAAR/NMS8PGVclwwRGAVDwLJcdv9lx1+KfAAg+qsR9jiN8/+NoSFvDAH5qCL+UqQDQO7hJ/KQit8OGBSQqhx2/1XFX4poAMgt+kQeZSlyW2JIiFPOpwCpx7/K8LcFPQDkEn1ij6YNuw0yIIQr5yEgVXXEXwp0AEg5/MQeMep3u2UwQBNS3/3XJZgBIMXoE3ukjsEgDJwCpKOu3b8UwACQSviJPbCs+98DAwHKkvLuv874Sw0PAGvXLDT510+E4APF9fr3wlAALKs7/lIAJwAxIfpAeTglKBcPA8SrifhLDAADEXygPp3/3hgG0EuKx/9NxV9iAFiB4ANh4HRgPJwCxKXJ+EsMAEQfiED73ymDQL5S3P03LcsBgOgDcWIQQCqa3v1LGQ0ARL9ZG9775uXfv/XsxgZXghTwfIH+UnwYILXdfwjxlzIYAAh/MZ2BTunvkhg4UsepAGISSvylRAcAor+s7tiGaNTrgIEhTgwCwGgKDQBmtlvSPkmzkg5J+rS793wXHzP7SuvrtQ4XuUafwJevyHXKkBAuBoG0HgZI6fg/pN2/VGAAMLNpSfslfdTdj5rZtyTdLekbPS77fkkbSl/lALmEn9CHZdD/D4aDMFQ5CJy3NTozs16b5s9qbe+9ELBCaPGXip0A3CrphLsfbX28X9Jn1DUAmNkVkv5K0sck3VXmIntJNfyEPn79/h8yGDTj2htfKm0IWJTp4S3v08HNu+SSTNKeuSO65+QzmpaX8ncAdSkyAGyV1PkTe15tfa7bX0ra7+5zZlbG2lZJMfohBP/c4qzm5q/U5pk3tG76YtPLSRaDQXPKOg1Yiv9OXZxavus8uHmnJOm+k09P9L2xWirH/yHu/qXiTwLsHG1X1d3Mdkr6bUl/MeibmNleSXvbH181O1voL08l/CHEvtOim75w/A49evr2y5/75NXf1QPbHte0sZupS6/bBUNBNSY5DThva3Rw864V8Zeki1MzOrh5p+469SwPB2CVUOMvFRsAjkva3vHxVkknui5zm6Qdkv65tfufNrNXJL3H3X/RvpC779PSkwklSf9q48aBlYk9/KEFv9sXjt+hx05/WBd8eRB77PSHJUmf236gqWVBq287DATlGXcIODOzvu8hv8t0Zma93nnxl5MtDqhRkQHgsKRrzGyHux+TdK+kJzsv4O5flfTV9sdmtuDu1427qFjDH3rwO51bnNWjp29fEX9JOu9X6NHTt+uz1zzJwwEB4ZSgXOMMAZvmz64+/mwxuTbNn518YbgsheP/kHf/kjQ17ALuvijpfklPmNnLks5KesTMbjGzg2Uu5tobX4ou/hve++blXzGZm79yoq+jeZ23vdhufyEY9f5mrS9oz9wRzV6aX/H52Uvz2jP3XHDH/7u2HB9+IVQm9PhLBZ8D4O7f09IRf6fDkvb0ufxI7wEQY/Rjt3nmjYm+jvDwdsvjGeU04J6Tz0haeuKfy2Ry7Zl77vLngZg0+k6AM1dciCb+KUS/07rpi/rk1d/VY6c/rPN+xeXPr7ULuuvqpzj+jxzDQDWm5brv5NO669SzvA9AhWI//o9h9y8l+lbAZUkt+t0e2Pa4JK14FcBdVz91+fNIA8PAcKM+J2CtL/CEP0SPAaBL6tHvNG2uz20/oM9e8yTvA5CJ9u2bQWC1Mt8wCPmKZfcvMQBcllP4u62bvqjt03NNLwM14lSgN4aA5sV8/B9T/KXMB4Ccow+0cSoA5GnoywBTxMumgNX4d7Eklicmp4jdf72yGgC4gwOG498JQwDykMVDALnfmRV19gN31vr3rT/E2w2HjIcGgGJi3P1LiQ8AhH9J3WEvapR1MSw0Z8N732QIiNCRk9uaXsJIYj3+jzX+UqIDQI7hDzXyZRn038dwUL0cTwN4RQBSl9QAkEP4Uw/9OPpdJwwG5eM0AFgW8+5fSmQASDn8BH983dcdA0E5cjwNQLViPf6PXdQDQGrhJ/bV6nX9MhSML4fTAB4GQD+x7/6liAeAVOJP9JvVef0zDIwuhyEA1WL335zoBoAUwk/0w8QwAKCIFHb/UmQDQMzxJ/pxaf//YhAYjlOA8MT2EsCYpBJ/KZIBINbwE/34cSpQDEMAxsHxf7OCHwBijD/hTxOnAoMxBCB1Ke3+pYAHgNjCT/TzwSCAkMVy/M/uv3lBDgAxxZ/w54tBYDVOAZCq1Hb/UoADQCzxJ/xoYxAARsPuPwxB/TjgGOJ/9gN3En/0xO1iSQz/joFRpLj7lwIaAGK40+AOHsNwG0GTYnj8P7bdf6rxlwIZAEKPP7t+jILbSjp4G2CkrPHnAIQcf+7IMS6eF4C6sfsvX8q7f6nhE4Cp9Zea/OsHIv4oQ663o5AHewBLgngIICQc96Ns3J5QNXb/5Ut99y8xAKzAHTUAIIf4SwwAlxF/AJ1ieQIgu3+MiwFAxB/Vy+02xrsBoi22+Oey+5cYALK7YwYwHLv/POUUfynzAYD4A0B1Ytv95ybbAYD4o27c5uLA7j9Pue3+pUwHAO6IAaBaMe3+c4y/lOEAQPyBavEEwOqFvvuPKf45y2oAIP4ABonh+D/0+Mcm192/lNkAAKBaMe/+Y4h/DGLa/eccfymjAYDdP5rGDwbCpELf/ccUf2Q0AACoFrv/ahH/cuW++5cyGQDY/QPoJ4b4o1zEf0kWAwDQtNSP/2Pd/ccSf3b/5SH+yxgAAEwk1vjHgvijKgwAQMVS3/3HKobdP/EvF7v/lZIfAHj8H01KPf6x7v6Jf36I/2rJDwBAU4h/mGKIfwxi2v0T/94YAIAKEP8wxRL/0Hf/xD8NDAAARkL8q0X8URcGAKBkKe/+iX+1iH+52P0PxgBQsnPzUzp+Zp3OzXPV5oj4h4f4l4P4p2dN0wtIxeIl6cFD79aBH14nl2SS7rzpFX32Az/SNLNA8lIOv0T8q0b8y5VD/H9w6f9M/D2ST1Ndd8wPHnq3Dhy5ThcWp3VxcVoXFqf1+JHr9OChd9fy96M5Kcf/rWc3Ev+KEf9yEf/ikh8A6nBufkoHfnidLixMr/j8+YVpHThyHQ8HJGr9oQPJxz9WxL8cxD9tlKkEP3vrCnmfr7kvfR1pSTn8EvGvQ+jxR5jK2v1LmTwHYP2hA5W+I+Cvbbgg6/M1s6WvIw2EP2zEvzzs/sNTZvwlTgBKsW7mku686RWtXbO44vNr1yzqzl2vaN3MpYZWhjIR/3D95PnriX+JiH94yo6/lMkJgFT9KcBnP/AjSdKBI9fJfWnnf8euVy5/HvFKPfxS/PGPBfEvH/EfXzYDQNWmp6Q//dCPtPe2H+tnb12hX9twgZ1/5Ah/+GKJP+GvBvGfTFYDQPsOvcqTgHUzl7RtU/o3ypTlEv5zmtHrulJX6Q2t03zTSxpJLOGXiH9ViP/kshoA2qp+OABxyiH8kvTLZ9+mv5n6dzpgH5TLZHLd6d/XH1/6tqb7vp4lHMS/XMQ/X4UGADPbLWmfpFlJhyR92t0XOr6+TdLDkrZIWpT03939L0pfbYkYAiDlE31p+bi/Hf8LNnv5a4/rg9KU9CeXnmxqeUPFFH6J+Fcll/hXvfuXCrwKwMymJe2X9Al3f5ekjZLu7rrYgqQH3P0GSTdLer+Z/V7Ziy1bTnf+WCn1N/Hp1Plufuc0syr+knTeZnXAPqhzmmliiUPFFP8jJ7cR/4oQ/3IVeRngrZJOuPvR1sf7JX288wLufsrdD7d+f1HSc5KuLXOhVcklAliSU/il1U/ye11Xyvu8a4XL9LqurGNZhcX08j4pjl2/RPxDVlf8pWIPAWyVdLzj41dbn+vJzK6S9DFJvzPZ0urDwwFpyyn4bf2e3X+V3pD1eZzf5LpKb1S5rJHEFH6J+FeJ+Fej6JMAO+8x+r3pnczsCklPSPqiu7/Y4+t7Je1tf/wvNvX9VrWr4xUCqE+O0ZeGv6xvneZ1p39fj+uDOt/xMMBav6g7/PtBvBqA8FeH+Ier7vhLxQaA45K2d3y8VdKJ7gu1nivwmKTD7v5gr2/k7vu09GRCSdKObdPBPeWY04B45Rp9abTX8//xpW9LU9IBLb8K4I7WqwCaFFv4pXjiH2P4JeJfNXMf3OBW2H8s6aPufszMDkj6X+7+ja7L7dfS6cB9PuybtuzYNu3P/83kjzmeW5zV3PyV2jzzhtZNX5z4+3ViGAgb0R9fSO8DEFv8Ywm/RPxDN0n837rw8ovuvmPcPz/0BMDdF83sfklPmFn7ZYCPmNktkj7v7nvM7DZJ90p6QdI/mZkkPeTuXxp3YUUsuukLx+/Qo6dvv/y5T179XT2w7XFNWzmHCzw0EJacg99W1rv3rdO8tupnpXyvccUWfon4Vy2X8EvN7fzbhp4AVGnSE4D/+uqdeuz0h3Xel3/c7lq7oLuufkqf215NKBgE6kf0l8T+tr2dCH+1Ygy/RPxHVfkJQKjOLc7q0dO364J3vZ7Zr9Cjp2/XZ695svSHAyROBOpA8Fci/M0j/tUj/vWLdgCYmx98cjA3f6W2T89V9vd3RophYDIEf7WUoi8R/roQ//CFEn8p4gFg88zg1ysP+3qZugPGQNAfsR+M8IeB8NeH+Den0QHg0tkpvfXsRm1475sj/9l10xf1yau/2/c5AFUc/xeV+0BA5EeTWvSleMMvEf86Ef9mBXECMO4Q8MC2xyVpxasA7rr6qcufD0W/IMY8GBD5yaQYfYnw14nwxyPE+EsNvwrg3VfN+P/9D//y8sfjDAFSte8D0LS6hwTCXp1Uoy8R/roR/3hUGf+kXgXQvoMcdRBYN32x0if8NYkgx43ohyvG8EvEPyah7vzbghoA2sZ9SAAIQcrRl+IPvxRn/GMOv0T8QxTkACCNfxoANCH16EuEv0kxxz+38EtxxF8KeABoYxBAqIh+PAh/M4h/2IIfANoYBNC0HILfRvibFXv4pfziH1P426IZANoYBFAnoh+nWMMvxR//3MIvxRl/KcIBoI1BAFXIKfhSWtGXCH/TiH9coh0A2hgEMIncgi8R/dCkEH6J+Mco+gGgrfOOnGEA/eQY/DbCH54U4p9j+KX44y8lNAB0YhhAW87Bl9KLvkT4Q0L845bkANCJYSAfuce+LcXoS4Q/JIQ/DckPAJ0YBtJA6FdKNfhSGtGX0gm/RPxTktUA0Kk7IgwEYSL2vaUcfYnwhyjX8Etpxl/KeADo1is0DAX1IfTDEf14pBR+ifinigFgAIaCchH50aQe/DbCHy7CnzYGgBH1ixiDAYEvA9GPE+FPSw7xlxgASlMkfrEOCYS9OrkEX0ov+lJ64ZeIfy7xlxgAakVIkVPwpTSjLxH+VOUUf4kBAKhUbsGX0o2+RPhTlVv42xgAgBLlGHyJ6MeK+Ocbf6nhAWD+whX6yfPX69obX2pyGcBYco19W8rRlwh/6nIOf1sQJwCdd6QMAwhV7sGX0o++RPhzQPyXBDEAdGIYQAiI/TKiHz/Cv4z4LwtuAOjEMIA6EPuVcgh+G+HPB+FfLegBoFP3nTQDAUZF6Psj+ukh/suIf2/RDADdOB3AIMR+sJyC30b480MaeaZHAAAGc0lEQVT4B4t2AOjU686eoSB9RH40RD9thH8l4j9cEgNALzxkkA5CP54cgy/lFX2J8Hcj/MUlOwB06xcRBoPmEfhy5Br8NsIP4j+abAaAfhgMqkfgq5F78KX8oi8R/l4I/3iyHwD6KRKtnIcEol4vYr+M6KMT8R8fA8AExolgaEMDIQ8TwV8px+hLhH8Qwj85BoCaEVx0I/a95Rp9ifAPQvjLwwAA1Izg95Zz8NsI/2DEv1wMAEBFCP1wRH8J4R+M8FeDAQAoAbEvjugvIfrDEf5qMQAAIyL2oyH4KxH+Yoh/9RodAM4vMH8gXIR+fER/JaJfHOGvT+MF7r6T3bXleEMrQa4I/eQIfm+EvzjCX7/GB4BuDASoApEvF8Hvj+iPhvA3J7gBoFuvO26GAvRC5KtD8Icj/KMj/s0KfgDohaEgX0S+HgS/GKI/HsIfhigHgF76hYHBIC4Evn7EfjREf3yEPyzJDAD9DAsKA0I9CHs4CP7oiP5kCH+Ykh8AhikSJoaE1Qh6HIj9+Ij+5Ah/2LIfAIoYN3YhDw4EPE0EfzJEvxyEPw4MABUisqgSsS8H0S8P4Y8LAwAQAWJfLqJfLsIfJwYAICCEvjpEv3yEP24MAEADCH31CH51CH8aGACAChH6ehH9ahH+tDAAACUg9M0g+PUg/GliAABGQOibR/TrQfTTV2gAMLPdkvZJmpV0SNKn3X1h1MsAMSDyYSH49SL8+Rg6AJjZtKT9kj7q7kfN7FuS7pb0jVEu08v5xalJ1g6MjciHi+A3g/Dnp8gJwK2STrj70dbH+yV9RivjXuQyPfW7I775ql8VWBrQG4GPB8FvDtHPW5EBYKukzve0fbX1uVEvMxIGA/RD3ONG8JtH+CEVfxKgd/zexr2Mme2VtLfjUwsPnvrS/yu4hiWnRro0VnuHpF80vYjMcJ3Xj+u8flzn9fv1Sf5wkQHguKTtHR9vlXRijMvI3fdp6YmCkiQzO+buOwqvFhPjOq8f13n9uM7rx3VePzM7NsmfL/IsvMOSrjGz9v/YeyU9OcZlAABAIIYOAO6+KOl+SU+Y2cuSzkp6xMxuMbODgy5T3bIBAMAkCj0HwN2/J6n7aOewpD1DLjPMvuEXQcm4zuvHdV4/rvP6cZ3Xb6Lr3Nx9+KUAAEBSeCceAAAyxAAAAECGahkAzGy3mR0zs5fN7CEzW/XcgyKXQXHDrk8z22Zm/9vMXjSzF8zsvzS11lSMchs2s6+YGT8rY0IF71s2mtk3zezHZvYjM/uDJtaaioLX+UfM7IetX/9oZjc0sdYUmNmXzezEoPuLcftZ+QDQ8XMCPuHu75K0UUs/J2Cky6C4gtfngqQH3P0GSTdLer+Z/V69K03HKLdhM3u/pA01Li9JI1znX5T0grv/a0k3SPp2fatMywjX+dcl/b673yTpYUmfr22R6Xlc0m/2++Ik/azjBKDXzwn4+BiXQXFDr093P+Xuh1u/vyjpOUnX1rrKtBS6DZvZFZL+StKf1ri2VA29zs3sbZJ+V0tDgHzJ6VpXmZai99Uu6e2t318p3sN1bO7+jLu/NuAiY/ezjmP2Rn6WQOZGuj7N7CpJH5P0OxWvK2VFr/O/lLTf3efM+r2rNgoqcp3/hqTXJP2tmf1W6/J/5O6v1LLC9BS9nd8t6e/N7JyktyTdVsPacjV2P+t6EmApP0sAIyl0fbZ2pE9I+qK7v1j5qtI28Do3s52SflsFfkomCht2O5+RdJOk/+HuN0v6n5IeqmNhCRt2O5+W9OeSdrv7Nkl/LembNa0tV2P1s44BoLSfJYDCCl2frX+oj0k67O4P1rS2VBW5zm/T0ptl/bOZvSJp2sxeMbN31LPE5BS9b/m5u/9D6+O/09JzXjCeItf5eyS93d2fb338qKQPVb+0bI3dzzoGAH6WQP2KXp9fl/RLSf+5roUlbOh17u5fdfct7n6du18nabH1e36C2niKXOevSXrBzG5pfep2SUeFcRW5b/mppOvN7JrWxx+RxOlidcbuZ+UDAD9LoH5FrnMzu01LN5RbJf1T6+U6/6mxRUeuyHWOco1wnf+hpC+b2XOS/kTSp+pfbRoK3p+fkvRnkr5jZkckPSDpvqbWHDsz+5qZndDSieGJ1sel9JO3AgYAIEO8EyAAABliAAAAIEMMAAAAZIgBAACADDEAAACQIQYAAAAyxAAAAECGGAAAAMjQ/wdguWKh19/j9AAAAABJRU5ErkJggg==\n", + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAgAAAAF4CAYAAAAi4UHLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALiAAAC4gB5Y4pSQAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3X2MXfV95/HPd4YZ/EAwXce0CzZ4q2hJLAUIhXYRG1yvQhW5TRUiRbgirLIQ0m0cbYnSLmq06h+oWiWrVUCbOlUiGaIAlR0h0Coq1SY0FS5qpJXVBrAh6SIVsAPChgYn4KfxzHf/mHs8d2buw7n3noffw/slITHj6/Hhepj355x75465uwAAQF6m2j4AAADQPAYAAAAZYgAAAJAhBgAAABliAAAAkCEGAAAAGWIAAACQoVIDwMy+ZmZHzezcgNvsMLMXzOwlM3vQzC6o7jABAECVyl4B2C/p1/r9oplNS9or6ZPu/j5JF0m6Y/LDAwAAdSg1ANz9GXd/Y8BNbpB01N0Pd97eK+kTkx4cAACoR1WX6TdLOtL19qud9y1jZrsl7e561wfMZio6BACxWqP1bR9Co9ZOW9uH0Lg10wttH0JyXjtz4hfufvG4v7/Kx+m7f6hAz89ud98jaU/x9tTUrK+bvbLCQwAQm+umdrR9CI26/pK1bR9CY67b+Iu2DyFpdxx66Ogkv7+qAXBE0hVdb2+WNNGBAUhfTvHPJfxEPx5VDYCDki43s23u/oKkOyU9XtHHBpAg4p8Wwh+fst8G+A0zOyppuvPtgN8ws+vN7ElJcvd5SXdLeszMXpJ0UtLDtR01gKgR/zRct/EX5/9BfEpdAXD33+/zSzu7bvMDSduqOCgA6col/qmHH/HjxXoANIb4x43wp4UBAKARxD9ehD9NDAAAtSP+8SH66WMAAKhVDvFPKfwS8c8FAwBAbYh/XAh/XhgAAGpB/ONB+PPEAABQOeIfB8KfNwYAgEoR//ARfkglXwkQAMog/uEj/ihwBQBAJVKPP+FHargCAGBixD9sxB+9cAUAwESIf7gIPwbhCgCAsRH/cBF/DMMVAABjIf5hIvzlXXPZkbF+37Ovban4SNrBAAAwMuIfJuK/ZNy4V/WxYxgJDAAAIyH+4ck9/HXGflz9jimkYcAAAFAa8Q9PbvEPMfajWHn8bQ4CBgCAUoh/eHKJf+zRH6T7v63pMcAAADAU8Q9P6vFPOfr9ND0GGAAABiL+YUk5/DlGv58mxgADAEBfxD8sKcaf6A9X3EdVDwEGAICeiH9YUos/4R9d1VcFGAAAViH+YUkl/kS/OtdcdkQ6NNnHYAAAWIb4hyWF+BP+MDEAAJxH/MMSe/wJf9gYAAAkpR//mBB+NIGfBgggi/jHcvZP/NEUrgAAmSP+4Yg5/oQ/PgwAIGPEPxyxxp/wx4sBAGSK+IcjxvgT/vjxHAAgQ8Q/HMQfbeEKAJAZ4h+O2OJP+NPCAAAykUP4JeJfF+KfHh4CADJA/MNC/BECrgAAicsl/rGIKf6EP21cAQASllP8Yzj7J/4ICQMASBTxDwvxR2h4CABIEPHHOAh/XhgAQEJyCr8UT/xjOPsn/vnhIQAgEcQ/TMQfoWIAAAkg/mEi/ggZAwCIHPEPE/FH6HgOABCp3MIfE+KPGHAFAIhQrvGP4eyf+CMWXAEAIpJr+KU44h8D4t+c03aB3p5Zp0vmTmqNn2v7cFZhAACRIP7hC/3sn/g3Y16mb1327/XkpmvkkkzSzuPP6tOvPaNpeduHdx4DAAhczuGXiH9ViH9zFuN/tc5OLSX2yU1XS5Lueu3v2jqsVXgOABAw4k/8q0D8m3PaLtCTm67R2amZZe8/OzWjJzddrdMWznl3OEcC4Lzcwy8R/6oQ/2a9PbOu70V+l+ntmXX6lbM/b/SY+uEKABAY4o+qEP/mXTJ3Utbn10yuS+ZONno8g3AFAAgE4V/C2f/kiP9yV37wJxP9/leev6rU7db4Oe08/mznOQBLDwPMLsxp5/HngvpuAAYA0DLCvxzxn1zu8Z809sM+5rAx8OnXnpG0+MQ/l8nk2nn8ufPvD4W5t/ctCVNTs75u9srW/nygTYR/NeJfjRwHQB3RH2bYEKj7dQCu/t5TL7r7tnF/P1cAgIYR/t5iiX/ocot/G+Hv9Wf3GgNr/FwwT/jrhQEANITw9xdT/EM++88l/m1Gv5/imMo+VyAEDACgZoR/MOJfjRziH2L4V7rygz+JZgQwAIAaEP1yiH81Uo9/DOHvFsvVAAYAUCHCX15M8Ud7Yot/TEoNADPbIWmPpFlJByR91n35UxrN7KOSvtx586Sku9z9xQqPFQgS0R9dbPHn7L95KYQ/9IcDhr4SoJlNS9or6ZPu/j5JF0m6o8dNvynp99z9WknfknRfhccJBOW6qR3n/8FoiH91iH/4Qv5vKfNSwDdIOuruhztv75X0iR63c0kXd/59g6TXJz88IBxEf3LEH8OEHMxxhfrfVOYhgM2Sumfmq533rXSHpL8ys1OS3pV00+SHB7SL2FcntviHLrWz/1AjWZUQHw4o+yTA7pcLXPVzDjoPE/yJpB3u/ryZ3Snp25J+e8XtdkvavfSe6REPF6gf0a9ejPEP+eyf+McptBFQZgAckXRF19ubJR1dcZsPSbrY3Z/vvP2IpK+t/EDuvkeLTyaUtPhSwCMdLVADgl8v4l8t4o+qlBkAByVdbmbb3P0FSXdKenzFbX4q6Sozu9zdfyrpo5L4DgAEieA3h/gDy4V0FWDoAHD3eTO7W9JjZlZ8G+DDZna9pPvcfae7v25mfyzpe2Z2TtI7ku6q9ciBEoh9e2KMf+g4+0eVSj0HwN1/IGnlTxw6KGln120ekvRQdYcGlEfowxJr/EM++yf+qBqvBIjoEPuwEX8MQ/zDwABAkIh8nIh/PVI6+yf+4TwPgAGA1hD5tMQafzSH+IeFAYBaEPe8xBx/zv6RKwYARkbcUYg5/BLxbxJn/+FhAOA8wo5REH+URfzDxADIBHFHlWKPfwxSOvtHmBgAiSDwaEoK8efsvzmc/YeLARAJAo+2pRB+KY74c/afthC+BVBiAASFyCNUxL85xB9NYQC0gNAjFqmEX4oj/qnh8n/YGAA1IvSIWUrxjwVn/5M5pRm9pQ3aqBNaq7m2Dyd4DIAKEHqkJMXwc/aftnmZHpi6Vftsu1wmk2uXP617Fp7QtLztwwsWA2AEhB6pI/7tSe3sv8nL/0X8z9js+fft13ZpSvriwuONHUcZoTwBUGIADETwkYsUwy/FE3+M75RmVsVfkk7brPZpuz6n7/JwQB8MgA5ijxylGn4prvindvbfpLe0QS7r+Wsu01vaoM16s+Gj6i2ks38p4wFA8JGzlMMvxRX/FDV5+X+jTsj6PM5vcm3UicaOJTbZDACCDywi/mHh7H8yazWnXf609mu7Tnc9DLDGz+o2fzqYy/+hnf1LCQ8Agg8sl3r4Y0T8q3HPwhPSlLRPS98FcFvnuwBCEGL8pYQGAMEHVsst+rGd/aMa03J9ceFxfU7f5XUARhDtACD4QH+5hV+KL/6c/VdvreaCecJfIdSzfymiAUDwgeFyDL8UX/yRh5DjLwU+AIg+MFyu0S/EGH/O/tMXevylAAcA0QfKyT38UpzxB0IRxAAg+kA5RH9JrPHn7D99MZz9Sy0PgDVaT/yBIYj+arHGH+mLJf5SIFcAACwh+P3FHn7O/tMWU/wlBgDQOoJfDvFHqGILf4EBADSM4I8u9vjn5pXnr2r05wG0Kdb4SwwAoFbEfnIpxJ+z/zTFHH+JAQBUiuBXK4X4I02xx19iAAAjI/L1Syn8nP2nJYXwFxgAwAoEvl3EPw0pPg8gpfhLDABkhriHLaX4Ix2phb/AAEDUCHoaUgx/zmf/hdivAqQa/gIDAEEg5Pki/mmLcQSkHv4CAwCVIeIYRYrhR9xyCX+BAYCeiDnqlHL8OftfrQhrqFcCcgt/gQGQCYKOEKQcfon4DxPawwG5hr/AAEgAcUfoUg8/ymv7akDu0e/GAIgAgUescgo/Z/+jaXIIEP3eGAABIPBITU7hl4j/JKoeAsS+PAZAgwg9Updb+CXiXxXC3TwGQE2IPXKSY/iB2DEAKkDskavcw8/ZP2LGABgRsQcIv0T8ET8GwADEHlhC9JcQf6SAAdCF4APLEf3ViD9SkfUAIPjAakS/P+KPlGQ1AAg+0BvRH474IzXJDwCiD6xG8EdD/JGiJAcA0QeWI/jjI/5IVRIDgOADyxH8ahB/pCzaAUD0gUXEvh7EHyF79rUtE3+MqAYA0UfOCH1ziD9CVkX8pQgGANFHToh8+4g/QlZV/KVABwDRR2oIe/gIP0JXZfylwAYA4UcMiHl6iD9CVnX4C60PAKKPKhFnjIr4I2R1xV8qOQDMbIekPZJmJR2Q9Fl3P7fiNhdJ+rqkfydpQdL97v6NQR937bSNc8wIEOFFbAg/Qldn/KUSA8DMpiXtlfQxdz9sZt+RdIekh1bc9KuSDrn7fzQzk7Sp8qNFLYg3ckL4Ebq6w18ocwXgBklH3f1w5+29kj6vrgFgZu+R9DuSrpAkd3dJx6o9VJRF0IHeiD9C11T8pXIDYLOk7v9rXu28r9uvSnpD0p+b2a93bv+H7v5yFQeJ5Qg8MBrCj9A1Gf5C2ScBete/93rgfkbStZK+5O7/2cw+I+lBSf+h+0ZmtlvS7uLtdVPrRjvajBB5YHKEHzFoI/5SuQFwRJ1L+x2bJR3tcZt/cfe/7rz9l5L+58oP5O57tPhkQknSxpl/5StvkxMiD9SD8CMGbYW/UGYAHJR0uZltc/cXJN0p6fHuG7j7G2Z2yMyud/eDkm6RdLjHx8oWsQfqR/gRi7bjL5UYAO4+b2Z3S3rMzIpvA3zYzK6XdJ+77+zc9A8k7TWz9ZLelvSZug46dMQeaA7RR2xCiL9U8jkA7v4DSdtWvPugpJ1dt3lB0o3VHVociD3QDsKP2IQS/kLrrwQYE2IPtIvoI0ahhb/AABiA4APtI/qIWajxlxgAqxB9oH1EH7ELOfyF7AcAwQfaR/CRkhjiL2U6AIg+0D6ij9TEEv5CNgOA6APtIfZIWWzhLyQ9AIg+0A6Cj1zEGn8psQFA8IFmEXrkKubwF5IYAIQfqA+RB5akEP5CtAOA6APVIPBAOSnFX4psABB9YDiCDlQrtfAXohgAhL85xAMAFqUa/kKwA4DoT46YA8DoUg9/IbgBQPjLIe4AUL1c4i8FNAAI/2pEHgCakVP4C60OgDXTC9mHn8gDQLtyjL8U0BWAHBB7AAhHruEvMABqQuwBIEy5h7/AAKgIwQeAsBH+5RgAYyL4ABAP4r8aA6Akgg8A8SH8/TEABiD6ABCvlOP/D2+9Z+KPwQDoQvABIH4ph1+qJv4SA4DoA0AiUg+/VF38pUwHANEHgLSkHv8qw1/IZgAQfQBIT+rhl+qJv5TBACD8AJAm4j+ZJAcA0QeAdBH+aiQ1AAg/AKSN+Fcn+gFA9Kt35Qd/0sqf+8rzV7Xy5wKIQ+rxbyr8hWgHAOEfrK2IT6LKY2ZMAOlIPfxS8/GXIhwAhH9RjIFv0rj3D8MBCAvxr080AyDH8BP55o16nzMYgPqkHv+2wl8IfgDkEH5CHy8GA1C91MMvtR9/KeABkGr4iX3eGAzAYMS/OcENgNTCT/AxiXE+fxgNiFXq8Q8l/IVgBkAq4Sf4aBujATEi/s1rfQDEHn6CjxTwXRNoE/FvR6sDYM0F59r848dG9IFFPKcBkyD87Wr9CkAsiD4wuTL/HzES8kD828cAGIDoA80b9P8d4yANxD8MDIAViD4Qrl7/fzIK4pJy/GMJf4EBIKIPxIxREA/iH5asBwDhB9LU/f82YyAMxD882Q0Aog/kZeX/8wyC5hH/MGUzAAg/AGnpawFDoBmpxj/m8BeSHwCEH0AvDIH6Ef+wJTsACD+AMhgC9SD+4UtuABD+4dbf+E7jf+a7P7yo8T8TGMWVH/wJIwADpRR/KaEBkHv424j6KEY9PgYD2sDVgGqkePafWvylBAZATuEPPfJVKvPfykhAXbgaMD7iH49oB0DK4c8p9JMYdj8xEDAJRsDoUot/quEvRDcAUgs/sa9Pv/uWYYCyGAHlEf/4RDMAUgk/wW8fwwCjYATkJ4f4S5EMgNjjT/Tj0OvviVEAiREwTEpn/7nEXwp8AMQafoKfDkYBMBjxj1ewAyC2+BP9fDAK8sRVgNWIf9yCGwAxhZ/oo7Dyc4FBkCZGQJpyjL8U2ACIJf6EH8NwlQCpS+XsP9f4SyUHgJntkLRH0qykA5I+6+7n+tz2651fLz0uYgg/0cekuEqQBq4CEP9UDI20mU1L2ivpY+5+2My+I+kOSQ/1uO2HJa0f5QBCjz/hR10YBPE6bRfo7Zl1umTupNb0PhdKFvFPR5mz9BskHXX3w52390r6vFYMADO7UNKXJX1c0u1l/vCZC8+UP9IGEX20ofvzjjEQpnmZHpi6Vfuu2a4Fn5JJ2nn8WX36tWc0LW/78FBSCvE/+PapiT9GmQGwWdKRrrdf7bxvpT+VtNfdj5vZxAfWhlzDf2p+VsfnNmjTzAmtnT7b9uFAjIFQPTB1q/bZdp2xWanzZe7JTVdLku567e9aPLJmpHD2T/yXlH2cvnvarqq7mV0t6Tck/bdBH8TMdkvaXbx96bqpkn98vXIN/7ybvnLkNj1y7Jbz7/vUpd/XvVv2a9o4mwkFYyAMpzSzFP8uZ6dm9OSmq3X76z9M+uEA4h+GquIvlRsARyRd0fX2ZklHV9zmJknbJP1z5+x/2sxelvQhd/9ZcSN336PFJxNKkt6/cabVyuQa/sJXjtymR499RGd86Qvao8c+Ikn60hX72josDMAYaM9b2iBfff4jSXKZ3p5Zp185+/OGjwplEf/VypyCH5R0uZlt67x9p6THu2/g7n/h7pe5+1Z33yppvvPvP1Ogco//qflZPXLsFp32C5e9/7RfqEeO3aJT87N9fidCsf7Gd87/g/pt1AlZn8f5Ta5L5k42fETNif3sP/b4H3z7VOXxl0oMAHefl3S3pMfM7CVJJyU9bGbXm9mTlR9RzfiCuej43IaJfh1hYQzUb63mtMuf1hpf/jyZ2YU57Tz+XLKX/4l/u+oIf6HUcwDc/QdavMTf7aCknX1uH9QLDEmc8a+0aebERL+OcBWf6zxEUL17Fp6QpqR92r74cIBLO48/p0+/9kzbh4YeiP9gwYW6aoS/t7XTZ/WpS7+vR499ZNnDAGvsjG6/9Cm+GyABDIHqTcv1xYXH9Tl9V29pg9499MvJnvlLcZ/9E//hkh4AxH+we7fsl6Rl3wVw+6VPnX8/0sAQqN5azWmz3tQrvrHtQ0EPxL+cZAcA8R9u2lxfumKfvnD547wOQAYYAhhFrGf/xL+85AZA7uE/efOusX7fe7X4Yg+TPI953QG+dTAGDAEMQ/zb0WT8pcQGQE7xHzf0dRp0TIyD8DAEkBLiP7okBkDq4Q8x9qPq99/AMGgfQwDdYj37j1kb8ZcSGAApxj+F4Je18r+VQdCe9Te+wwjIXKzxj/nsv634S5EPgFTin1Pwh2EQtIsRgNgQ//FFOwBijz/RL6f7fmIMNIMRkKcYz/6J/2SiHAAxx5/wj48x0BxGwHCvPH9V24eQNeI/uagGQKzhJ/rVYwzUjycH5iO2s3/iX41oBkCM8Sf8zSjuZ4ZAPbgakDbi35yQ4i9FMgBiij/Rbw9XBerDCAAmE1r8pRI/DrhtxB/jOHnzLv4+gCE4+29GiPGXAr8CEEv8CU24eHigOlwFQJuIf/WCvwIQMs4y48HfUzViGeV1S+E7AGI6+yf+9Qh2AIT+hYagxIfBBiwi/vULPf5SoA8BhBx/AhI/HhaYDA8FAIPFEH8pwCsAxB9N4e8T44j98j9n//WKJf5SYFcAQo0/oUgXVwPGw1WAOBH/esUUfymgKwDEH23i7xkIB/FvRjADIEREIS/8fWOYmC//x3L2T/ybE8QACPHsnxjkib93ALlofQAQfwAx4Oy/fpz9N6v1ARAa4g8+B5AS4l+fmOMvtTwAptYttPnHr8IXfhT4XEC3mM/+Y0D828EVgA6+4GMlPicQu1jO/mOTQvwlBoAkvtADSE8s8Y/t7D+V+EsMAGAgxiG4/F8f4t+u7AcAX+AB9BNr/GM4+yf+7ct6ABB/AKkh/tVLMf5S5gMAwHhy+DkAMZ79xxB/hCPbAcDZP8ricwWoDmf/4chyAPAFHcAgnP3Xg/iHJcsBAAD9EP96EP/wMAAAjCSHx/9jQvyrl0P8pQwHAJf/gfGlHv8Yz/5RrVziL2U4AIBRrTuwr+1DQANijD9n/9XKKf4SAwBASamf/ceG+GNSWQ0ALv8D6CW2s3/iX73czv6lzAYAgPGkfPYfW/xjQPzjwAAABuDx/7TFGP8Yzv5jkmv8JQYAgCFSPfsn/vWI6ew/5/hLDACgL87+041/jIh/tXKPv8QAAHoi/mnHP7azf+JfLeK/iAEAYBXiHw7iXy3iv4QBAKyQ+9k/8ccoYoo/lmMAAF1yj3/KYox/6Gf/scWfs//lGAAVOzU3pSNvr9WpOe5axCfVs3/iD+K/2gVtH0Aq5hek+w+8X/t+tFUuySTtuvZlfeHmH2uaLRCF3M/+iX84Yoh/TGf/xL83BkBF7j/wfu17dqvOzE+ff9/+Z7dKkv7oN3/c0lGhDMKfZvgl4l8X4p+GrM5N6/pCf2puSvt+tFVnzk0ve//pc9Pa9+xWHg4IGPEn/iEh/tUi/oNRpgq8+e6F8j6/5r746wgP8Sf+ISH+1SL+wzEAKvDe9WdkfX7NbPHXERbiT/xDQvwxqn9Y+NuJP0Z2A6COL/xrZxa069qXteaC+WXvX3PBvHZd87LWzixU/mdifMSf+IeE+Fcv9bP/KuIv8STAynzh5sUn+u17dqvcF8/8b7vm5fPvRxiIP/EPCfGvHvEvz9z7PXpdv21bpv35Bza08mefvHlXLR/31NyU3nz3Qr13/RnO/ANC+Al/aIh/9XKL/7tnXnrR3beN+/G4AlCxtTML2nJJ2p+EMck9/FL/+J/SjN7SBm3UCa3VXMNHVQ3iXx/iH5Yqz/wL2Q6AdQf21XYVAGHIPf79wj8v0wNTt2qfbZfLZHLt8qd1z8ITmu77/SzhIf71If55KPUkQDPbYWYvmNlLZvagmV2w4te3mNnfmNmLZnbIzP6snsOtVu6BSNW6A/uy/7sddMm/iP8Zm9VZm9EZm9V+264Hpm5t8AgnQ/xRyCH+dZz9SyUGgJlNS9or6ZPu/j5JF0m6Y8XNzkm6190/IOk6SR82s9+t+mCBQQj/YvgHxf+UZs7Hv9tpm9U+265Tmqn7ECfyyvNXEf+axXT2T/wnU+YKwA2Sjrr74c7beyV9ovsG7v66ux/s/PtZSc9JurLKA61L7sFIBX+P5Z7o95Y2yPu8aoXL9JbaeVJuGbGGXyL+dSD+kyvzHIDNko50vf1q5309mdlGSR+X9FuTHVpzeD5AvAj/aM/w36gTsj6P85tcG3WiqsOqFPGvH/EPS93xl8o/CbD7K0a/F72TmV0o6TFJX3X3F3v8+m5Ju4u3f/mSvh+qcYyAuBD+RaN+e99azWmXP6392q7TXQ8DrPGzus2fDu67AQh/M4h/WJqIv1RuAByRdEXX25slHV15o85zBR6VdNDd7+/1gdx9j6Q9xdvbtkwH9ZRjRkD4CP+iSb6v/56FJ6QpaZ+Wvgvgts53AYSE+DcjpvjnoKn4SyVeCKgT9n+S9DF3f8HM9kn6P+7+0Irb7dXi1YG7vOSrC1X1QkCn5md1fG6DNs2c0NrpsxN/PEZAWIj+kipf0CfU1wGIOfwS8a9T6mf/o8a/9hcCcvd5M7tb0mNmNivpgKSHzex6Sfe5+04zu0nSnZIOSfpHM5OkB939f417YGXMu+krR27TI8duOf++T136fd27Zb+mbfyLC0VwGALtIvxL6nglv7Wa02a9WfnHnQTxb0Zs4ZeIfx2ifing//7qLj167CM67Us/bneNndHtlz6lL11RTTwYAc0j/Mul/DK+BcLfHOIfnnHjn+1LAZ+an9Ujx27RGV/x/cx+oR45dou+cPnjlTwcwNWAZhD91XIIv0T8m0T80S3aAXB8bvCVg+NzG3TF9PHK/jyeIFg9ot8b4Y8H8a9XDvFv49J/odUBsHCy1CsR97RpZvD3Kw/79XFwNWByRL8/wh8X4l8v4l+/1q8AvPvDi7T+xndG/n1rp8/qU5d+v+9zAKq4/N8PQ2A0RH+wXMIvpRH/mMIvEf9QtR1/KYABII0/Au7dsl+Sln0XwO2XPnX+/XXrDhtjYAnBHy6n6EtphF8i/k0g/s1p9bsA3r9xxv/vf/rX598eZwRI1b8OwKRyHANEvxzCH6fYwi8R/1BVGf9svwug29rps5U+4W9SqV8ZIPajI/zxIv7NIP7NC2oAjPtQQMhWxjK2QUDsJ5Nb+CXi36YYwy8R/7YENQCkNEdAt15BDWEUEPpqEf64xRZ+ifiHLMT4SwEOACn9EbBS2fiOMhQIevNyjL6UVvgl4t8k4t+uIAeAlN8IKIOoh4nwpyHG8EvEP2Qhx18KeABIjACEi+ing/A3j/iHIegBIC19oWUIoG25Rl9KM/wS8W8D8Q9H8AOgwNUAtIXwpyfW8Evxxj+H8EvxxF+KaABIjAA0J+foS+mGXyL+bSD+YYpqAEiMANQn9+hLhD9UsYZfIv4hi24ASDwvANUh+mlHX4o7/BLxj0GM8ZciHQAFhgDGQfQXEf7wEf/wxRp/KfIBUOBhAQxD9BelHn2J8IeA+MchiQEgcTUAqxH9JTmEXyL+ISD+8UhmABQYAvki+MvlEn2J8Icil/inIrkBUGAI5IHor0b44xN7/HMLfwpn/1LCA6DAEEgLwe8tp+hLhD8kxD9eyQ+AQnc4GAPxIPj95RZ9KZ3wS8Q/NimFv5DNAOjGVYFwEfzBcoy+RPhDRPzjl+UAKHBVoF3Evhyin44U4p9T+KV04y9lPgC6MQbqR/DLyzX6EuEPGfFPCwOgh5WhYhCMjtiPLufoS2mGXyL+sUo9/hIDoBQGQX+EfjK5R18i/DEg/mliAIyhV/Tz7FfaAAAG7UlEQVRSHwWEvjpEP93oS4Q/drnEX2IAVGZQIEMfB8S9fkR/EeGPB/FPHwOgAVUEduWIINphI/hLUo5+gfjHL7f4SwyAaBD88BH95Qh/fHIMv5Rn/KWWB8DcmQvb/OOBiRD81XKIvpRe+CXin6PWrwAUX0Sv/OBPWj4SYDiiv1ou0ZfSDL9E/HPV+gAoMAQQIoLfH+GPX67hl4i/FNAAKDAE0CaCP1hO0ZfSDb9E/BHgACgwBNAEgj9cbtGXCH+qCP9ywQ6AAkMAVSL45eQYfSnt8EvEH8sFPwAK3V+4GQMoi+CXl2v0JcKfOuLfWzQDoBtXBdALsR9dztGX0g+/RPyJf39RDoACVwXyRvDHk3v0JcKfC+I/WNQDoBtjIG3EfnwEf0kO4ZeIv0T8y0hmAHRjDMSN2E+O6C9H+PNB+MtLcgB0WxkTBkE4CH11CP5quUS/QPyJ/6iSHwAr9YoOo6BehL4eRL83wp8n4j+67AZAL1wlqAahrxfB7y+36EuEvxvxHw8DoId+Ict9GBD4ZhH84XIMv0T8C4R/MgyAEQwLYMwDgbi3j+CXk2v0JcLfjfhPjgFQoXEiWtVoIODxIfijIfwoEP9qMABaRrjzQfBHl3P0JcLfC/GvDgMAqAnBH0/u0S8Q/+UIf/UYAEBFCP74iP4Swr8a8a8HAwAYA7GfHNFfjvD3RvzrwwAASiD41SD6qxH+3gh//RgAwArEvlpEvzfC3x/xb0arA+D0uQvOf7G95rIjbR4KMkXs60H0+yP8gxH/5gRzBaD7CzFjAHUg9vUi+oMR/sEIf/OCGQDdGAOYFLGvH8Evh/APR/zbUWoAmNkOSXskzUo6IOmz7n5u1NuMY+UXcgYBViL2zSH65RH+4Qh/u4YOADOblrRX0sfc/bCZfUfSHZIeGuU2VeHqQL4IffMI/ugIfznEv31lrgDcIOmoux/uvL1X0ue1PO5lblM5rg6kidC3i+iPjuiPhviHocwA2Cypu6yvdt436m1qxyCIC6EPA8EfH+EfDeEPS9knAXrXv9u4tzGz3ZJ2d73r3B2HHvp/JY9hdIdq+8gx+yVJP2v7IDLDfd487vPmcZ83799M8pvLDIAjkq7oenuzpKNj3EbuvkeLTxSUJJnZC+6+rfTRYmLc583jPm8e93nzuM+bZ2YvTPL7p0rc5qCky82s+Iu9U9LjY9wGAAAEYugAcPd5SXdLeszMXpJ0UtLDZna9mT056Db1HTYAAJhEqecAuPsPJK28tHNQ0s4htxlmz/CboGLc583jPm8e93nzuM+bN9F9bu4+/FYAACApZZ4DAAAAEsMAAAAgQ40MADPbYWYvmNlLZvagma167kGZ26C8YfenmW0xs78xsxfN7JCZ/Vlbx5qKUT6HzezrZjbxz8rIXcmvLReZ2bfN7J/M7Mdm9vttHGsqSt7nHzWzH3X++Xsz+0Abx5oCM/uamR0d9PVi3H7WPgC6fk7AJ939fZIu0uLPCRjpNiiv5P15TtK97v4BSddJ+rCZ/W6zR5qOUT6HzezDktY3eHhJGuE+/6qkQ+7+byV9QNITzR1lWka4z78p6ffc/VpJ35J0X2MHmZ79kn6t3y9O0s8mrgD0+jkBnxjjNihv6P3p7q+7+8HOv5+V9JykKxs9yrSU+hw2swslfVnSHzV4bKkaep+b2Xsk/Y4WR4B80bFGjzItZb9Wu6SLO/++QdLrDRxbktz9GXd/Y8BNxu5nE5fZo/lZAgkZ6f40s42SPi7pt2o+rpSVvc//VNJedz9u1u9VtVFSmfv8VyW9IenPzezXO7f/Q3d/uZEjTE/Zz/M7JP2VmZ2S9K6kmxo4tlyN3c+mngRYyc8SwEhK3Z+dM9LHJH3V3V+s/ajSNvA+N7OrJf2Gav4pmZkZ9nk+I+laSf/b3a+T9F1JDzZxYAkb9nk+LelPJO1w9y2S/oekbzd0bLkaq59NDIDKfpYASit1f3b+R31U0kF3v7+hY0tVmfv8Ji2+WNY/m9nLkqbN7GUz+6VmDjE5Zb+2/Iu7/3Xn7b/U4nNeMJ4y9/mHJF3s7s933n5E0m/Wf2jZGrufTQwAfpZA88ren9+U9HNJ/7WpA0vY0Pvc3f/C3S9z963uvlXSfOff+Qlq4ylzn78h6ZCZXd951y2SDgvjKvO15aeSrjKzyztvf1QSVxfrM3Y/ax8A/CyB5pW5z83sJi1+otwg6R87367zX1o76MiVuc9RrRHu8z+Q9DUze07SFyV9pvmjTUPJr+evS/pjSd8zs2cl3SvprraOOXZm9g0zO6rFK4ZHO29X0k9eChgAgAzxSoAAAGSIAQAAQIYYAAAAZIgBAABAhhgAAABkiAEAAECGGAAAAGSIAQAAQIb+P0Jsb8J3DTtGAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] diff --git a/pyproject.toml b/pyproject.toml index 1e65886c..5d5bdad7 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -4,7 +4,7 @@ build-backend = "poetry.masonry.api" [tool.poetry] name = "chaospy" -version = "4.0-beta3" +version = "4.0.1" description = "Numerical tool for perfroming uncertainty quantification" license = "MIT" authors = ["Jonathan Feinberg"] diff --git a/tests/distributions/kernel/test_gaussian_kde.py b/tests/distributions/kernel/test_gaussian_kde.py new file mode 100644 index 00000000..e76e6b58 --- /dev/null +++ b/tests/distributions/kernel/test_gaussian_kde.py @@ -0,0 +1,63 @@ +"""Tests for Gaussian kernel density estimation.""" +import numpy +import chaospy + + +def test_gaussian_kde_1d_integration(): + """Make sure that 1D distribution integration is correct.""" + dist = chaospy.GaussianKDE([0, 1, 2]) + t = numpy.mgrid[-2.6:4.6:2e5j] + scale = numpy.ptp(t) + assert numpy.isclose(numpy.mean(dist.pdf(t)*scale), 1.) + assert numpy.isclose(numpy.mean(t*dist.pdf(t)*scale), 1.) + + +def test_gaussian_kde_2d_integration(): + """Make sure that 2D distribution integration is correct.""" + dist = chaospy.GaussianKDE([[0, 2], [2, 0]]) + samples = dist.sample(1e4) + assert numpy.allclose(numpy.mean(samples, axis=-1), 1, rtol=1e-1) + + +def test_gaussian_kde_rotation(): + """Make sure rotation does not affect mapping.""" + dist = chaospy.GaussianKDE([[0, 0, 2], [0, 2, 0], [2, 0, 0]], rotation=[0, 1, 2]) + grid = numpy.mgrid[0.01:0.99:2j, 0.01:0.99:2j, 0.01:0.99:2j].reshape(3, 8) + inverse = dist.inv(grid) + assert numpy.allclose(dist.fwd(inverse), grid) + assert numpy.allclose(dist.pdf(inverse), + [2.550e-05, 2.553e-05, 2.553e-05, 2.552e-05, + 2.525e-05, 2.522e-05, 2.522e-05, 2.519e-05]) + assert numpy.allclose( + inverse, + [[-1.38424, -1.38424, -1.38424, -1.38424, 3.19971, 3.19971, 3.19971, 3.19971], + [-1.31003, -1.31003, 3.31003, 3.31003, -1.48391, -1.48391, 1.48429, 1.48429], + [ 0.51571, 3.48391, -1.48391, 1.48413, -1.48391, 1.48429, -1.48391, 1.48429]], + rtol=1e-5, + ) + dist = chaospy.GaussianKDE([[0, 0, 2], [0, 2, 0], [2, 0, 0]], rotation=[2, 1, 0]) + inverse = dist.inv(grid) + assert numpy.allclose(dist.fwd(inverse), grid) + assert numpy.allclose(dist.pdf(inverse), + [2.550e-05, 2.525e-05, 2.553e-05, 2.522e-05, + 2.553e-05, 2.522e-05, 2.552e-05, 2.519e-05]) + assert numpy.allclose( + inverse, + [[ 0.51571, -1.48391, -1.48391, -1.48391, 3.48391, 1.48429, 1.48413, 1.48429], + [-1.31003, -1.48391, 3.31003, 1.48429, -1.31003, -1.48391, 3.31003, 1.48429], + [-1.38424, 3.19971, -1.38424, 3.19971, -1.38424, 3.19971, -1.38424, 3.19971]], + rtol=1e-5, + ) + dist = chaospy.GaussianKDE([[0, 0, 2], [0, 2, 0], [2, 0, 0]], rotation=[0, 2, 1]) + inverse = dist.inv(grid) + assert numpy.allclose(dist.fwd(inverse), grid) + assert numpy.allclose(dist.pdf(inverse), + [2.550e-05, 2.553e-05, 2.553e-05, 2.552e-05, + 2.525e-05, 2.522e-05, 2.522e-05, 2.519e-05]) + assert numpy.allclose( + inverse, + [[-1.38424, -1.38424, -1.38424, -1.38424, 3.19971, 3.19971, 3.19971, 3.19971], + [ 0.51571, -1.48391, 3.48391, 1.48413, -1.48391, -1.48391, 1.48429, 1.48429], + [-1.31003, 3.31003, -1.31003, 3.31003, -1.48391, 1.48429, -1.48391, 1.48429]], + rtol=1e-5, + )