diff --git a/.github/workflows/tests.yml b/.github/workflows/tests.yml index 30bce949d..eb9bbd942 100644 --- a/.github/workflows/tests.yml +++ b/.github/workflows/tests.yml @@ -57,7 +57,7 @@ jobs: run: | python3 -m pip install --upgrade pip python3 -m pip install -e".[test]" - python3 -m pip install "jax[cuda12_local]" -f https://storage.googleapis.com/jax-releases/jax_cuda_releases.html + python3 -m pip install "jax[cuda12]" - name: Run nvidia-smi run: | diff --git a/docs/experimental/index.rst b/docs/experimental/index.rst new file mode 100644 index 000000000..27ec09d3b --- /dev/null +++ b/docs/experimental/index.rst @@ -0,0 +1,11 @@ +ott.experimental +================ +.. module:: ott.experimental + +The :mod:`ott.experimental` module groups experimental code that might be useful +for users but whose API we expect to change significantly in coming months. + +.. toctree:: + :maxdepth: 2 + + mmsinkhorn diff --git a/docs/experimental/mmsinkhorn.rst b/docs/experimental/mmsinkhorn.rst new file mode 100644 index 000000000..19098fa3b --- /dev/null +++ b/docs/experimental/mmsinkhorn.rst @@ -0,0 +1,16 @@ +ott.experimental.mmsinkhorn +=========================== +.. module:: ott.experimental.mmsinkhorn +.. currentmodule:: ott.experimental + +Solvers for multimarginal entropic OT problems, defined using :math:`k` point +clouds of variable sizes in dimension :math:`d`, as proposed in +:cite:`benamou:15`, and presented in :cite:`piran:24` (Algorithm 1). + +Multimarginal Sinkhorn +---------------------- +.. autosummary:: + :toctree: _autosummary + + mmsinkhorn.MMSinkhorn + mmsinkhorn.MMSinkhornOutput diff --git a/docs/index.rst b/docs/index.rst index 3d23091d8..1e7f671d5 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -117,6 +117,7 @@ Packages solvers/index initializers/index neural/index + experimental/index tools math utils diff --git a/docs/references.bib b/docs/references.bib index e996ca942..4296a702d 100644 --- a/docs/references.bib +++ b/docs/references.bib @@ -25,6 +25,16 @@ @article{cuturi:15 year = {2016}, } +@inproceedings{piran:24, + author = {Piran, Zoe and Klein, Michal and Thornton, James and Cuturi, Marco}, + publisher = {PMLR}, + booktitle = {Proceedings of The 41st International Conference on Machine Learning}, + month = jun, + series = {Proceedings of Machine Learning Research}, + title = {Contrasting Multiple Representations with the Multi-Marginal Matching Gap}, + year = {2024}, +} + @inproceedings{peyre:16, author = {Peyré, Gabriel and Cuturi, Marco and Solomon, Justin}, editor = {Balcan, Maria Florina and Weinberger, Kilian Q.}, @@ -356,6 +366,16 @@ @article{solomon:15 year = {2015}, } +@article{gangbo:98, + author = {Gangbo, Wilfrid and Święch, Andrzej}, + journal = {Communications on Pure and Applied Mathematics}, + number = {1}, + pages = {23--45}, + title = {Optimal maps for the multidimensional Monge-Kantorovich problem}, + volume = {51}, + year = {1998}, +} + @inproceedings{genevay:18, author = {Genevay, Aude and Peyre, Gabriel and Cuturi, Marco}, editor = {Storkey, Amos and Perez-Cruz, Fernando}, diff --git a/docs/spelling/technical.txt b/docs/spelling/technical.txt index f7997b48c..ccb8b390f 100644 --- a/docs/spelling/technical.txt +++ b/docs/spelling/technical.txt @@ -88,6 +88,7 @@ linearized logit macOS methylation +multimarginal neuroimaging normed numerics @@ -100,6 +101,7 @@ parameterization parameterizing piecewise pluripotent +polymatching polynomials positivity postfix diff --git a/docs/tutorials/Sinkhorn_Barycenters.ipynb b/docs/tutorials/Sinkhorn_Barycenters.ipynb index 1e3817275..85457104f 100644 --- a/docs/tutorials/Sinkhorn_Barycenters.ipynb +++ b/docs/tutorials/Sinkhorn_Barycenters.ipynb @@ -6,7 +6,9 @@ "id": "C4nZjbMHWcm_" }, "source": [ - "# Sinkhorn Barycenters" + "# Sinkhorn Barycenters\n", + "\n", + "This tutorial covers the computation of Wasserstein barycenters using regularized OT, based on {cite}`cuturi:14` and {cite}`benamou:15`." ] }, { @@ -56,20 +58,11 @@ }, { "cell_type": "code", - "execution_count": 3, + "execution_count": null, "metadata": { "id": "ysURew0UKhHE" }, - "outputs": [ - { - "name": "stderr", - "output_type": "stream", - "text": [ - "/home/michal/mambaforge/envs/ott/lib/python3.10/site-packages/nilearn/datasets/struct.py:850: UserWarning: `legacy_format` will default to `False` in release 0.11. Dataset fetchers will then return pandas dataframes by default instead of recarrays.\n", - " warnings.warn(_LEGACY_FORMAT_MSG)\n" - ] - } - ], + "outputs": [], "source": [ "n_subjects = 4\n", "dataset_files = datasets.fetch_oasis_vbm(\n", diff --git a/docs/tutorials/index.rst b/docs/tutorials/index.rst index 63429a468..dbe6f8d51 100644 --- a/docs/tutorials/index.rst +++ b/docs/tutorials/index.rst @@ -18,6 +18,7 @@ Linear Optimal Transport OTT_&_POT Hessians LRSinkhorn + mmsink sinkhorn_divergence_gradient_flow sparse_monge_displacements diff --git a/docs/tutorials/mmsink.ipynb b/docs/tutorials/mmsink.ipynb new file mode 100644 index 000000000..676f1dd8b --- /dev/null +++ b/docs/tutorials/mmsink.ipynb @@ -0,0 +1,156 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Multimarginal OT\n", + "\n", + "We consider in this tutorial the resolution of the multimarginal OT problem, as first formulated in {cite}`gangbo:98` and solved numerically in {cite}`benamou:15` using entropic regularization." + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [], + "source": [ + "import jax\n", + "import jax.numpy as jnp\n", + "\n", + "import matplotlib.patches as ptc\n", + "import matplotlib.pyplot as plt\n", + "\n", + "from ott.experimental import mmsinkhorn" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We sample 4 small and uniform point clouds, each of size 6, in dimension 2, and solve the regularized multimarginal OT problem using the {class}`~ott.experimental.mmsinkhorn.MMSinkhorn` solver. By default, the squared-Euclidean distance is used to compare the pairs of points." + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [], + "source": [ + "n_s, d = [6] * 4, 2\n", + "\n", + "rngs = jax.random.split(jax.random.PRNGKey(0), len(n_s))\n", + "x_s = [jax.random.uniform(rng, (n, d)) for rng, n in zip(rngs, n_s)]\n", + "a_s = None\n", + "\n", + "out = mmsinkhorn.MMSinkhorn()(x_s=x_s, a_s=a_s, epsilon=1e-2)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We now extract the biggest mass transfers across 4-tuples of points (one from each point cloud) using the {class}`multimarginal OT output object `." + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [], + "source": [ + "top_k = n_s[0] * 2\n", + "val, idx = jax.lax.top_k(out.tensor.ravel(), top_k)\n", + "indices = jnp.unravel_index(idx, n_s)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "A small function to help plot 4-tuples." + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [], + "source": [ + "def ccworder(A: jnp.ndarray) -> jnp.ndarray:\n", + " # https://stackoverflow.com/questions/5040412/how-to-draw-the-largest-polygon-from-a-set-of-points\n", + " A = A - jnp.mean(A, 0)[None]\n", + " return jnp.argsort(jnp.arctan2(A[:, 1], A[:, 0]))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We can now plot some elements of the multimarginal OT tensor, by representing 4-tuples as polygons.\n", + "\n", + "Because the marginal distributions are uniform, and the number of points in each point cloud is the same (here 6), we expect that the OT tensor will be close, numerically, to a polymatching tensor of size $6^4$. \n", + "\n", + "We list the 4-tuples with the largest 12 transport values, each displayed as a quadrilateral linking those 4 points, filled with a transparency that is proportional to the transported mass. We do observe that only 6 of these 4-tuples stand out, the remaining 6 being barely visible due to their low transparency." + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiMAAAGdCAYAAADAAnMpAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy88F64QAAAACXBIWXMAAA9hAAAPYQGoP6dpAACQiUlEQVR4nO3deXxbZ5k3/N99jnZZlhd5TZzYTuLse5s0TdONlAxLWYZCpzAtU7YHpsPLkHlmaIelA8zQDsM2DxT6UmDgeYFpGei+pA1pQ7eUtNkTO0kTJ/G+25ItWcs5537/kOXItiRrP+dI1/fzyQeiSNZt1da5dN3XfV2Mc85BCCGEEKISQe0FEEIIIaS4UTBCCCGEEFVRMEIIIYQQVVEwQgghhBBVUTBCCCGEEFVRMEIIIYQQVVEwQgghhBBVUTBCCCGEEFUZ1F5AMhRFQU9PDxwOBxhjai+HEEIIIUngnGN8fBz19fUQhPj5D10EIz09PWhoaFB7GYQQQghJQ2dnJxYuXBj333URjDgcDgDhb6a0tFTl1RBCCCEkGR6PBw0NDdPX8Xh0EYxEtmZKS0spGCGEEEJ0Zr4SCypgJYQQQoiqKBghhBBCiKooGCGEEEKIqigYIYQQQoiqKBghhBAyr4ngBHwhn9rLIAWKghFCCCEJhZQQbnvmNvz1s38NWZHVXg4pQBSMxCErMj6797N46PhDai+FEEJU9ZvW3+Ci5yLeHnsbvzv7O7WXQwoQBSNxPHruUbzW8xp+dPRHOD92Xu3lEEKIKvq9/Xjg6APTf//BoR9geHJYxRWRQkTBSAxj/jF8763vAQAYGL75xjfBOVd5VYQQkn//8dZ/IKSEpv8ekAP4/qHvq7giUohSDkZefvll3HzzzaivrwdjDI8//vi8j9m/fz82bdoEs9mMpUuX4pe//GUaS82f/zz8n/BJ4UItmcs41H8Iey7uUXlVhBCSXwd7D+L5i89D5pfrRGQu44nzT+DowFH1FkYKTsrBiNfrxfr16/HAAw/Mf2cAFy5cwHve8x7ccMMNOHr0KP7+7/8en/rUp/D888+nvNh8ODl0En94+w9QuDJ9GwPD/QfvhzfkVXFlhBCSPyElhG+88Q0IbO5lQmQivnHgG1TMSrIm5dk073rXu/Cud70r6fs/+OCDaGpqwne/+10AwMqVK/Hqq6/i+9//Pnbt2pXq0+eUrMj4+oGvQ2DCjE8CHBxj/jE8eOxB/MMV/6DiCgkhJD9+0/obXPJcivlvMpeni1lvW3FbnldGClHOa0YOHDiAnTt3zrht165dOHDgQNzHBAIBeDyeGX/y4dFzj+L0yOkZgUiEAgX/t/X/UjErIaTgzS5ajYeKWUm25DwY6evrQ01NzYzbampq4PF4MDk5GfMx9913H5xO5/SfhoaGXC9zRtFqPFTMSggpBvs798Mv+6f/zjmDHKgCl80z7ueTfHi56+U8r44UIk2eprnnnnvgdrun/3R2dub8OR8+8/CcmhDOw38iIsWsRwaO5Hw9hBCilncsfgdsBhu4YoLkbUJoZDskzzqE3JvBFROA8Iczh9GBGxfdqPJqSSFIuWYkVbW1tejv759xW39/P0pLS2G1WmM+xmw2w2w2x/y3XFlVuQocMzMejAGSdxEEaw8YUwAosBosaHDkPlNDCCFq8fttuKHiH/DoqdcAsOnbuWxFyL0RRuchQJDwD1f8A5xmp3oLJQUj55mRbdu2Yd++fTNu27t3L7Zt25brp07JtQuvxY4FOyAyccbtTAxBHl8FrhgAbsBn130eZeZKKApt1RBCCkdQUnCscwz/34GL+P2hLlQb1qHCUgnG2Iz7cakEsmczVpSvwQeXfVCl1ZJCk3IwMjExgaNHj+Lo0aMAwkd3jx49io6ODgDhLZY77rhj+v6f/exn0d7ejn/6p3/C6dOn8eMf/xi/+93v8MUvfjE730EW3bP1njm/eIJpAFwugTKxDg32Zty6/COQFY6grMAfkhGUFMgUmBBCdGpoIoAXT/fjoVfa8dKZAQx7gwAAgQm4vuH6mDVySsiB7aVfgKLM+SdC0pJyMPLWW29h48aN2LhxIwBg9+7d2LhxI772ta8BAHp7e6cDEwBoamrCM888g71792L9+vX47ne/i5/97GeaO9YLAA2OBnx67afBotKSTJAhmAahSBbsKPsC/CEGUbj87wrnCEUFJpKsUIErIUTTFIXjbP84fvdWJ379xiUc73IjJM+NLGrttVhesXxGrxHGGFZVrgKkCjx9vJeyxCQrGNfBldPj8cDpdMLtdqO0tDSnz+WX/Lj5sZsx4BuAgvAvJ5OdaBI+jJuX3AyjKODGFdVYUGZFMMYvbwRjgMgYBMYgCCzu/QghJJaJgITHj3QjKCWXfjAZBHxw4wLYzfFLAScCEo53jeFUtwfeoJTU1/WFfPhN228QUkJgYDCJJnxs5cdgMVgAAMuqHXj32to5WWVCgOSv3zkvYNUbi8GCL1/1ZXz+xc9P32Y2BfCXS29CKASEZAV7W/txzTIXWmocCMmxt2k4ByTOAXBABkSBhYMTCkwIIUk41e3GVx4/CYbwh5tEOAc4gJYaB7Y0Vcz5945hH451jeHCkBdKip8/bUYbttZtxavdr4KDY1v9tulABADeHhjH3laGd66uTenrEhKNgpEYrm+4HjsW7MBrPa9B4Qo+v/HzuKJ8IV4/PwQgvDXz8tlBTAQkbFpUDpHxhFkSAJAVDnkqMBFYeKtHpMCEEBLHlY0VWFJlnwogEt9XYMCSqhJcsbh8+raAJKO1x4PjXW6M+oIZrWWNaw3aRtogMhErK1bO+ffWXg+MBgE3LK/O6HlI8dJknxEtuGfrPRCYgMbSRty28jYsrS6BUZz5ch2+NIpX3x4CY4DFKCYdXETXmQQkmepMCCFzCALDP7xz+byBCAAoHPiHm1ogCAyD4wH8sbUfP3vlAv50djDjQAQIF7N+aNmH8MGlH4y7HXOscwyvnxvK+LlIcaLMSBwNjgb87J0/Q42tBkbBCAhAc5UdZ/rGZ9zvdJ8HvqCEG1dUwygKSWVJokW2cySFg7GprAlt5xBCAPzF6tp5syMCA5pcdiyutOGRNzvQ6/bHvmOGDML8l4uDF0dgMgi4onHuVhEhiVBmJIHNNZux0LFw+u8ramMX33SM+PDMiV74QzIEgaWUJYnGOWYcG45Xj0IIKQ7JZEcUDqxbWIYXWvtzFoik4tVzQzjeNab2MojOUDCSgiqHGa6S2J1hB8cDePJYD9yTIQCAURRgEjN7eWVl5rFhWeG0nUNIkYlkR2J9vmEAKuwmLKqw5X1dibx4egBtvfkZcEoKAwUjKVpe64j7b57JEJ461oOB8fCnE0FgMBsECFk48hapMwlICtWZEFJEEmVHOIBtzZXznrZRwwun+nFuYELtZRCdoGAkRUuq5hayRvOHZDx7vBcdwz4A4QZBJoOQ8DGp4hyQFI6AdHk7hxoPEVK4YmVHIlmRpdUlqq0rEQ6O50704tKwd/47k6JHwUiKTAYBS6oS//JLCsfetn6c7rucphSzmCWZjdrTE1LYYmVHtJwViZA5x9PHe9E9Nqn2UojGUTCShhV18bdqIjjnePXtIRy6NDJ9Wy6yJLNRe3pCCtOuVTVwlZjCTdCg7axItJCs4Imj3RjwqF9cS7SLgpE0uErMqHLELmSd7UjHGF4+OzhjGyWXWZJoytSR4eg6E9rOIUSfzg5MYGtTJab6Oms+KxItKCl47Eg3hicCai+FaBQFI2mKd8w3lrP943ihtW/GIKp8ZEmiRepMoo8NU2BCiD4oCscb7cNYWl2CCrtJN1mRaJMhGY8d6YbbF1J7KUSDKBhJU3OVHSZD8i9f1+gknjneC9+s4VT5ypLMRnUmhOhHa68H7skQGAM+snkhPrJ5oW6yItEmAhL+cLgLE4HkhvSR4kHBSJqMooCl8xSyzjY0EcBTx3rmfDLId5ZkNmpPT4h2yVNZkQiLSYTFJKq4osx4/CE8ergLk0FZ7aUQDaFgJAOJeo7EM+6X8OTxHvTHKOZSK0sSLfrYcECi7RxC1Ha8a6zgMgkj3iAePdIFf4gCEhJGwUgGKlMoZI0WCMl47kQvLg7NPX+vdpYkGrWnJ0RdIVnBWxdH1V5GTgyOB/Dk0Z4ZtXSkeKl/xdO5lXXJF7JGkxSOfacHcKrHHfPftZAlmY3a0xOSX8c6x+ANFlZWJFqPexJPHeuhDzmEgpFMNbvsMKdQyBqNc44D54dx8MJIzH+PZEkMGpzgS+3pCcmtgCTjrUuFmRWJFhk0StvBxY2CkQwZxPk7ss7neNcYXjozEPfTgUEUYDYImq2ep/b0hGTfkY6xoqmpaB+cwAutffSBpohRMJIFK9Lcqol2fmACz5/qQ1CKvX/KGIPZIGoySzIbHRsmJDP+kIzDHYWfFYl2um8cL54eUHsZRCUUjGRBhd2E6lJLxl+nZ2wSTx/vgTdB5bzWsySzUXt6QlJ36NJo3A8mhexEtxsvnx1UexlEBRSMZMnKNI75xjLiDeKpYz0Y9Qbj3kdPWZJo1J6ekPn5ghKOdo6pvQzVHO4YndFXhRQHCkaypCmDQtbZJgISnj7eg1534kmXesuSRKP29ITEdvDCSNEfd32jfbjotqmKHQUjWWIQBSytzk52BAACkoI9J/vQPjiR8H56zZLMRnUmhADj/hBOdMU+7l9sXj47iJPd9FoUCwpGsmhFlrZqImSF48XTA0n9Quo5SzIbtacnxerghRHI9LM+bV/bAM72j6u9DJIHFIxkUbndhFpn5oWss73RPpzUHmqhZEmiUXt6UizckyGc6vGovQxN4eBJZYiJ/lEwkmXZzo5EnOx248XT/UltXxhEASaxMLIk0ag9PSlkb7QPQ6GsyBwK53j2RC86R3xqL4XkEAUjWdZYaYfZmJuJmu2DXjx3sjepRkiCEM6SiAWUJZmN2tOTQjHiDeJ0L21HxCMpHE8e60Gfe+6AUVIYKBjJMoMoYFl1Zh1ZE+lz+/HM8d6kp3gaCzRLMhu1pyd69kb7MDjo5zWRkKzgsSPdGBwPqL0UkgMUjORArrZqIkZ94V4kIwl6kUQrhixJNGpPT/RkcDxARZpJCkgyHjvSlbAPE9EnCkZyoMxmQp3TmtPn8E71IukeS9yLJFokS1Js6Ngw0bLXzw+pvQRd8QVl/OFwFzz+kNpLIVlUfFemPFme4+wIAAQlBc+f7MO5geQrzQWBwWIsnizJbNSenmhJn9uPC0NetZehOxMBCY8e6ko4OoPoCwUjOdLkyl0hazSFc+w/M4BjKbaPLtYsSTRqT0/URlmR9I1NhvDoke6imWxc6Ir7apRDosDQksNC1tnevDiC188NpfQpv9izJNGoPT3Jt65RHzrouGpGhicCeOxId1EOFSw0FIzk0Ira0rw+X2uvB39sG4CU4lwLypLMRXUmJNdeP0/D4LKh3+PHE0e7U37fI9pCV6AcctqMqC/LbSHrbJeGvXjuZF/KqUvKksRH7elJtl0a9qInheJzklj32CSePt5LHxp0jIKRHMtHIets/R4/njrWk1a1OWVJEqP29CQbKCuSfReHvdhzso8+KOgUXXVyrLHSDkseCllnc0+G8NSxHgxNpN4giLIkyaH29CQd5wcn0O+hTqK58PbAOPa29qu9DJIGCkZyTBQYWmrynx0BgMmgjGeOpz/TwSgKMFKWJGnUnp4kg7IiudXa68FLZwbUXgZJEV1p8kCNrZqIkKxgb2t/2h0eRYHBbBAgFHo/+Syj9vQkljN94xhOI1tJUnOscwyvn6Nj03pCwUgeOK1GLMhzIWs0hXO8fHYQRzpG03o8YwwmA2VJ0kXt6QkAcM7xRjtlRfLl4MURvHVxRO1lkCTR1SVPVtTl95hvLIcujeLVt1PrRRKNsiTZQceGi1NrrwejPpqpkk+vnhtKuSEkUQcFI3myuMIGqyn/hayzne7zYG9rf9pn8ilLkl3Unr44KArHn9vpU7oaXjozgNYej9rLIPOgK0qeCALDcpUKWWfrGPHhmRO9GbVRpixJ9lF7+sJ1ssdNg91UtLe1H+cGaDKyllEwkkctKhayzjY4HsCTx3rgnkz/DZKyJLlD7ekLhyQrOHiBsiJq4uB47kQfLg3TUEKtoqtIHpVajFhYrl4h62yeqV4kA+OZ9TygLEnuUZ2Jfh3vdmOCpsuqTuYcTx/vRTd1vtUkCkbyLN/zaubjD8l49ngvOoYzG9hFWZL8ofb0+hGSFbxJWRHNCMkKnjjajQFqOqc5dOXIs0UVNthMBrWXMYOkcOxt68fpvsyLvChLkl/Unl7bjnSMYZJG3GtKUFLw2JFu6veiMRSM5JkgMCyvLVF7GXNwzvHq20M4dCnzT3GUJVEHtafXloAk49Cl9Hr7kNyaDMl49HA33D4qKtYKbX1ELxItNQ4c7XRrMrV+pGMM3oCMa5a6IGQ4m0YUGAQmICgr0OC3WvBkhUMGR0gGBMam/nuEg0WSBcPngUD8ExqHOgMIjE59+hbNgK08TwsjyfAGJfzhcBc+cmUDSsx0KVQb/RdQgcMS7sjaNZpZnUaunO0fhy8o4R0razLObjDGYDaIkGQFEn1KV43CORQ5/PozBohTwQkFJmkaPg/8cFPcf57kJhyR/wIz3mK3fJYCEo3x+EN49HAXPry5QRN9oIoZ5dFVsrJOO8d8Y+kancQzx3vhC2bnFIBBFGA2CKBrn/qoPX0WJMiIAMBbSgtCsz/ryVSjoEUj3iAePdKVUd8lkjkKRlSyqMIGu8ZTg0MTATx1rCdr+6qRLIkhw+0fkl10bDi7vNyMY3yJ2ssgKRgcD+CJo90IpdmZmmSOghGVMMbQopGOrImM+yU8ebwH/Vk8CkdZEu2i9vSZO6isgARK+etNr9uPJ4/2UDCuEgpGVLS81qGLPftASMZzJ3pxcSh73QspS6J91J4+dR5uxUnepPYySJo6R8OjMujnPP8oGFFRidmABg11ZE1EUjj2nc7+wCnKkugDtadPjpvbUQ6agaJn7YMTeKG1jzKCeUbBiMq01pE1Ec45Xj8/lPU5G5Ql0Z9YdSb05g00CEP4a8M+fEz8Izazs7CDOn3q0em+cbx4ekDtZRQVbVdQFoGGCitKzAZdza443jUGb1DCtcuqIGYxgDCIAkSBU18SnaFjw3NVMQ+qxJO4hp9EB69GG1+M8wIDtdjSjxPdbhhFAde2VKm9lKJAwYjKGGNoqXXgsM46NZ4fmMBkUMbOlTUwGbKXYKO+JPrGOSBN1ZowNtVsjbGMG+jpFWPAYjaAxRhAcFMJ3hZr0NY7rtkeQ2Smwx2jMBkEXNVcqfZSCl5aV5EHHngAjY2NsFgs2Lp1Kw4ePJjw/j/4wQ+wfPlyWK1WNDQ04Itf/CL8fkpfRiyv0Uch62w9Y5N4+ngPvDnI6hhEASaRakn0jNrTz2QSGVbXO3HL5oX4xDVNuHqJCxV2k9rLIvN4o30Yhzv09WFRj1IORh555BHs3r0b9957Lw4fPoz169dj165dGBiIvb/229/+FnfffTfuvfdetLW14ec//zkeeeQR/PM//3PGiy8UdrMBiypsai8jLSPeIJ461oNRbzDrX1sQwlmSbG4FEfXIysxjw8VcZ1JqMWJLUwXu2NaIv7pyEdYvLIPFSMeBterls4M42e1WexkFLeVg5Hvf+x4+/elP484778SqVavw4IMPwmaz4Re/+EXM+7/++uvYvn07PvrRj6KxsRHvfOc7cdttt82bTSk2K2q133MknomAhKeP96DPnZtsl5GyJAUn0s8k+tiwrgITcxq/r3EeU+u04IYV1fjMjmbcvL4eS6tLINIPu+bsaxvA2X46KZUrKdWMBINBHDp0CPfcc8/0bYIgYOfOnThw4EDMx1x99dX49a9/jYMHD2LLli1ob2/Hs88+i9tvvz3u8wQCAQQCl1snezzZPU6qRQvLrXBYDBj366eQNVpAUvDcyV5c11KF5qrsTyUWBAazIBZ9qr8QRdeZAOEBi5qvM6lcAnz+8Lxt4aeZHeHHJCAIDEuqSrCkqgT+kIyz/eNo6/WgN0dBPkkNB8eek30wCCwn73HFLqVgZGhoCLIso6amZsbtNTU1OH36dMzHfPSjH8XQ0BCuueYacM4hSRI++9nPJtymue+++/D1r389laXpXqQjq55HjssKx4unB+ALylizwJmT5zCKAkQWrkMghSkybRhR04Y1uVU3T3CRCYtRxLqFZVi3sAxjviBaez043TsOj5/O46hJ4RzPnujF+zcsQINOt9a1Kud9Rvbv349vfetb+PGPf4zDhw/j0UcfxTPPPINvfvObcR9zzz33wO12T//p7OzM9TI1YXmtA0IBpGffaB/GG+3DOfv6gsBgMVItSTGg9vRAmc2Eq5e48IlrmnDL5oVYXe/M6gk2khpJ4XjyWA963ZNqL6WgpJQZcblcEEUR/f39M27v7+9HbW1tzMd89atfxe23345PfepTAIC1a9fC6/XiM5/5DL785S9DEOb+UpnNZpjN5lSWVhBspnAh68Xh7LVdV8vJbjd8QQnXtVTnLGigLElxUTiHMtUJNtLPRND6dk6WLSy3YWG5DTcsr0L7kBdtvR5cGvZBKbIATW0hWcHjR3pwy+aFqHIU37UqF1IKr00mEzZv3ox9+/ZN36YoCvbt24dt27bFfIzP55sTcIhiuGq82D7hJGNFnX4LWWdrH/TiuZO9CEi5G81NWZLiVOzt6Q2igJYaB96/YQE+taMJ17ZU0UUxzwKSjMeOdOXkJGExSjnXt3v3bjz00EP41a9+hba2Nnzuc5+D1+vFnXfeCQC44447ZhS43nzzzfjJT36Chx9+GBcuXMDevXvx1a9+FTfffPN0UEIuW1huQ6nFqPYysqbP7cfTx3pz3mE2cuKGFKdibk9vMxmwaVE5PrZ1Mf76qsXYvLgcJWbqZ5kPvqCMPxzuolqeLEj5J/bWW2/F4OAgvva1r6Gvrw8bNmzAnj17potaOzo6ZmRCvvKVr4Axhq985Svo7u5GVVUVbr75Zvzbv/1b9r6LAtNS68BbF7M7/0VNo75wL5Jdq2tz2uRJEBgsdOKm6BVze3pXiRk7llXhmqUudIz40NbrwflBL0K0lZkzEwEJjx7qwoevaICdgsC0Ma6Djw8ejwdOpxNutxulpfoZLJeuyaCM/z7YUXD7wCaDgHesrMGCstxPKlYUqiUhMxVre/qgpODtgXGc7h1H1+gkOArrfUUrKkvM+PDmhdS8bpZkr9+U19Ygq0nE4srCOzYWlBQ8f7IP5wYmcv5ckVqSQjidRLKjWNvTmwwCVtc78aHNC/GJaxqxfSm1oc+F4YkAHjvSjaBEH4LSQcGIRi3XcUfWRBTOsf/MAI51juXl+UwGAUaqJSExFGN7eofFiCsbw23ob9uyCOsbymClT/JZ0+/x44mj3ZAoK5syepfWqIXlNpRaC6eQdbY3L47g9XNDeXnzFwUGs0GgLAmJS/ft6dNQU2rBDcur8WlqQ59V3WOTePp4b1Fk3bKJghENW15TmNmRiNZeD/a1DeTlUwRjjLIkJCmRY8MBqTiODUfa0L93XT0+fW0zblxRjTqnRe1l6drFYS/2nOwr+IA2m+idWcNaagqjI2siF4e9eO5kH/yh3PUiiUZZEpKqWMeGC1WkDf2tVy7C31zdiK1NlQXVaiCf3h4Yx97W/vnvSABQMKJpVpOIRlfhFbLO1u/x46ljPXk7q09ZEpKuYmpPX2YzYduSSnzimiZ8+IoGrKE29Clr7fXgpTMDai9DF+hor8b0uifxnv/zKrxTTcI455AT/BcSAHxo80Isrdb/FEmrScSu1bVwleSvkyTnHCGZF9wxapJfxdKeXpIVakOfhisbK7B9qUvtZagi2es3dWjRGKfVCIWH96uTIQMosRTGf8bJoIxnjvfiHSursbA8PxmhcJaETZ+sICQdnAMS58DUtGFRKMx+JpE29C01DviCEk73jaOt14PB8YDaS9O0Ny+OwGQQcGVjhdpL0SzKuWmMzWTAXdcvRTJvYYwBLTUlqC0tnGKzkKzghVP9ONs/ntfnpVoSkk3F0J5+dhv6KxZXUBv6BF47N5S3lgZ6RMGIBn3sqkVw2uYvGuMcuLalKg8ryi+Fc7x8dhBHOkbz+rxUS0JyoRiODbtKzLhmmQufvKYJf7lxIVbUltLvUQwvnRlAa49H7WVoEoWxGhTJjnzr2ba4jZsZA5ZVF1ZWZLZDl0bhDcjYvrQyr3NFRIFBYALVkpCsi2znSAovyPb0jDEsqrRhUaUNQUnBuYEJtPV6qA19lL2t/TAZGJZWF3brhlRR6KpR82VHCjUrMtvpPg/2tvbnvaMhZUlIrhV6e3qTQcCq+lJqQz8LB8dzJ/pwadir9lI0hd5pNSpR7Ugh1ook0jHiwzMnevPWiyRapJaESklIrhVye3pqQz+TzDmePt6L7rFJtZeiGRSMaFi87EixZEWiDY4H8OSxHrgn89OLJBpjDGaDCEOBpNKJ9hVynUl0G/r3ra/HsmpHUbahD8kKnjjajQGPX+2laAL1GdG4h15un1M70lJTgo9c0aDamtRkMYp45+oaVDvUyQpxHk6ra/+3hhSqdI8Nf/PpVrx6bijp+79rdS3+/qaWVJeXFn9Ixtv94fqSHndxZQssRhEf3rwQlXnsr5RP1GekQHzsqkV4YP85jPkuZwSKLSsSzR+S8ezxXty4ogaLKvPfnTaSJZFkBVIB7e8T/ZAVDnmqn4nAWDg4SSIw6R6bxJm+5I/Mb2woy2CVqbEYRaxd6MTahU64fSG09npwus+jSiY03/whGY8e7sZHrmhI6hRloaJtGo2Lrh0RGcPmRWWoLtAIOlmSwrG3rR+n+9Q7ImcQBaolIapLpT39/3PjsqS/rsCAv71+abaWmRKnzYhtSypx5/bLbejNhsKuL/EGJfzhcBcmpjpvFyMKRnQgUjsic46/v6mFjpsivF3y6ttDOHRpRLU1UC0J0RJl6shwdJ1J9LThVfWleOfqmnmzKKLAcMvmhapkHmdbUGbFzlU1+PSOJrxnbR2aXPaCbUzo8Yfw6OEuTAbzX6ivBbRNowM2kwH//qF1ONs3jh3LqtAx7MNEQEoqNVvojnSMwRuQcc1Sl2q9GgyiAFGgWhKiHfHa0//9O1rwwqnEk2Q55/i7G5LPouSDQRSwrMaBZVNt6M/0jaOtdxwD44VV/DniDeLRI1340KaFsBTZaSPKjOjErtW1+Pw7wm8QzVUlBdWPIFNn+8fxQmufqrNlKEtCtCzSz6S5yo6dK6sR78dUS1mReGwmAzYuKsdHty7C7QXYhn5wPIAnjnYX3awsOk2jQ31uPx4/0gXGkitcKxauEjPeuboGNpO6b0yKwhFSKEtCtKmt14MP/vj1mP8mMGD//75B08FILJxzdI5MorXXg/ODE7m/kPtGATmF4YCiGbCVp/QUDeU2fGDjAt2/x9NpmgJW67TAYTXC7Qvp/gc1m4YmAnjqWA/+YnWdqlXpgsBgFujEDdGmlXWleMfKarzYNjCjZYAoMHxo0wLdBSLAzDb0IVmZPiackzb0vlHg4IOpP27LZ1MKSDpHw80e37u2rmDGBSRC2zQ61VRpB2Og7ZpZxv0Snjzeg34NNBIyiAJMIp24IdrSOzaJ9Qudcy7RWqwVSYdRvNyG/pM7mnDNUhcqs9mGPpWMSIaPax+cwPOn+gqm4V0iFIzoVFOVHSJjkJTC6cyYLYGQjOdO9OLikPqzHwQhXEtCGSyituGJAPa29uOPbf2wGkW01JRM/xsDNF8rko4SswFXNFbg9m2N+OiWRdigwzb0Z/rHsa9tQO1l5Bxt0+hUbakFdosBbl8IssJhEOliF01SOPadHsC25kqsqle/zsgoChAZ1ZKQ/JsISDjWOYYLQ94ZH1yuXebC2f4JAAAHCiIrkkh1qQXVpRZcu6wKF4e9aOsdR/vgBGQd/EKe7HHDZBAKuuElBSM6xRjD4go7Tk6OQeYcIudgtB8wA+ccr58fwkRAwpamCrWXM11LUmjTWYk2BSQZJ7vdONM3HvPnrabUgpaaEpztn8DmRWUFlxWJRxAYmqtK0FxVoqs29Ic7RmEyCLiquVLtpeQEBSM61uiyoa3XDVnilB1J4HjXGLxBCdcuq9LEdkkkSxIssqN7JD9kheN0nwcnuz0ISokbaF3fUgWvX8JVzZUYmgjAVWTdnWe3oW/r86CtN4M29ByIOWo9S95oH4bJIGDTotRO5ugBBSM6Vu+0wmIyQFJClB2Zx/mBCUwGZexcWQOTQf1SKUFgsFCWhGRZ++AEjnaOwZtkW/Eqhxl/s70RANA54iu6YCSa02bEVc2VuKq5Ej1jk2jr9eBs/wQC8wR0AMA5w3lehxLmRy1y2xX65bODMIkC1ixw5vR58k39d2WSNkFgWFxhC3/a53SyZj49Y5N4+nhP0m/U+WCcOnFDSCZ6xybx9PFevHZuKO2f784RX5ZXpT/nBsZxrHMMg+PhLNHWpgosqSqByBgGPAH0u/0YngjOeEyIG9DKF2MITnRxF2Se+9/nfW0DKQ091APKjOhck8uOs/3jYIxBVjhEgbIjiYx4g3jqWA92ra5FeTaP+2WAsiQkXSPeIA5fGkVvFuodBsYD8IckWIzFeVk42z+Od37/5XnvV4ch3Dr1Enm5BWf4QgQR7mskQUQPr0QDG8zlUsHB8fypPhjFcP1LIaCPZDq3oMwKs0GYroWgJlvzmwhIePp4D/rc6vciiUZZEpIsb0DCa+eG8OyJ3qwEIkC44LtrVFu/E/m0pKoEiyps85Z8VLBx3C7uRR2GcY4vmA5EInpRiRDPfUCncI5nT/QWTEaL3vl0ThAYFlXaw7MmWLgVOfUdmV9AUvDcyV60D06ovZQZBIHBYqS+JCS2gCTj0KURPHG0G+2DE1n/Xe8cUb83j1pEgWH3TS1J9Wt9my9ALyqxjp3HSnYJLrghIFyQroChi7tyu9gpksLx5LGerAWkaqJgpAA0VtrCc2oYZUdSISscL54ewMlut9pLmYOyJCSarHCc6nHjiaM9aO3x5Gw7r2tU/xe1TNy8vj6cHUnwWcDNbXhDWQUAYAxwMh+WCj3YzM5iCeuBE14MoAx+np9t4JCs4PEjPRgcT7MzrEbQu10BWFhug8kgTM8vUBQOhbIjSXujfRhvtA+rvYw5IlkSgWqAitqFwQk8eawHhy+NIhCa/2RHJvwhWROjFNQynR1J8PbpZbH7sYiMo4q5sVLowEZ2Dlb4UQlP1B1yd1IpIMl47EgXRr3B+e+sUcVZqVRgRIGhodyG84MTYIyB83DfEYH6jiTtZLcbvqCE61qqNbdFYjIIkBVedCPFi12vexKHLo1h1JvfT7xdo5OoKbXk9Tm15Ob19fje3rPoHPXFDErGWSn8mz4NC4t/asmEcMuRm1aHt9Bb3Qac9ZXAF8xdMOkLyvjD4S585MoGlFrUGxSaLgpGCkSTy47zgxMQBQZJ5uHsCONFMe0xW9oHvZgM9mHnqmqYDdqaXyEKDAITEJIp61XoRr1BHO4YQ8+YOoWJHSNebF5ceE21khXJjvz9I0fn/BsDsHFRGSylydWEvOKx4sNXNKAawLUKx6URH9p6PWgfnMjJdvpEQMIfDnXhI1c0wG7W1+WdtmkKxMJya7izp8CmOwBS7Ujqet2TePpYLyY01IskgjEGk0GAkWpJCpI3IOH1c8N45kSvaoEIAAxPBDEZ1N7Pfz7Fqx0RBZZS99PusUlcmBrYKQgMTS473r22Dp++thk7V9agvsyazWUDANyTITx6uAv+HG/pZRu9qxUIgyigoSL8gx0pZOU8nCEhqRn1hXuRjGh0/1UUGMwGgWpJCkRQUnC4YwRPHOvB+cFx1U/Dcc7RMVLchayxakciWRGrKbWs6avnhub8NzUbRKxZ4MRHrmjAnVc34armSpRZs7e1MuwN4rEj3QhK+tnapWCkgDRW2gFgRs0DZUfS453qRdI9ps03ZcqS6J+scLT2uvH40W6c6vZA1lBNUNdoYfSuyEQkOxKRalYkYngigNZeT9x/j7Sh/5vtTfjIFQ1Yu8CZlW3ifo8fTxzthqShn6tE6J2sgDRU2GAQGBhj07UikWJWkrqgpOD5k304N6CtXiTRKEuiT5ETMocu5v6ETDq6xyZVz9CoLZIdiUgnKxJx4PxwUkFBfZkV71hZg89c24z3rK1Dc1VJRr/b3VNjAvRwDaBgpIAYRQELysORvBj1A6yHH0StUjjH/jMDONY5pvZS4qIsiX70uf145ngvXj03hAl/mpNh8yBQ5Ed8I25eX4+qEjMMaWZFIiYCEo51jSV9f1FgWFbjwPvW1+PTO5px/fLqtE84XRz24rmTvZoPLvVVbkvm1eSy49KwF4LAwJTwMd9IdkRrR1b15M2LI/AGJGxbUqnZ2T904ka7Rr1BHOkYQ7eKhamp6hjxodaZ/QJLPREFhk9f24QzfeNpZ0UiDl4Yxep6JyzG1L6O1SRiQ0MZNjSUYcQbRFuvB229npSK7M8NTOCF1n7sWl2b6rLzhoKRArNoaoqvrHAIDJCnrkmSokBggmYvpHrQ2uuBLyjj+uVVMGg0CxHOkjDqS6IR3oCEY51utA9lv3V7rnWOTmJLk9qrUF+J2YiGitiNzlIRkGS8dXEU1yxLv1V8hd2E7UtduHpJJbpGJ9Ha68G5gYmkftfbej0wGQTcsLw67efPJQpGCozJIGBBmRUdI77poAQAwMPbNQZqhJaRcMpTxk2ralL+hJNPlCVRV1BScLLHjdN945oqTE3FyEQA3kAIdrP+GmhlkyeL22lHO0exvsEJR4ZNyRhjaKiwoaHChhtXKDg3MIG2Xg86RybBE0zXOdY5BpMoYPtSFzjnmDx8GKO//W9MHjkCLkkwLliAsltuQel73g3Bkt/Gd4zrIFz3eDxwOp1wu90oLS1Vezma93b/OP50NjzCOiQrl4/3MsAkUnYkG5xWI3atqdVFp0NJVuhUVZ4oCseZ/nGc7Hbrrs9DLDuWVWFFXfG+54ZkBQ+8dC6rX3NVXSnemaPtkomAhDN9HrT2jmN4In7n3m31NtT+x1fhO/AGIIqAPPWzyhjAOcSyMiz88Y9h27Qx4zUle/3WZq6ZZKShwjZdgT2jToRTMWu2uCdDeOpYD4YS/MJrhUEUYDYICYd/kcxdHPLiyWM9eOviSEEEIgAKZjx9usb92W/+1jZPoJCJErMBmxdX4ParFuOjWxdh46Jy2GbVunBZxnP/+X9x5Exv+AY56md1KjchezzouPNO+Ftbc7LOWCgYKUAWo4j6snCKTWBsRiZEnipoJZmbDMp45nivLnoyMMZgNogwUBFz1vW5/Xj2RC9eeXsQ4xo+IZOObvdkUTdO9Exm/78nB8er54ay/nVnq3ZYcF1LFT51TTPev2EBWmocMAgMgbY2SH19eKuqBe2ldbEfrCjgoRD6vvmvOV9nBAUjBarJVTL9/yk7kjshWcELp/pxtn9c7aUkhbIk2TPmC+LFtgHsbe3L2SddtYUkBX1FfMQ3F5kRALgw5M1bQ8XoNvSf2tGE9a88gSr/GADgQN0adJTEKWhVFEweOQL/2bP5WWdenoXk3eLKy1s1AsP0vBogHIxQdiR7FM7x8tlBHOkYVXspSaEsSWZ8QQkHzg/j6eO9ujqqm66OYa/aS1BNNotXZ3v17cGcfe14xIF+NJ46iHdeehPvb38Fa4fO43jVUvTYKmM/QBAwsf9PeVkbnaYpUBajiDqnBd1jk2CMQWQMclQAIikcRjpZk1WHLo3CG5Cxfal2e5FEM4gCRIEjKCsxR6WTmYKSglO9brT16veETDq6RrU5EiEfcrnt1uv249zAOJZWO3L2HLPJE5cDy5KQH+uG27FuuB2j5hJwzPjMGiYIUMbzk/WlzEgBa3TZp///7IZnCmVHcuJ0nwd7W/t1Mw+CsiTzUxSOtl4PnjjajZNd7qIKRIDw4MhCq4VJlmcyt9OLXzs3nNf3YdEZ+zRLeWBibiACAIoS9zHZRsFIAWustE9/QmezClkBGqKXKx0jPjxzoldXJyqoliS2Qjwhk47OIp3im8ttGiAc6J3sjj9EL9uMtbUwr1gBCEle+hUFJTe+I7eLmkLBSAGzmkTUlJqn/x4rO0INsXJjcDyAJ4/1wJ2DavxcoSzJZX1uP5470VeQJ2TS0TlSfHUjisLhDeQ+AH2jfTiv3ZIrbv9rQEni+UQRti1bYG7OTxteCkYKXNPsrZpZ1xk6WZM7nqleJAPj+jqNYBCFqeZ4aq8k/8Z8Qbx0OnxCZmhCX//dcqnH7S+694rxgJSwm2m2eIMSjnSM5fx5Ipzvex9s27Ylzo4IAgSLBbX3fi1v66JgpMA1Vtpn/F1kMbIjRfYmk0/+kIxnj/eiY1hfpy4EobiyJL6ghDemTsjooW9Mvkmygh53cW3V5DMj9talEUwG87MNyIxGNPz4ATh27gzfIEY1RZsKUAzV1Vj8m1/DvGRJXtYE0Gmagmc3G1BTapkeBz5jXs0USeEwFclFRw2SwrG3rR/bl1ZiRa2+WmsbRAEC4wgphXniplhPyKSjY9iLhvLMB8bpRa6LV6MFJQV/vjCM6/M0xE6wWrHw//wn/GfOYvTh/8bkkaPgwSCMDQ0ou+VDcNxwA5ghv+EBBSNFoMllnw5GGGMQBDYjG8J5ODsiUECSM5xzvPr2ELwBCZsXV6i9nJQIAoNZEBGSlYJJ1SsKx9mBcZzoKowZMvlQbEd8810rdKLLjY2LyuG05m/elWV5C+ruvTdvz5cIbdMUgegjvsDcrRqATtbky5GOMbx8dlCXW2PGqVoSvbs0dULmzQvFfUImVZ7JENyTQbWXkTeeHHVfjUfmHAfO575NvFbp/52FzKvEbECV4/KpGkGYe8yXc14wn3q17mz/OF5o7ctrBX22CAKDxSjOOZmlB/2e8AmZl+mETNo6hosnO6LGz8jpvnHdFbxnCwUjRaJpdnYkxsWEgpH86RqdxDPHe+EL5vfTV7boKUvi9oXw0ukBvHCKTshkqpiO+OZiSF4yXsvDED0tSuvd5IEHHkBjYyMsFgu2bt2KgwcPJrz/2NgY7rrrLtTV1cFsNqOlpQXPPvtsWgsm6Zm9VRPrgy1lR/JraCKAp471wO3T56d0rWdJJoMy3mgfxlPHe+iETJb0efy66S6cqYmAOh8ULg37dHf6LhtSDkYeeeQR7N69G/feey8OHz6M9evXY9euXRgYGIh5/2AwiJtuugkXL17E73//e5w5cwYPPfQQFixYkPHiSfJKLUa4Si5v1UQKWWeTkmmGQ7Jm3C/hyeM90wXGeqS1LElIVnC0cxSPH+3G2/3jRTv2YDKkwBeUk/rjl5L7vZcVXhSFrL6gpGod3atFmB1J+TTN9773PXz605/GnXfeCQB48MEH8cwzz+AXv/gF7r777jn3/8UvfoGRkRG8/vrrMBrDVcKNjY2ZrZqkpdFlx1DUqHNx1qkaAAAPv+Fo9dNuIQqEZDx3ohfXL6+ek8HSC0FgsKh84kZRON4eDJ+QyVfPBq26OOTDbw52pPSYT1y9GHVl1nnv1znq0+3PabLyeaw3loFxP870jWN5bf6G6KktpY8zwWAQhw4dws5IsxQAgiBg586dOHDgQMzHPPnkk9i2bRvuuusu1NTUYM2aNfjWt74FWY7/ZhEIBODxeGb8IZmbXTcixJhXA4SzI8X6aVItksKx7/QAWnv0/bOuVpakYzh8QuZge/6aR2lZVakJqQzlNokM5VGZ00SKITOihQLn188P6fLUXbpSetcYGhqCLMuoqamZcXtNTQ36+vpiPqa9vR2///3vIcsynn32WXz1q1/Fd7/7Xfzrv/5r3Oe577774HQ6p/80NDSkskwSh9NqRIXdNOO2mBkQTsWsauCc4/XzQzh4YUTtpWQkn7UkAx4/9pzsw5/O0gmZaHaTAVc2Jd/P5uolLlgMyV0OJvwhjHoD899Rx3I9IC8Z7skQjne71V5G3uT8I4yiKKiursZPf/pTbN68Gbfeeiu+/OUv48EHH4z7mHvuuQdut3v6T2dnZ66XWTRiFrLGuGbInFN2RCXHu8bw0pkB3QeEucySuH0h7D8zgOdP9WGwSI9Czueq5oqksiMmkWFzY3lKX7ujwKf45rvHSDx/bh9GMMl6Hr1LqWbE5XJBFEX09/fPuL2/vx+1tbUxH1NXVwej0Qgxqv/9ypUr0dfXh2AwCJPJNOcxZrMZZnNyKUOSmqZKOw5fGp3+O2MMImOQZwceU9kRQyq5XpI15wcmMBmUsXNlDUxJfmLVIkFgMDMBITlqQvRwOxCcSP6LmEqAymYA4RMyx7rGcH5woqhS2OmIZEfeaE+caUslKxLRMeLD+oayDFanbWod651tMiTjrUsjuHqJS+2l5FxKwYjJZMLmzZuxb98+fOADHwAQznzs27cPf/d3fxfzMdu3b8dvf/tbKIoCYWoIz9mzZ1FXVxczECG5VW43ocxmwpjvcifFWPNqgEghK49ZV0Jyr2dsEk8f78Gu1bWwm/U7uYExBpMh/DMWGjwH/GRbyl8j9JnXcSpQgbbe8aI5WpoNVzVX4M0LI5DjxG3pZEWAcIFlUFJ0HSgnMq6RzAgQ7tq8fmGZrt8DkpHyT9Lu3bvx0EMP4Ve/+hXa2trwuc99Dl6vd/p0zR133IF77rln+v6f+9znMDIygi984Qs4e/YsnnnmGXzrW9/CXXfdlb3vgqSkyTVz2FW8Y74AtYlX24g3iKeO9WDUq/823KLAYJYmICD5YELhDGeVBXjixABOdLkpEEmCwjkCIRkT/hBCIRmr6+MPZ0wnKwKETy4Vcu8WLdSMRITk8BC9QpdyqHXrrbdicHAQX/va19DX14cNGzZgz54900WtHR0d0xkQAGhoaMDzzz+PL37xi1i3bh0WLFiAL3zhC/jSl76Uve+CpKSx0o4jHWMzbhMYg4K5gYeicHDKjqhqIiDh6eM9uGlVLWqdFrWXkxHGGExMhswVhOZ5++lQqnGEL4GH2wGJguKQrCAoKZBkBSFZQUjm4f9VFIQUDknmCEl8TsBWURI7A51uViSic8SH5qqStB+vVQFJ1lydxsluDzYtKkeZzYRQTw+koSEwsxmmpiYIBbLDwLgOqhQ9Hg+cTifcbjdKS/U1gl2rfvdW55x90YAkI0Y8AkFgMGqoqVWxEgWG61qq9H0B6DkK/PQ6AADnQAgilFkJ2gHuxGFlGQZ52eUbt90FlNbnb515IiscQUlGSA4HEcGoIENSFIRkTAUf6ReUc85xbnAcF4ZnFvpe31KF7Usr0167zWTAx65anPbjtWpoIoBfv3FJ7WXMxIGFvedwxXO/xuSRI9M3C85SlN/6V6i443YYXNqsK0n2+l3Ym1AkriaXHcc6x2bcFrOQFeHsiCJwCJQdUZWscLx4egC+oIw1C5xqLydjjAEmXM6SuLkNR5Vl6OBVai8tIwrnCEnh7EUkkJgOMhRlKvAI30fJ4WdBzjk4whmQd6yowc9fuzT9WcMoABsXl2X09X1BCYPjflQ59J2tm00rxavTOMf4iy/haGsrajovIDrkUNweDP/85xh79FEs/tUvYV6yRLVlZoqCkSLVWBkjGIlTyAqEL4QCnazRhDfahzERkHBVc/qfarVEZBxQJLyqrMUo127WZ/Y2SXBGBoNDUmJvk6ghEogYRIYmVwksJhGbF5fhrUtjAIBtS1wwiwI4z2wLtmt0suCCES0VrwKA7+BBBFpbAQBHXEtxU8dbM+8gy5BHR9Fx5yfQ/MzTEB367NpKwUiRqnKY4bAYZvziRQpZYx2ZVBQOhfG4ha4kv052u+ELSriupbogWveLAsdN7BCelbfCDxMUDnCw6T+5cnmbRIEUncGYqsWQsrBNkk+RIAQIByKNU4EIAFzZWI5Dl8YgMmDT4nCtiMIBAekHJB0jPmxclH7diRZpqXiVB4PwHb68LTNgLUe33YUF3lmza2QZ0uAg3I89hoo77sjzKrODgpEi1lhpx4lZHf7EOIWsQPhkjakALnyFon3Qi8lgH3auqobZIM7/AI0rZZN4h3gYe+QtABMQLmAK/ywqAsAZC98yT1CgcI7g1DaJPGObJFzoma9tknyLDkQEgWFxhR120+WfC5vJgPeur4NJZDNO0GQSkAyOB+APSbAYC+dSoqXMiP/MGWDW6JQjVctQ7x2KGaKP/H+/Rvntt+vywEHh/ASRlDW55gYjgsDAFBbzDZ9zDkWh7IiW9Lon8fSxXuxaU4uSAuhDUMU8uF44iheVTTNCYoExBJXLp0iCkoyAHAk65KgAgxfd5OnoIAQIv1aLym0oscz9eWipLolZ+8V5uIYnnefuHJnEshp9bg3EoqWaEWlwMPwfJur92G0uwYXSOjR7emfemXOEOjvB/X4w6/wDD7VG/+9eJG3VpRaUmA2YCMz8JCAKDFKcLkmUHdGeUV+4F8mu1bVzZg9pXZCLmOQm+GGGj5vhgxmTMMHEAzimLEEIBoRggNQdADePxf06kQtyASU65jU7CAHC162FZVY4bcbUvhaQ9geNrlFfQQUjWsqMIE5g7TXGDza4pKH1p4CCkSK3uNKOUz2zsiMJ3o8451OdWSkg0RLvVC+SnStrUJ/EGPhcUhQOX1DCZEiGLyjDH5IxGVQwGZIx2TeJSWkL/DBjkhshIf72kp350cMjZwcSRxmMsXDaeurHstCDk3iBSJ3DgvI4fUWiHxtrIBVHeIsr1VNzhTTFV1Y4vEHtXMxFx9yjsCUhH1YNX4h5f6GkBILdHvPftI6CkSLX5JobjCQqZAVAwYhGBSUFz5/qw45lVVhand1TKZzzqaBChi8khwONgIyAFA40/JHbQlLihlE+A8CTO5bcwAYR5AYMoQwQUvukX6jBSawgBAgHIlUlZlRl2BSPc4CnWD/iD8no9/hRU6r/UzVam/xsXrkCvjffnHHblX2nIcb6KRBFlH34w2CCPntCUTBS5GpKzbCZDPAF527VxAtGKDuiXbLCsf/MALwBKalBZkFJgS8oXQ4mguGAYzIUzmj4Q0r475KUnQu53QVc/yVASm4EfTPnCI4yeJTMLnTRwUmkHkpvwUnMQISF461ymwl1WcqIpVPQ2jHsLZBgRDtZEQAQS0thampC8OJFQFGwaLwf9b4YreEZAxhD+W1/lfc1ZgsFI0WOMYbGShtaez0zbhcYA2OxC1kBQFIUiIL+T3AUIoVzvPL2ALpHfVi9oBT+kDIdZPilcHARDj4kdWYP2ZPvFMkAtJRytPZ44Atk50IRucjqJTiJlw2JBCKlVgMWlCcfiCTzLaYakHSNTeLKpFegXW4NFa9GON5xI8Z+/wew0WFsHjgz9w5TgciCb/87TIsW5X+BWULBCEGjyz4nGAESF7KC03ZNPnE+1ftCVqJah/Pp46rT80okPt247kSXB6+8PYRlNSW6PgElCgwraktxstuNoCTP/4AUxQ1OeHIX7lyaLxApsRiwqMKecp1HMs3OUjlhMzQewGRQgtWk70uK1jIjAMAsFpR9+Baseu0Z2NsVcMYAQQj/B1IUmJcvR/U//m+UbN+u9lIzou+fHJIVdU4LrCYRk8GZb/SiwCApADjgCY3j5c6XcW3DtSg1hivnJUWBwARdnmnXistDzy73w5DkqMBDmerymeaguBFvEK0941hR54BBxx10jQaGFXUOnOrxQM5xh9Pp4GRWvQnyGJzEDUIAgIUzl1ajmFYgkvQakFpBa8fIJJbX6vtUjZYankWrqizFTf/6JfB/+izG//hHSIODYGYLbFdshnXdOrWXlxUUjBAwFm6QdLovRnZkal7Nnzr/hPaxdnBw3Nz83vA/TmVH9HyRy4Vw0y358sCzqM6eUnSQISl52RYY94dwssuNVfWlMBn1WdwGAFaTiJYaB870evLarCxeMWyugpNkAhGzQUBjpQ2GHGe8Uilo7Rzx6j4Y0WJmBABuXFEdzm6WlaHsllvUXk5OUDBCAACNLlvsYERgOD92Ae1j7QCA9rF2XPJcwuLS8LROmXOIGc630IOE2ySRbIYyc5tESyZDMk50u7GirhR2s35rfUqtBjRXl+Bc/7hqa8hVcJIwCJl6PoExGIVwnZchg0nasQ/3xpZs/Uj3mD/jWTdq01LDs4hVdaWqH9fPBwpGCACg3mmF2SgiEJq5VaNAxksdL04XszLGsK/zRXx81R0Qmaj77Mjlse1zZ5NIWdgm0ZKgpOBUtxvLax0pN8XSksoSE4KSHR3DXrWXAiDzY8TzBiFTX1tgDCJjWOyyw5Rp+/9UohEkF5AEJRl9Hj/qnPq8cHLO5zSAVJvFKGLHMn1PsU4WBSMEQGSWhQ1nZ33iPNj7JjwhN8DDb36cc4wHxnF44BCurNkCIFLIqp1PRIrCEZTn2SaZ+jctnp7IJVnhaOv1YElVCapKzWovJ211ZRb4QwoGPNpruJVKcJJKICJMBSJWY+aZrXR+7JMpaO0Y9uk2GPEGZc3NKtq+xAWrSb+ZzFRQMEKmNbnsM4IRd9CNAz2vg7HIxNLL70R/7jmI5RUro4pZOYw5zI7wqeFnkXHt0QFFaEY2Q5vbJFrCOXBuYAIBScHCCn1eOACgqcqKkCxj1BtUeykJxQpOFM6RzI9p+NQmm6rrss0YfJdvybSM7xr1YSsq87eoLNLaFk1tqQVrFybXILAQUDBCptWXWWE2CAhMddB8sePFyxN8GQf45TchBeGi1kgxq6Jw8DSyI/Nuk0QFGSS7Okd8CEoKmqttSClnrxkMS6tL0NY7jgmNnoKIhzEGcZ7MyXQgAoaGMmvMwXf5Nt8JmxFvEN5ACHaz/rYBtVS8ysBw48pqtZeRV+r/dBPNEAWGhgobzg1MoN3djnOj56L+VQFwuWCOcz6nmDWSHVGUqNMkSnSQEb1VUpzbJFrT7/EjKCto0WkvEkFgWF7rwKluN/yh7Pcgyab4rdzndoeNTqUscJqzXuMTbz5Nco9NfMKmY2QSK+v0F4xo6VjvugYnqh3672ibCgpGyAxNLjvODUzgpc6XZnRgZSz8iYjLRoCL4BDBYMAL5w7jHYsqLwcXivqNokhqRr1BnOrxYGVdqS4LkQ0iw4q6UpzqdiOU4x4k6UqqNgSRLZnLvU5qS8xwlVrC/a00FLknKmjtGvFhZd3cAW9ap5W5NHaTAVcv0edWVyb023SA5MSCMitMBgE2gy3Gv3JwyQklVA0eqoQSKgOTnRj1hTARkBCQFEgaesMkyZvwSzjR5YY/pM2L+XzMRgHLa0shaGxIWKQ+JNnfiuhApKrEjJoya/gUjcBgFAUYhPD/z7RYPBu/pQqPyuRE6XZPxp1rpWWeSW1s0+xoccGc6WkpHdLWby5RnUEU0FBuw85FO+e8YzHGIZqHwYyj07e1VCybcR+F85hvUET7/CEZJ7vcmNDQ3nkq7BYRy6pLNHGqi0/9HqTymxAdiJRbjTF7S7CpUzWGqOAkMkdKDbF+1UOSgl4NnnKajxYyIw3lNqyo1V9WKRsoGCFzNLrsqLJVYVPNpllvcuF3HtEwAdE0gMWl9bAb546ql2mjRrdCcrgXyZhX/TfmdJTZjWh02VVdw3RBagqPiQ5EHGYDFlbEykzGelw4UxIJTsQUgpNsfGiIFLTO1jnsy/hr55tH5SBcZAw3riiuotVoFIyQORrKrTCKArYv2A6LeLmIKvr9zWZm2LV8LRwxKvy1tr9NUqNw4HSfBwOegNpLSUt1qRn1KUyxzZZUt2QiogMRu0nE4sr0580IGQQn6Yr1+945qq/MiD8kq15vtHlxOcrtJlXXoCYKRsgcBlHAwnIrzKIZNyy6Yda/ht90diy8Fg6TFcuqSmJWfSuUHdE1zoHzAxPoHNHfJ1wAaKiwwZXH0wipbskACE/ejQpErEYRTZX2rE7CjhucILuF5nxW/ciYL6iJbY9kqX2SptRixJamClXXoDY6TUNianTZcWHIi9WVq3F04Ch6vb3hdvACQ52tFisrVgAIv5E2lFtRYhZxadgHOWr8uoLkJ34SbeoamURQ4miusmmiFiMVzVU2BCUFnsncNUVLKwgBwoEILgciJlEIByIZzJtJhsDY9LHhSNATCSQyDU6iT9hII6M48p2nsOjUQfBQCMaGBpT95V/Cvv1qMI0VGQPqF69ev7wqo1lDhYCCERLToorwRFBJ4Xjn4nfiV6d+Ff4HznHDohvn3L/cZoLFKODCkA+TU/0eUhk/TrRrwONHUJLRUuvI6qf2XGOMoaW2BK09HvhyMHMkW4GIQWBodtlhNOT3YsR5uLdQdE+TcAF6+t+bEpIw8cLz8J86hbN+N6q6jwMA/G1tGH/uORgXL0bDT34Mc3Nz1r6PbFAzi9NcVYLmqrm1d8WmuEMxEpdRFLCgPFxEFylmBYCN1RvgsrhiPsZqNGBFjQMVtvC+Z7ziNqI/Y74QWrs9COmsE64oMKyoLc18sFyUdGtDAMwJRETG0FxlhzkL82ayYXpbRxTmbOvMhysKxp54HL7WNgDAoMkBOfJIOfwBJdTVhYu3fRTBjo4cfQfpUat41SgKuH55cQzCmw8FIySuJtfliv5rFlyDaxdei+0Ltyd8jCAwNLnsWFRum26URgrDREDCya4x+IP66kViNDCsqHNkZQsko+2MqYF3kUBEYAyNLjusRnUS1Ml8H9HBSaTHSbxsp/9UG4LtFwDOoQCQBQGD1rKZd5JlKBMT6Pvmv2a6/KxSKzOypakCpRb9davNBQpGSFyLKi4X05lEE7bWbYVZTK7au8phRkt1CYyiQAFJAfGHFJzo1l8vEqtJREuNI+1tw4yyIcB0IDL9V8awuMKOErN+dsojPU5mN2CLfF++I4dxucU8AwfQZ4/RSVSW4X31VQS7uvK19HmpUTNSYTdh86LyvD+vVlEwQuIyGYQ5jZdSeSsvMRuxss4Bu1mkRmgFRJrqRTIyoe1pubOVWg1ork5tbz7jIASYE4gAQEOZFaVWlQORDH8no4MTwTsB3tcDgctTYQjAwdBri3OxZQzje/+Y0fNnkxqZkRuWV+tyHlSuUDBCEmqa1UAq1RMVBkHA0qoS1DiLa+hToVM4cLZ/HH1uffUiqSwxoaEyuYZi2ThhEisQWeC0aqKfRDY/Hih+f7gWBoDAOUSuQOAKfAYL/LGyqYIAeWwsiytIX0hWpovu82VFrQMNSTa2KxYUjJCEFlXY5ryZphqQMMZQ67RgaXV2eygQdXEOXBicwCWdddusL7OiujR+cJyVbAgQMxCpLbXA5TBn+pU1h5nmBhwMgFGRYZZjZNAUBUKJup1yI8bzvOVoMgjYsYyKVmejYIQkZDGKqC+b+cad1rY7BxwWI1bVO2Aza+PkAMmOntFJvN0/oautuKYqW8zsRFayIQj/jswORFx2E2oSBEF6JpaVQSwvx+yN3PLAeOytXUVBybXX5mNp88r3Fs3VS1yw66hWKF8oGCHzmj3rI93chqxwmA0iVtaVFuSnw2I2NB5AW884ZN1Ma2VYWl0C+9RJhqxlQzCzq2pEWZzBd2rLVgDJANivuGLO7RX+8bl3FgRY1q+HZfnyrDx3pvJZvFrtsGD9Qmfenk9PKBgh82qstM94c82kE6ekhBuhNbnsaHLZQLs2hcM9GcKpbjdCkj4CEkFgWFHrgNkgZK1+IlYg4jAb0FChzQ622fwvZV2/HoaaGoBdvqyUBzwz7yQIYAYDar/y5Sw+c2bylRlhCA/C0+LPgRZQMELmZTGKqMtSAaqi8OlPYy6HBSvqHDAb6cewUHgDMk50jWEymN+CwHQonIMxYFltKYxZ6EESKxDJdPCdnjCjERW3/RWM9XVTNwio9E8FI+EXB4LdhoafPQTr2rXqLXSWfM2lWbOgFLVUyB8XbVyRpDRW2tEzFp7Emen7qqRwGMWpKaVmI1bXO9E+OIExn34Ga5H4ApKCk11urKgrhUPt46txRNqeA4DFKGBZjQOn+8ahKOk1dIsViFiMQtYH32mdYLWi4vbbEbx4EcKhN2AbPQkuyzDW16P8Ix9G6c3vg6iRwtWIfHRftRpFbF8au3M1CdPmOwXRnCaXHQfah8PD8jL8WorCoQiX59aIAsOyGgd6xibRrbPR4yQ2SeFo7XFjaY0DlSXqH2ONiA5CopVYDFhaVYK3B8ZTrqOIFYiYRAHNlSU5H3yXKc6RfhFYHIwxmJua0LxlLVr+UzvbMfHk4zTNNctcsGik5b9Wafs3hWiG1SSipjRcdMqmJ3+mL1ahY32ZFS21JTAYiueTZCFTOHC2bxy9Y361lwLOOWQldiASUWY3YlFFap/aYwUiag2+05oqHRSpKwrPeTfh+jIrVtdT0ep8ivu3haSksfLyG3Vyo7PiUxQes02802rCqtpSlFgoaVcoLg55cXFIvV4kCudI9pBPjdOMOmdyp15iBSICC2cRtTL4bj65LDWu1kEwMhGUkJ0zVLEJLFy0SuZHwQhJWnQ31mzU40lxJsCajSJW1DpQXar9NzOSnN6xSZztS30LJBPJZENiaai0obIkQaEhixeIhE+J2Uw6CqRz9N9DFBhcJdr//fVM5rZObUNDmS5eBy2gYIQkzW42THeuzMZGCuccSpyPrIwxLK60o7mqOE4iFIPhiSBae8Yh5aEXSSrZkFiaq2yxp6lO7VDODkQYY1hUYUOJmSawAkCF3ayLuSu5rBcpMRtwVXOMQYEkJgpGSEqaprZqsnVWfr4LU2WJGavqHbCY9JH2Jol5JkM41eVGMEe9SNLNhszGGMPSWsfMLEecQAQAFpZZ4bTqLxDJVViohy0aILeZketaqmAq8rqhVNArRVLSVDUVjGTp60UuHolYTQasqitFuYZOZZD0+YIyTuagF0mm2ZDZDALD8tpSmAzi9JyZWIFIvdOKCg0MvtOSKp1sseYqM9JYaceyGkdOvnahomCEpKTEbECVw5yVmpGIZFqIiwLD0qoSNFRYs/rcRB2RXiTZaMWtZCkbEovRwLC81gGjIXZmrsZh1sWpkURyUcejm8xIDhqeGQSGG5ZT0WqqKBghKWty2bPa0jiZ7EhErdOK5TWOoj82WQgkhaOtx42h8RhTXZMwPU8mhyUoDAwOixEt1Y45tUuVdhNqkzx5U0xsJgMcseptNCgXmZErGivgtOnj+9cSHZV9E1UNnwcC4aFXi0MKDnp8YHKCC4HBBFgrkv7yssKT7lTpsBqxqs6B9iEvxvM45Ipkn8KBt/vHEZBsWFCe/IU910EIEA5EIj+TTqsRTa4SnB8M/w44LUYs0ODgu3Rku++ZXrIiQPbn0pRZjbiyMfn3PXIZBSNkfsPngR9umv6rE0ClvB19Sjl4ouTa1X+XdEASyY4kG5CYDCKW1zjQNTqJPrf6TbVIZjqGfQhKCpqqbEh0aeRZrguJJzoQiahymBGQZLh9QSyq1ObgOy3QS72ILyhl/WTXDSuqi6r9fzZRrpvMLzB3DHgj65v/cVJq6XdJUVLav2aMoaHChqXVxTX/o1D1uf040zsR97h3tgtU4xHY3EAkYt1CJ25cVVNQx82znWGq0klfjWxv0SyrdmBxpbbm7ugJBSMkLc2sDyzbBwN5csWss5XbzVhV54DVTMd/9W7EO9WLJKohXraO6yZDYCxuoFFqMWBbcyWubq7EgjJb7hejQ4LAdFPQm81jvUZRwHXLq7L29YoRBSMkLU7mRQWbyPrXlTlPq7rfYjJgVW2ppoaykfSM+0M42eVGMKTkLRsCJA5EbCYB25dWwmwUwRjD9mWVqLDr46KbT2U2EwwaHw4Ykc1pvVc1V6LETFUPmdDHTw3RpCbWk/0vmmZ2BAh/KmuuKkFjpY2O/+qcLyjhWNcYJvzZ7UUST6JAxGwQsH2pC9aoBmhGUcANK6p0c2okkWzGetU62aIBsnest7LEjI0NZVn5WsWMghGStiWsFyZk/5x+OCWf/ltkVakFK2odMBnpx1uPOA+PLgtKClp73HD7cjs/JFEgYhQZti+pjBl0WIwirl9epZuheHFlcf9LL8WrQPZqRt6xoloXre+1jvJKJG0VwgRuZ39EkBsQROSPESEYEOQGBBpMCJVVICjLCIQUBGWOoKQgKCsISQoCkoKgJMfMhMgKh0FM/xe8xGLE6rpStA95c34xI9kRCUKiyQrH6T4Pml0lObnQJQpERIHhqqZKlCXorlpqNeK6ZS7sOz2QdkZPbVnNjOikXgTITs3IqrpS1BfIEW+1UTBCMmJGCDITYIYEMyQAU8dsGYByI1BXOu/XkBWOkMQRlBUEJHk6YBEYg6RMBTCSEhXUhAOZgBQOapQ4n+wMooCWGge6x3zoGaXjv1oVKwiZ+e/A+cEJBCUFCyqy98YvxmnvDoSn8m5pLE8qAKoqtWD7kkq8cm44r1OJtcZsFFFm00/NVqaZEYtRxI5lVLSaLRSMkIwwBlh4ED4kGLk+D1FgEE0MFgiI/pE0igIsSaTAQ3IkyxL+EwloAlNBzZoFTnSN+nCs043JkAxFViAp4RMa+ZggS+KbLxCJ1jk61YukOnEvkmQkCkQAYGNDGepS+MTbUGnHpqCMQ5dGM1qXnrl0VC8SeZ/IxPYlLlhpgGfWUDBCMiYyDjMPIoDsfioKyQqMojBvDxGjKMAoCkDC98JK/MXqWvzp7CAGJwJRt3NIcrg1uSxzSJxDngpWJEWBJE8FLTKfuo1DVpSp25S8nfQoNKkEIdH6x/0IygqW1ZSktU/PwCCwxFOn19Q70OhKvV/EirpSeEMSTvfM7cujdZzzjJu4VeuoXiTT4tXaUgvWLJg/60uSR8EIyQoTkyFzGRKy+0khKClZ+/RRYjHiXWvq8OcLIzjT75m6lcEgIlyfksbBCEXhkDkuBzBRQYskK9PZFzk6sIm6rRiz+ukGIhGjviBaez1YUVs6XVfE/QH4T57E5KmTUMa9YCYjzEuWwLppE8TyMgDJBSItNSVoqU3/IrN5UQV8AQUdw960v4YastESvkZH9SKZbNEwMNy4opo68GYZBSMkaywIwgsLeBYnXUiKAlmZPzuSLEFg2LakElUOE95oH854m0YQGAQARjG9gElWOGQFkKa3jpRZAQ2HzPnlf58R4Ogrksk0CIk24ZdwstuNFXWlYO1n4H7iKXBp6gIz9STS8BC8f34D9iu2wLHzHRDExFszjZU2rFngzHhtVy+pRCAko99TPHVKjDFUOdLfqs23TIpX1zU4UV2qn+9VL9IKRh544AH8x3/8B/r6+rB+/Xr88Ic/xJYtW+Z93MMPP4zbbrsN73//+/H444+n89REDWZHUneL1I9MRvZLzCVZefqAJMNmym7cvLTagUq7CS+dGczJGPFkiQKDKIRn7aQq3Jk0aotJia6FuRzUzPi3qK2ofMYy2QxEIvwhGceOnEH1/udgk6S5x0KmvkHvwYNgigznu98V92vVl1mwcVFZVtYlCgw7WlzYe2oA7sn0JhLrjdNqhElHk7TTzYzYTQZcvaQyy6shQBrByCOPPILdu3fjwQcfxNatW/GDH/wAu3btwpkzZ1BdXR33cRcvXsT//t//Gzt27MhowUQFlUuAzx+OOaNmNgMAk8QRNNiBiqasPH2kPiPbnR3L7Wa8d109Xj03hI4RfaXVgfCn0ZlbTKkFNEokKzNdA3M5aJmxpTSVlVGi62iSnJqbiyDk8hdXMH6yFeOl9Wj09MER8sW4DyCCY/Ktt2Bbvw7GBQvm3KXKYcaVjRVZTbubDSJuWFGFF071wxfU/mRpnuE+jZ6O9ALp14zsaHHBnMYHBzI/xlM8i7Z161ZceeWV+NGPfgQAUBQFDQ0N+PznP4+777475mNkWca1116LT3ziE3jllVcwNjaWUmbE4/HA6XTC7XajtJSKhvTAF5Syuo0gCizr2ZFoJ7rcONI5GveYMJmNQ1IAJU6NjKSETzZFBzWRwt9QlrIyUv8AJk8cn/77gvFBVAY8039nHBAioZDAYF23Fs73vX/G1yi3GXHNMle4ADoHRiYC+OPpAYQkJSdfP1sSDQdMxvZlLqzIoNYm3x55swO9KU77bii34UObF+ZoRYUr2et3Su/uwWAQhw4dwj333DN9myAI2LlzJw4cOBD3cd/4xjdQXV2NT37yk3jllVfmfZ5AIIBA4PKJB4/Hk+De8+Oc42evXMDAePI/fO/fsCAr+8fFymIQ4QvJWeu7ELmI5eqisXahE9UOE/50dgi+kPY/yaqPwSAAENicM1RKElmT6S2mqWyMHF0jMytTIyvhraXLgU248Dc0NIDwx/nwk3U7qiCJBtT4RmYGIuFFwd92ekYw4rCEU+65+pkCgIoSM3YsdWH/2cG404gLgZ7awAOpb9OIjOGGFfEz/yRzKQUjQ0NDkGUZNTU1M26vqanB6dOnYz7m1Vdfxc9//nMcPXo06ee577778PWvfz2VpSXEOfD9P57FZFCet6unMjUbpdRipGAkA4LAYDYI8IeyN1skKOUuGAGAGqcV711Xh/1nB1MKXEkYT2Go3YwtpjSmUsgKx/Bb/fCPXoIkiJAhQGYiFEFAUDDAKs9Nw/NgEJFzI1ajgO1LKvPSyr2uzIqtjeU40D6S8+dKVyZhktEgoDxBl1qtkRWOiUBqwcimxeWo0NH3qEc5rTgaHx/H7bffjoceeggulyvpx91zzz1wu93Tfzo7OzNahyAwfGpHMxgDQjJP+EdWOMwGAbdtXZTRc5Ko/h9ZovBwdiSXbGYD/mJ1LVbVUSCainxO1wXC23ZmiwlmRYI95Icz6ENlwIPqybGYgQgAMJMJAINJZLhmqQu2PE5Zba52YO1CDf9MZZDBrLSbdXXMdSLFrEipxYitTRU5Wg2JSOm30eVyQRRF9Pf3z7i9v78ftbW1c+5//vx5XLx4ETfffPP0bYoSvpgYDAacOXMGS5YsmfM4s9kMszm7ab9Pbm/Cz15phy+Y+JO6wIC/ubpRV90EtcxsECArPGu1GLnOjgDh4HVLUwWqHWa8dn4QIblw0+uZSiUbkm2WZcsRPHEqubpLgcG8cnl48N1SFxzW/E/bXbewDJNBGecGJvL+3PPJ5D9hVYEXr163vCrrxfNkrpReYZPJhM2bN2Pfvn3TtymKgn379mHbtm1z7r9ixQqcOHECR48enf7zvve9DzfccAOOHj2KhoaGzL+DJDltRnxqRzPmq9EyigI+fW1zfhZVBBhjSbV0T5bCw7Nq8qHRZcd719WjXEfzNvIp39kQINK4jEFkDPaVyyHarMmdAlE4HJuvwNamClW3FK5srEB9mU2158+FGh11XgVSC0aaXHYsqcpOiwKSWMrh3u7du/HQQw/hV7/6Fdra2vC5z30OXq8Xd955JwDgjjvumC5wtVgsWLNmzYw/ZWVlcDgcWLNmDUym/L4pfHJ7U8ILI2VFckMUWFaPwwVlJW8DyZxWE96ztg7NLnpDiggXnyZ3tDcbogMQUQj/f8YYIBrgvPm9kTslZNu8CduvWad6sypBYLhmWSUq7IXzHlOjswZgnsnktmmMooDrl1PRar6kHIzceuut+M53voOvfe1r2LBhA44ePYo9e/ZMF7V2dHSgt7c36wvNhvmyI5QVyR2TQYBByE6qk/PwhN98MYgCrm2pwpbGyrjj5otFPrMhDDECkFnMy1eg7EMfAhMN4YAk+i5Tv+i2LVtw3Wf+CvXl2hj1Hr7IVaHEkv+tokTSCfAdFmNWM5/5MJ5kZuTKxgo4VdjOK1Yp9xlRQzb7jLh9IWy7f9+c2hEG4JPXNOFL71oRfk9jbOp/E8+yIMnjnMMbzM5xX8YY7CYx7/9tBsf9+NPZwZSr8fUumeO62cDAwBhSDvq434/J48cxeeIElImp2TRLl8K2eTPWrV6M5XXa64Hh8YXwQls/Alk8cZaJSNCXiuaqEt0def39oS50jcZokBelwm7Cx7YuztoYimKWkz4jhSCSHfnRi2/P+IQnCsBHrmyIatQ1852XsZmDtihQSR1jDBaDgMksvPlyzhGQlLx/KqtyWPDedfV4+ewgetyTeX1uNUQ6qOYyEEk3AJnxNSwW2LZsgW3WWIpl1XZNBiIAUGoz4rplLuw7PaC7OUMReiteBZLLjNywvJoCkTwryhLh2bUjAgM+eU0zJJnjdG/sBmuchz8dSlPNt4KygoCkwB+SEZBkBCUFIVmZ7kKpKDxvdQ16Ysjicd9QHmtHolmMIm5aVYN1C8ry/tz5FNmSycVLHKsOJNsWV1qxdmFZ1r9uNlWVWrB9SaVuP9TorQ08MP/R3hW1DjRUFFaRsR4UZTAyu3YkUiuyqr4UlSVm7GvrT2mqYyRQiXSRnC9YKfZAxWIUs/apI6BSm23GGDYtLseNK6p1NSAsGbksUI1ZiJoDdU4zNi0qz8nXzraGSjs2ZWlIXyZS/e8tCkx3xf4TAQlygm/UZBCwY1lVHldEIgrrXTQF0dmR6BM0DRU2XLOsCi+/PYQT3W7wLKRPZwcrlFUJt4vPxoUoNBXcqWVRhR03r6tHZYGcjshFgWokC2IQhJwGIBFVDjO2NOkr27CirhQr6pKbjq0VlSVmCDrbyphvi+bqJS7Y89gMj1xWtMGI02bE565bghKzYc4JGqfViPetr8fwRABPHu/BqC93Y8BTzarIBRKoRNrFZ0M+T9bE4rAY8e61dVhapa+LSbRsZ0OiA5BcbcPEUmYLd8vU437/5sUVWFRpV+35U/1Pr8d6kUTHeqsdFqzXcpfcAle0wQgA/N2NS3HgnhtjphpNBgF/saYOTa4SPHm0B4cvjaryCTxWsFIoWZVstYsPTX3vahIFhmuWuXD1EpfuLoTZyobkow4kkRKziO1LKnW9bXb1kkpUO/TRt0OP9SLxMiMMDDeuqNZVNq3Q6Pe3NgsYY3DMc9Z/S1MFdq6qQWuvB08c7cbgeCDh/fMpnayK1iaHmg1CVi5a+erKOp+WGgfevaYOpRrrIRFLNrIhcRuS5ZnVKGD7UldeBt/lkigwXLvcBadVhS6xKf4g6DEYidd9dc2CUtQ69REEFqqiDkaStaSqBB/cuAAyB5450YuD7SOQVN4aSEa8YEVLWZVstYuPjJrXgsoSM96zrg4Ly7VbkZ9pNmQ6CFExAIkwTg2+K5S9frNBxA0rqmA15TewSuXHwWY2aK5pWzLGY5yksRpFbF+a/CBXkhsUjCSpssSMWzYvRJ3TglO9Hjx2pBs9Y/rtM5FuViUXwYoosKyk1gOSNppHAeELys6VNdjYUA4tZX6VDLIhatWBJGIQGLYvUWfwXS7ZzQZc36LdAW3VOjtFE+GJEYxcs8yluy6yhUibP+kaZTGKuHldPdYsKMVEQMbzp/rx2rkhBDV0EcyW+YKVbAcqZoOYcbt4WeGay1itbyjDTStrVX+zS3dLRosBSITAGLY2VaCipDAHGVaUmLFjaaUmT6xU6Ww4XsTslg31ZVasrqeiVS2gYCRFgsCwY1kVrl9eBVEAzvZP4LEj3egYSdxeuNDkIqtiMQoZp/vVPlkTS32ZFTevq0NViTp70qluyahdiJoMxoArG8tQU+D7/PXlNmxpzF+/lGQ/UOitXmTfpX040n8Coaj3B4GFi1aJNlAwkqaVdaV43/oFsJlE+IIK9rUNYP/pAQRC2rsY5luyW0CzA5VIu/hMyFPPpzV2sxHvWlOLFbX5a02eajYkmcF0WrFxoRMLNFyTk01Lqh1YszA/PzfJ/KgIOmt2dm70HHb/aTe+8Md7EJIvZ0Y2NJTp6vsodBSMZKDWacGHNi+c/pRwYdiHPxzuxIXBCZVXpm2JWuvLUwPZMtkC0srJmtkEgeGq5kpcu6wKhhyn3pPNhszehtFyABKxqtaBxqoStZeRV+sXlqNZI99zmc2k2VqW2Tjn+OYb3wQADHt9ODRwCABQYjbgquZKNZdGZtHHT5SGlZgNeP+GerTUhN8oAhLH/rND+GNb/5zJwGR+nId7vIBfrgEJyQqCkoJApFZlagsoXqCicG1mRyKaq0rw3rX1OTm+mUw2RMt1IPNZWm3HinptDr7Lta1NFahzWtVehq6KV/dc3IPDA4ehcAVcseDowFGM+cdwXUuVrvvRFCL6r5EFBlHAO1bWYNuSyumTE50jk3jsSDfO9MUevEcSs8Yp+OSchz/1zxOo+AKSppvAldlNeO+6OjRmseNmomyIHupA5rO4wop1Gh98l0uCwLCjxYXyHI4eSOZXRS/Fq96QF/cfvB/hGesAl8P1RW+NPIel1drIMpHLCuNgvkZsaChDpd2Eva39CEjhi+Tr50fQPuTFNUuq4LDSy50sJjCYDWJKx3Uj4+7BAQUcvoAEY9SnH8Yuj6tnwFRdBFTbmjCKAq5fXo2T3W4c7hiFMjEISCk01TOYAbsLCo+dCYn+XvWw/ZJIrdOMjToZfJdLRlHADcur8PypPngDiafP5opeilcfPPYgxvxjmHpXAFcsUCDhovQM/thxNW5afJPKKyTRGNfix8ZZPB4PnE4n3G43Sku1n6J1+0LYc6oXI97LxVIGITzldXWdA5pqPKFx/qAMSUlvy4UxBrtJTOr1nh2oMMam/jc/F/L+i6fx8s/+EV6e2hu9ct0/gdtmNmyKfB96zH7EUlliwvYllbqpU8gHty+IF1oHst5WgDGWsJ7JbBTx11ctzupz5sK50XP40FMfgsIvv3cER6+EYBqG0X4BLqsLT3/wadiMxVEEraZkr9/0250DTpsRH9y4EI2Vl3/QJYXj4IURPHW8F6Pe3A3eKzQWY/rt4jnnSR/1Vbu1fo3Jj/cKB1DLRpK6v8IBmTPwUDiTouc6kEScViOubqZAZDanzYTrWnIwB2mez6ZVOqgXiRStRrZnIhjjEG0XwcExNDmEh048pNIKSSz0G54jJoOAd62tw+bFM1PLQxNBPHW8B8c6x8A10r5c0zJsFx+SecozN2LJR2t9qxDCLvEQVgsXE65D5gwcDAwcAoOu60ASiQy+M1KhYUzVpRZcnecTIdU6qBd54dILODxwGDK/nDXiXIBoawdj4Q8nHBz/dfK/0OnpVGuZZBb6Lc+xLU0VeOfqGhjFyxcKWQEOd4zhyWM9GPZqZ/CeVgkZtIvnnCOQ46O+2Wytzxhwpfg2bhCOwYSZNQEKD/eBEMAhQoHIuOb7gaTLYggPvrPkeT6L3ixy2bFxUVnWvt584bMe6kX6vf1zbmNMgWCamXWUuYxh/3C+lkXmQcFIHkQG7TksMwtYR3whPHWsF29eHIGs4aOoWmAyiGmnpEOykpXsSLoSBSuXAxUOiQuQOYPCGRYJA3iP+AYq2MTU4wEGPhWA8IIuOzKKDDuWFc7gu1xbVe/E8lpHzp+HMYYqh/Y73t7ScguqrFVztmmiiUzEtQuuxYbqDflbGEmIgpE8iQzaqy+b+cvMOXCy24PHj/ag3+NXaXX6YDGIaWcBcp0dSdd0oMIBCSJCMCAIAwIwwswkXCOcgIEpMDAF4lRxbSEzCAxXL6ksuMF3ubZ5cTkWVWTnmHi87UWn1aiL3hw2ow1f3vrl6VM0sTDGcPfWu/O4KjIf7f9kFZDoQXuzefwSnj3RhzfODyOk0Qun2sLHfdP7kQ3Jii5rdIIwTJ2MAYSprIg49UeYqhlhUX/0TGAMW5srUKmDIkmtYYzh6qWVOc1c6GGLJuLGRTfiqrqrILK523wMDJ9Z+xk0OBpUWBmJh4KRPJs9aG+2tr5xPHa0G92jk/lfnA4YRAHGNE9WaHGI3nz8PH6X1sixYyHqjyhEilmhu0Bl8+Iy1JRqfxtAq0SB4doWF0ozzCrFC9n10uwMCP9efOWqr8y5XYCAGlsN7lxzpwqrIolQMKKS6EF7s3kDMl5o7ccrZwc1O2dFTWZDesd99Zgd8SO9lvFsqrBVL4HKhoVONFRQz4dMWYwiblheBWsOCn/11AYeABaXLsYn1nxiRu2IAgVfueorsBgo6NUaCkZUFBm0VxUn/Xlu0Is/HOrExSFvnlemcRkc99Vq7UgsDByTLPsXAK0FKitrHWim9txZU2Ix4rqWqrR7s8QqGTEaBJTbsz9LKdc+tfZTcFldUxOpRexYsAPXNVyn9rJIDBSMqKzEbMAHogbtzeaXOF46M4gX2/oxSYP3pglT7eJTJSlKVpuV5ZIAjskUO7JmKt+BypJqO1YW6eC7XKosMWPH0koIWWqK5iox6/IIeXQxK2MM92y9R+0lkTjo7JwGRAbtVZaY8Ub7cMxPJpdGJtHj7sLWpgosq8n9MT49MBqEqeOyqWU7gpKivf4V5rnBqAAFk5gnGDHkL1iJXIxmXJKm/hI5gRH9ozvfaeqGcivWF/Hgu1yrL7fhysZy/Lk9ua6+icTL3urBjYtuxEdaPoLlFcupaFXDaDaNxnSO+PBCa3/CWpH6Mgu2L3WhhPowAJzDG5RT73ZqFCFqrcX4SDsQmJj+q0kEnnl7Eu5AnO/NYAbsVXlaXPpiBSo1DjOuas7eJ3cS39GuUZzqSn56eKz5NDetqsGiLE6YJsUj2es3Xc00pqHChls2LcRzJ3sx6gvFvE/PmB+PH+nG5kXlWFnsg/cYg8UgYDKU2hZWUFZg1VowUtE846+CUcRkbydg0U+dSyyzMyoVdhO2L3XBIAqXAxV+OViZnr5MsmLDwnJMBhS0D07Mf2dgKqU18z2lmk45kRzT2LsxAcKD9v5y08xBe7OFZI43Lozg2ZN9cMcJWoqFmMZxX1nhkDR81FdgDJKsFNxpqlKrAdcsvTz4brpGRbhcmxI5vm0Qwp/Qxaj6lSIOuzOytakCdU5rWo91WIwZzYciJBkUjGhUZNDepnnmTvR7AnjiWDeOd7l1d2w1m8zG1NvFa/lCLzDAl2K2R+vsJhE7lrpgSrLwmAKV7BEEhmuWuVBun7/2Y/a7iJ6anRH9omBE47Y2V+KmVTMH7c0mK8ChS6N48ngPRrzBPK5OW1JtF69w7WZHBMbgL6DTU2aDgGuWuWA1ZWdnmAKV1JkMAq5fnvrMHz01OyP6RcGIDiytLsEHYgzam23EG8JTx3pw6NKobo6vZlM67eK1mh0RBAZfgQQjxqmuww5LfubNUKASn81kwA3Lq+bNTkUXhFNmhOQDBSM64YozaG82hQPHu9x44mg3BscDaT0X59q8QCcj1XbxCueamwUU6S5bCMGIQWDYvtSFMps2Bt+lEqgwFOZgQqfNhGuXuRJua0ZCEVEUUJnE1g4hmaJgREcSDdqbbWxSwjMnenGwfSSlrYj/PPyf+Nizfw1J1m9RbKrt4oOyMn9TjDyKXCMmQ5K6C8lQZPCdSyefrGcHKgZRKNhApcZpwbbmypj/xgE8d+E57O/aj0q7iY5fk7ygYERnIoP2rmuJPWgvGufAqV4PHjvSjZ6x+QfvtY204tetv8bZ0bP4zenfZGnFKkixXTznXFND9FiBZEY2Ly5L+wSH1hRioLLYZcfGGAXy7aPncGLoJP7cexBuqSP/CyNFiYIRnVpVH3/Q3mwTARnPn+rHa28PIRTnosu5gm/9+VsQWPhH4v899lP0+/qzuuZ8SrVdfEjmmsmOTGdGdByMrF/oxOIiaZKl50BlVb0Ty2svd3QOKRL2db4EhnB28TdnfwJZ0e/PIdEPCkZ0bL5Be7OdHZjAo4e70Dnim/NvT55/Cm3DpyHz8BuPzGV8763vZXW9+WY0CDAIyf2IayU7Ep75Er5cpdrITStW1pXSyIIpeghUNi8ux6KKcOD4Zv9BTITGwcGhcI6OyZP4n7P/o8KqSLGhYETnIoP2liU59dQXVPDHtgHsPz2AQCh88fUEPPjB4e/PGLUtcxn7Ovbhjd43crLufLEYhaSP+2ohOxJd66LHzEizy47VNPguKVoJVBhj2LakAoJhAn/uPThdvMqEAAQxiO8f+j5G/JnPtyEkEQpGCoBBFLBzVQ2uaq5IujP8hWEf/nC4E+2DE3jg6APwhnyY3YRbgID7D96v62LWSLv4ZHDOEVD5ZE0kGAlIMiSdHc9uKLfGrEEgqct3oCIKDK+P/hxMuFxbxgzh9vEBOYDvH/p+hs9ASGIUjBSQjYvK8Z61dTAlefENSBy/O3oKjxxuhSzN7WGiQEH3eLe+i1mRWrv4kKyo2slWr/UiNaVmXNlYocsx83qTi0Dl1e5X8Obg6xDsbQALAWBgYni4nsxlPH7ucRwdOJrD74oUOwpGCkxk0F55En0dOBTs6/gjIFVC8myCHKiNcR+u+2JWILV28WrVjkTXi+jpJE2F3YhtNIFXE9IJVPySH//x1nfAwMDEAAwlZwAmgxnGp7+uwAR848A3qJiV5AwFIwXIaTPig5sWYHGCQXsA0D3eg37vABTOwbkI2bsEIc9qcGnmccyQEsIfzv4hl0vOi2TbxauVHZlRL6KT4tVSiwHbl7imB98R7YoXqOzt2IMBXx84ZAAcgmECBtvbEIze6ccqXMHbY29jf9d+1dZPChu9gxQos0HEu9bUJhy057K6YDGYZ6RxuVSGoGcjJO8ycDl8SoeD48raK3O74DxIpV28GrUjeitetZtE7FjmgpkmuuraEufSqawcwBgHE2SIlkEwNvN3QGACmpxNKq2SFDoKRgoYYyzhoD2LwYJrF1w7Z0onA4Psr0bIfQWUySW4YcE7CyIYAZJvFy8pSt7n+0T/F9L6Nk22B98R9ayvXod3N70HoiAATAJjypxCeIEJ+Piqj6PZ2azOIknBo2CkCEQG7ZXEGLS3yrUKtfaaOe3TGQufcuX+BixmH8WRjjHNTrhNVbLt4vM9RC+65kLL2zRGMb+D70huCYzhH7d8AWaDGPM0HgNDubkc/2v9/8r/4kjRoGCkSLhKzPhwjEF7DALesWjnjCmd0//GgG31V8Ei2vHWxRE8/GYnWns8qp42yYok28VLigI5TwHY7OBoMqjNuTSiwHD1Eu0MviPpYwwwiQJMBgFVtip8YdMXYt6Pg+PuLXfDbiyOjrpEHRSMFJHIoL3ZTamqbdVYV7V+xgVRYAzl5jJsqtkIUQgXvk0GZbx2bgiPHOrE+YGJfC8/qwSBJXUEOl8na2YfRNHiNg1jwNamiqQ7/hJtYgwwigLMBnFGNu7W5bdiiXPJ9EgIABCZiCtqrsCuxl1qLJUUEQpGiowgMFzbMnfQ3tX1V8MkmqbrFhTO8Y7F74DIwhmESEACAOOTEl48PYA/xGktrxcmw/zHfWWF5yU7Eh0Ics7hD2lvS2zTonLUlxXG4LtiFQlCYv3cGwQDvrrtq1D45Z89Do6vXvVV6h9Dco6CkSK1qr4UN6+vnx60F13MKoChpXwZGhyLZjwmOiABgJGJIPac7MPTx3swMO7P5/Kzxmqc/7hvPk7WRH9C9YcUKBoZ2hexboETTS5K0+uVQQhvTc4XfG+u2Yz3Nr8XIhMhYKpotYyKVknuUTBSxOqc1hmD9qaLWQUR1y68LuZjZgckANA75scTR3qwt7UPY75gztedVUm0i1c4z2nx7ux6EZ/G6kWW1zrQUkuD7/QoEoSk0gfmH674B5hEE8otVLRK8ofxWJWLGuPxeOB0OuF2u1FaSkO4sk2SFew/M4i3BybglyYxEfLCZXUlfIys8LhFry3VDmxaXIYSHZ22CIRkhBIEHAJjsJlzc4zVMNWAKqJr1If9ZwZz8lypanLZsXlxudrLICkSp7qtpru9cmzwGKwGK1rKW7K8MlJskr1+U5MAMj1or7LEhD9fGIHFMH9dgCgwyArmBCScA2f6x3FucByr6pzYuKhMF02xzEYRssLjbo8onCMkKTAm2TQtFbMvGFppeLagzJqwaR7RHoExGMX0g5CI9VXrs7QiQpJD2zRk2sZF5Xj3muQH7YkCi9uvQ1aAE91uPPxmJ450jOqiR4l1nqApGI6+sv68WjxJU+0wY0sTDb7TC4GFuwubDAL9NyO6RMEImWFRZfKD9oBw4WWiBmJBScFbF0d10aMk3C4+fkDCOc/6Ud/o4XgRagcj5TYjti2pTHqwIFFPdK8QCkKInlEwQuZIdtBexHwBCQDd9CgxGgQYhPi/FiGZZzU7Eut186vYfdVhMeCapa6kWuYT9cTrFUKIXtE7DokpMmhvY5I1A4LAknpT1EOPEosxfrv4bGdHYj2PWpkRm0nAjqU0+E7rEvUKIUSvKBghcTHGcFWCQXuzCSy5gATQeI+SedrFZzM7EuvlUuNor0lkuGZZVc5ODJHMJdsrhBA9SisYeeCBB9DY2AiLxYKtW7fi4MGDce/70EMPYceOHSgvL0d5eTl27tyZ8P5EexIN2pstlYAE0G6PkkTt4jnnWWmEFqteRFGy87VTYZwKREp1dBS7mKTTK4QQvUn5p/uRRx7B7t27ce+99+Lw4cNYv349du3ahYGBgZj3379/P2677Ta89NJLOHDgABoaGvDOd74T3d3dGS+e5I+rxIxbNs0dtBeLwFjKn94uDvnw+0NdePnMICb8oXSXmVWJ2sWHZCXjYtyYWzR5rhcRGMO2Zhcq7Ka8Pi+ZnyiET8hQEEKKQcpNz7Zu3Yorr7wSP/rRjwAAiqKgoaEBn//853H33XfP+3hZllFeXo4f/ehHuOOOO5J6Tmp6ph2KwvHquSGc6vHMe1/OOeQ0LtiiAM30KOEKhy8kx2zwZhSFjNZnFIU5wc7AuB8vnOpP+2umgjHgqqZKLCineTNakq1eIYRoQbLX75RC7mAwiEOHDmHnzp2Xv4AgYOfOnThw4EBSX8Pn8yEUCqGioiLufQKBADwez4w/RBsig/aubXHFrHeIxtLIkADa6lHCpj6dxpJpdiTWK+MP5u973byojAIRDREYo2O6pGilFIwMDQ1BlmXU1NTMuL2mpgZ9fX1JfY0vfelLqK+vnxHQzHbffffB6XRO/2loaEhlmSQPVtc78b4N9bAaE/8ITQckaby3aqVHiUEU4h51zaS+I1ZtjS+Un+LVNfWlaHSV5OW5SGLRvULomC4pVnndjLz//vvx8MMP47HHHoPFEr/24J577oHb7Z7+09nZmcdVkmTNHrQXD2MMIksvIAG00aPEbIh93FdSFChpBEnxjg7n41hvS40DK+pou1Nt1CuEkMtSCkZcLhdEUUR//8w97f7+ftTW1iZ87He+8x3cf//9eOGFF7Bu3bqE9zWbzSgtLZ3xh2iTw2LEBzbUY1l14k/ZmQYkgMo9ShIc9w2mkR2Jd+3x5zgYaaq0Y91CZ06fg8yPeoUQMlNKwYjJZMLmzZuxb9++6dsURcG+ffuwbdu2uI/79re/jW9+85vYs2cPrrjiivRXSzQpMmjvquYKJNrqzkZAAqjXo0SI0y5eUhTIKda1qJEZqXdasGlxWc6+Ppkf9QohJLaUt2l2796Nhx56CL/61a/Q1taGz33uc/B6vbjzzjsBAHfccQfuueee6fv/+7//O7761a/iF7/4BRobG9HX14e+vj5MTGi3JThJTzKD9rIVkAAze5SM5qlHSbx28al2ZY2Xls/V0d4qhxlbmyupMFIlIvUKISShlNst3nrrrRgcHMTXvvY19PX1YcOGDdizZ890UWtHRweEqDfrn/zkJwgGg7jllltmfJ17770X//Iv/5LZ6onmLKq04UObFmDPyT6M+mL3C2GMQQQggwNZqEm9OOTDpWEfWqod2LS4DCU5bt5lMQrwBvmM476ywiHLCsQkLjaJ5vhM5qD7apnNiG3NNPhODaLAYBDomC4h80m5z4gaqM+I/gQkGfvaBnBpOH5tB+ccCkfMHh7pylePEllWMDkriyEwllQ7dYPAYn5CDskKHnkzu8XaJWYRNyyvVr1fS7GhXiGEhOWkzwghyUpm0F7k2G8237Dz1aNEjHHcV+E8qeeL9/1mu17EahSwY1kVBSJ5RL1CCEkPBSMkZy4P2qtOOGgv2wEJkJ8eJeYYhYjJnKyJt1symcVgxCgy7FhWBTsNvssL6hVCSGYoGCE5t7TagfdvSDxoLxcBCZD7HiUWgzhj3QrnCCUISGINx4vI1rRegxAOREqtNPgu16hXCCHZQcEIyYsqx/yD9nIVkAAze5R0jHiz9nVjtYsPygoQpw4mYfFqFk7SCFPZKBp8l3vUK4SQ7KFghOSN1STi5nX1WF0fv4gplwEJEO5R8vzJfjx9vAf9nuz0KJndLp5zHveob+KTNJkFI4wBW5rKUeucf7IySR/1CiEk+ygYIXmVzKC9XAckQLhHyZNHs9ejZHa7+JDMY2ZHEl2/Mi1g3dRQhoXltoy+BomPeoUQkjtU3UZUsbreiQq7Cc+f7MNkaG4WQRQYFCVcg5FLWetRMtUuPlL3EcmOmKI6tiaqFwEy26ZZXVeKpioafJcL1CuEkNyjEJ+oJjJoz1USu75BEFjCbY1s4Rw40z+O373ViTfODyOQZlAwu1387OzIfN9Luts0LTUlWJlg64ukR2DheiCjSMd0Cck1CkaIqhwWIz64cUHcQXv5CkiA7PQoiW4XzzlHIOpkzXzfRTrbNI0VNqxbWJby40h81CuEkPyjYISoLjJob2ucQXuCwPJ6bDLTHiUW4+WLWEhWph+fKKjyh+SUt6TqnBZsWlye0mNIfNQrhBD1UDBCNGNTgkF7AstvQAJk0KOEMViivofIyZpE6091i8ZVYsLWpgq6aGYJ9QohRF0UjBBNiQzaK7PNLSRVIyAB0utREt0uPiQr8w4ETGVar9NqxPYlLjrVkQVGUaBjuoRoAL2bEc0ps5nwl5sWYHHl3GOqwtQ8GzWk2qMkul28pCSuP0k2M1JiFnHtMheMMbJHJHnUK4QQbaF3NKJJiQbtMRUDEuByj5IXTs3foyTSLl5WOOQEtSfJBCMWAw2+y5Q41TGXskqEaAv1GSGaFRm05yoxYf+ZwfBR2ah/EwUkvMDn2qVhHzpGfFhWXYLNi8tj9iiJtIvnCBfGWk2xA4n55tIYxXCzOBp8lx7qFUKIttE7G9G8pdUOOK0m7DnVhwn/5Yv2dEDC+bw1GbnCOXC2fwLnByewss6JTYvK5mQuTAYxPEBPViDJLOan8kQ1IwaBYftSFw2+S4PAGIwiBSGEaB3lKokuRAbt1c2au8IYg8jY/E08ckxWgJNTPUoOXxqdMblXYJfbxcebWeOPs00jMIarllTCVWLOyboLVfQxXQpECNE+CkaIblhNIt63fu6gPa0EJEB4K+bQpVE8/GbHdI8SxsKfzC1GEbLCYzZTi9fw7MrGctSW0uC7ZEWCEDqmS4i+0DYN0ZXIoL0KuwmvnRtCpGSEMQYRgAz1tmyi+UMKXjs3hGNdY9i+pBIttaUQBQaTQUBAUmAQBSjBIILnz0P2TcI7wCCUV8z4GhsbytBQQYPvkmUUBTodQ4hOUTBCdGnNAicqS2YO2tNaQAIAEwEJ+04P4miXG1ubKrC40o7A4AB6fvsbjP/uEShuN/yiCQMNm2Gsr4Nt61ZY163DqvpSLInTIp/MZBBi1+EQQvSDcZ7jsahZ4PF44HQ64Xa7UVpKA8HIZeP+EPac7MPQxOUjtpxzVYtao80+hlwVHEfdv/8zbAM9sIQCAIAxUwleq18b3mPgHMvXNGLnN78EZqDPColQEEKI9iV7/abfZKJrkUF7S6OyCJEaEi0ULkYvgQeDaPvlf+OP9iYcrFqOYbMDAOA3TE0t5hy1E0NoeO73GPg//0eF1eoD9QohpPDQbzPRPYMo4KZZg/YiGQm1A5LoZ/e3tUIZHwc4R5+9EnsXXYmDNSsxZrIDAConx7Bq5CLAOcZ++9+QRkdVWbNWRYIQo0gnZAgpNBSMkIKxaVE53rWmdsagPVUDEoYZz+07fGTGPwvgaC+tw0sNmzBqsmP5yCUIU3tLXJLgfurpvC5XqwRGQQghhY6CEVJQFlfa5wzaUysgYdF5EUWGPDIS4z483DmNMfxx8ZVoq1iMEBMBQUDgzOk8rlZ7qFcIIcWDghFScGIN2lMjIIl+tnh14gxARWAcIpcREgxorWjE841bcK60DrKUuEV8oaJeIYQUHwpGSEGKNWgv3wFJ9FMx0QDBao19PwAi5xC5AgaOgGjCMddSPFe+Em8PjMcNZAqRkYIQQooSBSOkYEUG7d20qhqGqYtbPgOS2c9jWbduZoQy+/64HJQIigJ51Vq8dHoQjx7uRsewN8erVZdBCHeopaZlhBQnCkZIwVta7cAHNi5AiSXct0MUGIQcBySxAh7r+nVgwvy/cowxWJoaYa6ogMAYRrxB7DnVjyeP9qDP7c/FclUTCULomC4hxY3eAUhRmD1oT8hxQBLrK4uOUpS+973h7Ei852YMYlkZSt/1F9PHkw2iAJEx9Hv8ePJYD54/2YcRbzD243WCeoUQQqJRB1ZSVBSF45VzQ2jt8Uz/XcnBr0Ci7aBQVycm/vQyQn194RumOq8yQYB51So4rrsOzBJ7OJ7CORSFAwxYVl2CzYvL4bAYY95Xi0SBwaCB/i+EkPxI9vpNwQgpSie73dOD9qYv8FmUzCd+aWAAoc5OKJIE0W6DeelSMEvsItfZFB4OogQGrKpzYuOiMliMYqbLzhmBMRhFCkIIKTbJXr9p+AUpSmsWOFFhN+GFU1OD9gRkLSBJ9oJrqK6Gobo6recQWHibiXOOkz1unO7zYN1CJ9YuKJvR9E1tjAFGQaDTMYSQhLTzrkVIntWXWfGhzQtRWWIKX9yzdMHM52WXMQaDIIBz4EiHG4+82YGT3W7IWc70pL4u6hVCCEkeBSOkqEUP2staQKLCtTdS7BqSOd5oH8Hv3upUrUcJ9QohhKSKghFS9IxRg/ZEgWXc60LNS3AkKJkMyvjTmaG89iihXiGEkHRRzQghUzYtKkel3YQ/tg0gEJLT2+pgydeM5BJjDCID3JMhvNA6gJpSM65srECtM8YpneF2IDiR3Bc2lQCVzTNuMkwdPyaEkHTRaRpCZhn1BrHnVB9GvUHInAMp/IZEMhNapCgcDRVWbF5cjsoSc/jG4XbgJ9tS+0KfOwBUNtMxXULIvJK9ftPHGUJmKbdHBu3ZITKW0r6Llq/LgsDQPebHE8d68KczAxj3h5LPiER/ndAEzAYBRpGm6RJCsoOCEUJiMBtEvHttLTYuKk8pINHDpZmBoX3Ihz8c7sIbXX74eXJN0wQoMCMEE/ULIYRkGQUjhMTBGMO2JZW4aVUNzKIwf6ShkXqRZHHO0DYk4X/k63BEWYogj900jYHDBAkmJms680MI0S8qYCVkHstqHHDajNhzsg/uydB0DYk0OADfseMIdXaCyzKMJXbY1qyGZeUqMJNJ3UWnQGIGHOVL0aY0YIPQjuWsAyLjYOAwQIbINF9WRgjROQpGCElCtcOCWzYvxPMn+9A5NAHPCy8gcPbs9FwZAAhNTGC8rw8Tf3oZzvfdDFNjk8qrTo4REmwIwC74MQQnRL4Aq3AJRkFRe2mEkCJBwQghSbKZDLh5XR0e++I3MHRpOHxj1GE0NpUy4aEQxh59DGW33ALTokVqLBWMAVajCKtJhN1kgM10+f9bTSLsZhE2kwG2IQ+Mh59TZY2EEBJBwQghKZjc/xLW7P0dzM56HK5eDpldLruaUU7BOcb3voDKT35y9r9kxCAw2MwibKapYGLG/878/0nVr4hUBEIIUR8FI4SkYOTXvwFEEcvcPSgLevFK/Vr4RfN0ViSaPOZGsKMzqeyIxSjMyFpYjYap7MXlAMNqEmE2aHcyLyGEpIuCEUKSpAQC8L3xxvTfqybd2HXpTby8YAPcJvuc+wsMMLSfRc2aZXOCisjWic1sgM1Ic1wIIcWNghFCkqR45854sUsB3NTxJi6U1sGkSLBIQVjkAKxSECbG4VxdhvpNC1VYLSGE6AcFI4QkSSgpmXF6JsLAFSxzd899gMEAUevjC8yO/DyGEEISoGCEkCQJJhPs12yH9/UDgCzP/wBJguOmnblfWCYqlwCfPwwExpO7v9kRfgwhhGQRBSOEpKD8Yx+D95VX57+jIMDUuBjWzZtzv6hMUXBBCFEZtYMnJAUl110Hx7velXginiCAiSLq/u3fdNUenhBC1ELBCCEpYIxhwbf/Hc4P/WX4BjHqqK0Q/nUSS0vR8POfwbZxoworJIQQ/WGcc80PnvB4PHA6nXC73SjVekEgKRqBCxcw9vAj8P75z+B+Pwy1tSj74Afg+Iu/gGA2q708QghRXbLXb6oZISRN5qYm1Nxzt9rLIIQQ3aNtGkIIIYSoioIRQgghhKiKghFCCCGEqIqCEUIIIYSoioIRQgghhKgqrWDkgQceQGNjIywWC7Zu3YqDBw8mvP///M//YMWKFbBYLFi7di2effbZtBZLCCGEkMKTcjDyyCOPYPfu3bj33ntx+PBhrF+/Hrt27cLAwEDM+7/++uu47bbb8MlPfhJHjhzBBz7wAXzgAx/AyZMnM148IYQQQvQv5aZnW7duxZVXXokf/ehHAABFUdDQ0IDPf/7zuPvuuT0Xbr31Vni9Xjz99NPTt1111VXYsGEDHnzwwaSek5qeEUIIIfqTk6ZnwWAQhw4dwj333DN9myAI2LlzJw4cOBDzMQcOHMDu3btn3LZr1y48/vjjcZ8nEAggEAhM/93tdgMIf1OEEEII0YfIdXu+vEdKwcjQ0BBkWUZNTc2M22tqanD69OmYj+nr64t5/76+vrjPc9999+HrX//6nNsbGhpSWS4hhBBCNGB8fBxOpzPuv2uyHfw999wzI5uiKApGRkZQWVmZ0RRUj8eDhoYGdHZ20nZPntBrnn/0mucfvebqoNc9/1J9zTnnGB8fR319fcL7pRSMuFwuiKKI/v7+Gbf39/ejtrY25mNqa2tTuj8AmM1mmGcNGisrK0tlqQmVlpbSD26e0Wuef/Sa5x+95uqg1z3/UnnNE2VEIlI6TWMymbB582bs27dv+jZFUbBv3z5s27Yt5mO2bds24/4AsHfv3rj3J4QQQkhxSXmbZvfu3fj4xz+OK664Alu2bMEPfvADeL1e3HnnnQCAO+64AwsWLMB9990HAPjCF76A6667Dt/97nfxnve8Bw8//DDeeust/PSnP83ud0IIIYQQXUo5GLn11lsxODiIr33ta+jr68OGDRuwZ8+e6SLVjo4OCMLlhMvVV1+N3/72t/jKV76Cf/7nf8ayZcvw+OOPY82aNdn7LpJkNptx7733ztkCIrlDr3n+0Wuef/Saq4Ne9/zL1Wuecp8RQgghhJBsotk0hBBCCFEVBSOEEEIIURUFI4QQQghRFQUjhBBCCFFVQQUjDzzwABobG2GxWLB161YcPHgw4f3/53/+BytWrIDFYsHatWvx7LPP5mmlhSWV1/2hhx7Cjh07UF5ejvLycuzcuXPe/05krlR/1iMefvhhMMbwgQ98ILcLLECpvuZjY2O46667UFdXB7PZjJaWFnqPSVGqr/kPfvADLF++HFarFQ0NDfjiF78Iv9+fp9Xq38svv4ybb74Z9fX1YIwlnCEXsX//fmzatAlmsxlLly7FL3/5y/SenBeIhx9+mJtMJv6LX/yCnzp1in/605/mZWVlvL+/P+b9X3vtNS6KIv/2t7/NW1tb+Ve+8hVuNBr5iRMn8rxyfUv1df/oRz/KH3jgAX7kyBHe1tbG/+Zv/oY7nU7e1dWV55XrV6qvecSFCxf4ggUL+I4dO/j73//+/Cy2QKT6mgcCAX7FFVfwd7/73fzVV1/lFy5c4Pv37+dHjx7N88r1K9XX/De/+Q03m838N7/5Db9w4QJ//vnneV1dHf/iF7+Y55Xr17PPPsu//OUv80cffZQD4I899ljC+7e3t3ObzcZ3797NW1tb+Q9/+EMuiiLfs2dPys9dMMHIli1b+F133TX9d1mWeX19Pb/vvvti3v8jH/kIf8973jPjtq1bt/L/9b/+V07XWWhSfd1nkySJOxwO/qtf/SpXSyw46bzmkiTxq6++mv/sZz/jH//4xykYSVGqr/lPfvIT3tzczIPBYL6WWHBSfc3vuusufuONN864bffu3Xz79u05XWehSiYY+ad/+ie+evXqGbfdeuutfNeuXSk/X0Fs0wSDQRw6dAg7d+6cvk0QBOzcuRMHDhyI+ZgDBw7MuD8A7Nq1K+79yVzpvO6z+Xw+hEIhVFRU5GqZBSXd1/wb3/gGqqur8clPfjIfyywo6bzmTz75JLZt24a77roLNTU1WLNmDb71rW9BluV8LVvX0nnNr776ahw6dGh6K6e9vR3PPvss3v3ud+dlzcUom9dRTU7tTdXQ0BBkWZ7uAhtRU1OD06dPx3xMX19fzPv39fXlbJ2FJp3XfbYvfelLqK+vn/MDTWJL5zV/9dVX8fOf/xxHjx7NwwoLTzqveXt7O1588UV87GMfw7PPPotz587hb//2bxEKhXDvvffmY9m6ls5r/tGPfhRDQ0O45pprwDmHJEn47Gc/i3/+53/Ox5KLUrzrqMfjweTkJKxWa9JfqyAyI0Sf7r//fjz88MN47LHHYLFY1F5OQRofH8ftt9+Ohx56CC6XS+3lFA1FUVBdXY2f/vSn2Lx5M2699VZ8+ctfxoMPPqj20grW/v378a1vfQs//vGPcfjwYTz66KN45pln8M1vflPtpZEkFERmxOVyQRRF9Pf3z7i9v78ftbW1MR9TW1ub0v3JXOm87hHf+c53cP/99+OPf/wj1q1bl8tlFpRUX/Pz58/j4sWLuPnmm6dvUxQFAGAwGHDmzBksWbIkt4vWuXR+zuvq6mA0GiGK4vRtK1euRF9fH4LBIEwmU07XrHfpvOZf/epXcfvtt+NTn/oUAGDt2rXwer34zGc+gy9/+cszZqaR7Ih3HS0tLU0pKwIUSGbEZDJh8+bN2Ldv3/RtiqJg37592LZtW8zHbNu2bcb9AWDv3r1x70/mSud1B4Bvf/vb+OY3v4k9e/bgiiuuyMdSC0aqr/mKFStw4sQJHD16dPrP+973Ptxwww04evQoGhoa8rl8XUrn53z79u04d+7cdOAHAGfPnkVdXR0FIklI5zX3+XxzAo5IMMhpBFtOZPU6mnLJq0Y9/PDD3Gw281/+8pe8tbWVf+Yzn+FlZWW8r6+Pc8757bffzu++++7p+7/22mvcYDDw73znO7ytrY3fe++9dLQ3Dam+7vfffz83mUz897//Pe/t7Z3+Mz4+rta3oDupvuaz0Wma1KX6mnd0dHCHw8H/7u/+jp85c4Y//fTTvLq6mv/rv/6rWt+C7qT6mt97773c4XDw//7v/+bt7e38hRde4EuWLOEf+chH1PoWdGd8fJwfOXKEHzlyhAPg3/ve9/iRI0f4pUuXOOec33333fz222+fvn/kaO8//uM/8ra2Nv7AAw/Q0V7OOf/hD3/IFy1axE0mE9+yZQt/4403pv/tuuuu4x//+Mdn3P93v/sdb2lp4SaTia9evZo/88wzeV5xYUjldV+8eDEHMOfPvffem/+F61iqP+vRKBhJT6qv+euvv863bt3KzWYzb25u5v/2b//GJUnK86r1LZXXPBQK8X/5l3/hS5Ys4RaLhTc0NPC//du/5aOjo/lfuE699NJLMd+fI6/zxz/+cX7dddfNecyGDRu4yWTizc3N/L/+67/Sem7GOeWvCCGEEKKegqgZIYQQQoh+UTBCCCGEEFVRMEIIIYQQVVEwQgghhBBVUTBCCCGEEFVRMEIIIYQQVVEwQgghhBBVUTBCCCGEEFVRMEIIIYQQVVEwQgghhBBVUTBCCCGEEFVRMEIIIYQQVf3/5CjnSq5XDegAAAAASUVORK5CYII=", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "# Setttings for plot\n", + "markers = \"vsdo\"\n", + "\n", + "fig, ax = plt.subplots()\n", + "for i in range(len(n_s)):\n", + " ax.scatter(\n", + " out.x_s[i][:, 0],\n", + " out.x_s[i][:, 1],\n", + " s=400 * out.a_s[i],\n", + " marker=markers[i],\n", + " )\n", + "\n", + "for j in range(top_k):\n", + " points = [out.x_s[i][indices[i][j], :] for i in range(len(n_s))]\n", + " # reorder to ensure polygons have maximal area\n", + " points = [points[i] for i in ccworder(jnp.array(points))]\n", + " points = ptc.Polygon(points, fill=True, alpha=0.5 * float(val[j] * n_s[0]))\n", + " ax.add_patch(points)" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "ott", + "language": "python", + "name": "ott" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.10.6" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} diff --git a/src/ott/__init__.py b/src/ott/__init__.py index c40402511..b542553bf 100644 --- a/src/ott/__init__.py +++ b/src/ott/__init__.py @@ -15,6 +15,7 @@ from . import ( datasets, + experimental, geometry, initializers, math, diff --git a/src/ott/experimental/__init__.py b/src/ott/experimental/__init__.py new file mode 100644 index 000000000..241d46d35 --- /dev/null +++ b/src/ott/experimental/__init__.py @@ -0,0 +1,14 @@ +# Copyright OTT-JAX +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +from . import mmsinkhorn diff --git a/src/ott/experimental/mmsinkhorn.py b/src/ott/experimental/mmsinkhorn.py new file mode 100644 index 000000000..a0e5e49db --- /dev/null +++ b/src/ott/experimental/mmsinkhorn.py @@ -0,0 +1,424 @@ +# Copyright OTT-JAX +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +from typing import Any, NamedTuple, Optional, Tuple, Union + +import jax +import jax.numpy as jnp +import jax.tree_util as jtu +import numpy as np + +from ott.geometry import costs, pointcloud +from ott.math import fixed_point_loop +from ott.math import utils as mu + +__all__ = ["MMSinkhornOutput", "MMSinkhorn"] + + +class MMSinkhornState(NamedTuple): + potentials: Tuple[jnp.ndarray, ...] + errors: jnp.ndarray + + def solution_error( + self, + cost_t: jnp.ndarray, + a_s: Tuple[jnp.ndarray, ...], + epsilon: float, + norm_error: float = 1.0 + ) -> float: + coupl_tensor = coupling_tensor(self.potentials, cost_t, epsilon) + marginals = tensor_marginals(coupl_tensor) + errors = jnp.array([ + jnp.sum(jnp.abs(a - marginal) ** norm_error) ** (1.0 / norm_error) + for a, marginal in zip(a_s, marginals) + ]) + return jnp.sum(errors) + + def set(self, **kwargs: Any) -> "MMSinkhornState": + """Return a copy of self, with potential overwrites.""" + return self._replace(**kwargs) + + +class MMSinkhornOutput(NamedTuple): + r"""Output of the MMSinkhorn solver used on :math:`k` point clouds. + + This class contains both solutions and problem definition of a regularized + MM-OT problem involving :math:`k` weighted point clouds of varying sizes, + along with methods and properties that can use or describe the solution. + + Args: + potentials: Tuple of :math:`k` optimal dual variables, vectors of sizes + equal to the number of points in each of the :math:`k` point clouds. + errors: Vector of errors, along iterations. This vector is of size + ``max_iterations // inner_iterations`` where those were the parameters + passed on to the :class:`~ott.experimental.mmsinkhorn.MMSinkhorn` solver. + Follows the conventions used in + :attr:`~ott.solvers.linear.sinkhorn.SinkhornOutput.errors` + x_s: Tuple of :math:`k` point clouds, ``x_s[i]`` is a matrix of size + :math:`n_i \times d` where `d` is common to all point clouds. + a_s: Tuple of :math:`k` probability vectors, each of size :math:`n_i`. + cost_fns: Cost function, or a tuple of :math:`k(k-1)/2` such instances. + epsilon: Entropic regularization used to solve the multimarginal Sinkhorn + problem. + ent_reg_cost: The regularized optimal transport cost, the linear + contribution (dot product between optimal tensor and cost) minus entropy + times ``epsilon``. + threshold: Convergence threshold used to control the termination of the + algorithm. + converged: Whether the output corresponds to a solution whose error is + below the convergence threshold. + inner_iterations: Number of iterations that were run between two + computations of errors. + """ + potentials: Tuple[jnp.ndarray, ...] + errors: jnp.ndarray + x_s: Optional[jnp.ndarray] = None + a_s: Optional[Tuple[jnp.ndarray, ...]] = None + cost_fns: Optional[Union[costs.CostFn, Tuple[costs.CostFn, ...]]] = None + epsilon: Optional[float] = None + ent_reg_cost: Optional[jnp.ndarray] = None + threshold: Optional[jnp.ndarray] = None + converged: Optional[bool] = None + inner_iterations: Optional[int] = None + + def set(self, **kwargs: Any) -> "MMSinkhornOutput": + """Return a copy of self, with potential overwrites.""" + return self._replace(**kwargs) + + @property + def n_iters(self) -> int: # noqa: D102 + """Total number of iterations that were needed to terminate.""" + return jnp.sum(self.errors != -1) * self.inner_iterations + + @property + def cost_t(self) -> jnp.ndarray: + """Cost tensor.""" + return cost_tensor(self.x_s, self.cost_fns) + + @property + def tensor(self) -> jnp.ndarray: + """Transport tensor.""" + return jnp.exp( + -remove_tensor_sum(self.cost_t, self.potentials) / self.epsilon + ) + + @property + def marginals(self) -> Tuple[jnp.ndarray, ...]: + """:math:`k` marginal probability weight vectors.""" + return tensor_marginals(self.tensor) + + def marginal(self, k: int) -> jnp.ndarray: + """Return the marginal probability weight vector at slice :math:`k`.""" + return tensor_marginal(self.tensor, k) + + @property + def transport_mass(self) -> float: + """Sum of transport tensor.""" + return jnp.sum(self.tensor) + + +def cost_tensor( + x_s: Tuple[jnp.ndarray, ...], cost_fns: Union[costs.CostFn, + Tuple[costs.CostFn, ...]] +) -> jnp.ndarray: + r"""Create a cost tensor from a tuple of :math:`k` :math:`d`-dim point clouds. + + Args: + x_s: Tuple of :math:`k` point clouds, each described as a + :math:`n_i \times d` matrix of batched vectors. + cost_fns: Either a single :ott:`ott.geometry.costs.CostFn` object, or a + tuple of :math:`k (k-1)/2` of them. Current implementation only works for + symmetric and definite cost functions (i.e. such that + :math:`c(x, y) = c(y, x)` and :math:`c(x, x) = 0`). + """ + + def c_fn_pair(i: int, j: int) -> costs.CostFn: + if isinstance(cost_fns, costs.CostFn): + return cost_fns + return cost_fns[i * k - (i * (i + 1)) // 2 + j - i - 1] + + k = len(x_s) # TODO(cuturi) padded version + ns = [x.shape[0] for x in x_s] + cost_t = jnp.zeros(ns) + + for i in range(k): + for j in range(i + 1, k): + cost_m = pointcloud.PointCloud( + x_s[i], x_s[j], cost_fn=c_fn_pair(i, j) + ).cost_matrix + axis = list(range(i)) + list(range(i + 1, j)) + list(range(j + 1, k)) + cost_t += jnp.expand_dims(cost_m, axis=axis) + return cost_t + + +def remove_tensor_sum( + c: jnp.ndarray, u: Tuple[jnp.ndarray, ...] +) -> jnp.ndarray: + r"""Remove the tensor sum of :math:`k` vectors to tensor of :math:`k` dims. + + Args: + c: :math:`n_1 \times \cdots n_k` tensor. + u: Tuple of :math:`k` vectors, each of size :math:`n_i`. + + Return: + Tensor :math:`c - u[0] \oplus u[1] \oplus ... \oplus u[n]`. + """ + k = len(u) + for i in range(k): + c -= jnp.expand_dims(u[i], axis=list(range(i)) + list(range(i + 1, k))) + return c + + +def tensor_marginals(coupling: jnp.ndarray) -> Tuple[jnp.ndarray, ...]: + return tuple(tensor_marginal(coupling, ix) for ix in range(coupling.ndim)) + + +def tensor_marginal(coupling: jnp.ndarray, slice_index: int) -> jnp.ndarray: + k = coupling.ndim + axis = list(range(slice_index)) + list(range(slice_index + 1, k)) + return coupling.sum(axis=axis) + + +@jtu.register_pytree_node_class +class MMSinkhorn: + r"""Multimarginal Sinkhorn solver, aligns :math:`k \,d`-dim point clouds. + + This solver implements the entropic multimarginal solver presented in + :cite:`benamou:15` and described in :cite:`piran:24`, Algorithm 1. + The current implementation follows largely the template of the + :class:`~ott.solvers.linear.sinkhorn.Sinkhorn` solver, with a much reduced + set of hyperparameters, controlling the number of iterations and convergence + threshold, along with the application of the :cite:`danskin:67` theorem to + instantiate the OT cost. The iterations are done by default in log-space. + + Args: + threshold: tolerance used to stop the Sinkhorn iterations. This is + typically the deviation between a target marginal and the marginal of the + current primal solution. + norm_error: power used to define p-norm of error for marginal/target. + inner_iterations: the Sinkhorn error is not recomputed at each + iteration but every ``inner_iterations`` instead. + min_iterations: the minimum number of Sinkhorn iterations carried + out before the error is computed and monitored. + max_iterations: the maximum number of Sinkhorn iterations. If + ``max_iterations`` is equal to ``min_iterations``, Sinkhorn iterations are + run by default using a :func:`~jax.lax.scan` loop rather than a custom, + unroll-able :func:`~jax.lax.while_loop` that monitors convergence. + In that case the error is not monitored and the ``converged`` + flag will return :obj:`False` as a consequence. + use_danskin: when :obj:`True`, it is assumed the entropy regularized cost + is evaluated using optimal potentials that are frozen, i.e. whose + gradients have been stopped. This is useful when carrying out first order + differentiation, and is only valid mathematically when the algorithm has + converged with a low tolerance. + """ + + def __init__( + self, + threshold: float = 1e-3, + norm_error: float = 1.0, + inner_iterations: int = 10, + min_iterations: int = 0, + max_iterations: int = 2000, + use_danskin: bool = True, + ): + self.threshold = threshold + self.inner_iterations = inner_iterations + self.min_iterations = min_iterations + self.max_iterations = max_iterations + self.norm_error = norm_error + self.use_danskin = use_danskin + + def __call__( + self, + x_s: Tuple[jnp.ndarray, ...], + a_s: Optional[Tuple[jnp.ndarray, ...]] = None, + cost_fns: Optional[Union[costs.CostFn, Tuple[costs.CostFn, ...]]] = None, + epsilon: Optional[float] = None + ) -> MMSinkhornOutput: + r"""Solve multimarginal OT for :math:`k` :math:`d`-dim point clouds. + + Takes :math:`k` weighted :math:`d`-dim point clouds and computes their + multimarginal optimal transport tensor. The :math:`d` dimensional point + clouds are stored in ``x_s``, along with :math:`k` probability vectors, + stored in ``a_s``, as well as a :class:`~ott.geometry.costs.CostFn` + instance (or :math:`k(k-1)/2` of them, one for each pair of point clouds + ``x_s[i]`` and ``x_s[j]``, ``i MMSinkhornState: + """Return the initial state of the loop.""" + errors = -jnp.ones((self.outer_iterations, 1)) + potentials = tuple(jnp.zeros(n) for n in n_s) + return MMSinkhornState(potentials=potentials, errors=errors) + + def _converged(self, state: MMSinkhornState, iteration: int) -> bool: + err = state.errors[iteration // self.inner_iterations - 1, 0] + return jnp.logical_and(iteration > 0, err < self.threshold) + + def _diverged(self, state: MMSinkhornState, iteration: int) -> bool: + err = state.errors[iteration // self.inner_iterations - 1, 0] + return jnp.logical_not(jnp.isfinite(err)) + + def _continue(self, state: MMSinkhornState, iteration: int) -> bool: + """Continue while not(converged) and not(diverged).""" + return jnp.logical_and( + jnp.logical_not(self._diverged(state, iteration)), + jnp.logical_not(self._converged(state, iteration)) + ) + + @property + def outer_iterations(self) -> int: + """Upper bound on number of times inner_iterations are carried out. + + This integer can be used to set constant array sizes to track the algorithm + progress, notably errors. + """ + return np.ceil(self.max_iterations / self.inner_iterations).astype(int) + + def tree_flatten(self): # noqa: D102 + aux = vars(self).copy() + aux.pop("threshold") + return [self.threshold], aux + + @classmethod + def tree_unflatten(cls, aux_data, children): # noqa: D102 + return cls(**aux_data, threshold=children[0]) + + +def run( + const: Tuple[jnp.ndarray, Tuple[jnp.ndarray, ...], float], + solver: MMSinkhorn, state: MMSinkhornState +) -> MMSinkhornOutput: + + def cond_fn( + iteration: int, const: Tuple[jnp.ndarray, Tuple[jnp.ndarray, ...], float], + state: MMSinkhornState + ) -> bool: + del const + return solver._continue(state, iteration) + + def body_fn( + iteration: int, const: Tuple[jnp.ndarray, Tuple[jnp.ndarray, ...], float], + state: MMSinkhornState, compute_error: bool + ) -> MMSinkhornState: + cost_t, a_s, epsilon = const + k = len(a_s) + + def one_slice(potentials: Tuple[jnp.ndarray, ...], l: int, a: jnp.ndarray): + pot = potentials[l] + axis = list(range(l)) + list(range(l + 1, k)) + app_lse = mu.softmin( + remove_tensor_sum(cost_t, potentials), epsilon, axis=axis + ) + pot += epsilon * jnp.log(a) + jnp.where(jnp.isfinite(app_lse), app_lse, 0) + return potentials[:l] + (pot,) + potentials[l + 1:] + + potentials = state.potentials + for l in range(k): + potentials = one_slice(potentials, l, a_s[l]) + + state = state.set(potentials=potentials) + err = jax.lax.cond( + jnp.logical_or( + iteration == solver.max_iterations - 1, + jnp.logical_and(compute_error, iteration >= solver.min_iterations) + ), + lambda state, c, a, e: state.solution_error(c, a, e, solver.norm_error), + lambda *_: jnp.inf, state, cost_t, a_s, epsilon + ) + errors = state.errors.at[iteration // solver.inner_iterations, :].set(err) + return state.set(errors=errors) + + fix_point = fixed_point_loop.fixpoint_iter_backprop + state = fix_point( + cond_fn, body_fn, solver.min_iterations, solver.max_iterations, + solver.inner_iterations, const, state + ) + converged = jnp.logical_and( + jnp.logical_not(jnp.any(jnp.isnan(state.errors))), state.errors[-1, 0] + < solver.threshold + ) + + out = MMSinkhornOutput( + potentials=state.potentials, + errors=state.errors, + threshold=solver.threshold, + converged=converged, + inner_iterations=solver.inner_iterations + ) + + # Compute cost + if solver.use_danskin: + potentials = [jax.lax.stop_gradient(pot) for pot in out.potentials] + else: + potentials = out.potentials + + cost_t, a_s, epsilon = const + ent_reg_cost = 0.0 + for potential, a in zip(potentials, a_s): + pot = jnp.where(jnp.isfinite(potential), potential, 0) + ent_reg_cost += jnp.sum(pot * a) + + ent_reg_cost += epsilon * ( + 1 - jnp.sum(coupling_tensor(potentials, cost_t, epsilon)) + ) + return out.set(ent_reg_cost=ent_reg_cost) + + +def coupling_tensor( + potentials: Tuple[jnp.ndarray], cost_t: jnp.ndarray, epsilon: float +) -> jnp.ndarray: + return jnp.exp(-remove_tensor_sum(cost_t, potentials) / epsilon) diff --git a/src/ott/math/utils.py b/src/ott/math/utils.py index 6e1aeb7fa..888be6f78 100644 --- a/src/ott/math/utils.py +++ b/src/ott/math/utils.py @@ -182,7 +182,9 @@ def logsumexp_jvp(axis, keepdims, return_sign, primals, tangents): @functools.partial(jax.custom_vjp, nondiff_argnums=(2,)) def softmin( - x: jnp.ndarray, gamma: float, axis: Optional[int] = None + x: jnp.ndarray, + gamma: float, + axis: Optional[Union[int, Sequence[int]]] = None ) -> jnp.ndarray: r"""Soft-min operator. diff --git a/src/ott/solvers/linear/sinkhorn.py b/src/ott/solvers/linear/sinkhorn.py index 08a5ac54a..6083170c5 100644 --- a/src/ott/solvers/linear/sinkhorn.py +++ b/src/ott/solvers/linear/sinkhorn.py @@ -93,7 +93,7 @@ def recenter( """Re-center dual potentials. If the ``ot_prob`` is balanced, the ``f`` potential is zero-centered. - Otherwise, use prop. 2 of :cite:`sejourne:22` re-center the potentials iff + Otherwise, use Prop. 2 of :cite:`sejourne:22` re-center the potentials iff ``tau_a < 1`` and ``tau_b < 1``. Args: @@ -296,8 +296,8 @@ class SinkhornOutput(NamedTuple): (without having to materialize it when not needed). Args: - f: dual variables vector of size ``ot.prob.shape[0]`` returned by Sinkhorn - g: dual variables vector of size ``ot.prob.shape[1]`` returned by Sinkhorn + potentials: list of optimal dual variables, two vector of size + ``ot.prob.shape[0]`` and ``ot.prob.shape[1]`` returned by Sinkhorn errors: vector or errors, along iterations. This vector is of size ``max_iterations // inner_iterations`` where those were the parameters passed on to the :class:`~ott.solvers.linear.sinkhorn.Sinkhorn` solver. diff --git a/tests/experimental/mmsinkhorn_test.py b/tests/experimental/mmsinkhorn_test.py new file mode 100644 index 000000000..932318013 --- /dev/null +++ b/tests/experimental/mmsinkhorn_test.py @@ -0,0 +1,133 @@ +# Copyright OTT-JAX +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +import pytest + +import jax +import jax.numpy as jnp +import numpy as np + +from ott.experimental import mmsinkhorn +from ott.geometry import costs, pointcloud +from ott.solvers import linear + + +class TestMMSinkhorn: + + @pytest.mark.fast.with_args( + a_none=[True, False], b_none=[True, False], only_fast=0 + ) + def test_match_2sinkhorn(self, a_none: bool, b_none: bool, rng: jax.Array): + """Test consistency of MMSinkhorn for 2 margins vs regular Sinkhorn.""" + n, m, d = 5, 10, 7 + rngs = jax.random.split(rng, 5) + x = jax.random.normal(rngs[0], (n, d)) + y = jax.random.normal(rngs[1], (m, d)) + 1 + if a_none: + a = None + else: + a = jax.random.uniform(rngs[2], (n,)) + a = a.at[0].set(0.0) + a /= jnp.sum(a) + + if b_none: + b = None + else: + b = jax.random.uniform(rngs[3], (m,)) + b.at[2].set(0.0) + b /= jnp.sum(b) + cost_fn = costs.PNormP(1.8) + geom = pointcloud.PointCloud(x, y, cost_fn=cost_fn) + out = linear.solve(geom, a=a, b=b, threshold=1e-5) + + ab = None if a is None and b is None else [a, b] + solver = jax.jit(mmsinkhorn.MMSinkhorn(threshold=1e-5)) + out_ms = solver([x, y], ab, cost_fns=cost_fn) + assert out.converged + assert out_ms.converged + + for axis in [0, 1]: + f = out.potentials[axis] + f_ms = out_ms.potentials[axis] + f -= jnp.mean(f) + f_ms -= jnp.mean(f_ms) + np.testing.assert_allclose(f, f_ms, rtol=1e-3, atol=1e-3) + np.testing.assert_allclose( + out.ot_prob.geom.epsilon, out_ms.epsilon, rtol=1e-5, atol=1e-5 + ) + np.testing.assert_allclose(out.matrix, out_ms.tensor, rtol=1e-2, atol=1e-3) + np.testing.assert_allclose( + out.ent_reg_cost, out_ms.ent_reg_cost, rtol=1e-6, atol=1e-6 + ) + + @pytest.mark.fast.with_args( + a_s_none=[True, False], costs_none=[True, False], only_fast=0 + ) + def test_mm_sinkhorn(self, a_s_none: bool, costs_none: bool, rng: jax.Array): + """Test correctness of MMSinkhorn for 4 marginals.""" + n_s, d = [13, 5, 10, 3], 7 + + rngs = jax.random.split(rng, len(n_s)) + x_s = [jax.random.normal(rng, (n, d)) for rng, n in zip(rngs, n_s)] + + if a_s_none: + a_s = None + else: + a_s = [jax.random.uniform(rng, (n,)) for rng, n in zip(rngs, n_s)] + a_s = [a / jnp.sum(a) for a in a_s] + + if costs_none: + cost_fns = None + else: + cost_fns = [costs.PNormP(1.5) for _ in range(3)] + cost_fns += [costs.PNormP(1.1) for _ in range(3)] + + out_ms = jax.jit(mmsinkhorn.MMSinkhorn(norm_error=1.1) + )(x_s, a_s, cost_fns=cost_fns) + assert out_ms.converged + np.testing.assert_array_equal(out_ms.tensor.shape, n_s) + for i in range(len(n_s)): + np.testing.assert_allclose( + out_ms.marginals[i], out_ms.a_s[i], rtol=1e-4, atol=1e-4 + ) + + def test_mm_sinkhorn_diff(self, rng: jax.Array): + """Test differentiability (Danskin) of MMSinkhorn's ent_reg_cost.""" + n_s, d = [13, 5, 7, 3], 2 + + rngs = jax.random.split(rng, 2 * len(n_s) + 1) + x_s = [ + jax.random.normal(rng, (n, d)) for rng, n in zip(rngs[:len(n_s)], n_s) + ] + + deltas = [ + jax.random.normal(rng, (n, d)) for rng, n in zip(rngs[len(n_s):], n_s) + ] + eps = 1e-3 + x_s_p = [x + eps * delta for x, delta in zip(x_s, deltas)] + x_s_m = [x - eps * delta for x, delta in zip(x_s, deltas)] + + solver = mmsinkhorn.MMSinkhorn(threshold=1e-5) + ent_reg = jax.jit(lambda x_s: solver(x_s).ent_reg_cost) + out_p = ent_reg(x_s_p) + out_m = ent_reg(x_s_m) + ent_g = jax.grad(ent_reg) + g_s = ent_g(x_s) + first_order = 0 + for g, delta in zip(g_s, deltas): + first_order += jnp.sum(g * delta) + + np.testing.assert_allclose((out_p - out_m) / (2 * eps), + first_order, + rtol=1e-3, + atol=1e-3) diff --git a/tests/solvers/linear/sinkhorn_test.py b/tests/solvers/linear/sinkhorn_test.py index 5d3fc7751..693b2e62f 100644 --- a/tests/solvers/linear/sinkhorn_test.py +++ b/tests/solvers/linear/sinkhorn_test.py @@ -476,7 +476,7 @@ def test_restart(self, lse_mode: bool): assert num_iter_restarted == 1 @pytest.mark.cpu() - @pytest.mark.limit_memory("50 MB") + @pytest.mark.limit_memory("80 MB") @pytest.mark.fast() def test_sinkhorn_online_memory_jit(self): # test that full matrix is not materialized.