Created
July 8, 2023 05:57
-
-
Save smsharma/9b781aa9310342993bc102ddf9af9852 to your computer and use it in GitHub Desktop.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| { | |
| "cells": [ | |
| { | |
| "cell_type": "code", | |
| "execution_count": 370, | |
| "metadata": {}, | |
| "outputs": [], | |
| "source": [ | |
| "import jax.numpy as jnp\n", | |
| "import jax\n", | |
| "from jax import random\n", | |
| "\n", | |
| "from flax import linen as nn\n", | |
| "from einops import rearrange\n", | |
| "import optax\n", | |
| "import tensorflow as tf\n", | |
| "import numpy as np\n", | |
| "from tqdm import tqdm, trange\n", | |
| "import matplotlib.pyplot as plt\n", | |
| "\n", | |
| "import tensorflow_probability.substrates.jax as tfp" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 371, | |
| "metadata": {}, | |
| "outputs": [], | |
| "source": [ | |
| "from jetnet.datasets import JetNet\n", | |
| "\n", | |
| "particle_data, jet_data = JetNet.getData(jet_type='q', data_dir=\"../data/\", num_particles=30)" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 372, | |
| "metadata": {}, | |
| "outputs": [], | |
| "source": [ | |
| "n_feat = 3\n", | |
| "n_particles = 15\n", | |
| "\n", | |
| "x = particle_data[:, :n_particles, :n_feat]\n", | |
| "\n", | |
| "x_std = x.std(axis=(0,1))\n", | |
| "x_mean = x.mean(axis=(0,1))\n", | |
| "\n", | |
| "# Normalize\n", | |
| "x = (x - x_mean) / x_std\n", | |
| "\n", | |
| "x_flattened = rearrange(x, 'b p c -> b (p c)')" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 373, | |
| "metadata": {}, | |
| "outputs": [], | |
| "source": [ | |
| "class MLP(nn.Module):\n", | |
| " \"\"\" A simple MLP in Flax. This is the score function.\n", | |
| " \"\"\"\n", | |
| " hidden_dim: int = 128\n", | |
| " out_dim: int = 2\n", | |
| " n_layers: int = 4\n", | |
| "\n", | |
| " @nn.compact\n", | |
| " def __call__(self, x):\n", | |
| " for _ in range(self.n_layers):\n", | |
| " x = nn.Dense(features=self.hidden_dim)(x)\n", | |
| " x = nn.gelu(x)\n", | |
| " x = nn.Dense(features=self.out_dim)(x)\n", | |
| " return x\n" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 374, | |
| "metadata": {}, | |
| "outputs": [], | |
| "source": [ | |
| "import functools\n", | |
| "\n", | |
| "import ott\n", | |
| "from ott.geometry import costs, geometry, pointcloud\n", | |
| "from ott.problems.linear import linear_problem\n", | |
| "\n", | |
| "from ott.solvers.linear import acceleration, sinkhorn\n", | |
| "from ott.tools.sinkhorn_divergence import sinkhorn_divergence" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 447, | |
| "metadata": {}, | |
| "outputs": [], | |
| "source": [ | |
| "class VAE(nn.Module):\n", | |
| " \"\"\" A simple variational auto-encoder module.\n", | |
| " \"\"\"\n", | |
| " num_latents: int = 4\n", | |
| " num_out: int = 2\n", | |
| "\n", | |
| " def setup(self):\n", | |
| " self.encoder = MLP(out_dim=self.num_latents * 2)\n", | |
| " self.decoder = MLP(out_dim=self.num_out)\n", | |
| "\n", | |
| " def __call__(self, x, beta, z_rng):\n", | |
| "\n", | |
| " # Concatenate x and beta\n", | |
| " x = jnp.concatenate([x, beta], axis=-1)\n", | |
| "\n", | |
| " # Get variational parameters from encoder\n", | |
| " enc = self.encoder(x) # Shape (batch_size, num_latents * 2)\n", | |
| " enc = rearrange(enc, 'b (n c) -> b n c', c=2) # Reshape to (batch_size, num_latents, 2)\n", | |
| " mu, logvar = enc[:, :, 0], enc[:, :, 1]\n", | |
| "\n", | |
| " # Sample from variational distrib. of latents\n", | |
| " z = tfp.distributions.Normal(loc=mu, scale=jnp.exp(0.5 * logvar)).sample(seed=z_rng)\n", | |
| "\n", | |
| " # Concatenate z and beta\n", | |
| " z = jnp.concatenate([z, beta], axis=-1)\n", | |
| " \n", | |
| " # Decode\n", | |
| " recon_x = self.decoder(z)\n", | |
| "\n", | |
| " return recon_x, mu, logvar\n", | |
| "\n", | |
| "@jax.vmap\n", | |
| "def kl_divergence(mu, logvar):\n", | |
| " \"\"\" KL-divergence between latent variational distribution and unit Normal prior\n", | |
| " \"\"\"\n", | |
| " prior_latent = tfp.distributions.Normal(loc=0., scale=1.) # Prior\n", | |
| " q_latent = tfp.distributions.Normal(loc=mu, scale=jnp.exp(0.5 * logvar)) # Variational latent distrib.\n", | |
| "\n", | |
| " return tfp.distributions.kl_divergence(q_latent, prior_latent)\n", | |
| "\n", | |
| "@jax.vmap\n", | |
| "def reco_gap(pred, true, beta=0.01):\n", | |
| " \"\"\" Gaussian MSE\n", | |
| " \"\"\"\n", | |
| " log_prob = tfp.distributions.Normal(loc=pred, scale=beta).log_prob(true)\n", | |
| " return -log_prob\n", | |
| "\n", | |
| "# @jax.jit\n", | |
| "# def get_sinkhorn(x1, x2):\n", | |
| "# geom = pointcloud.PointCloud(x1, x2)\n", | |
| "# ot_prob = linear_problem.LinearProblem(geom)\n", | |
| "# solver = sinkhorn.Sinkhorn()\n", | |
| "# ot = solver(ot_prob)\n", | |
| "# return ot.reg_ot_cost\n", | |
| "\n", | |
| "# @jax.vmap\n", | |
| "# def reco_gap(pred, true, beta=1.):\n", | |
| "# \"\"\" Sinkhorn reco\n", | |
| "# \"\"\"\n", | |
| "\n", | |
| "# pred_unflatten = rearrange(pred, '(n c) -> n c', c=n_feat)\n", | |
| "# true_unflatten = rearrange(true, '(n c) -> n c', c=n_feat)\n", | |
| "\n", | |
| "# sinkhorn_div = 1 / beta * get_sinkhorn(pred_unflatten, true_unflatten)\n", | |
| "# return sinkhorn_div" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 448, | |
| "metadata": {}, | |
| "outputs": [], | |
| "source": [ | |
| "num_latents = 64\n", | |
| "num_out = x_flattened.shape[-1]\n", | |
| "\n", | |
| "vae = VAE(num_latents=num_latents, num_out=num_out)\n", | |
| "key = jax.random.PRNGKey(42)\n", | |
| "key, z_key = random.split(key)\n", | |
| "_, params = vae.init_with_output(key, x_flattened[:16], jnp.ones((16, 1)), z_key)\n" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 449, | |
| "metadata": {}, | |
| "outputs": [ | |
| { | |
| "name": "stdout", | |
| "output_type": "stream", | |
| "text": [ | |
| "4602346500.0\n" | |
| ] | |
| } | |
| ], | |
| "source": [ | |
| "@jax.jit\n", | |
| "def loss_fn(params, x_batch, log_beta_batch, z_rng):\n", | |
| "\n", | |
| " beta_batch = jnp.power(10., log_beta_batch)\n", | |
| " \n", | |
| " recon_x, mean, logvar = vae.apply(params, x_batch, beta_batch, z_rng)\n", | |
| "\n", | |
| " reco_loss = reco_gap(recon_x, x_batch, beta_batch).mean(-1)\n", | |
| " kld_loss = kl_divergence(mean, logvar).mean(-1)\n", | |
| "\n", | |
| " loss = reco_loss + kld_loss\n", | |
| " return loss.mean()\n", | |
| "\n", | |
| "print(loss_fn(params, x_flattened[:128], jnp.ones((128, 1)) * -5., key))" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 450, | |
| "metadata": {}, | |
| "outputs": [ | |
| { | |
| "name": "stderr", | |
| "output_type": "stream", | |
| "text": [ | |
| "100%|██████████| 5000/5000 [01:13<00:00, 67.73it/s, val=2330.8193] \n" | |
| ] | |
| } | |
| ], | |
| "source": [ | |
| "n_steps = 5000\n", | |
| "n_batch = 128\n", | |
| "\n", | |
| "lr_schedule = optax.warmup_cosine_decay_schedule(\n", | |
| " init_value=0.0,\n", | |
| " peak_value=1e-3,\n", | |
| " warmup_steps=int(0.1 * n_steps),\n", | |
| " decay_steps=n_steps - int(0.1 * n_steps),\n", | |
| " end_value=0.0\n", | |
| ")\n", | |
| "\n", | |
| "opt = optax.adam(learning_rate=lr_schedule)\n", | |
| "opt_state = opt.init(params)\n", | |
| "\n", | |
| "with trange(n_steps) as steps:\n", | |
| " for step in steps:\n", | |
| "\n", | |
| " # Draw a random batches from x\n", | |
| " key, subkey = jax.random.split(key)\n", | |
| " idx = jax.random.choice(key, x.shape[0], shape=(n_batch,))\n", | |
| " \n", | |
| " x_batch = x_flattened[idx]\n", | |
| "\n", | |
| " # Draw random batch of betas\n", | |
| " key, subkey = jax.random.split(key)\n", | |
| " log_beta_batch = jax.random.uniform(key, shape=(n_batch, 1), minval=-3., maxval=1.)\n", | |
| "\n", | |
| " # Get loss and update\n", | |
| " loss, grads = jax.value_and_grad(loss_fn)(params, x_batch, log_beta_batch, key)\n", | |
| " updates, opt_state = opt.update(grads, opt_state, params)\n", | |
| "\n", | |
| " params = optax.apply_updates(params, updates)\n", | |
| "\n", | |
| " steps.set_postfix(val=loss)" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 531, | |
| "metadata": {}, | |
| "outputs": [ | |
| { | |
| "data": { | |
| "text/plain": [ | |
| "Text(0, 0.5, '$D$')" | |
| ] | |
| }, | |
| "execution_count": 531, | |
| "metadata": {}, | |
| "output_type": "execute_result" | |
| }, | |
| { | |
| "data": { | |
| "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjIAAAGwCAYAAACzXI8XAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy89olMNAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA9zUlEQVR4nO3deXiU5b3/8c/MZJmEJBMSCDOBgBHZYmRfTJFTKyixlKJytCLU9UCl+HNv1R4VUY+obT3WVkFpC3IQsXqkilY8IIoigbBFNo0QgmwJ0cRMwjIhzDy/P4BIyMKEJPPMTN6v68pVM8+d4RsfdT69n/v+3hbDMAwBAACEIKvZBQAAAJwrggwAAAhZBBkAABCyCDIAACBkEWQAAEDIIsgAAICQRZABAAAhK8LsAlqbz+fTgQMHFB8fL4vFYnY5AADAD4ZhqLKyUqmpqbJaG553Cfsgc+DAAaWlpZldBgAAOAd79+5Vly5dGrwe9kEmPj5e0om/EQkJCSZXAwAA/FFRUaG0tLSaz/GGhH2QOfU4KSEhgSADAECIOduyEBb7AgCAkEWQAQAAIYsgAwAAQhZBBgAAhCyCDAAACFkEGQAAELIIMgAAIGQRZAAAQMgiyAAAgJAV9p19A8nrM5RbWKaSSo9S4u0amp4km5WDKgEAaC0EmRaydGuRZizZriK3p+Y1l8Ou6WMzlJ3pMrEyAADCF4+WWsDSrUWaumBjrRAjScVuj6Yu2KilW4tMqgwAgPBGkGkmr8/QjCXbZdRz7dRrM5Zsl9dX3wgAANAcBJlmyi0sqzMTczpDUpHbo9zCssAVBQBAG0GQaabiioZDzOlKKv0bBwAA/EeQaYalW4v0xHvb/BqbEm9v5WoAAGh72LV0jk4t8D3byheLJKfjxFZsAADQspiROQeNLfA93akOMtPHZtBPBgCAVkCQOQdnW+B7SlK7KM2aNJA+MgAAtBKCzDlYvr3Yr3EPj+lDiAEAoBURZJrI6zO0OG+/X2OdjphWrgYAgLaNINNEuYVlKjtcfdZxye2iWOALAEArI8g0kb+Plcb1T2WBLwAArYwg0wRen6FF6/f6NfbyDGcrVwMAAAgyTfCXFTt0uMp71nE8VgIAIDAIMn7y+gzN/Xy3X2N5rAQAQGAQZPyUW1im8qNnX+QrSY6YyFauBgAASAQZvzXl0MfXc/fI6ztb318AANBcBBk/NeXQx+KKKq3ZVdqK1QAAAIkg47eh6UlyOfwPM9Ne26ilW4tasSIAAECQ8ZPNatH0sRl+jy8/Wq2pCwgzAAC0JoJME2RnuvTSDQPVlA1JM5ZsZ70MAACthCDTRD/t69JfJgzwa6whqcjtUW5hWesWBQBAG0WQOQc/7Zuq2ZMGKtHPbdZN2fEEAAD8R5A5R9mZLr14w0C/xjZlxxMAAPAfQaYZLu6eLJfDrsaWzLgcdvl8ht7J26+cglLWywAA0IIshmGE9SdrRUWFHA6H3G63EhISWvz9l24t0tQFGyWdWBNzpphIq45W+2q+dznsmj42Q9mZrhavBQCAcOHv5zczMs2UnenSrEkD5Tyjx4w98sTf2tNDjCQVuz1sywYAoIVEmF1AOMjOdOnyDKdyC8tUUulRh7ho3fuPPHmqq+qMNSRZdGJb9uUZTg6XBACgGZiRaSE2q0VZ3ZM1rn9nWS0WHayoG2JOYVs2AAAtgyDTCvzdbs22bAAAmocg0wr83W7NtmwAAJqHINMKTh0webZt2UPTkwJWEwAA4Ygg0wpOP2CyoTBzYWoCC30BAGgmgkwraWhbdmLsiWMNln9Zovk5u02oDACA8MH261Z05rbslPgTj5NmryzQ7z/M12PvblPnxBhd2iulzhhmawAAODs6+5rAMAw9+L9b9Mb6vYqyWRUfE6HSQ8dqrtP9FwDQ1tHZN4hZLBY9eXWm+rjidczrqxViJLr/AgDgL9ODzP79+zVp0iQlJycrJiZGF110kdavX19z/eabb5bFYqn1lZ2dbWLFLcNqsajs8LF6r52aIpuxZDuHTAIA0AhT18h8//33Gj58uH7yk5/ogw8+UMeOHbVjxw61b9++1rjs7GzNnTu35vvo6OhAl9ricgvL/O7+m9U9OXCFAQAQQkwNMs8884zS0tJqhZT09PQ646Kjo+V0Ov16z6qqKlVV/RAQKioqml9oK6D7LwAAzWfqo6V3331XgwcP1rXXXquUlBQNGDBAc+bMqTPuk08+UUpKinr16qWpU6eqtLS0wfecOXOmHA5HzVdaWlpr/grnjO6/AAA0n6m7luz2Ex/S9957r6699lqtW7dOd911l2bPnq2bbrpJkrRo0SLFxsYqPT1dBQUF+t3vfqe4uDjl5OTIZrPVec/6ZmTS0tKCateSJHl9hi55ZoWK3R41dANcDrtWPXAZW7EBAG2Ov7uWTA0yUVFRGjx4sFavXl3z2p133ql169YpJyen3p/ZtWuXunfvruXLl2vkyJFn/TOCcfv1KUu3Fmnqgo2SVG+YefhnffQfl5wf2KIAAAgCIbH92uVyKSMjo9Zrffr00Z49exr8mfPPP18dOnTQzp07W7u8VtdQ998o24nbMnfVbtbIAADQCFMX+w4fPlz5+fm1Xvv666/VrVu3Bn9m3759Ki0tlcsVHs3i6uv+26tTvP599mrt+u6wpszfoEVTLpY9su5jNAAA2jpTZ2TuuecerVmzRk899ZR27typhQsX6pVXXtG0adMkSYcOHdJvfvMbrVmzRrt379ZHH32kcePG6YILLtDo0aPNLL1F2awWZXVP1rj+nZXVPVlJcVH6281D5IiJVN7ect3/5hcK8wbMAACcE1ODzJAhQ7R48WK9/vrryszM1BNPPKHnn39eEydOlCTZbDZt3rxZP//5z9WzZ0/ddtttGjRokD777LOw6CXTmPQO7TR70iBF2ix6b3ORnl++w+ySAAAIOpy1FOT+sX6vfvvWZknSn67vr5/1TeWASQBA2AuJXUuBEOpBRpJmfvClXl65SxFWixwxkSo9zAGTAIDwFhK7luCfB0b3Vr8uDh33GbVCjMQBkwCAto0gEwIMScUV9W/D5oBJAEBbRpAJAU05YBIAgLaEIBMCOGASAID6EWRCAAdMAgBQP4JMCBianiSXw67GNlm7HCe2YgMA0JYQZEKAzWrR9LEnzqRqKMzcd3lP+skAANocgkyIaOiAyVPh5Z95B3Tc6zOjNAAATENDvBDj9Rm1Ovs6YiL177NX68gxryaPSNd/jsk4+5sAABDk/P38NvX0azTdqQMmT/eHa/vp169t1JzPCpXZ2aFx/TubVB0AAIHFo6Uw8NOLXPr1pd0lSQ/872ZtO+A2uSIAAAKDIBMm7ruily7t1VGeap+mzN+gsjOOMgAAIBwRZMKEzWrRn34xQN2SY7W//KjuWLiRxb8AgLBHkAkjjthIvfLLwYqNsml1Qame/uArs0sCAKBVEWTCTC9nvJ67rp8k6a+rCrV40z6TKwIAoPUQZMJQdqZLd/zkAknSg/+7RV/sLVdOQaneyduvnIJSTskGAIQNtl+HqXsu76ltB9z6OP9bXf3S5zo9u7gcdk0fm6HsTJd5BQIA0AKYkQlTNqtFP+uXKkk6cwKm2O3R1AUbtXRrkQmVAQDQcggyYcrrM/SHD/PrvXYq18xYsp3HTACAkEaQCVO5hWUqcnsavG5IKnJ7lFtYFriiAABoYQSZMFVS2XCIOZdxAAAEI4JMmEqJt599UBPGAQAQjAgyYWpoepJcDrssjYxxOewamp4UsJoAAGhpBJkwZbNaNH1shiQ1GGZ+3j9VNmtjUQcAgOBGkAlj2ZkuzZo0UE5H7cdH7aJtkqTX1+7Rvu+PmFEaAAAtwmIYRljvv62oqJDD4ZDb7VZCQoLZ5ZjC6zOUW1imkkqPUuLt6p+WqOvnrNEXe8s1sGui3vhVliJtZFoAQPDw9/ObT682wGa1KKt7ssb176ys7smKibLpLxMGKN4eoY17yvXcsq/NLhEAgHNCkGmj0pJi9cz4vpKkWZ8UaOXX35pcEQAATUeQacN+epFLky7uKkm69408lVTQUwYAEFoIMm3cw2My1NsZr9LDx3TXojyOLAAAhBSCTBtnj7TpxYkDFRtlU86uUr348U6zSwIAwG8EGah7xzg9eVWmJOn55V9r7a5SkysCAMA/BBlIkq4Z2EXjB3aRz5DuXLRJZYePmV0SAABnRZBBjcfHXajzO7bTwYoq3f/mF/KxXgYAEOQIMqjRLjpCL94wUFERVq34qkR/W1VodkkAADSKIINa+rgS9OjPTpzR9MzSr7Thm++VU1Cqd/L2K6eglF1NAICgEmF2AQg+E4d1VU5Bqd7fUqTrZq+W97Ts4nLYNX1shrIzXeYVCADASczIoA6LxaKRfVIkqVaIkaRit0dTF2zU0q1FJlQGAEBtBBnU4fUZ+v2H+fVeO5VrZizZzmMmAIDpCDKoI7ewTEXuho8rMCQVuT3KLSwLXFEAANSDIIM6Sir9O3PJ33EAALQWggzqSIm3t+g4AABaC0EGdQxNT5LLYZelkTEuh11D05MCVhMAAPUhyKAOm9Wi6WNP9JJpKMzccdkFslkbizoAALQ+ggzqlZ3p0qxJA+V01H58FGk7EV7e2XSAXUsAANPREA8Nys506fIMp3ILy1RS6VFKvF3OBLvG/mWVcneX6S8rduquUT3MLhMA0IYxI4NG2awWZXVP1rj+nZXVPVnpHdvpyasyJUl/+uhrrdvNFmwAgHkIMmiyqwZ01jUDO8tnSHcvypP7SLXZJQEA2iiCDM7J4+MydV5yrPaXH9XvFm+RYbBeBgAQeAQZnJO46Aj96foBirBa9P6WIr2xbq/ZJQEA2iCCDM5Zv7RE/WZ0L0knzl7aWXLI5IoAAG0NQQbNMnnE+RrRo4OOVnv1/17fJE+11+ySAABtCEEGzWK1WvTHa/spqV2Uviyq0DNLvzK7JABAG0KQQbOlJNj1x2v7SZLmfr5bH39VYnJFAIC2giCDFvGT3im6Zfh5kqT73/xCJRWcjA0AaH0EGbSYB6/srT6uBJUePqZ7//GFfBxhAABoZaYHmf3792vSpElKTk5WTEyMLrroIq1fv77mumEYevTRR+VyuRQTE6NRo0Zpx44dJlaMhkRH2PTnCQMUE2nTqp3fac5nu8wuCQAQ5kwNMt9//72GDx+uyMhIffDBB9q+fbv++Mc/qn379jVjnn32Wb3wwguaPXu21q5dq3bt2mn06NHyeHh0EYwuSImrOTn79x/m64u95eYWBAAIaxbDxJasDz74oD7//HN99tln9V43DEOpqam67777dP/990uS3G63OnXqpHnz5un6668/659RUVEhh8Mht9uthISEFq0f9TMMQ9MWbtS/thSrW3Ks3r9zhOKiOZ8UAOA/fz+/TZ2ReffddzV48GBde+21SklJ0YABAzRnzpya64WFhSouLtaoUaNqXnM4HBo2bJhycnLqfc+qqipVVFTU+kJgWSwWzby6rzonxuib0iN65J9blFNQqnfy9iunoFRe1s4AAFqIqUFm165dmjVrlnr06KEPP/xQU6dO1Z133qlXX31VklRcXCxJ6tSpU62f69SpU821M82cOVMOh6PmKy0trXV/CdTLERupP13fXxZJizcd0IQ5a3TXojxNmLNGlzyzQku3FpldIgAgDJgaZHw+nwYOHKinnnpKAwYM0JQpUzR58mTNnj37nN/zoYcektvtrvnau5czgMzy3aEq1Tf3Uuz2aOqCjYQZAECzmRpkXC6XMjIyar3Wp08f7dmzR5LkdDolSQcPHqw15uDBgzXXzhQdHa2EhIRaXwg8r8/QjCXb6712KtzMWLKdx0wAgGYxNcgMHz5c+fn5tV77+uuv1a1bN0lSenq6nE6nPvroo5rrFRUVWrt2rbKysgJaK5omt7BMRe6Gd5YZkorcHuUWlgWuKABA2DF1K8k999yjH/3oR3rqqad03XXXKTc3V6+88opeeeUVSScWjd5999168skn1aNHD6Wnp+uRRx5RamqqrrrqKjNLx1mUVPq3Pd7fcQAA1MfUIDNkyBAtXrxYDz30kB5//HGlp6fr+eef18SJE2vG/Pa3v9Xhw4c1ZcoUlZeX65JLLtHSpUtlt9tNrBxnkxLv3/3xdxwAAPUxtY9MINBHxhxen6FLnlmhYren3gW/kuRMsOvzBy+TzWoJaG0AgOAXEn1kEL5sVktNh9+GYkpqol1kGABAcxBk0GqyM12aNWmgnI7aj4+S46Jks0ob95Trb6sKTaoOABAO6BuPVpWd6dLlGU7lFpappNKjlHi7hqYn6bW13+jRd7Zp5gdf6aLODg07P9nsUgEAIYgZGbQ6m9WirO7JGte/s7K6J8tmteiXF3fTVf1T5fUZuuP1TSqpYPcSAKDpCDIwhcVi0VPXXKReneL1bWWVfv3aRlV7fWaXBQAIMQQZmCY2KkKzfzlI8dERWv/N95r5r6/MLgkAEGIIMjBVeod2+uN1/SRJf/+8UO9tPmByRQCAUEKQgemuuNCpqZd2lyT99q3N2nGw0uSKAAChgiCDoHDf5T31o+7JOnLMq9sXbNChquNmlwQACAEEGQSFCJtVL0wYIGeCXQXfHtZv3/pCYd50GgDQAggyCBod4qL10qSBirRZ9K8txTTLAwCcFUEGQWVg1/Z65GcnjjaY+cFXWrur1OSKAADBjCCDoHN6s7xpCzfpIM3yAAANIMgg6JxqltfbGa/vDlVpGs3yAAANIMggKMVGRWjWJJrlAQAaR5BB0DqzWd67efuVU1Cqd07+r9fHriYAaOssRpjvca2oqJDD4ZDb7VZCQoLZ5eAcPLP0K836pEAWSaf/w+py2DV9bIayM11mlQYAaCX+fn4zI4Ogl5l64h/gMxN3sdujqQs2aunWosAXBQAICgQZBDWvz9CT739Z77VTwWbGku08ZgKANoogg6CWW1imInfD268NSUVuj3ILywJXFAAgaBBkENRKKv3rIePvOABAeCHIIKilxNtbdBwAILwQZBDUhqYnyeWwy9LImJT4aA1NTwpYTQCA4EGQQVCzWS2aPvbE2UuNhZmKo9WBKQgAEFQIMgh62ZkuzZo0UE5H7cdHnRKilRQbpZLKKk35n/XyVHtNqhAAYBYa4iFkeH2GcgvLVFLpUUq8XUPTk1Tw7SGNn7ValZ7j+nm/VD3/i/6yWhubuwEAhAIa4iHs2KwWZXVP1rj+nZXVPVk2q0U9O8Vr9qRBirBa9O4XB/Tcsq/NLhMAEEAEGYS84Rd00FPXXCRJ+svHO/WPdXtNrggAECgEGYSF6wan6Y6fXCBJ+t3iLVq14zuTKwIABAJBBmHjvit66uf9UnXcZ2jqgg36+mCl2SUBAFoZQQZhw2Kx6PfX9tXQ85JUWXVct8xdR8dfAAhzBBmElegIm17+5SCld2in/eVH9R+vrteRY8fNLgsA0EoIMgg77dtFae7NQ5TULkqb97l15+t5nI4NAGGKIIOwdF6Hdppz4yBFRVi1/MuDevL97WaXBABoBQQZhK1B3ZL03HX9JElzP9+teZ8XmlwRAKClEWQQ1n7WN1W/ze4lSXr8ve1avv2gyRUBAFoSQQZhb+qPu+v6IWnyGdL/e32Ttuxzm10SAKCFEGQQ9iwWi564KlMjenTQ0Wqvbn11nfaUHVFOQaneyduvnIJSFgMDQIji0Ei0GRWeal03O0dfFVcqwmrR8dPCi8th1/SxGcrOdJlYIQDgFA6NBM6QYI/UjVndJKlWiJGkYrdHUxds1NKtRWaUBgA4RwQZtBlen6E/r9hZ77VTsWbGku08ZgKAEEKQQZuRW1imInfDRxYYkorcHuUWlgWuKABAsxBk0Gb4e+4S5zMBQOggyKDNSIm3t+g4AID5CDJoM4amJ8nlsMvSyJjkdlEamp4UsJoAAM1DkEGbYbNaNH1shiQ1GGYOVx1XfnFl4IoCADQLQQZtSnamS7MmDZTTUfvxkdNh1wUd4+Q57tNNc3O1p/SISRUCAJqChnhok7w+Q7mFZSqp9Cgl3q6h6Uk6VHVcv3j5RMO885Jj9dbUH6lDXLTZpQJAm0RDPKARNqtFWd2TNa5/Z2V1T5bNapEjJlKv3jpUXdrHaHfpEd0yd50OVR03u1QAQCMIMsBpOiXYNf/WoUpqF6Ut+9361f+sV9Vxr9llAQAaQJABznB+xzjNu2WI2kXZ9PnOUt37jy/o9gsAQYogA9Sjb5dEvfzLwYq0WfT+5iLNWLJNYb6cDABCUpOCzEcffaSLL75Ydrtd8fHxGjJkiJ555hlVVrJdFeHnkh4d9Nx1/WWxSPNzvmnwnCYAgHn8DjJr167VlVdeqejoaD388MN65JFH1LdvX/3hD39QZmamNm/e3Jp1AqYY2y9Vj429UJL03LKvtXDtHpMrAgCczu/t1+PHj5fVatWbb75Z6/UjR47oV7/6lT755BNt2bJFiYmJrVHnOWP7NVrCH/8vX39esVNWi/TSxIHKznSZXRIAhLUW336dk5OjO+64o87rsbGxevXVV9WlSxfNnj373KoFgty9l/fUhKFp8hnSnYvytGZXqdklAQDUhCDz7bffKj09vf43sVp111136f3332+xwoBgYrFY9MS4TF2R0UnHjvs0+dX12nbAbXZZANDm+R1kvF6v7PaGTwUeNGiQ8vPzm/SHP/bYY7JYLLW+evfuXXP90ksvrXP99ttvb9KfAbSUCJtVL0wYoKHpSaqsOq6b/r6OowwAwGRN2rU0f/58rV27Vh6Pp861hIQElZeXN7mACy+8UEVFRTVfq1atqnV98uTJta4/++yzTf4zgJZij7Rpzo2D1dsZr+8OVemXf1+rbyur5PUZyiko1Tt5+5VTUErfGQAIkAh/B44YMUJPPPGEKisrFRERoV69emnQoEEaOHCgBg0apE6dOsnrbXoH1IiICDmdzgavx8bGNnodCDRHTKTm3zpU42ev1jelR3TNS5/rmNengxVVNWNcDrumj81gUTAAtDK/Z2RWrlwpt9ut/Px8zZ8/X1deeaX27dunxx57TCNGjFCvXr3OqYAdO3YoNTVV559/viZOnKg9e2pvb33ttdfUoUMHZWZm6qGHHtKRI41P5VdVVamioqLWF9DSUhLsmn/rMMVFR2jv90drhRhJKnZ7NHXBRi3dWmRShQDQNvg9I3NKjx491KNHD11//fU1rxUWFmr9+vXatGlTk95r2LBhmjdvnnr16qWioiLNmDFDI0aM0NatWxUfH68bbrhB3bp1U2pqqjZv3qwHHnhA+fn5evvttxt8z5kzZ2rGjBlN/bWAJuuaFKvoSKsOVdW9ZkiySJqxZLsuz3DKZrUEujwAaBP87iMTCOXl5erWrZuee+453XbbbXWur1ixQiNHjtTOnTvVvXv3et+jqqpKVVU/fLJUVFQoLS2NPjJocTkFpZowZ81Zx70++WJldU8OQEUAED787SPT5BmZ1pSYmKiePXtq5876W8EPGzZMkhoNMtHR0YqOjm61GoFTSirrLnpvzjgAQNMF1aGRhw4dUkFBgVyu+hdI5uXlSVKD14FASolvuB3BuYwDADSdqUHm/vvv18qVK7V7926tXr1aV199tWw2myZMmKCCggI98cQT2rBhg3bv3q13331XN954o/7t3/5Nffv2NbNsQJI0ND1JLoddja1+cTnsGpqeFLCaAKCtMTXI7Nu3TxMmTFCvXr103XXXKTk5WWvWrFHHjh0VFRWl5cuX64orrlDv3r113333afz48VqyZImZJQM1bFaLpo/NkKQGw8wvhqSx0BcAWlFQLfZtDRwaida2dGuRZizZriL3D2th7BFWeY77FBcdoUVTLlZmZ4eJFQJA6PH385sgA7QAr89QbmGZSio9Som3q28Xh257dZ3W7CpTh7govXX7j3Reh3ZmlwkAIYMgcxJBBmap8FTrFy+v0ZdFFeqaFKu3pmax8BcA/OTv53dQ7VoCwkmCPVKv3jpEXZNitafsiG76+zpVeKrNLgsAwgpBBmhFKfF2/c9tQ9UhLkpfFlVo8qvr5alu+plkAID6EWSAVtYtuZ3m3TJUcdERWltYprsX5XE6NgC0EIIMEACZnR165cZBirJZtXRbsR7+51aF+fI0AAgIggwQID/q3kHPX99fFov0eu4e/feyr80uCQBCHkEGCKCfXuTSE+MyJUkvrNipV1fvNrcgAAhxBBkgwCZd3E33jOopSXpsyTa9t/mAyRUBQOgiyAAmuHPkBboxq5sMQ7rnjTyt2vGd2SUBQEgiyAAmsFgsmj72Qo25yKVqr6Ff/c96bd5XbnZZABByCDKASWxWi577RT8NvyBZh495dcvcddpZckg5BaV6J2+/cgpK2aYNAGfBEQWAyQ5VHdf1r+Ro6/4K2SyS97R/I10Ou6aPzVB2psu8AgHABBxRAISIuOgI3Zh1nqTaIUaSit0eTV2wUUu3FgW+MAAIAQQZwGRen9FgT5lTuWbGku08ZgKAehBkAJPlFpapyO1p8LohqcjtUW5hWeCKAoAQQZABTFZS2XCIOZdxANCWEGQAk6XE21t0HAC0JQQZwGRD05PkcthlaWSMIyZSQ9OTAlYTAIQKggxgMpvVouljMySpwTBTcbRaK74qCVxRABAiCDJAEMjOdGnWpIFyOmo/PnI57Mo6P1mGpDsWbtSGb1jwCwCnoyEeEES8PkO5hWUqqfQoJd6uoelJ8hmGpsxfr4/zv1VibKTeuj1LF6TEm10qALQqfz+/CTJACDhy7LgmzFmrL/aWq3NijP536o/qzN4AQDihsy8QRmKjIjT35iE6v0M77S8/qpvn5sp9tNrssgDAdAQZIEQktYvSq7cOVcf4aH1VXKkp89fLU+01uywAMBVBBgghaUmxmnfLEMVFR2htYZnu/UceRxcAaNMIMkCIuTDVoVd+OUiRNov+taVYM5ZsU5gvdQOABhFkgBD0ows66Lnr+kuS5ud8o5c+KTC3IAAwCUEGCFFj+6Xq0Z+daKT3+w/z9eb6vSZXBACBR5ABQtitl6TrVz8+X5L04Ntb9DHdfwG0MQQZIMQ9MLq3rhnQWV6foV+/tlGb9nwvr89QTkGp3snbr5yCUhYEAwhbEWYXAKB5rFaLnvn3vvru8DF9+vW3mvTXtYqNsunbQ8dqxrgcdk0fm6HsTJeJlQJAy2NGBggDkTarZk0cqG7JsTp8zFsrxEhSsdujqQs2aunWIpMqBIDWQZABwoQ90qajx+pvkHfqwdKMJdt5zAQgrBBkgDBx4rDJqgavG5KK3B7lFnKCNoDwQZABwkRJpadFxwFAKCDIAGEiJd6/07D9HQcAoYAgA4SJoelJcjnssjQyxuWwa2h6UsBqAoDWRpABwoTNatH0sSc6/TYUZq7I6CSbtbGoAwChhSADhJHsTJdmTRoop6P246N20TZJ0oK1e/RxPt1/AYQPixHmx+ZWVFTI4XDI7XYrISHB7HKAgPD6jJO7mDxKibdrcLf2+u3/btbiTfsVE2nT61MuVv+0RLPLBIAG+fv5TZAB2ohjx3267dV1+mzHd0pqF6W3bs/S+R3jzC4LAOrl7+c3j5aANiIqwqrZkwapbxeHyg4f041/z1VJBVuxAYQ2ggzQhrSLjtDfbx6i85Jjte/7o7pp7jpVeKrNLgsAzhlBBmhjOsRFa/6tw9QhLlpfFlVoyvz1qjpe/9EGABDsCDJAG9Q1OVbzbhmiuOgIrdlVpnvf+EI+zmACEIIIMkAbldnZodmTBinSZtH7W4r0+HvbFeZr/wGEIYIM0IZd0qOD/nhdf0nSvNW7NWtlgbkFAUATRZhdAABz/bxfqr6rrNLj723Xs0vz1TEuWtcM7FKrD83Q9CQ6AgMISvSRASBJmvnBl3p55S5ZLZIjJlLfH/lhN5PLYdf0sRnKznSZWCGAtoQ+MgCa5MHs3hqWniSfoVohRpKK3R5NXbBRS7cWmVQdANSPIANAkuQzpG9Kj9R77dS07Ywl2+VldxOAIEKQASBJyi0sU3EjnX4NSUVuj3ILywJXFACcBUEGgCSppNK/4wr8HQcAgUCQASBJSom3t+g4AAgEggwASdLQ9CS5HHY1tsna6TixFRsAggVBBoAkyWa1aPrYDElqMMx0iIuSL7w7NgAIMaYGmccee0wWi6XWV+/evWuuezweTZs2TcnJyYqLi9P48eN18OBBEysGwlt2pkuzJg2U01H78VFyuyhFWC3aur9C97/JuUwAgofpnX0vvPBCLV++vOb7iIgfSrrnnnv0/vvv680335TD4dAdd9yha665Rp9//rkZpQJtQnamS5dnOOt09v306281ef56vZN3QPH2CD0xLlMWC91+AZjL9CATEREhp9NZ53W3262//e1vWrhwoS677DJJ0ty5c9WnTx+tWbNGF198caBLBdoMm9WirO7JtV77Se8UPfeL/rpr0SYtWLNHjphI/WZ07wbeAQACw/Q1Mjt27FBqaqrOP/98TZw4UXv27JEkbdiwQdXV1Ro1alTN2N69e6tr167Kyclp8P2qqqpUUVFR6wtAy/h5v1Q9eVWmJOnFjwv0yqccMgnAXKYGmWHDhmnevHlaunSpZs2apcLCQo0YMUKVlZUqLi5WVFSUEhMTa/1Mp06dVFxc3OB7zpw5Uw6Ho+YrLS2tlX8LoG2ZOKybfpvdS5L01L++0qLcPSZXBKAtM/XR0pVXXlnz13379tWwYcPUrVs3/eMf/1BMTMw5vedDDz2ke++9t+b7iooKwgzQwn596QVyH63Wyyt36aHFWxRvj9SYvhwoCSDwTF8jc7rExET17NlTO3fu1OWXX65jx46pvLy81qzMwYMH611Tc0p0dLSio6MDUC3Qtj2Y3VsVR4/r9dw9uvuNTWoXbdOIHh3rLBK2WVkQDKD1BFWQOXTokAoKCvTLX/5SgwYNUmRkpD766CONHz9ekpSfn689e/YoKyvL5EoBWCwWPXlVpio91Xpvc5Emz1+vBHukSg8fqxnjctg1fWyGsjOZrQHQOkxdI3P//fdr5cqV2r17t1avXq2rr75aNptNEyZMkMPh0G233aZ7771XH3/8sTZs2KBbbrlFWVlZ7FgCgoTNatFz1/XXhakJqvYatUKMJBW7PZq6YKOWbi0yqUIA4c7UGZl9+/ZpwoQJKi0tVceOHXXJJZdozZo16tixoyTpv//7v2W1WjV+/HhVVVVp9OjReumll8wsGcAZbFaLSg9V1XvN0IkuwTOWbNflGU4eMwFocRbDCO9+4xUVFXI4HHK73UpISDC7HCDs5BSUasKcNWcd9/rki+v0pgGAhvj7+W16HxkAoa2k0tOi4wCgKQgyAJolJd5+9kFNGAcATUGQAdAsQ9OT5HLYGzwxW5I6JURraHpSwGoC0HYQZAA0i81q0fSxGZLUcJgx1OCCYABoDoIMgGbLznRp1qSBcjpqPz5KiY9WUmyUDlZW6fo5a1gnA6DFsWsJQIvx+ow6nX33f39U17+SowNuj7p3bKdFU7LUMZ7u2wAa5+/nN0EGQKvbU3pEv3glR0Vuj3qkxGnh5IsJMwAaxfZrAEGja3KsFk25WM4Eu3aUHNINc9boO9bMAGgBBBkAAdEtuZ0WTblYnRKia8IMC4ABNBdBBkDAnNfhxBqZlPhofX3wkCb+da1KD1XJ6zOUU1Cqd/L2K6egVF5fWD/xBtCCWCMDIOAKvj2kCa+sUUlllTon2lXtNVRS+cPsDKdmA2CNDICg1b3jiQW/CfYI7S/31AoxEqdmA/AfQQaAKdI7tFNURP3/CTo1TTxjyXYeMwFoFEEGgClyC8v03aFjDV43JBW5PcotLAtcUQBCDkEGgCk4NRtASyDIADAFp2YDaAkEGQCm8OfU7A5xUZyaDaBRBBkApvDn1OxKz3GtLSwNXFEAQg5BBoBpGjo125lgV29nnKqO+3Tz3HX6cFuxSRUCCHY0xANguvpOzT7u8+mu1/O0dFuxrBbpmfF9de3gNLNLBRAgnH59EkEGCF3HvT499PYWvblhnyTp4TF99B8jzje5KgCB4O/nd0QAawKAJomwWfXsv/dVYmyk5nxWqCff/1Luo9W69/Ke8hmqM4tjsza2dBhAOCLIAAhqFotFv/tpHyXGRun3H+brzyt2avM+t/KLK1Vc8UOPGc5nAtomFvsCCHoWi0XTfnKBnrwqU5K08utva4UYifOZgLaKIAMgZEwY2lWJMZH1XuN8JqBtIsgACBm5hWUqP1rd4HXOZwLaHoIMgJDB+UwAzkSQARAyOJ8JwJkIMgBChj/nM8VF2zS4W/uA1QTAXAQZACHDn/OZDlV5NfW1Dar0NLyWBkD4IMgACCkNnc/kcth1y/DzFBVh1fIvSzR+1mrtKT1iUpUAAoUjCgCEpPrOZ7JZLfpib7kmz1+vksoqtY+N1EsTBymre3KD4wEEJ85aOokgA7Q9Bys8mjJ/vb7Y51aE1aJ/H9xFK/O/VZGbTsBAqPD385tHSwDCTqcEu974VZbG9U/VcZ+hRbl7a4UYiU7AQLggyAAIS/ZIm/54bT/FR9d/pBydgIHwQJABELbW7f5elVXHG7xOJ2Ag9BFkAIQtOgED4Y8gAyBs0QkYCH8EGQBhy59OwJK089tKhfkGTiBsEWQAhC1/OgFL0iP/3KbbXl3PIyYgBBFkAIS1xjoBv3TDQD08po+iIqxa8VWJsp//TB9uK64Z4/UZyiko1Tt5+5VTUMruJiAI0RAPQJvQWGff/OJK3f1Gnr4sqpAkXTe4i7K6d9CzS7+iiR5gEjr7nkSQAeCPquNePbfsa73y6S419F/FU4+nZk0aSJgBWhmdfQGgCaIjbHroyj5aeNswNXQEE030gOBDkAGA01ksaiyj0EQPCC4EGQA4DU30gNBCkAGA0/jbHK+hM5wABBZBBgBO428Tvd+89YUWrt3DWhnAZAQZADhNY030Tn3fKT5apYer9bvFWzTmhc+0asd3NWPoPQMEFtuvAaAeS7cWacaS7fX2kbmsdyctWPON/vTRDrmPVkuSRvZO0YieHfTyyl30ngFaAH1kTiLIADhXjTXRk6TvDx/Tnz7aof9Z802DMy/0ngHODUHmJIIMgNaWX1ypn/9llaqO++q9bpHkdNi16oHLagUhAA2jIR4ABEjZ4WMNhhiJ3jNAayLIAEAz+dtTZtsBdytXArQ9NEIAgGbyt/fMk+9/qbWFZZo84nwNOa+9LJYfHjOdbT0OgPoRZACgmU71nil2e9TQosPoCKuqjvu0bPtBLdt+UP26OPQfI87XlZlOLf/yYIM7pFggDDSOxb4A0AKWbi3S1AUbJalWmDl919IFKfH626pCvb1xX82amqTYKJUdOVbn/djthLYu5Bb7Pv3007JYLLr77rtrXrv00ktlsVhqfd1+++3mFQkADcjOdGnWpIFyOmo/ZnI67DVh5IKUOM285iKtfvAy3T2qh5JiI+sNMRInbQP+CopHS+vWrdPLL7+svn371rk2efJkPf744zXfx8bGBrI0APBbdqZLl2c4z7rWJTkuWneP6qkBXRN109/XNfh+p+92yuqe3MrVA6HJ9CBz6NAhTZw4UXPmzNGTTz5Z53psbKycTqcJlQFA09msFr9DR/mRar/GfbitSP3TEhUTZatzjUXCaOtMDzLTpk3TmDFjNGrUqHqDzGuvvaYFCxbI6XRq7NixeuSRRxqdlamqqlJVVVXN9xUVFa1SNwA0l7+7neat/kZvrt+n0Rc6NbZ/qi65oIMibdZGj1FgXQ3aClODzKJFi7Rx40atW1f/1OoNN9ygbt26KTU1VZs3b9YDDzyg/Px8vf322w2+58yZMzVjxozWKhkAWow/u53iom1yxERqf7lHb2/ar7c37VdSuyhd1NmhlV9/W2d8sdujqQs2skgYbYZpu5b27t2rwYMHa9myZTVrYy699FL1799fzz//fL0/s2LFCo0cOVI7d+5U9+7d6x1T34xMWloau5YABCV/djuNvtCpjXvK9W7efr23uUilh+tfIHz6z3IkAkJd0J+19M9//lNXX321bLYfnvl6vV5ZLBZZrVZVVVXVuiZJhw8fVlxcnJYuXarRo0f79eew/RpAsGvKI6LjXp/+uqpQT3/w1Vnf9/XJFze4Xoe1NQh2/n5+m/ZoaeTIkdqyZUut12655Rb17t1bDzzwQJ0QI0l5eXmSJJeL6VIA4cPf3U6SFGGzyuXwb23Ng29v1jUDuujHvTqqb2eHrCffj7U1CCdB1RDv9EdLBQUFWrhwoX76058qOTlZmzdv1j333KMuXbpo5cqVfr8nMzIAwk1OQakmzFnTpJ9JahelET06qH1spOat/qbOdRrwIdgE/YzM2URFRWn58uV6/vnndfjwYaWlpWn8+PF6+OGHzS4NAEx1tkXCFkkd46N196ge+mzHd1q14zuVHT6md/IONPiexsmfm7Fkuy7PcNY7G8TjKASjoJqRaQ3MyAAIR/4sEj41s1Lt9WnTnnK9tuYbvfNFw2HmlDk3DtLlGbX7d/E4CoEW9It9A4UgAyBcNTVcvJO3X3ctyvPrvXs74zXkvCQNSU/S4arj+t3bW+rM/pztcRQzOGgOgsxJBBkA4awpYeFc1tacTUNbvZnBQXOF/BoZAMDZNeVIBH/W1jgddi3+9XBt2vO9cneX6eP8Eu3+7kiD73nqPKj/21asKy86EVBOPfY6889orFkfszc4V8zIAEAb0pS1NVLTHkelOuzq28WhzwtKVek5Xu+Y+mZwmL1Bffz9/LYGsCYAgMmyM12aNWmgnGf0onE67PXOlPh7HpRF0gG3R0u3HWwwxEi1T/SWfghWp4cY6YfZm6Vbi2q97vUZyiko1Tt5+5VTUCqvL6z/vzj8wIwMALRB/j7K8foMXfLMirM+jlp697/py6IKLVy7R+/6sTNqULdEXZHh1Muf7lJZA0cunDl7w8xN28Ji35MIMgDQPE15HNUaC4pfn3yx3EeP1bvupqFHYqy5CX0s9gUAtIhTj6POnA1x1jMb4s+J3kntonTL8PO0fPtBfbHPfdY//88rvtbW/RX1vl99jfz8mbkh6IQPZmQAAH7x98Pf3xmclp698XfmRhJBJwTwaOkkggwABJ6/syKNrb+RpMTYSA3smqgVX3171j8zuV2kDlV5VXXcV+91y8n3+/5Idb3XpLMHHX8P90TzEWROIsgAgDn8mdXwZ/bGERPV4utu6nO2oGOcvF5+2vXTwxmzOC2LIHMSQQYAgtvZZm/82TmVkhCt8QM766VPdgWs7lN/tiRN+bd0vftFUZ3f4ZExGWrfLqpWuJFE4PEDQeYkggwABL+zzWYE08xNcyTGRkpSnVmdMwPPoG7tteGb7xv8vi2EH4LMSQQZAAgPzZ25CSVWi3R6r78zv09qF6Wr+qdqZJ9OkiF9d7hKHeKi6/x1yaEqlR2qUmJslMqPHKvzv0ntouR0xARlMCLInESQAYDw0dyZG0dspNxHqkM+6LS0eHuEftQ9WTGRNp0ZCwxD+u5Qlb6t9OhQlVeVnmM6esyQ94z3WHjzMP2od4cWq4kgcxJBBgDalsZmbiQRdFrZ7qfHtMj7EGROIsgAQNvT2MzNuQQdNE1LhBmCzEkEGQDAmZoadNqf3JZ9ahv2KWd+jxNa4jETQeYkggwAoKnqCzrLthfXO5Pz834uvfJpoSRCzemaOyvDWUsAAJwjm9WirO7JtV7LznQ12Nl3QNf2dUIOAoMgAwCAn+oLOFL9Ief7w1V64v0v6zyiMlS7jwyahyADAEALqC/kjM501ZnBkXTWwHO2PjLBbuHNwwL2ZxFkAABoJQ3N4Jwt8DTU2Xf59mItztuvssPBPaPTkv1kzobFvgAAhJDTFyLX183Xn86+e8oO6+1N+1XpObOtXfPRR6aFEWQAAKjrVCAqrvDou8oqlR2pUlG5J+Q6+/JoCQCANqihx16hxmp2AQAAAOeKIAMAAEIWQQYAAIQsggwAAAhZBBkAABCyCDIAACBkEWQAAEDIIsgAAICQRZABAAAhK+w7+55qtVxRUWFyJQAAwF+nPrfPdpJS2AeZyspKSVJaWprJlQAAgKaqrKyUw+Fo8HrYHxrp8/l04MABxcfHy2Kx+PUzFRUVSktL0969ezloMshwb4Ib9yd4cW+CG/enLsMwVFlZqdTUVFmtDa+ECfsZGavVqi5dupzTzyYkJPAPVJDi3gQ37k/w4t4EN+5PbY3NxJzCYl8AABCyCDIAACBkEWTqER0drenTpys6OtrsUnAG7k1w4/4EL+5NcOP+nLuwX+wLAADCFzMyAAAgZBFkAABAyCLIAACAkEWQAQAAIYsgc4YXX3xR5513nux2u4YNG6bc3FyzS2oTPv30U40dO1apqamyWCz65z//Weu6YRh69NFH5XK5FBMTo1GjRmnHjh21xpSVlWnixIlKSEhQYmKibrvtNh06dCiAv0V4mjlzpoYMGaL4+HilpKToqquuUn5+fq0xHo9H06ZNU3JysuLi4jR+/HgdPHiw1pg9e/ZozJgxio2NVUpKin7zm9/o+PHjgfxVws6sWbPUt2/fmiZqWVlZ+uCDD2quc1+Cx9NPPy2LxaK777675jXuT8sgyJzmjTfe0L333qvp06dr48aN6tevn0aPHq2SkhKzSwt7hw8fVr9+/fTiiy/We/3ZZ5/VCy+8oNmzZ2vt2rVq166dRo8eLY/HUzNm4sSJ2rZtm5YtW6b33ntPn376qaZMmRKoXyFsrVy5UtOmTdOaNWu0bNkyVVdX64orrtDhw4drxtxzzz1asmSJ3nzzTa1cuVIHDhzQNddcU3Pd6/VqzJgxOnbsmFavXq1XX31V8+bN06OPPmrGrxQ2unTpoqefflobNmzQ+vXrddlll2ncuHHatm2bJO5LsFi3bp1efvll9e3bt9br3J8WYqDG0KFDjWnTptV87/V6jdTUVGPmzJkmVtX2SDIWL15c873P5zOcTqfx+9//vua18vJyIzo62nj99dcNwzCM7du3G5KMdevW1Yz54IMPDIvFYuzfvz9gtbcFJSUlhiRj5cqVhmGcuBeRkZHGm2++WTPmyy+/NCQZOTk5hmEYxr/+9S/DarUaxcXFNWNmzZplJCQkGFVVVYH9BcJc+/btjb/+9a/clyBRWVlp9OjRw1i2bJnx4x//2LjrrrsMw+Dfm5bEjMxJx44d04YNGzRq1Kia16xWq0aNGqWcnBwTK0NhYaGKi4tr3RuHw6Fhw4bV3JucnBwlJiZq8ODBNWNGjRolq9WqtWvXBrzmcOZ2uyVJSUlJkqQNGzaourq61v3p3bu3unbtWuv+XHTRRerUqVPNmNGjR6uioqJm9gDN4/V6tWjRIh0+fFhZWVnclyAxbdo0jRkzptZ9kPj3piWF/aGR/vruu+/k9Xpr/QMjSZ06ddJXX31lUlWQpOLiYkmq996culZcXKyUlJRa1yMiIpSUlFQzBs3n8/l09913a/jw4crMzJR04u99VFSUEhMTa4098/7Ud/9OXcO527Jli7KysuTxeBQXF6fFixcrIyNDeXl53BeTLVq0SBs3btS6devqXOPfm5ZDkAHgt2nTpmnr1q1atWqV2aXgpF69eikvL09ut1tvvfWWbrrpJq1cudLsstq8vXv36q677tKyZctkt9vNLies8WjppA4dOshms9VZMX7w4EE5nU6TqoKkmr//jd0bp9NZZ1H28ePHVVZWxv1rIXfccYfee+89ffzxx+rSpUvN606nU8eOHVN5eXmt8Wfen/ru36lrOHdRUVG64IILNGjQIM2cOVP9+vXTn/70J+6LyTZs2KCSkhINHDhQERERioiI0MqVK/XCCy8oIiJCnTp14v60EILMSVFRURo0aJA++uijmtd8Pp8++ugjZWVlmVgZ0tPT5XQ6a92biooKrV27tubeZGVlqby8XBs2bKgZs2LFCvl8Pg0bNizgNYcTwzB0xx13aPHixVqxYoXS09NrXR80aJAiIyNr3Z/8/Hzt2bOn1v3ZsmVLrbC5bNkyJSQkKCMjIzC/SBvh8/lUVVXFfTHZyJEjtWXLFuXl5dV8DR48WBMnTqz5a+5PCzF7tXEwWbRokREdHW3MmzfP2L59uzFlyhQjMTGx1opxtI7Kykpj06ZNxqZNmwxJxnPPPWds2rTJ+OabbwzDMIynn37aSExMNN555x1j8+bNxrhx44z09HTj6NGjNe+RnZ1tDBgwwFi7dq2xatUqo0ePHsaECRPM+pXCxtSpUw2Hw2F88sknRlFRUc3XkSNHasbcfvvtRteuXY0VK1YY69evN7KysoysrKya68ePHzcyMzONK664wsjLyzOWLl1qdOzY0XjooYfM+JXCxoMPPmisXLnSKCwsNDZv3mw8+OCDhsViMf7v//7PMAzuS7A5fdeSYXB/WgpB5gx//vOfja5duxpRUVHG0KFDjTVr1phdUpvw8ccfG5LqfN10002GYZzYgv3II48YnTp1MqKjo42RI0ca+fn5td6jtLTUmDBhghEXF2ckJCQYt9xyi1FZWWnCbxNe6rsvkoy5c+fWjDl69Kjx61//2mjfvr0RGxtrXH311UZRUVGt99m9e7dx5ZVXGjExMUaHDh2M++67z6iurg7wbxNebr31VqNbt25GVFSU0bFjR2PkyJE1IcYwuC/B5swgw/1pGRbDMAxz5oIAAACahzUyAAAgZBFkAABAyCLIAACAkEWQAQAAIYsgAwAAQhZBBgAAhCyCDAAACFkEGQAAELIIMgAAIGQRZACEpB//+MeyWCyyWCyKiopSnz59tHDhQrPLAhBgBBkAIccwDG3atEl/+MMfVFRUpPz8fGVnZ+vGG29UYWGh2eUBCCCCDICQs2PHDlVWVio7O1tOp1Pp6em67bbb5PV6lZ+fb3Z5AAKIIAMg5GzYsEHt27dXRkaGJGnfvn36z//8T0VHR6tv374mVwcgkCLMLgAAmmrjxo1yu92Kj4+X1+uVx+NRTEyMZs+erdTUVLPLAxBAFsMwDLOLAICmGDlypC688ELdeeedKi8v1/3336/hw4frv/7rv8wuDUCAEWQAhJz27dtr1qxZuv766yVJ27dvV9++fbVz506dd9555hYHIKBYIwMgpOzatUvl5eXKzMyseS0jI0Pdu3dn+zXQBhFkAISUDRs2KDIyUj179qz1+siRI7V48WKTqgJgFoIMgJCyceNG9ejRQ1FRUbVeHzVqlDZs2KB9+/aZVBkAM7BGBgAAhCxmZAAAQMgiyAAAgJBFkAEAACGLIAMAAEIWQQYAAIQsggwAAAhZBBkAABCyCDIAACBkEWQAAEDIIsgAAICQRZABAAAh6/8DNG0aclx+RRkAAAAASUVORK5CYII=", | |
| "text/plain": [ | |
| "<Figure size 640x480 with 1 Axes>" | |
| ] | |
| }, | |
| "metadata": {}, | |
| "output_type": "display_data" | |
| } | |
| ], | |
| "source": [ | |
| "@jax.jit\n", | |
| "def eval(params, z):\n", | |
| " def eval_model(vae):\n", | |
| " return vae.decoder(z)\n", | |
| " return nn.apply(eval_model, vae)(params)\n", | |
| "\n", | |
| "@jax.jit\n", | |
| "def encode_decode(params, x, beta, key):\n", | |
| " x = jnp.concatenate([x, beta], axis=-1)\n", | |
| " def encode_decode_model(vae):\n", | |
| " enc = vae.encoder(x)\n", | |
| " enc = rearrange(enc, 'b (n c) -> b n c', c=2)\n", | |
| " mu, logvar = enc[:, :, 0], enc[:, :, 1]\n", | |
| " z = tfp.distributions.Normal(loc=mu, scale=jnp.exp(0.5 * logvar)).sample(seed=key)\n", | |
| " z = jnp.concatenate([z, beta], axis=-1)\n", | |
| " return vae.decoder(z)\n", | |
| " return nn.apply(encode_decode_model, vae)(params)\n", | |
| "\n", | |
| "def get_RD(params, x, beta, key):\n", | |
| " x = jnp.concatenate([x, beta], axis=-1)\n", | |
| " def encode_decode_model(vae):\n", | |
| " enc = vae.encoder(x)\n", | |
| " enc = rearrange(enc, 'b (n c) -> b n c', c=2)\n", | |
| " mu, logvar = enc[:, :, 0], enc[:, :, 1]\n", | |
| " q_latent = tfp.distributions.Normal(loc=mu, scale=jnp.exp(0.5 * logvar))\n", | |
| " prior_latent = tfp.distributions.Normal(loc=0., scale=1.)\n", | |
| " R = tfp.distributions.kl_divergence(q_latent, prior_latent)\n", | |
| " z = q_latent.sample(seed=key)\n", | |
| " z = jnp.concatenate([z, beta], axis=-1)\n", | |
| " return vae.decoder(z), R\n", | |
| " \n", | |
| " reco, R = nn.apply(encode_decode_model, vae)(params)\n", | |
| " D = reco_gap(reco, x[..., :-1], jnp.ones((x.shape[0], 1)))\n", | |
| " return R.sum(-1).mean(), D.sum(-1).mean()\n", | |
| "\n", | |
| "ii = 20\n", | |
| "\n", | |
| "log_beta_ary = jnp.linspace(-3., 1., 100)\n", | |
| "\n", | |
| "R_list = []\n", | |
| "D_list = []\n", | |
| "\n", | |
| "n_agg = 16\n", | |
| "\n", | |
| "for log_beta in log_beta_ary:\n", | |
| " beta = jnp.ones((n_agg, 1)) * 10 ** log_beta\n", | |
| " R, D = get_RD(params, x_flattened[ii:ii + n_agg], beta, key)\n", | |
| " R_list.append(R)\n", | |
| " D_list.append(D)\n", | |
| "\n", | |
| "plt.plot(R_list, D_list, 'o-')\n", | |
| "plt.xlabel(r'$R$')\n", | |
| "plt.ylabel(r'$D$')" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 532, | |
| "metadata": {}, | |
| "outputs": [], | |
| "source": [ | |
| "ii = 4\n", | |
| "beta = jnp.ones((1, 1)) * 1.e-3\n", | |
| "x_samples = encode_decode(params, x_flattened[ii:ii + 1], beta, key)[0]" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 533, | |
| "metadata": {}, | |
| "outputs": [ | |
| { | |
| "data": { | |
| "text/plain": [ | |
| "<matplotlib.collections.PathCollection at 0x349f36400>" | |
| ] | |
| }, | |
| "execution_count": 533, | |
| "metadata": {}, | |
| "output_type": "execute_result" | |
| }, | |
| { | |
| "data": { | |
| "image/png": "iVBORw0KGgoAAAANSUhEUgAAAi8AAAGfCAYAAACX9jKsAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy89olMNAAAACXBIWXMAAA9hAAAPYQGoP6dpAABCxklEQVR4nO3deXxU9b3/8deZycqSBQkJgbDLpggYIMYNrVEQW6G1CJQKIgVri13ibQtdpOrtjbfSXn9VW6vXpVUq6HUtWihGccEIGECQTUGWACYsIQkkkGRmvr8/BoKRLDMwZyYneT8fj/PAnPmeOZ85j8i8+Z7v+X4tY4xBRERExCFckS5AREREJBgKLyIiIuIoCi8iIiLiKAovIiIi4igKLyIiIuIoCi8iIiLiKAovIiIi4igKLyIiIuIoCi8iIiLiKAovIiIi4ihR4TjJI488wgMPPEBxcTFDhw7loYceYtSoUc0et2jRIqZMmcL48eN55ZVXAjqXz+dj//79dOzYEcuyzrFyERERCQdjDEePHiU9PR2Xq5m+FWOzRYsWmZiYGPPkk0+aTZs2mVmzZpmkpCRTUlLS5HE7d+403bp1M1dccYUZP358wOcrKioygDZt2rRp06bNgVtRUVGz3/WWMfYuzJiVlcXIkSN5+OGHAX/PSEZGBnfeeSdz585t8Biv18uVV17JbbfdxnvvvUdZWVnAPS/l5eUkJSVRVFREQkJCqD6GiIiI2KiiooKMjAzKyspITExssq2tt41qamooLCxk3rx5dftcLhc5OTkUFBQ0ety9995Lly5dmDlzJu+9916T56iurqa6urru56NHjwKQkJCg8CIiIuIwgQz5sHXA7qFDh/B6vaSmptbbn5qaSnFxcYPHvP/++zzxxBM8/vjjAZ0jLy+PxMTEui0jI+Oc6xYREZGWq0U9bXT06FFuueUWHn/8cTp37hzQMfPmzaO8vLxuKyoqsrlKERERiSRbbxt17twZt9tNSUlJvf0lJSWkpaWd0X7Hjh3s2rWLb3zjG3X7fD6fv9CoKLZt20bfvn3rHRMbG0tsbKwN1YuIiEhLZGvPS0xMDJmZmeTn59ft8/l85Ofnk52dfUb7gQMHsnHjRtavX1+33XjjjVx99dWsX79et4RERETE/nlecnNzmT59OiNGjGDUqFE8+OCDVFZWMmPGDACmTZtGt27dyMvLIy4ujgsvvLDe8UlJSQBn7BcREZG2yfbwMmnSJA4ePMjdd99NcXExw4YNY+nSpXWDePfs2dP8ZDQiIiIiJ9k+z0u4VVRUkJiYSHl5uR6VFhERcYhgvr/DsjyAiEhY+HywcwVsehmqSiGmA/S/DgZ+A6JiIl2diISIwouItA5ffAzPT4cjO8EVBT4PWG7YsAjanQfj/wwDxka6ShEJAQ02ERHnK94IT46Fsj3+n30e/5/G6/+zqhSemwxb34hMfSISUgovIuJ8r/wAPNWnw8oZTg7te/l2qD0etrJExB4KLyLibPsKoXhDE8HlFAPVFf7xMCLiaAovIuJsW9/wj3EJiAu2LrG1HBGxn8KLiDhbdQXQ/Cq0fj6oOmJnNSISBgovIuJscYnUjWlpjuWC9ufZWo6I2E/hRUScbdCNp58uao7xwaDx9tYjIrZTeBERZ+t6EXTL9M/p0hTLBfHJMPjG8NQlIrZReBER55vwKMS0bzzAWC7AgpuegKjYsJYmIqGn8CIizpfSH773JnQZ5P/ZFeUPMqeeQuqYBre8BP2uiVyNIhIyWh5ARFqHlAHw/fdh70ew6SWoOnxybaOx/tDiaua2kog4hsKLiLQelgUZI/2biLRaum0kIiIijqLwIiIiIo6i8CIiIiKOovAiIiIijqLwIiIiIo6i8CIiIiKOokelRSSsyqtq2X7wGD5j6NGpHakJcZEuSUQcRuFFRMLi05Kj/PWdHbz68X48Xv8q0BYwekAKs6/ow6X9Oke2QBFxDIUXEbHdim0HmP33QrzG4PWZuv0GeO+zQ6zYdpC7vz6Y2y7vHbkiRcQxNOZFRGy1/cAxbn+mkFqvr15wOeXUvnuXbObfm4rDXZ6IOJDCi4jY6smVO/H4DGfGlvpcFvwp/7Ow1CQizqbwIiK2qaz28GLh3gZ7XL7KZ+CT/RVs2l8ehspExMkUXkTENnuPHKfa4wvqmE9LjtpUjYi0FgovImIbn2m+x+WMY4LLOiLSBim8iIhtuiXHE+W2gjqmd0p7m6oRkdZC4UVEbJMQF82NF6XjdjUfYCygX5cODM9Isr0uEXE2hRcRsdXMKwKbu8UAc67uh2UF11MjIm2PwouI2OqC9ET+ePNQXBY0dAfpVKfMD6/uy4Th3cJbnIg4kmbYFRHbjR/WjfSkeB7K/4z3PjtUb86XQV0T+P7ovnxjaHrE6hMRZ1F4EZGwGNmrE3+fmUVRaRVbvqg4uTBjewanJ0S6NBFxGIUXEQmrjE7tyOjULtJliIiDacyLiIiIOIrCi4iIiDhKWMLLI488Qq9evYiLiyMrK4vVq1c32vall15ixIgRJCUl0b59e4YNG8YzzzwTjjJFRETEAWwPL4sXLyY3N5f58+ezdu1ahg4dypgxYzhw4ECD7Tt16sSvfvUrCgoK2LBhAzNmzGDGjBksW7bM7lJFRETEASxjzmLxkSBkZWUxcuRIHn74YQB8Ph8ZGRnceeedzJ07N6D3uPjii7nhhhu47777mm1bUVFBYmIi5eXlJCToKQYREREnCOb729ael5qaGgoLC8nJyTl9QpeLnJwcCgoKmj3eGEN+fj7btm3jyiuvbLBNdXU1FRUV9TYRERFpvWwNL4cOHcLr9ZKamlpvf2pqKsXFxY0eV15eTocOHYiJieGGG27goYce4tprr22wbV5eHomJiXVbRkZGSD+DiIiItCwt8mmjjh07sn79etasWcPvfvc7cnNzWbFiRYNt582bR3l5ed1WVFQU3mJFREQkrGydpK5z58643W5KSkrq7S8pKSEtLa3R41wuF/369QNg2LBhbNmyhby8PK666qoz2sbGxhIbGxvSukVERKTlsrXnJSYmhszMTPLz8+v2+Xw+8vPzyc7ODvh9fD4f1dXVdpQoIiIiDmP78gC5ublMnz6dESNGMGrUKB588EEqKyuZMWMGANOmTaNbt27k5eUB/jEsI0aMoG/fvlRXV/PGG2/wzDPP8Je//MXuUkVERMQBbA8vkyZN4uDBg9x9990UFxczbNgwli5dWjeId8+ePbhcpzuAKisr+cEPfsDevXuJj49n4MCBPPvss0yaNMnuUkVERMQBbJ/nJdw0z4uIiIjztJh5XkRERERCTeFFREREHEXhRURERBxF4UVEREQcReFFREREHEXhRURERBxF4UVEREQcReFFREREHEXhRURERBxF4UVEREQcReFFREREHEXhRURERBxF4UVEREQcReFFREREHEXhRURERBxF4UVEREQcReFFREREHEXhRURERBxF4UVEREQcReFFREREHEXhRURERBxF4UVEREQcReFFREREHEXhRURERBxF4UVEREQcReFFREREHEXhRURERBwlKtIFiIhIK+aphk+XQdkesFyQ0h/6XA0ud6QrEwdTeBERkdCrPQHvPgBr/hdOlIHlBgwYHySkQ/adkPV9cOkGgARP4UVEREKrphKe+SbsXeMPKwDGe/r1iv2wbB7s/Qhuely9MBI0RV4REQmtV35YP7g0ZtNLsCIvPDVJq6LwIiIioVP6OWx+ufngAoCBgkeg+pjtZUnrovAiIiKh89FTJ8e3BKi2Cja+YF890iopvIiISOgUrao/vqU5rijYV2hfPdIqKbyIiEjoeKqDa29M8MdIm6fwIiIiodMxzT+fS6AsCzp0sa8eaZUUXkREJHSGTAxwsO5JPg8M+bZ99UirFJbw8sgjj9CrVy/i4uLIyspi9erVjbZ9/PHHueKKK0hOTiY5OZmcnJwm24uISAsy6EaI7wRYzbe13NB1KKQPt70saV1sDy+LFy8mNzeX+fPns3btWoYOHcqYMWM4cOBAg+1XrFjBlClTePvttykoKCAjI4PrrruOffv22V2qiEho+IIYsNraRMXANx89+UMTAcZygzsabnwoLGVJ62IZY4ydJ8jKymLkyJE8/PDDAPh8PjIyMrjzzjuZO3dus8d7vV6Sk5N5+OGHmTZtWrPtKyoqSExMpLy8nISEhHOuX0SkWcbA7pWw+nH4dCl4TkBULPS7FkbNgt6j/WM72pItS+DF28Bb478+dVyAD+KT4TsvQMbISFUoLUww39+2Lg9QU1NDYWEh8+bNq9vncrnIycmhoKAgoPeoqqqitraWTp06Nfh6dXU11dWnR6pXVFScW9EiIsGoPQ4vzoKt//Q/9uvz+Pd7quHTf8HWJf4QM/FpiO0Q0VLPRXlVLS8UFvH6hi8oraqhfUwUV/ZPYWpWDzI6tTvzgEFfh9ytsH6hf+6X8iJ/b0vnfjByln+cS0z78H8QaRVsDS+HDh3C6/WSmppab39qaipbt24N6D1+8YtfkJ6eTk5OToOv5+Xlcc8995xzrSISRsb417Up/tj/310GQ89Lndc74fPBCzPgs2Unf/Z85fWTt4925MPi78LU/wO3PX/tVtV4eG39fl77eD8Hj1YTG+1mRM9kpmb14PzUjuf03s98uJv7/rmZWq+PL/ehbCuu4K/v7GDqJT347TcuIMr9lZEI7TrBpXf6N5EQatELM95///0sWrSIFStWEBcX12CbefPmkZubW/dzRUUFGRkZ4SpRRIK1/U3492/gwGZOj4kw0Kkv5MyHweMjWV1wtr3h711pjvHB52/DppfhookhL+PFwr3c/eonVNZ4sazTd2m2fFHB0x/s4msDu/Dg5GEkxEUH/d5PrdzJPf/c3OBr3pPnWfjhHo6e8PDgpGFYTgug4ki2Dtjt3LkzbrebkpKSevtLSkpIS0tr8tgFCxZw//338+9//5uLLrqo0XaxsbEkJCTU20SkhfrkRXj223Bgy8kd5uSGf02c56fBR09Gqrrgrf5r4FPhWy5/+xB79sPd3PXCx1TW+Ht5vjy8xOvz//DOtgNM/msBldWeht6iUUWlVdy7pOHg8mUGeHX9fpZ+UhzU+4ucLVvDS0xMDJmZmeTn59ft8/l85Ofnk52d3ehxv//977nvvvtYunQpI0aMsLNEEQmXoyXw0u0nf2joOYGT+17PhcM7wlXV2TtRDjvfDXwqfOPzr7R8tKT5tgHacfAYv3n1k2bbeQ1sLT7Kfy8N7Hb9Kf9YvSfgLwm3BU9/sCuo9xc5W7Y/Kp2bm8vjjz/O3/72N7Zs2cIdd9xBZWUlM2bMAGDatGn1BvT+93//N7/5zW948skn6dWrF8XFxRQXF3PsmFYdFXG0dX8/+UXf3AOOLmf0vhw/cpbHlYashGc/3I0rwNs0PgPPf1TE0RO1Ab//i4V7624NNcdrYNXOUkoqTgT8/iJny/bwMmnSJBYsWMDdd9/NsGHDWL9+PUuXLq0bxLtnzx6++OKLuvZ/+ctfqKmp4dvf/jZdu3at2xYsWGB3qSJip09eCmzmVeOFjf9nfz3nKvosn5SJbuDJnLNQ6/WxeE1R3a2hQFTX+liy4YvmG550uLIm6LoOHtU6RWK/sAzYnTNnDnPmzGnwtRUrVtT7edeuXfYXJCLhd6Is8LbV5baVETLtO0OnPlC6k+Z7k05KSIfE7iE5/ZHKGqpqgpsML8ptsetwZcDto91WUOEIIDZKq86I/fRbJiLhEd/wXE0Niku2r45QsSzI+n4Q7V0w6nZwBTjA1y5BZJHhGcm4gnh4KCEuih7nhaZnSaQpCi8iEh5DJga22rDlhotutr+eUBg6xd+T0twTR5YbOqTCxc3PEh6o5PYxxEcHF4Q8XtPwhHKNmH5pTwLteHFbFlOyehAbFeFwJm2CwouIhMfwW/xr2QSyYN+IGbaXExJxCTD9NeiY1niAsdz+W0zTXvVP2hYi0W4X387sjjuIrpFot4tvDE0PuH3OoFT6p3Zo9hwuC+Jj3EzL7hXwe4ucC4UXEQmP9ufBt58Gl6uRHhjLv034MyT3Cm9t56JTH7j9Xbj8p/71er4sLtE/u+zt70HKgJCf+pbsnvgC7BpxWxY3ZXYjMT7wieqi3C7+dtsouifFN3r7yO2yiIt28/SMkXRLig/4vUXOhe0LM4abFmYUaeF2rYQ3fwt7V9ffn3YRXDMfzm94KRBH8NTA/nVQXQGxHaHrMIhueHbwUHni/Z3c18xEcm6XRe/O7XnpB5ee1Sy7ZVU1/PXdz1n44W4qTpye6C7abTFhWDfuuKovfVKcu26TtAzBfH8rvIhIZJRshuINJ9c2GgTpwyJdkWM9++Fu7luymRqP/1H0U3+pu13+p4Uu7Xsej3znYpLbx5zTeao9XtbsPMKRqhrax7q5uEcySe3O7T1FTlF4UXgRkTam4kQtLxXu5dX1+zl4rJq4KDcX90xiWnYvLuyWGOnyRJql8KLwIiIi4ijBfH9rwK6IiIg4isKLiIiIOIrCi4iIiDiKwouIiIg4isKLiIiIOIrCi4iIiDiKwouIiIg4isKLiIiIOIrCi4iIiDiKwouIiIg4isKLiIiIOIrCi4iIiDiKwouIiIg4isKLiIiIOIrCi4iIiDiKwouIiIg4SlSkCxARZ/H5DCt3HOLdTw9yvNZLasc4JgzvRkandpEuTUTaCIUXEQnYml2l3PX8x+wprSLKZQFgDPxx+aeMvTCN//72RSTERUe4ShFp7XTbSEQCsmZXKd95/EP2HqkCwOMzeHwGrzEY4N+bSpjy2IdU1XgiW6iItHoKLyLSLJ/PkLt4PV6fwWcabuM1hi1fVPDEezvDW5yItDkKLyLSrPe3H6LoyPFGg8spPgN/L9iNx+sLT2Ei0iYpvIhIs9759GDdGJfmHDxWzaclx2yuSETaMoUXEWnW8Vqvre1FRIKh8CIizUrtGIdp5pZRvfYJsfYVIyJtnsKLiDRrwvB0fAGkF5cFI3ol0z1Zc76IiH0UXkSkWT3Pa8+1g1NxW02Pe/EZuP3KvmGqSkTaKoUXEQnIAxOHMiCtAw2N2z21765r+3Pt4NTwFiYibY7Ci4gEJDE+mhe+fyk/uuZ8zmsfU++1i3sk89gtmdx5zfkRqk5E2hLLmGCG4bV8FRUVJCYmUl5eTkJCQqTLEWmVar0+thUf5UStl9SEOK1rJCLnLJjv77D0vDzyyCP06tWLuLg4srKyWL16daNtN23axE033USvXr2wLIsHH3wwHCWKSBCi3S4u7JbIiF6dFFxEJOxsDy+LFy8mNzeX+fPns3btWoYOHcqYMWM4cOBAg+2rqqro06cP999/P2lpaXaXJyIiIg5je3j54x//yKxZs5gxYwaDBw/m0UcfpV27djz55JMNth85ciQPPPAAkydPJjZWc0WIiIhIfbaGl5qaGgoLC8nJyTl9QpeLnJwcCgoKQnKO6upqKioq6m0iIiLSetkaXg4dOoTX6yU1tf6jk6mpqRQXF4fkHHl5eSQmJtZtGRkZIXlfEZG2zOP1caSyhooTtbSy5zqkFYiKdAHnat68eeTm5tb9XFFRoQAjInIWjDGsKyrj2YLd/HPDfmq9/tDSqV0M372kB1OyetA1MT7CVYrYHF46d+6M2+2mpKSk3v6SkpKQDcaNjY3V2BgRkXN0otbLf7zwMUs2fIHbZeH1ne5tKa2q4eG3t/PI2zu4Z/wFfPeSnhGsVMTm20YxMTFkZmaSn59ft8/n85Gfn092dradpxYRkQB5vD7uWFjIGxu/AKgXXE7xGfAaw69f+YRnPtwd7hJF6rH9tlFubi7Tp09nxIgRjBo1igcffJDKykpmzJgBwLRp0+jWrRt5eXmAf5Dv5s2b6/573759rF+/ng4dOtCvXz+7yxURaXOeW72HFVsPEujIlvmvfsJV/VPOnOPHUw2uaHBp8naxl+3hZdKkSRw8eJC7776b4uJihg0bxtKlS+sG8e7ZswfXl37R9+/fz/Dhw+t+XrBgAQsWLGD06NGsWLHC7nJFRNoUYwxPrtwV1DEWsHDVHuZePxAObYePnoB1C6G63P9q+jAYdTtc8E2IjrOhamnrtDyAiEgbtmZXKRMfDX7qioS4KNZdvxv30p8DLjDe0y9aLjA+6NwfbnkFEruFrF5pvVrc8gAiItIybT9w7KyOu7bmLdz/+hkYUz+4gD+4AJTugL99HU6Un2OVIvUpvIiItGG1Xh9WkMfEUMtvop9pfoyMzwtHdsGaJ86uOJFGKLyIiLRh57WPDXig7inXu1aRZFUGFnqMD1Y/5g8yIiGi8CIi0oaNHpBCfLQ7uGPcG/ESxDFHv4DD24OsTKRxCi8iIm1Yh9goJo7ojtsV+M2jeE7gwhfciWrObmyNSEMUXkRE2rjbR/elQ2wUgeaXMtMBb7BfH/HJwRcm0giFFxGRNq5bUjwLv5dFu5jAbgUt9Y0iikDHsLigy2BI7n32BYp8hcKLiIjQPjaKQGf9etc3hL2+znhNIF01Phg1G6xgn2kSaZzCi4iI8D/LP+WEJ7BxLAYXd9XegcFqOsBYLuh1BQz/boiqFPFTeBERaeMOHavm9Y1fNLggY2NWmUFMr/0Fx4j3z1P35QenXSdXnhlwPXxnMbijQ1yxtHW2r20kIiIt2+I1RZzNSjErfUPIqn6EG6M+5I7ED+kdfQSiYqHHpTByJnQdakO1IgovIiJt3uqdpQTR6VLPCWJ53jOaj6xxvPWTq0Jal0hjFF5ERNq40sqaJl934aMDx3HjpYo4qok5o03FiVq7yhM5g8KLiEgbFxt95vDHC63P+aZ7JcNdnzHY2k2cdTqc7DfnsdbXj/d9Q3jNeylVxBHj1hBKCR+FFxGRNq57UjzrrDK8xnCVax13Rb3AENcuao2LKHxnPOWcbh2mi+sIN7hW8ZuoZ3jO+zVWdLwtMsVLm6SoLCLSxn3z4u60N8dYEP0Xno55gMHWbgCirTODyylRJ19rb1Uzw72Mx4/eAdvfDGPV0pYpvIiItHFXpJzgX/G/5puulQC4reBG77otH3E1pfDsTVDwiB0litSj20YiIm3ZsQO4/jaOrhzCZQW52OKXWObksct+Ca5oyJodogJbBmMMa/cc4R+r9rCt5BhRLouRvZL5TlZPenduH+ny2hzLnM3D/S1YRUUFiYmJlJeXk5CQEOlyRERapL1Hqni5cC9j1t9B32PrcAe8VlEgLJiVD90yQ/iekVNxopYfPLuW97cfwu2y6ibzc1sWXmO47bLe/OqGQUGtzC1nCub7Wz0vIiJtyOb9FTzw762s2HaQm10r6B/9UehPYrngpdnw/ZUQHRf69w+jWq+PmU+voXD3EYB6sxB7T/7b/6mVO3FZ8OuvD45IjW2RwouISFtQfYxP33qajwreY7QPMlyp3BX1Aj4DIe8wMF44vAPW/g2ybg/xm9vkiw2w7yPweeG8ftB7NLhcLNtUzJpdR5o81ABPvL+Tadm96HFeu/DU28YpvIiItGbGwHt/wPfuHzi/topelhvcEI3X/oWeV/215a8ovbsAls2D/etO7rAAA4nd4ap5PL2qDy6LZmcgtiz4x+o9zL1+oN0VC3raSESk9TIG/vVzeOs+XJ4qLAtiLC8xVhiCCwZKd0DRKrtPdPY+Ww5/+zp88fGXdp5MKeV74dUfcsW+JwNaOsFnYN2epntoJHQUXkREWqvP34bVj0Xu/JYb9hRE7vxNOVEOz0/33yYyjT9l9WP3C2Ra2wJ6y2BW5ZZzo/AiItJarXoMXO7I1rB/fWTP35iPF0FtFXU9LY3w4ObWqGXNvp3bZTEgrWOIipPmKLyIiLRGnmr4bKm/ZyFSjBcOBtZrEXYb/y+gZlF4ud61GhdNz4Hj9RmmjOoRisokAAovIiKtUfVR/5iXSPOciHQFDTt+mOZ6XU6Jsnx0dFU3+rrLgmsHp3Jht8QQFSfNUXgREWmNYjrgf3ImwqJiI11Bw+KTA25qLDedEpOA+g9OnZqU7rJ+nfl/k4eFsDhpjh6VFhFpjaLjoF8O7HjLf/smEiw3dO4fmXM354Jvwd5Cmu19cUVhDRjHG9+8mlfX7+OZD3fz+cFK3C4Y3iOZ6dm9uHpgF82uG2YKLyIirVXW7bB9eWRrSB8W2fM3Zth3IP9e/9igpgKMzwOjZhMf42byqB5M1riWFkHhJQC7D1eyeE0ROw9V4nZZDO2exLczu5PcPibSpYmINK5fDlw83T/TbSQYL2RkRebczYlPhpuegOdvOZldGgkwl+dC7yvCWZkEQAszNuF4jZefv/gx//z4C9wuC5/PYFn+X3G3y+JHXzufO7/WD6slzx4pIm2bzwcr/gtW/gnjrcFj/EMdo2yeYddgURHXlX+OfoORvTu33MeId7wNS+fBwS3+NZmw/KGrfRcY/XMY+b2WPUNwKxLM97fCSyNqvT6mP7maDz8/3OTsindc1ZdfjNV00CLSwh0vgw3P8/mm1azbU0qy9xBXuz627XvZZyz+y/tdnvBcjwEyeyYz/dJejL0gjZioFvasiDGwdw3s/ch/m6jz+dDvWnDr5kQ4KbyEILw8v6aIn7+4IaC2//7plfRPbaH/qhAR+Yoaj49/b9rPwNdvol/N1pC/v8e42Gc6c13N76nGf3v91PpAKR1jeWL6CC7qnhTy84qzBfP93cLib8tgjOGplTsD+heJ22Xx7Ie77S9KRCREYqJcfH1od/r94gNI7hPy93fhI7f2jrrgAqcXNiw9Vs3ERwt4/7NDIT+vtB0KLw0oq6plS/HRgOZ38voMb24psb8oEZFQc7thxuvQsevJ8R6hkef5DoVmQIOveQ3UeH187+9r2LS/PGTnlLZF4aUBVbXBzYlwoqbpaaNFRFqshHS4ban/T+vs10HyGn9X9e9rJ/G49+tNtjXGf+vqx4vW08pGLkiYKLw0ILldNO4AR7FZQJeEFjqDpIhIIJJ7wex3YdA3Tu4IbhSvx7gopz2zanL5s3d8QMf4DGw/cIyfv7iBkooWuoSAtFhhCS+PPPIIvXr1Ii4ujqysLFavXt1k+xdeeIGBAwcSFxfHkCFDeOONN8JRZp12MVGMvTAt4BkTJ47IsLkiERGbtT8Pbv4b3Pz307Piuhp/2sZnLHzG4oSJZpH3ar5W/QeW+0YEfdoXPtpLdl4+c/6xluJyhRgJjO3hZfHixeTm5jJ//nzWrl3L0KFDGTNmDAcOHGiw/QcffMCUKVOYOXMm69atY8KECUyYMIFPPvnE7lLrue3y3nibekYa/+j5djFuvn1x9zBVJSJis8Hj4YerYMZSGH4LpF4Arvq3kw6bjqzwDeVezy2Mqv4zv/bMpIyzf+LSZ+BfnxRzzR9X8H+Fe3UrSZpl+6PSWVlZjBw5kocffhgAn89HRkYGd955J3Pnzj2j/aRJk6isrGTJkiV1+y655BKGDRvGo48+ekb76upqqqtPr/ZZUVFBRkZGSCape2rlTu755+a6R/y+zO2yiHJZPDVjJJf27dz0G/l8UFsF0e3ApTt1IuIwnmr/PDE+D69tLuNHr+6y5TQW/klAcwZ14U9ThtMuRvOstCUt5lHpmpoaCgsLycnJOX1Cl4ucnBwKCgoaPKagoKBee4AxY8Y02j4vL4/ExMS6LSMjdLdwZlzWm8enjWBw1/oX0QJG90/hxTsubTq4FK2BF78H/9kF8rrBf6bA4mmwa2XIahQRsV1ULHRMhcRufFETZ9tpTv0b8a2tB5j6+CqOnqi17VzibLbG2kOHDuH1eklNTa23PzU1la1bG54Yqbi4uMH2xcXFDbafN28eubm5dT+f6nkJlWsHp3Lt4FQ2769g12H/2kYXpCfQPbld0we+9wf/ol+uKP+MjeD/c9vrsOVV/3oZ19ytaadFxDGqPV4WrS6y/Tw+Ax/vLeO2p9fwzMws4qLP/ikoaZ0c3ycXGxtLbKz9T/sMTk9gcHqAt6HWP+cPLnA6uJxy6uf3/+h/NHHUrNAVKSJioweWbmPX4cqwnMtnoHD3Ef7rjS3cO/7CsJxTnMPW8NK5c2fcbjclJfUncSspKSEtLa3BY9LS0oJq3+L4fPD27wJru+J+yLwV3NG2liQicq4+2lXKE+/v/Mray4ZUjjDEtZP+1l7aW8dxYagmmiJfFzaa3uww6Xg5u54Tn4G/F+zm+gu7kt33vFB8DGklbA0vMTExZGZmkp+fz4QJEwD/gN38/HzmzJnT4DHZ2dnk5+fzk5/8pG7f8uXLyc7OtrPU0Nn1LpQH2K1adQg+XQaDmp7QSUTkq4wxrN5ZSsHnh6n2+EjtGMvXh6bTuUPoe6JP1Hr56eL1WJZ/grmB1h6+617OOPcqOlnHAP9cL76T88NYQHSUf7LPahPFat9AnvFeS77v4qCDjMuC3OfX82buaNrHOv5mgYSI7b8Jubm5TJ8+nREjRjBq1CgefPBBKisrmTFjBgDTpk2jW7du5OXlAfDjH/+Y0aNH84c//IEbbriBRYsW8dFHH/HYY4/ZXWpoHPqM02Pmm+Fyw6FP7a5IRFqZD3Yc4jevfMKOg/5xeC7AYwz/+foWvjm8G7+98YKQftH/84P1nF/+Abe6N5LjKqSn6yAeYxFlnf57LspqeKbxWMtDtmszV7g/4YBJ5DHP13nKOzbgEOMzUFxxgkVriph5ee+QfB5xPtvDy6RJkzh48CB33303xcXFDBs2jKVLl9YNyt2zZw+uLz0+fOmll/KPf/yDX//61/zyl7/k/PPP55VXXuHCCx1yz9MVxL8qTJDtRaTNe3vbAWY+vabun0den+HUgiYeY3hx7V62lRxl8exs4mPO4e+X2uPw8SLMqkeZeHArE0+usXhqco0vB5fmnAo2KZTzq6iFTHC/z09rf8hnJsA5sgz87YNd3HZZLyw95CCEYZ6XcAvmOXFb7F8Hj10VePvpS6D3FbaVIyKtx7FqD1m/e5OqGm+TfbsuC2Zd0Yd54wad3Yn2rIKXZ8OR3RjACqQnOQge48JgcY9nGs96rw34uIXfy+Kyfs3MqyWO1WLmeWmT0odD16EBrNBqQae+0OvysJQlIs73yrp9VDYTXMB/q+Ufq/ZwIshFZgFY9yw8OQbKigAT8uAC/p6YaMvLf0Y/xQ/crwR0jNtl8fK6fSGvRZxJ4cUO1/0n/nEvjXVvntw/5r80z4uIBOyfH+8PeMnEo9UeVm4/FNwJNr8Kr84BDJizCD5n4efRz3Ob+1/NtvP6DIW7j4ShInEChRc79L4SJj3rn5Xyq+HEcvknrrvpf2HA2MjUJyKOdPhYTVD9IEeqgpih9mgxvPpDAnrYIMR+HfUsw63Pmm2363AlVTWeZttJ66fwYpeB4yB3C+Tc67+NlNAdUofA134NuZthyLcjXaGIOExCfHDPWHQI5omjf/8GaqqCrCg0fFj8T/SfiaWmyXbGwJYvKsJUlbRkemjeTu06wWU/8m8iIudo7IVprCsqI5DHLGKjXFzWL8CJ3Y4dgE0vhu1W0VdFWT4yOMBPo17kfs+UJtuWVmq9I1F4ERFxjImZGTywbBu13qbTi9uyuCmzOx3jzpy9+1i1h1fX7+OTfeXUeg3pSfHc6l5KJ9PwPC3h4rYMt7qX8mfPjVTQvtF21Z7IBCxpWRReREQcIrl9DP/1zSH87P82NNrG7bJIT4zjP64bUG+/z2d46K3tPPrODo7XeolynXyOyECvqGWMd1u4IzDe5cti8PAt93s87W18PGC0W6MdRGNeREQcZeKIDP7f5GEkxPn/7el2Wbgs/58AI3sl8+IPLqVT+5i6Y4wx/OrljfzPm59y/OTj0x6f8U9wZwyDrd24iWzPi5/hVvcymho03DFO/+YW9byIiDjO+GHdGHNBGm9s/IKCHf61jdIS4/jWxd0YmHbm5F6vb/yC59Y0vuZaPNV2lhswlwW9rBJ6WAfYY1IbbDOogc8nbY/Ci4iIA8VFu/nWxd351sXNT7H/5Ps7cVn+yesaUklciKs7N0OsnQ2Gl7TEOJK/1KMkbZduG4mItGJFpVWs3VPWaHAB2Gx64TEt4+vAGLjIteOM/S4Lhmckhb8gaZFaxm+riIjYorjiRLNt1vv64moRY17883pmuzafsd9n/I+Ki4DCi4hIqxbI0zmve7MwLejrIN06fMa+pPhohRep03J+W0VEJOTO79KB+Gh3k22OkMA60zegye/CoSP1Z/p1WzD1kh7ERjX9OaTtUHgREWnF2sdGMXFE97pHqRvzT292i1knNtby1AUYC4iNdjMtu1dEa5KWReFFRKSVm3VFH+KiXTSVXz7x9Q5fQQHob/kf7TbAb2+8gNSElvVElESWwouISCuX0akdf79tFO1ionA30L3idllsNj0jPL9ufe2tE7hdFqP7pzAxs/nHwaVtUXgREWkDMnt2YtlPr2TmFb3p+KXVpi0LvjawC0/MGk0LuWsEQLTlo0vHWB749kVYLeV+lrQYmqRORKSN6JYUzy/HDSL32v7sKztOrddHascvTfxmuSO2svRXxcTGs2j2JXTR7SJpgMKLiEgbExftpm9KhzNfSOwGZXvCX1AD7r5lLF3Pa3x1aWnbdNtIRET8el4W6QoAMNHt6dprcKTLkBZM4UVERPy6Dot0BRjASh9Oi3luW1okhRcREfHrFfmeF4MLel8R6TKkhVN4ERERv7Qh0HU4RPC5I2MMJf0mRuz84gwKLyIiclrW7RChGV88xsWbvkwWbvZE5PziHAovIiJy2gXfhPYpYIX/6yHK8vG/nnE8++EevL6WNGWetDQKLyIiclp0HIx/BIwvrKf1GhfPea5ijRlIaVUNOw4eC+v5xVkUXkREpL7+Y+Ciyf5J68LAaywOkcDvPN+t27dhb3lYzi3OpPAiIiJnuv5+SMoAl70BxmcsfLj4ce0cjtEOgCiXxaclR209rzibwouIiJwpPhmmL4H2qbb1wPiMf2jwnNo7+dBXf1K6Y9UatCuNU3gREZGGJWXA95ZDck+8If668BgXHqKYXZvLMt+oM15vaPVrkVMUXkREpHGJ3eH293inww2Af3zKuTIGNpleXF+TR74vs8E2ye2iz/k80nopvIiISNNiO7B+6N1Mrf01+0xnwN9zEgyf8W/HTQx5nil8s+ZedphuDbb1+AwXdEs857Kl9dKq0iIi0qwh3RL5k3cwo73/w5WujUxzL+Nq13pcFtQaN268uL7SKVNr3EThxbJgp+nKU96xvOK9rG5gbnPnE2mMwouIiDQru+95xEW7OFEL7/iG8o5vKF04wnDXdi507WSotYOu1mFiqcWDmzI6sNHXm09Mbzb4+rDNZBDIsgOWBf27dKBrYpz9H0ocS+FFRESa1SE2iomZGfxj9enZbw+QzDLfSJb5RobsPMbAjMt6Y2nArjRBY15ERCQg0y/tiTH2TdtvAYnx0dw4LN22c0jrYFt4KS0tZerUqSQkJJCUlMTMmTM5dqzp6Z4fe+wxrrrqKhISErAsi7KyMrvKExGRIPXr0pEfXNUPuzpFDJD3rSG0i9FNAWmabeFl6tSpbNq0ieXLl7NkyRLeffddZs+e3eQxVVVVjB07ll/+8pd2lSUiIufgzmv60adze9xfHZ17jlwWjBuSxrghXUP6vtI6WcaGPsAtW7YwePBg1qxZw4gRIwBYunQp48aNY+/evaSnN90luGLFCq6++mqOHDlCUlJSk22rq6uprq6u+7miooKMjAzKy8tJSEg4588iIiL1FZVW8c0/r+RIVW1IVn92WXBht0Sem3UJ7WPV69JWVVRUkJiYGND3ty09LwUFBSQlJdUFF4CcnBxcLherVq0K6bny8vJITEys2zIyMkL6/iIiUl9Gp3a8dMdlpCbEnvF49NnI7JnMwu9lKbhIwGwJL8XFxXTp0qXevqioKDp16kRxcXFIzzVv3jzKy8vrtqKiopC+v4iInKnHee1Y+pMrmTjC/w/GYEOM22Xhdln8bMwAnpt1CR3jNKOuBC6o8DJ37lwsy2py27p1q121Nig2NpaEhIR6m4iI2C8hLpr/vukiFn4vi1G9OwH+UNJYjnFZp7exF6Txxo+u4IdX9yPKrQdfJThB9dHddddd3HrrrU226dOnD2lpaRw4cKDefo/HQ2lpKWlpaUEXKSIiLddl/TpzWb/ObD9wjNfW72NdURkb9pZTfry2rk1qQizDMpLI7JnMhOHd6NJRk9DJ2QsqvKSkpJCSktJsu+zsbMrKyigsLCQz07/o1ltvvYXP5yMrK+vsKhURkRatX5cO5F43AABjDNUeHx6fIcbtIiZKvSsSOrb8Ng0aNIixY8cya9YsVq9ezcqVK5kzZw6TJ0+ue9Jo3759DBw4kNWrV9cdV1xczPr169m+fTsAGzduZP369ZSWltpRpoiI2MSyLOKi3XSIjVJwkZCz7Tdq4cKFDBw4kGuuuYZx48Zx+eWX89hjj9W9Xltby7Zt26iqqqrb9+ijjzJ8+HBmzZoFwJVXXsnw4cN57bXX7CpTREREHMaWeV4iKZjnxEVERKRliPg8LyIiIiJ2UXgRERERR1F4EREREUdReBERERFHUXgRERERR1F4EREREUdReBERERFHUXgRERERR1F4EREREUdReBERERFHUXgRERERR1F4EREREUdReBERERFHUXgRERERR1F4EREREUdReBERERFHUXgRERERR1F4EREREUdReBERERFHUXgRERERR1F4EREREUdReBERERFHUXgRERERR1F4EREREUdReBERERFHUXgRERERR1F4EREREUdReBERERFHUXgRERERR1F4EREREUdReBERERFHUXgRERERR1F4EREREUdReBERERFHUXgRERERR1F4EREREUdReBERERFHsTW8lJaWMnXqVBISEkhKSmLmzJkcO3asyfZ33nknAwYMID4+nh49evCjH/2I8vJyO8sUERERB7E1vEydOpVNmzaxfPlylixZwrvvvsvs2bMbbb9//37279/PggUL+OSTT3j66adZunQpM2fOtLNMERERcRDLGGPseOMtW7YwePBg1qxZw4gRIwBYunQp48aNY+/evaSnpwf0Pi+88ALf/e53qaysJCoq6ozXq6urqa6urvu5oqKCjIwMysvLSUhICM2HEREREVtVVFSQmJgY0Pe3bT0vBQUFJCUl1QUXgJycHFwuF6tWrQr4fU59iIaCC0BeXh6JiYl1W0ZGxjnXLiIiIi2XbeGluLiYLl261NsXFRVFp06dKC4uDug9Dh06xH333dfkraZ58+ZRXl5etxUVFZ1T3SIiItKyBR1e5s6di2VZTW5bt24958IqKiq44YYbGDx4ML/97W8bbRcbG0tCQkK9TURERFqvhu/FNOGuu+7i1ltvbbJNnz59SEtL48CBA/X2ezweSktLSUtLa/L4o0ePMnbsWDp27MjLL79MdHR0sGWKiIhIKxV0eElJSSElJaXZdtnZ2ZSVlVFYWEhmZiYAb731Fj6fj6ysrEaPq6ioYMyYMcTGxvLaa68RFxcXbIkiIiLSitk25mXQoEGMHTuWWbNmsXr1alauXMmcOXOYPHly3ZNG+/btY+DAgaxevRrwB5frrruOyspKnnjiCSoqKiguLqa4uBiv12tXqSIiIuIgQfe8BGPhwoXMmTOHa665BpfLxU033cSf/vSnutdra2vZtm0bVVVVAKxdu7buSaR+/frVe6+dO3fSq1cvO8sVERERB7BtnpdICeY5cREREWkZWsQ8LyIiIiJ2UHgRERERR1F4EREREUdReBERET9PDfh8ka5CpFm2Pm0kIiIt3IGtsOZ/YcMiqD4KWJA2BEbNhgtvgph2ka5Q5AzqeRERaYuMgbfz4M9Z8NGTJ4MLgIGST+C1OfBQJhzcFtEyRRqi8CIi0ha983t4537/f5uvTAJqTt46OlYCT10PZXvCW5tIMxReRETamrI9sCKv+XbGCyfK4M17bC9JJBgKLyIibc1HT4EV4F//Pi9sfgWOHbS1JJFgKLyIiLQ1G58/81ZRU3we2Pa6ffWIBEnhRUSkrakqDa69yw1Vh+2pReQsKLyIiLQ1UbHBtTcGouLtqUXkLCi8iIi0NT0v8/emBMr4ICPLvnpEgqTwIiLS1oya5R+IGwjLBWkXQfdMe2sSCYLCi4hIW9N7NHQfAVYAvS/GB1fNs78mkSAovIiItDWWBVMWQ+fzG39k+lSwGXs/DBwXvtpEAqDwIiLSFrXvDDOXw+U/hbikM1/vkQXffREuuSPspYk0xzLGmEgXEUoVFRUkJiZSXl5OQkJCpMsREWn5PNXw+Qo4dsD/JFLXoZAyINJVSRsTzPe3VpUWEWnromKh/5hIVyESMN02EhEREUdReBERERFHUXgRERERR1F4EREREUdReBERERFHUXgRERERR1F4EREREUdReBERERFHUXgRERERR1F4EREREUdReBERERFHUXgRERERR1F4EREREUdReBERERFHUXgRERERR1F4EREREUdReBERERFHUXgRERERR7E1vJSWljJ16lQSEhJISkpi5syZHDt2rMljbr/9dvr27Ut8fDwpKSmMHz+erVu32lmmiIiIOIit4WXq1Kls2rSJ5cuXs2TJEt59911mz57d5DGZmZk89dRTbNmyhWXLlmGM4brrrsPr9dpZqoiIiDiEZYwxdrzxli1bGDx4MGvWrGHEiBEALF26lHHjxrF3717S09MDep8NGzYwdOhQtm/fTt++fc94vbq6murq6rqfKyoqyMjIoLy8nISEhNB8GBEREbFVRUUFiYmJAX1/29bzUlBQQFJSUl1wAcjJycHlcrFq1aqA3qOyspKnnnqK3r17k5GR0WCbvLw8EhMT67bG2omIiEjrYFt4KS4upkuXLvX2RUVF0alTJ4qLi5s89s9//jMdOnSgQ4cO/Otf/2L58uXExMQ02HbevHmUl5fXbUVFRSH7DCIiItLyBB1e5s6di2VZTW7nOsB26tSprFu3jnfeeYf+/ftz8803c+LEiQbbxsbGkpCQUG8TERGR1isq2APuuusubr311ibb9OnTh7S0NA4cOFBvv8fjobS0lLS0tCaPP3UL6Pzzz+eSSy4hOTmZl19+mSlTpgRbroiIiLQyQYeXlJQUUlJSmm2XnZ1NWVkZhYWFZGZmAvDWW2/h8/nIysoK+HzGGIwx9QblioiISNtl25iXQYMGMXbsWGbNmsXq1atZuXIlc+bMYfLkyXVPGu3bt4+BAweyevVqAD7//HPy8vIoLCxkz549fPDBB0ycOJH4+HjGjRtnV6kiIiLiILbO87Jw4UIGDhzINddcw7hx47j88st57LHH6l6vra1l27ZtVFVVARAXF8d7773HuHHj6NevH5MmTaJjx4588MEHZwz+FRERERvUnoAju6BsD3hrI11Ng2yb5yVSgnlOXERERE46vAM+/AusXwi1/k4F4pJgxG2QdTsHTBIvr9vH3iPHiXa7GNkrmZzBqUS7Q9MPEsz3t8KLiIhIW7fjbXhuMvg8/u1LjOWm0tWBm4/PY6vpgcuyAPD4DOe1j2H+jRdw49DAJp5tSouYpE5EREQc4PAOf3DxVJ8RXAAs4yXOc5S/R/8XHcwxPD6Dx+fv9zhcWcOPnlvH82vCO8eawouIiEhb9uGfwVcLNH4jJsry0YmjTHS/2+Dr817eyIGKhudjs4PCi4iISFtVexzWLQRfIIsfG6a5/93wK8awKIy9LwovIiIibVXFfvAcD6ipy4IM6wAWvjNe8xl4fcMXoa6u8VrCdiYRERFpWazgYoDBwmA1+Fr58fA9Vq3wIiIi0lYldIO4xICaeo3Fp6Y7NBJekttHh7Cwpim8iIiItFVRMZB5K1juZpu6MPzNO6bh1yz4Rggelw6UwouIiEhblvV9iEtoMsB4jIudJo1XvJc1+LrbZXHziAy7KjyDwouIiEhblpAO016D+KQzxsD4x7hAkenCd2t+yQli671undwWTBxK5w71X7OTwouIiEhb1/UimPMR5NwDiT1O7rSwugzEM+6PPDvsWQ66U7AsiHJZRLn8417Sk+J5bNoIxg/rFtZytTyAiIiI1Of1gGWB6/StpLKqGv758X72lh0nxu0is2cyV56fgsvV8ADeYAXz/R0VkjOKiIhI6+E+Mx4ktYvhluxe4a+lAbptJCIiIo6i8CIiIiKOovAiIiIijqLwIiIiIo6i8CIiIiKOovAiIiIijqLwIiIiIo6i8CIiIiKO0uomqTs1YXBFRUWEKxEREZFAnfreDmTi/1YXXo4ePQpARkb4VrcUERGR0Dh69CiJiYlNtml1axv5fD72799Px44dsazQrLfQllVUVJCRkUFRUZHWigoRXdPQ0zUNLV3P0NM1bZ4xhqNHj5Keno7L1fSollbX8+JyuejevXuky2h1EhIS9D9ciOmahp6uaWjpeoaermnTmutxOUUDdkVERMRRFF5ERETEURRepEmxsbHMnz+f2NjYSJfSauiahp6uaWjpeoaermlotboBuyIiItK6qedFREREHEXhRURERBxF4UVEREQcReFFREREHEXhRURERBxF4UXO8Lvf/Y5LL72Udu3akZSUFNAxt956K5Zl1dvGjh1rb6EOcTbX0xjD3XffTdeuXYmPjycnJ4fPPvvM3kIdpLS0lKlTp5KQkEBSUhIzZ87k2LFjTR5z1VVXnfE7+v3vfz9MFbc8jzzyCL169SIuLo6srCxWr17dZPsXXniBgQMHEhcXx5AhQ3jjjTfCVKlzBHNNn3766TN+H+Pi4sJYrbMpvMgZampqmDhxInfccUdQx40dO5YvvviibnvuuedsqtBZzuZ6/v73v+dPf/oTjz76KKtWraJ9+/aMGTOGEydO2Fipc0ydOpVNmzaxfPlylixZwrvvvsvs2bObPW7WrFn1fkd///vfh6Halmfx4sXk5uYyf/581q5dy9ChQxkzZgwHDhxosP0HH3zAlClTmDlzJuvWrWPChAlMmDCBTz75JMyVt1zBXlPwLxXw5d/H3bt3h7FihzMijXjqqadMYmJiQG2nT59uxo8fb2s9Thfo9fT5fCYtLc088MADdfvKyspMbGysee6552ys0Bk2b95sALNmzZq6ff/617+MZVlm3759jR43evRo8+Mf/zgMFbZ8o0aNMj/84Q/rfvZ6vSY9Pd3k5eU12P7mm282N9xwQ719WVlZ5vbbb7e1TicJ9poG8/ernEk9LxIyK1asoEuXLgwYMIA77riDw4cPR7okR9q5cyfFxcXk5OTU7UtMTCQrK4uCgoIIVtYyFBQUkJSUxIgRI+r25eTk4HK5WLVqVZPHLly4kM6dO3PhhRcyb948qqqq7C63xampqaGwsLDe75fL5SInJ6fR36+CgoJ67QHGjBmj38eTzuaaAhw7doyePXuSkZHB+PHj2bRpUzjKbRVa3arSEhljx47lW9/6Fr1792bHjh388pe/5Prrr6egoAC32x3p8hyluLgYgNTU1Hr7U1NT615ry4qLi+nSpUu9fVFRUXTq1KnJ6/Od73yHnj17kp6ezoYNG/jFL37Btm3beOmll+wuuUU5dOgQXq+3wd+vrVu3NnhMcXGxfh+bcDbXdMCAATz55JNcdNFFlJeXs2DBAi699FI2bdpE9+7dw1G2o6nnpY2YO3fuGYPDvro19j9ZICZPnsyNN97IkCFDmDBhAkuWLGHNmjWsWLEidB+iBbH7erZFdl/T2bNnM2bMGIYMGcLUqVP5+9//zssvv8yOHTtC+ClEApOdnc20adMYNmwYo0eP5qWXXiIlJYW//vWvkS7NEdTz0kbcdddd3HrrrU226dOnT8jO16dPHzp37sz27du55pprQva+LYWd1zMtLQ2AkpISunbtWre/pKSEYcOGndV7OkGg1zQtLe2MQZAej4fS0tK6axeIrKwsALZv307fvn2DrtepOnfujNvtpqSkpN7+kpKSRq9fWlpaUO3bmrO5pl8VHR3N8OHD2b59ux0ltjoKL21ESkoKKSkpYTvf3r17OXz4cL0v39bEzuvZu3dv0tLSyM/PrwsrFRUVrFq1KugnwJwk0GuanZ1NWVkZhYWFZGZmAvDWW2/h8/nqAkkg1q9fD9Bqf0cbExMTQ2ZmJvn5+UyYMAEAn89Hfn4+c+bMafCY7Oxs8vPz+clPflK3b/ny5WRnZ4eh4pbvbK7pV3m9XjZu3Mi4ceNsrLQVifSIYWl5du/ebdatW2fuuece06FDB7Nu3Tqzbt06c/To0bo2AwYMMC+99JIxxpijR4+a//iP/zAFBQVm586d5s033zQXX3yxOf/8882JEyci9TFajGCvpzHG3H///SYpKcm8+uqrZsOGDWb8+PGmd+/e5vjx45H4CC3O2LFjzfDhw82qVavM+++/b84//3wzZcqUutf37t1rBgwYYFatWmWMMWb79u3m3nvvNR999JHZuXOnefXVV02fPn3MlVdeGamPEFGLFi0ysbGx5umnnzabN282s2fPNklJSaa4uNgYY8wtt9xi5s6dW9d+5cqVJioqyixYsMBs2bLFzJ8/30RHR5uNGzdG6iO0OMFe03vuuccsW7bM7NixwxQWFprJkyebuLg4s2nTpkh9BEdReJEzTJ8+3QBnbG+//XZdG8A89dRTxhhjqqqqzHXXXWdSUlJMdHS06dmzp5k1a1bd/7RtXbDX0xj/49K/+c1vTGpqqomNjTXXXHON2bZtW/iLb6EOHz5spkyZYjp06GASEhLMjBkz6oXBnTt31rvGe/bsMVdeeaXp1KmTiY2NNf369TM/+9nPTHl5eYQ+QeQ99NBDpkePHiYmJsaMGjXKfPjhh3WvjR492kyfPr1e++eff97079/fxMTEmAsuuMC8/vrrYa645Qvmmv7kJz+pa5uammrGjRtn1q5dG4GqnckyxpiIdPmIiIiInAU9bSQiIiKOovAiIiIijqLwIiIiIo6i8CIiIiKOovAiIiIijqLwIiIiIo6i8CIiIiKOovAiIiIijqLwIiIiIo6i8CIiIiKOovAiIiIijvL/Aat1dHraVNSkAAAAAElFTkSuQmCC", | |
| "text/plain": [ | |
| "<Figure size 640x480 with 1 Axes>" | |
| ] | |
| }, | |
| "metadata": {}, | |
| "output_type": "display_data" | |
| } | |
| ], | |
| "source": [ | |
| "x_samples_unflatten = rearrange(x_samples, '(n c) -> n c', c=n_feat)\n", | |
| "\n", | |
| "# Plot samples, first two is xy and last is s\n", | |
| "plt.scatter(x_samples_unflatten[:, 0], x_samples_unflatten[:, 1], s=100 * np.exp(x_samples_unflatten[:, 2]))\n", | |
| "plt.scatter(x[ii, :, 0], x[ii, :, 1], s=100 * np.exp(x[ii, :, 2]))" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": null, | |
| "metadata": {}, | |
| "outputs": [], | |
| "source": [] | |
| } | |
| ], | |
| "metadata": { | |
| "kernelspec": { | |
| "display_name": "Python 3", | |
| "language": "python", | |
| "name": "python3" | |
| }, | |
| "language_info": { | |
| "codemirror_mode": { | |
| "name": "ipython", | |
| "version": 3 | |
| }, | |
| "file_extension": ".py", | |
| "mimetype": "text/x-python", | |
| "name": "python", | |
| "nbconvert_exporter": "python", | |
| "pygments_lexer": "ipython3", | |
| "version": "3.9.13" | |
| }, | |
| "orig_nbformat": 4 | |
| }, | |
| "nbformat": 4, | |
| "nbformat_minor": 2 | |
| } |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment