From 8f187cd4fd777ea6d6bc252da034dccb86585745 Mon Sep 17 00:00:00 2001 From: Valentin Pratz Date: Fri, 11 Oct 2024 16:13:18 +0000 Subject: [PATCH] Documentation: Consistency Model example notebook --- examples/ConsistencyModel.ipynb | 681 ++++++++++++++++++++++++++++++++ 1 file changed, 681 insertions(+) create mode 100644 examples/ConsistencyModel.ipynb diff --git a/examples/ConsistencyModel.ipynb b/examples/ConsistencyModel.ipynb new file mode 100644 index 000000000..00bdba4b4 --- /dev/null +++ b/examples/ConsistencyModel.ipynb @@ -0,0 +1,681 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "009b6adf", + "metadata": {}, + "source": [ + "# Consistency Models for Posterior Estimation" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "id": "d5f88a59", + "metadata": { + "ExecuteTime": { + "end_time": "2024-09-23T14:39:46.551814Z", + "start_time": "2024-09-23T14:39:46.032170Z" + } + }, + "outputs": [], + "source": [ + "import matplotlib.pyplot as plt\n", + "import numpy as np\n", + "import seaborn as sns\n", + "\n", + "# ensure the backend is set\n", + "import os\n", + "if \"KERAS_BACKEND\" not in os.environ:\n", + " # set this to \"torch\", \"tensorflow\", or \"jax\"\n", + " os.environ[\"KERAS_BACKEND\"] = \"jax\"\n", + "\n", + "import keras\n", + "\n", + "# for BayesFlow devs: this ensures that the latest dev version can be found\n", + "import sys\n", + "sys.path.append('../')\n", + "\n", + "import bayesflow as bf" + ] + }, + { + "cell_type": "markdown", + "id": "eadaf793-ab63-4f69-b962-178e343ca21b", + "metadata": {}, + "source": [ + "In this notebook, we use Consistency Models (CMs) as a plug-in replacement to obtain posterior samples with fewer sampling steps.\n", + "\n", + "CMs can be trained in two ways: First, they can be used to _distill_ an existing score-based diffusion model, thereby massively decreasing the sampling time at the expense of an additional training phase. Second, they can be trained from scratch using a procedure named _Consistency Training_. For now, we only support the latter." + ] + }, + { + "cell_type": "markdown", + "id": "6286c800-460a-4881-87d8-c3aca7aeec70", + "metadata": {}, + "source": [ + "## Background" + ] + }, + { + "cell_type": "markdown", + "id": "fdff817a-6321-4af0-9d41-7ec80097f93b", + "metadata": {}, + "source": [ + "Consistency Models [1] leverage some nice properties of score-based diffusion to enable few-step sampling. Score-based diffusion initially relied on a stochastic differential equation (SDE) for sampling, but there is also a ordinary (non-stochastic) differential equation (ODE) has the same _marginal_ distribution at each time step $t$ [2]. This means that even though SDE and ODE produce different paths from the noise distribution to the target distribution, the resulting distributions when looking at many paths at time $t$ is the same. The ODE is also called Probability Flow ODE." + ] + }, + { + "cell_type": "markdown", + "id": "4a2e996d-355a-4fab-8347-728e563c6014", + "metadata": {}, + "source": [ + "CMs now leverage the fact that there is no randomness in the ODE formulation. That means, if you start at a certain point in the latent space, you will always take the same path and always end up at the same point in the data space. The same is true for every point on the path: if you integrate to get to time $t=0$, you will end up at the same point as well. In short: for each path, there is exactly one corresponding point in latent space (at $t=T$) and one corresponding point in data space (at $t=0$). The goal of CMs is now the following: each point at a time point $t$ belongs to exactly one path, and we want to predict where this path will end up at $t=0$. The function that does this is called the _consistency function_ $f$. If we have the correct function for all $t\\in(0,T]$, we can just sample from the latent distribution ($t=T$) and use $f$ to directly map to the corresponding point at $t=0$, which is in the target distribution. So for sampling from the target distribution, we avoid any integration and only need one evaluation of the consistency function. In practice, the one-step sampling does not work very well. Instead, we leverage a multi-step sampling method where we call $f$ multiple times. Please check out the [1] for more background on this sampling procedure." + ] + }, + { + "cell_type": "markdown", + "id": "25023294-3096-4ebc-83a6-a372208e0504", + "metadata": {}, + "source": [ + "When only reading the above you might wonder why we also learn the mapping to $t=0$ of all intermediate time steps, and not only for $T=0$. The main answer is that for efficient training, we do not want to actually compute the two associated points explicitly. Doing so would require to do a precise integration at training time, which is often not feasible as it is too computationally costly. Learning all time steps opens up the possibility for a different training approach where we can avoid this." + ] + }, + { + "cell_type": "markdown", + "id": "9e6eb121-f89f-4268-9a0d-6fa733c758ff", + "metadata": {}, + "source": [ + "The details of this become a bit more complicated, and we advise you to take a look at [1] if you are interested in a more thorough and mathematical discussion. Here we will give a rough description of the underlying concepts." + ] + }, + { + "cell_type": "markdown", + "id": "9b201899-b946-4432-bcac-106b9d580d32", + "metadata": {}, + "source": [ + "First, we know that at $t=0$, it holds that $f(x,t=0)=x$, as $x$ is part of the path that ends at $x$. This _boundary condition_ serves as an \"anchor\" for our training, this is the information that the network knows at the start of the training procedure (we encode it with a time-dependent skip-connection, so the network is forced to be the identity function at $t=0$)." + ] + }, + { + "cell_type": "markdown", + "id": "a11f7ac7-9c11-49c1-b29c-3f0d44c4bf49", + "metadata": {}, + "source": [ + "For training, we now somehow have to propagate this information to the rest of the part. The basic idea for this is simple. We just take a point $x_1$ closer to the data distribution (smaller time $t_1$) and intergrate for a small time step $dt$ to a point $x_2$ on the same path that is closer to the latent distribution (larger time $t_2=t_1+dt$). As we know that for $t=0$ our network provides the correct output for our path, we want to propagate the information from smaller times to larger times. Our training goal is to move the output of $f(x_2, t=t_2)$ towards the output of $f(x_1, t=t_1)$. How to choose $x_1$, $t_1$ and $dt$ is an empirical question, see the [1] for some reasoning on what works." + ] + }, + { + "cell_type": "markdown", + "id": "73ff5ed5-dbd0-423b-b4f4-8d209947ed87", + "metadata": {}, + "source": [ + "In the case of _distillation_, we start with a trained score-based diffusion model. We can use it to integrate the Probability Flow ODE to get from $x_1$ to $x_2$. If we do not have such a model, it seems as if we were stuck. We do not know which points lie on the same path, so we do not know which outputs to make similar. Fortunately, it turns out that there is an _unbiased approximator_ that, if averaged over many samples (check out the paper for the exact description), will also give us the correct score. If we use this approximator instead of the score model, and use only a single Euler step to move along the path, we get an algorithm similar to the one described for distillation. It is called Consistency Training (CT) and allows us to train a consistency model using only samples from the data distribution. The algorithm for this was improved a lot in [3], and we have incorporated those improvements into our implementation." + ] + }, + { + "cell_type": "markdown", + "id": "283e8fea-9a36-4f8a-80f5-19e2fe98bd09", + "metadata": {}, + "source": [ + "As we have made several approximations to get there, how we have to choose our hyperparameters unfortunately becomes quite unintuitive. We have to rely on empirical observations to see what works. This was done in [4], we encourage you to use the values provided there as starting points. If you happen to find hyperparameters that work significantly better, please let us know. This will help others to find the correct region in the hyperparameter space." + ] + }, + { + "cell_type": "markdown", + "id": "8da3535f-0354-40a4-991f-845c33ff75b8", + "metadata": {}, + "source": [ + "To make this work for simulation-based inference, we can make the whole process conditional, so we can produce conditional distributions as well. Below, you can see a conceptual visualization of posterior estimation with CMs." + ] + }, + { + "cell_type": "markdown", + "id": "6c105be4-c490-4100-9502-60dd7405cfb4", + "metadata": {}, + "source": [ + "![Visualization of the way consistency models map from the path to the end point in the data distribution. Depicts the concepts described in the main text.](https://arxiv.org/html/2312.05440v2/extracted/5435837/figures/cmpe_main.png)" + ] + }, + { + "cell_type": "markdown", + "id": "baafb8fd-5b14-4ddf-a6dd-8272f706deaa", + "metadata": {}, + "source": [ + "### References\n", + "\n", + "[1] Song, Y., Dhariwal, P., Chen, M., & Sutskever, I. (2023). Consistency Models. *arXiv preprint*. [https://doi.org/10.48550/arXiv.2303.01469](https://doi.org/10.48550/arXiv.2303.01469)\n", + "\n", + "[2] Song, Y., Sohl-Dickstein, J., Kingma, D. P., Kumar, A., Ermon, S., & Poole, B. (2021). Score-Based Generative Modeling through Stochastic Differential Equations. In _International Conference on Learning Representations_. [https://openreview.net/forum?id=PxTIG12RRHS](https://openreview.net/forum?id=PxTIG12RRHS)\n", + "\n", + "[3] Song, Y., & Dhariwal, P. (2023). Improved Techniques for Training Consistency Models. *arXiv preprint*. [https://doi.org/10.48550/arXiv.2310.14189](https://doi.org/10.48550/arXiv.2310.14189)\n", + "\n", + "[4] Schmitt, M., Pratz, V., Köthe, U., Bürkner, P.-C., & Radev, S. T. (2024). Consistency Models for Scalable and Fast Simulation-Based Inference. *arXiv preprint*. [https://doi.org/10.48550/arXiv.2312.05440](https://doi.org/10.48550/arXiv.2312.05440)" + ] + }, + { + "cell_type": "markdown", + "id": "c63b26ba", + "metadata": {}, + "source": [ + "## Simulator: Two Moons" + ] + }, + { + "cell_type": "markdown", + "id": "9525ffd7", + "metadata": {}, + "source": [ + "We will use the Concistency Model as a plug-in replacement for Flow Matching. Refer to the tutorial \"Two moons toy example with flow matching\" for more details on this simulator." + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "4b89c861527c13b8", + "metadata": { + "ExecuteTime": { + "end_time": "2024-09-23T14:39:46.747091Z", + "start_time": "2024-09-23T14:39:46.744830Z" + } + }, + "outputs": [], + "source": [ + "simulator = bf.benchmarks.simulators.TwoMoons()" + ] + }, + { + "cell_type": "markdown", + "id": "f6e1eb5777c59eba", + "metadata": {}, + "source": [ + "We generate some data to see what the simulator does:" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "id": "e6218e61d529e357", + "metadata": { + "ExecuteTime": { + "end_time": "2024-09-23T14:39:46.798575Z", + "start_time": "2024-09-23T14:39:46.790581Z" + } + }, + "outputs": [], + "source": [ + "# generate 64 random draws from the joint distribution p(r, alpha, theta, x)\n", + "sample_data = simulator.sample((64,))" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "id": "46174ccb0167026c", + "metadata": { + "ExecuteTime": { + "end_time": "2024-09-23T14:39:46.854911Z", + "start_time": "2024-09-23T14:39:46.852129Z" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Type of sample_data:\n", + "\t \n", + "Keys of sample_data:\n", + "\t dict_keys(['parameters', 'observables'])\n", + "Types of sample_data values:\n", + "\t {'parameters': , 'observables': }\n", + "Shapes of sample_data values:\n", + "\t {'parameters': (64, 2), 'observables': (64, 2)}\n" + ] + } + ], + "source": [ + "print(\"Type of sample_data:\\n\\t\", type(sample_data))\n", + "print(\"Keys of sample_data:\\n\\t\", sample_data.keys())\n", + "print(\"Types of sample_data values:\\n\\t\", {k: type(v) for k, v in sample_data.items()})\n", + "print(\"Shapes of sample_data values:\\n\\t\", {k: v.shape for k, v in sample_data.items()})" + ] + }, + { + "cell_type": "markdown", + "id": "fee88fcfd7a373b0", + "metadata": {}, + "source": [ + "## Data Adapter\n", + "\n", + "The next step is to tell BayesFlow how to deal with the simulated variables. You may also think of this as informing BayesFlow about the data flow, i.e., which variables go into which network.\n", + "\n", + "For this example, we want to learn the posterior distribution $p(\\theta | x)$, so we **infer** $\\theta$, **conditioning** on $x$. In the output from the last command, we see that the simulator provides $\\theta$ as `\"parameters\"` and $x$ as `\"observables\"`." + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "id": "c9637c576d4ad4e5", + "metadata": { + "ExecuteTime": { + "end_time": "2024-09-23T14:39:46.905081Z", + "start_time": "2024-09-23T14:39:46.903091Z" + } + }, + "outputs": [], + "source": [ + "data_adapter = bf.ContinuousApproximator.build_data_adapter(\n", + " inference_variables=[\"parameters\"],\n", + " inference_conditions=[\"observables\"],\n", + ")" + ] + }, + { + "cell_type": "markdown", + "id": "254e287b2bccdad", + "metadata": {}, + "source": [ + "## Dataset\n", + "For this example, we will sample our training data ahead of time and use offline training with a `bf.datasets.OfflineDataset`." + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "id": "39cb5a1c9824246f", + "metadata": { + "ExecuteTime": { + "end_time": "2024-09-23T14:39:46.950573Z", + "start_time": "2024-09-23T14:39:46.948624Z" + } + }, + "outputs": [], + "source": [ + "batch_size = 64\n", + "num_training_batches = 512\n", + "num_validation_batches = 128" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "id": "9dee7252ef99affa", + "metadata": { + "ExecuteTime": { + "end_time": "2024-09-23T14:39:53.268860Z", + "start_time": "2024-09-23T14:39:46.994697Z" + } + }, + "outputs": [], + "source": [ + "training_samples = simulator.sample((num_training_batches * batch_size,))\n", + "validation_samples = simulator.sample((num_validation_batches * batch_size,))" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "id": "51045bbed88cb5c2", + "metadata": { + "ExecuteTime": { + "end_time": "2024-09-23T14:39:53.281170Z", + "start_time": "2024-09-23T14:39:53.275921Z" + } + }, + "outputs": [], + "source": [ + "training_dataset = bf.datasets.OfflineDataset(training_samples, batch_size=batch_size, data_adapter=data_adapter)\n", + "validation_dataset = bf.datasets.OfflineDataset(validation_samples, batch_size=batch_size, data_adapter=data_adapter)" + ] + }, + { + "cell_type": "markdown", + "id": "2d4c6eb0", + "metadata": {}, + "source": [ + "## Training a neural network to approximate all posteriors\n", + "\n", + "The next step is to set up the neural network that will approximate the posterior $p(\\theta|x)$.\n", + "\n", + "Consistency models use _scheduling functions_ to adjust some of the hyperparameters, for example the time discretization during training. Consequently, we have to specify the total number of training steps (_gradient updates_) before the start of the training.\n", + "For offline training with a given number of epochs, we can calculate it as below:" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "id": "516ac3c4-b66f-4cf0-a443-b00705a6ace5", + "metadata": {}, + "outputs": [], + "source": [ + "epochs = 30\n", + "total_steps = epochs * num_training_batches" + ] + }, + { + "cell_type": "markdown", + "id": "2c7980ec-5623-43e3-847e-16a5b6eb8777", + "metadata": {}, + "source": [ + "Apart from the usual parameters like learning rate and batch size, CMs come with a number of different hyperparameters. Unfortunately, they can heavily interact, so they can be hard to tune. The main hyperparameters are:\n", + "\n", + "- Maximum time `max_time`: This also serves as the standard deviation of the latent distribution. You can experiment with this, values from 10-200 seem to work well. In any case, it should be larger than the standard deviation of the target distribution.\n", + "- Minimum/maximum number of discretization steps during training `s0`/`s1`: The effect of those is hard to grasp. 10 works well for `s0`. Intuitively, increasing `s1` along with the number of epochs should lead to better result, but in practice we sometimes observe a breakdown for high values of `s1`. This seems to be problem-dependent, so just try it out.\n", + "- `sigma2` modifies the time-dependency of the skip connection. Its effect on the training is unclear, we recommend leaving it at 1.0 or setting it to the approximate variance of the target distribution.\n", + "- Smallest time value `eps` ($t=\\epsilon$ is used instead of $t=0$ for numerical reasons): No large effect, as long as it is kept small enough. Probably not worth tuning." + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "id": "09206e6f", + "metadata": { + "ExecuteTime": { + "end_time": "2024-09-23T14:39:53.339590Z", + "start_time": "2024-09-23T14:39:53.319852Z" + } + }, + "outputs": [], + "source": [ + "inference_network = bf.networks.ConsistencyModel(\n", + " subnet=\"mlp\",\n", + " subnet_kwargs=dict(\n", + " depth=6,\n", + " width=256,\n", + " ),\n", + " total_steps = total_steps\n", + ")" + ] + }, + { + "cell_type": "markdown", + "id": "851e522f", + "metadata": {}, + "source": [ + "This inference network is just a general CM architecture, not yet adapted to the specific inference task at hand (i.e., posterior appproximation). To achieve this adaptation, we combine the network with our data adapter, which together form an `approximator`. In this case, we need a `ContinuousApproximator` since the target we want to approximate is the posterior of the *continuous* parameter vector $\\theta$." + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "id": "96ca6ffa", + "metadata": { + "ExecuteTime": { + "end_time": "2024-09-23T14:39:53.371691Z", + "start_time": "2024-09-23T14:39:53.369375Z" + } + }, + "outputs": [], + "source": [ + "approximator = bf.ContinuousApproximator(\n", + " inference_network=inference_network,\n", + " data_adapter=data_adapter,\n", + ")" + ] + }, + { + "cell_type": "markdown", + "id": "566264eadc76c2c", + "metadata": {}, + "source": [ + "### Optimizer and Learning Rate\n", + "\n", + "We use an Adam optimizer with [cosine decay](https://keras.io/api/optimizers/learning_rate_schedules/cosine_decay/) to decrease the learning rate towards zero over the training time. " + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "id": "e8d7e053", + "metadata": { + "ExecuteTime": { + "end_time": "2024-09-23T14:39:53.433012Z", + "start_time": "2024-09-23T14:39:53.415903Z" + } + }, + "outputs": [], + "source": [ + "initial_learning_rate = 5e-4\n", + "scheduled_lr = keras.optimizers.schedules.CosineDecay(\n", + " initial_learning_rate,\n", + " total_steps,\n", + ")\n", + "\n", + "optimizer = keras.optimizers.Adam(learning_rate=scheduled_lr)" + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "id": "51808fcd560489ac", + "metadata": { + "ExecuteTime": { + "end_time": "2024-09-23T14:39:53.476089Z", + "start_time": "2024-09-23T14:39:53.466001Z" + } + }, + "outputs": [], + "source": [ + "approximator.compile(optimizer=optimizer)" + ] + }, + { + "cell_type": "markdown", + "id": "708b1303", + "metadata": {}, + "source": [ + "### Training\n", + "\n", + "We are ready to train our deep posterior approximator on the two moons example. We pass the dataset object to the `fit` method and watch as Bayesflow trains." + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "id": "0f496bda", + "metadata": { + "ExecuteTime": { + "end_time": "2024-09-23T14:42:36.067393Z", + "start_time": "2024-09-23T14:39:53.513436Z" + } + }, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "INFO:bayesflow:Fitting on dataset instance of OfflineDataset.\n", + "INFO:bayesflow:Building on a test batch.\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Epoch 1/30\n", + "\u001b[1m512/512\u001b[0m \u001b[32m━━━━━━━━━━━━━━━━━━━━\u001b[0m\u001b[37m\u001b[0m \u001b[1m5s\u001b[0m 8ms/step - loss: 0.3707 - loss/inference_loss: 0.3707 - val_loss: 0.3482 - val_loss/inference_loss: 0.3482\n", + "Epoch 2/30\n", + "\u001b[1m512/512\u001b[0m \u001b[32m━━━━━━━━━━━━━━━━━━━━\u001b[0m\u001b[37m\u001b[0m \u001b[1m4s\u001b[0m 8ms/step - loss: 0.3199 - loss/inference_loss: 0.3199 - val_loss: 0.3060 - val_loss/inference_loss: 0.3060\n", + "Epoch 3/30\n", + "\u001b[1m512/512\u001b[0m \u001b[32m━━━━━━━━━━━━━━━━━━━━\u001b[0m\u001b[37m\u001b[0m \u001b[1m4s\u001b[0m 8ms/step - loss: 0.3031 - loss/inference_loss: 0.3031 - val_loss: 0.2138 - val_loss/inference_loss: 0.2138\n", + "Epoch 4/30\n", + "\u001b[1m512/512\u001b[0m \u001b[32m━━━━━━━━━━━━━━━━━━━━\u001b[0m\u001b[37m\u001b[0m \u001b[1m4s\u001b[0m 8ms/step - loss: 0.2939 - loss/inference_loss: 0.2939 - val_loss: 0.2869 - val_loss/inference_loss: 0.2869\n", + "Epoch 5/30\n", + "\u001b[1m512/512\u001b[0m \u001b[32m━━━━━━━━━━━━━━━━━━━━\u001b[0m\u001b[37m\u001b[0m \u001b[1m4s\u001b[0m 8ms/step - loss: 0.2921 - loss/inference_loss: 0.2921 - val_loss: 0.3400 - val_loss/inference_loss: 0.3400\n", + "Epoch 6/30\n", + "\u001b[1m512/512\u001b[0m \u001b[32m━━━━━━━━━━━━━━━━━━━━\u001b[0m\u001b[37m\u001b[0m \u001b[1m4s\u001b[0m 8ms/step - loss: 0.2898 - loss/inference_loss: 0.2898 - val_loss: 0.2698 - val_loss/inference_loss: 0.2698\n", + "Epoch 7/30\n", + "\u001b[1m512/512\u001b[0m \u001b[32m━━━━━━━━━━━━━━━━━━━━\u001b[0m\u001b[37m\u001b[0m \u001b[1m4s\u001b[0m 8ms/step - loss: 0.2862 - loss/inference_loss: 0.2862 - val_loss: 0.2233 - val_loss/inference_loss: 0.2233\n", + "Epoch 8/30\n", + "\u001b[1m512/512\u001b[0m \u001b[32m━━━━━━━━━━━━━━━━━━━━\u001b[0m\u001b[37m\u001b[0m \u001b[1m4s\u001b[0m 8ms/step - loss: 0.2798 - loss/inference_loss: 0.2798 - val_loss: 0.2343 - val_loss/inference_loss: 0.2343\n", + "Epoch 9/30\n", + "\u001b[1m512/512\u001b[0m \u001b[32m━━━━━━━━━━━━━━━━━━━━\u001b[0m\u001b[37m\u001b[0m \u001b[1m4s\u001b[0m 8ms/step - loss: 0.2768 - loss/inference_loss: 0.2768 - val_loss: 0.2516 - val_loss/inference_loss: 0.2516\n", + "Epoch 10/30\n", + "\u001b[1m512/512\u001b[0m \u001b[32m━━━━━━━━━━━━━━━━━━━━\u001b[0m\u001b[37m\u001b[0m \u001b[1m4s\u001b[0m 8ms/step - loss: 0.2802 - loss/inference_loss: 0.2802 - val_loss: 0.2083 - val_loss/inference_loss: 0.2083\n", + "Epoch 11/30\n", + "\u001b[1m512/512\u001b[0m \u001b[32m━━━━━━━━━━━━━━━━━━━━\u001b[0m\u001b[37m\u001b[0m \u001b[1m4s\u001b[0m 8ms/step - loss: 0.2744 - loss/inference_loss: 0.2744 - val_loss: 0.2741 - val_loss/inference_loss: 0.2741\n", + "Epoch 12/30\n", + "\u001b[1m512/512\u001b[0m \u001b[32m━━━━━━━━━━━━━━━━━━━━\u001b[0m\u001b[37m\u001b[0m \u001b[1m4s\u001b[0m 8ms/step - loss: 0.2732 - loss/inference_loss: 0.2732 - val_loss: 0.2556 - val_loss/inference_loss: 0.2556\n", + "Epoch 13/30\n", + "\u001b[1m512/512\u001b[0m \u001b[32m━━━━━━━━━━━━━━━━━━━━\u001b[0m\u001b[37m\u001b[0m \u001b[1m4s\u001b[0m 8ms/step - loss: 0.2697 - loss/inference_loss: 0.2697 - val_loss: 0.3025 - val_loss/inference_loss: 0.3025\n", + "Epoch 14/30\n", + "\u001b[1m512/512\u001b[0m \u001b[32m━━━━━━━━━━━━━━━━━━━━\u001b[0m\u001b[37m\u001b[0m \u001b[1m4s\u001b[0m 8ms/step - loss: 0.2691 - loss/inference_loss: 0.2691 - val_loss: 0.2102 - val_loss/inference_loss: 0.2102\n", + "Epoch 15/30\n", + "\u001b[1m512/512\u001b[0m \u001b[32m━━━━━━━━━━━━━━━━━━━━\u001b[0m\u001b[37m\u001b[0m \u001b[1m4s\u001b[0m 8ms/step - loss: 0.2672 - loss/inference_loss: 0.2672 - val_loss: 0.2959 - val_loss/inference_loss: 0.2959\n", + "Epoch 16/30\n", + "\u001b[1m512/512\u001b[0m \u001b[32m━━━━━━━━━━━━━━━━━━━━\u001b[0m\u001b[37m\u001b[0m \u001b[1m4s\u001b[0m 8ms/step - loss: 0.2624 - loss/inference_loss: 0.2624 - val_loss: 0.2902 - val_loss/inference_loss: 0.2902\n", + "Epoch 17/30\n", + "\u001b[1m512/512\u001b[0m \u001b[32m━━━━━━━━━━━━━━━━━━━━\u001b[0m\u001b[37m\u001b[0m \u001b[1m4s\u001b[0m 8ms/step - loss: 0.2592 - loss/inference_loss: 0.2592 - val_loss: 0.2846 - val_loss/inference_loss: 0.2846\n", + "Epoch 18/30\n", + "\u001b[1m512/512\u001b[0m \u001b[32m━━━━━━━━━━━━━━━━━━━━\u001b[0m\u001b[37m\u001b[0m \u001b[1m4s\u001b[0m 8ms/step - loss: 0.2604 - loss/inference_loss: 0.2604 - val_loss: 0.3606 - val_loss/inference_loss: 0.3606\n", + "Epoch 19/30\n", + "\u001b[1m512/512\u001b[0m \u001b[32m━━━━━━━━━━━━━━━━━━━━\u001b[0m\u001b[37m\u001b[0m \u001b[1m4s\u001b[0m 8ms/step - loss: 0.2564 - loss/inference_loss: 0.2564 - val_loss: 0.2005 - val_loss/inference_loss: 0.2005\n", + "Epoch 20/30\n", + "\u001b[1m512/512\u001b[0m \u001b[32m━━━━━━━━━━━━━━━━━━━━\u001b[0m\u001b[37m\u001b[0m \u001b[1m4s\u001b[0m 8ms/step - loss: 0.2608 - loss/inference_loss: 0.2608 - val_loss: 0.2988 - val_loss/inference_loss: 0.2988\n", + "Epoch 21/30\n", + "\u001b[1m512/512\u001b[0m \u001b[32m━━━━━━━━━━━━━━━━━━━━\u001b[0m\u001b[37m\u001b[0m \u001b[1m4s\u001b[0m 8ms/step - loss: 0.2532 - loss/inference_loss: 0.2532 - val_loss: 0.2182 - val_loss/inference_loss: 0.2182\n", + "Epoch 22/30\n", + "\u001b[1m512/512\u001b[0m \u001b[32m━━━━━━━━━━━━━━━━━━━━\u001b[0m\u001b[37m\u001b[0m \u001b[1m4s\u001b[0m 8ms/step - loss: 0.2553 - loss/inference_loss: 0.2553 - val_loss: 0.3509 - val_loss/inference_loss: 0.3509\n", + "Epoch 23/30\n", + "\u001b[1m512/512\u001b[0m \u001b[32m━━━━━━━━━━━━━━━━━━━━\u001b[0m\u001b[37m\u001b[0m \u001b[1m4s\u001b[0m 8ms/step - loss: 0.2544 - loss/inference_loss: 0.2544 - val_loss: 0.2090 - val_loss/inference_loss: 0.2090\n", + "Epoch 24/30\n", + "\u001b[1m512/512\u001b[0m \u001b[32m━━━━━━━━━━━━━━━━━━━━\u001b[0m\u001b[37m\u001b[0m \u001b[1m4s\u001b[0m 8ms/step - loss: 0.2503 - loss/inference_loss: 0.2503 - val_loss: 0.2571 - val_loss/inference_loss: 0.2571\n", + "Epoch 25/30\n", + "\u001b[1m512/512\u001b[0m \u001b[32m━━━━━━━━━━━━━━━━━━━━\u001b[0m\u001b[37m\u001b[0m \u001b[1m4s\u001b[0m 8ms/step - loss: 0.2532 - loss/inference_loss: 0.2532 - val_loss: 0.3162 - val_loss/inference_loss: 0.3162\n", + "Epoch 26/30\n", + "\u001b[1m512/512\u001b[0m \u001b[32m━━━━━━━━━━━━━━━━━━━━\u001b[0m\u001b[37m\u001b[0m \u001b[1m4s\u001b[0m 8ms/step - loss: 0.2545 - loss/inference_loss: 0.2545 - val_loss: 0.1789 - val_loss/inference_loss: 0.1789\n", + "Epoch 27/30\n", + "\u001b[1m512/512\u001b[0m \u001b[32m━━━━━━━━━━━━━━━━━━━━\u001b[0m\u001b[37m\u001b[0m \u001b[1m4s\u001b[0m 8ms/step - loss: 0.2511 - loss/inference_loss: 0.2511 - val_loss: 0.2579 - val_loss/inference_loss: 0.2579\n", + "Epoch 28/30\n", + "\u001b[1m512/512\u001b[0m \u001b[32m━━━━━━━━━━━━━━━━━━━━\u001b[0m\u001b[37m\u001b[0m \u001b[1m4s\u001b[0m 8ms/step - loss: 0.2481 - loss/inference_loss: 0.2481 - val_loss: 0.2809 - val_loss/inference_loss: 0.2809\n", + "Epoch 29/30\n", + "\u001b[1m512/512\u001b[0m \u001b[32m━━━━━━━━━━━━━━━━━━━━\u001b[0m\u001b[37m\u001b[0m \u001b[1m4s\u001b[0m 8ms/step - loss: 0.2512 - loss/inference_loss: 0.2512 - val_loss: 0.2382 - val_loss/inference_loss: 0.2382\n", + "Epoch 30/30\n", + "\u001b[1m512/512\u001b[0m \u001b[32m━━━━━━━━━━━━━━━━━━━━\u001b[0m\u001b[37m\u001b[0m \u001b[1m4s\u001b[0m 8ms/step - loss: 0.2498 - loss/inference_loss: 0.2498 - val_loss: 0.2691 - val_loss/inference_loss: 0.2691\n" + ] + } + ], + "source": [ + "history = approximator.fit(\n", + " epochs=epochs,\n", + " dataset=training_dataset,\n", + " validation_data=validation_dataset,\n", + ")" + ] + }, + { + "cell_type": "markdown", + "id": "b90a6062", + "metadata": {}, + "source": [ + "## Validation" + ] + }, + { + "cell_type": "markdown", + "id": "ca62b21d", + "metadata": {}, + "source": [ + "### Two Moons Posterior\n", + "\n", + "The two moons posterior at point $x = (0, 0)$ should resemble two crescent shapes. Below, we plot the corresponding posterior samples and posterior density. \n", + "\n", + "These results suggest that our flow matching setup can approximate the expected analytical posterior well. You can achieve an even better fit if you use online training, more epochs, or better hyperparameters." + ] + }, + { + "cell_type": "code", + "execution_count": 16, + "id": "8562caeb", + "metadata": { + "ExecuteTime": { + "end_time": "2024-09-23T14:42:38.584554Z", + "start_time": "2024-09-23T14:42:36.076923Z" + } + }, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAg0AAAIQCAYAAAACUF3RAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAABkkUlEQVR4nO3de3xT9cE/8M/Jpen9Qpu2lLbcii1SKIKjXCcVAW9UcZsyXPFB3WQ8Pg7ZRZzPfqIvJ3NTBs82dNujiFMBN29FebgMinItCkWoQG25BSil6S1N0yZNcs7vj5KYtGl70iZpUz7v14sX5PSck+/5Nno++d6OIEmSBCIiIqJuKPq6AERERBQcGBqIiIhIFoYGIiIikoWhgYiIiGRhaCAiIiJZGBqIiIhIFoYGIiIikoWhgYiIiGRhaCAiIiJZGBqIiIhIFoYGIiIikoWhgTr15ptvQhAE5x+VSoXU1FQsXrwYly9f9vn7HThwACtXrkRDQ4PPzw18ez3nz5/3y/kHgv5cR/7+fHjSm/roz3XZU88//zxuvPFGiKLo3NbU1IRly5YhJSUFoaGhGD9+PDZt2tSj88s91+uvv44hQ4bAZDL1+FqohySiTqxfv14CIK1fv146ePCgtHv3bmnlypWSRqORhg8fLjU1Nfn0/f7whz9IAKRz58759LwO1dXV0sGDByWz2eyX8w8Ejt+5v34HveHvz4cnvfnMDLTP2+XLl6WIiAjpn//8p9v22bNnS7GxsdJrr70m7d69W3r00UclANI777zj9XvIPZfVapVGjRol/b//9/96dU3kPYYG6pTjBvLFF1+4bf/Nb34jAZDefvttn76fv24KJpPJp+fz93n7EkNDm4H4u+2tX/3qV9KQIUMku93u3Pbpp59KAKR3333Xbd/Zs2dLKSkpks1mk31+b8/18ssvSzExMfxdBRi7J8hrkydPBgBcuHDBuW3fvn2YNWsWoqKiEB4ejqlTp+LTTz91O06v1+MnP/kJ0tLSoNFooNVqMW3aNPz73//GypUr8ctf/hIAMHz4cGeXyJ49e5zHl5eXY+HChUhMTIRGo8Ho0aPxl7/8xe09Vq5cCUEQcPToUXz/+99HXFwcRo4cCaDz5mI5Ze/qvJ50da0OFRUVWLx4MUaNGoXw8HAMGTIE8+bNw4kTJzy+9/Hjx/GDH/wAMTExGDRoEJYvXw6bzYaysjLcfvvtiIqKwrBhw/D73//e4/ElJSW47777EB0djZiYGPzoRz+CXq/v9Bq8qXe51+yJnHqQ8/nwpLe/W0+fmY8//hjjxo2DRqPBiBEjsHbtWuc5XHk61rHf119/jR/+8IeIiYlBUlISHn74YRgMhi6v5cqVK4iMjMSCBQvctn/yySdQq9V45plnujy+N1pbW/H6669j4cKFUCi+vW18+OGHiIyMxA9+8AO3/RcvXozKykoUFxfLfg9vz/Xggw+isbGxx10h1DOqvi4ABZ+KigoAgFarBQB89tlnmD17NsaNG4fXX38dGo0G69atw7x587Bx40Y88MADAICCggIcPXoUv/3tb3HDDTegoaEBR48eRW1tLR599FHU1dXhT3/6Ez744AMMHjwYAHDjjTcCAE6ePImpU6ciPT0dr7zyCpKTk7F9+3Y88cQTqKmpwbPPPutWxvvuuw8LFizAkiVLuuz3lFt2b8/b1bU6VFZWIj4+Hr/73e+g1WpRV1eHDRs2IDc3FyUlJcjMzHQ75/33348f/ehHeOyxx7Bz5078/ve/h9Vqxb///W8sXboUv/jFL/Duu+/iqaeeQkZGBu677z634+fPn4/7778fS5Yswddff43f/OY3OHnyJIqLi6FWqz1ehzf1LueaPZFTD919Pjzxx+9227ZtuO+++/Dd734Xmzdvhs1mw8svv4yrV692eY3tfe9738MDDzyARx55BCdOnMDTTz8NAHjjjTc6PWbw4MH41a9+5QxQEydOxJ49e/CDH/wAP/3pT/Hb3/7W43GSJMFut8sql0rl+ZZQXFyM2tpa5OXluW0vLS3F6NGjOxw3btw458+nTp0q6729PVdycjKysrLw6aef4uGHH5b1HuQDfd3UQf2Xo6n60KFDktVqlYxGo/TJJ59IWq1WioqKkqqqqiRJkqTJkydLiYmJktFodB5rs9mk7OxsKTU1VRJFUZIkSYqMjJSWLVvW6ft11fw8d+5cKTU1VTIYDG7bH3/8cSk0NFSqq6uTJEmSnn32WQmAx75OT03vcsve1Xk96e5aPbHZbFJra6s0atQo6cknn3Rud7z3K6+84rb/+PHjJQDSBx984NxmtVolrVYr3XfffR2Odz2nJEnSO++806GbqX0dya33nl6zJ53Vg7fdE7743bavj+985ztSWlqaZLFYnPsYjUYpPj5eav+/U0+fN8d7/f73v3fbd+nSpVJoaKizTJ0xmUxSSkqKNGvWLOnw4cNSVFSUtHjx4i6PKyoqkgDI+tNZ3b700ksSAOd/8w6jRo2S5s6d22H/yspKCYD04osvdnk9vT3Xgw8+KCUlJcl+D+o9dk9QtyZPngy1Wo2oqCjcfffdSE5Oxv/93/8hKSkJJpMJxcXF+P73v4/IyEjnMUqlEgUFBbh06RLKysoAAJMmTcKbb76JF154AYcOHYLVapX1/mazGbt27cL8+fMRHh4Om83m/HPnnXfCbDbj0KFDbsd873vf6/a83pTdm/MC8q7VZrPhxRdfxI033oiQkBCoVCqEhISgvLwcp06d6rD/3Xff7fZ69OjREAQBd9xxh3ObSqVCRkaGW9eRw4MPPuj2+v7774dKpUJRUZHHa/C23nv6+/W2HuTwx+/WZDLhyy+/xL333ouQkBDn9sjISMybN8+r8uXn57u9HjduHMxmM6qrq7s8Ljw8HC+88AJ27dqFvLw83HHHHfj73//eoWvE1cSJE/HFF1/I+pOSkuLxHJWVlRAEAQkJCR1+1tV7d/Uzb/f39LPExERUV1fDZrN59T7UcwwN1K233noLX3zxBUpKSlBZWYnjx49j2rRpAID6+npIkuRsLnbl+B+Qo3l68+bNeOihh/C///u/mDJlCgYNGoRFixahqqqqy/evra2FzWbDn/70J6jVarc/d955JwCgpqbG7RhP5WnPm7J7c15A3rUuX74cv/nNb3Dvvfdiy5YtKC4uxhdffIGcnBy0tLR0OOegQYPcXoeEhCA8PByhoaEdtpvN5g7HJycnu71WqVSIj4/vtPvA23rv6e/X23qQwx+/W8c5k5KSOvzM07auxMfHu73WaDQAIOt6b7jhBgBtN9E333wTSqWyy/0jIyMxfvx4WX9cw5CrlpYWqNXqDu/V2eenrq4OQMfPbFd6cq7Q0FBIkuTx807+wTEN1K3Ro0fj5ptv9vizuLg4KBQKXLlypcPPKisrAcD57SQhIQFr1qzBmjVroNPpUFhYiBUrVqC6uhrbtm3r9P3j4uKc3xD/8z//0+M+w4cPd3st5xuON2X35ryO47q71rfffhuLFi3Ciy++6HZsTU0NYmNjZb2PN6qqqjBkyBDna5vNhtra2g43MAdv672nv19/1IM/frdxcXEQBMHj+IXugpGvHDt2DHfffTemTZuG/fv344033uj0d+Pw2WefdRiL0Jlz585h2LBhHbYnJCSgtbUVJpMJERERzu1jx47Fxo0bYbPZ3MYiOAaxZmdny3rfnp6rrq4OGo3GrTWJ/IstDdQrERERyM3NxQcffOD2LUkURbz99ttITU11fjNylZ6ejscffxyzZ8/G0aNHAXT+bSs8PBx5eXkoKSnBuHHjcPPNN3f409mNzx9l95anawXablKOa3b49NNP/bJwFgC88847bq/fe+892Gw2zJw50+P+van3zq7ZE7n14M23cX/8biMiInDzzTfjo48+Qmtrq3N7U1MTPvnkE6/O1RNlZWWYO3cupkyZgqKiItxzzz1YuXJlt7MufNE9kZWVBQA4c+aM2/b58+ejqakJ77//vtv2DRs2ICUlBbm5ubKvryfnOnv2bJeDYcn32NJAvbZq1SrMnj0beXl5+MUvfoGQkBCsW7cOpaWl2LhxIwRBgMFgQF5eHhYuXIisrCxERUXhiy++cI5GB9q+aQDA2rVr8dBDD0GtViMzMxNRUVFYu3Ytpk+fjhkzZuCnP/0phg0bBqPRiIqKCmzZsgW7d+/2W9m9JedagbYxCm+++SaysrIwbtw4HDlyBH/4wx+Qmprao2vpzgcffACVSoXZs2c7Z0/k5OTg/vvv7/QYufUu95o9kVsPXX0+PPHH7/b555/HXXfdhblz5+JnP/sZ7HY7/vCHPyAyMtLZjO4P58+fx2233YbMzEy8//77UKvV+N3vfofs7Gy8+OKLeOmllzo9NioqqtOWQrkcwfLQoUPO2QwAcMcdd2D27Nn46U9/isbGRmRkZGDjxo3Ytm0b3n777Q7dGYIg4JZbbvE4Vdbbc4miiMOHD+ORRx7p1bWRl/p0GCb1a50t7uTJ3r17pVtvvVWKiIiQwsLCpMmTJ0tbtmxx/txsNktLliyRxo0bJ0VHR0thYWFSZmam9Oyzz7otzvL0009LKSkpkkKhkABIRUVFzp+dO3dOevjhh6UhQ4ZIarVa0mq10tSpU6UXXnjBuY9jdLper+/0etqPEO+u7N2dtz2511pfXy898sgjUmJiohQeHi5Nnz5d2rt3r3TLLbdIt9xyS7fv/dBDD0kREREd3v+WW26RxowZ0+H4I0eOSPPmzZMiIyOlqKgo6Yc//KF09erVbutITr3LvWZP5NaDJHX9+fCkt79bT/Xx4YcfSmPHjpVCQkKk9PR06Xe/+530xBNPSHFxcd0e29l7dbWoVmVlpTRy5EhpwoQJHWax/PjHP5Y0Gk1AFryaMWOGdOedd3bYbjQapSeeeEJKTk6WQkJCpHHjxkkbN270uB8AacGCBZ2+h9xzSZIk7dq1y/m5psARJEmSAh9ViChQVq5cieeeew56vd7j6HfqHavVivHjx2PIkCHYsWNHXxfHb95//3088MADuHDhgtvYGLm2bt2Ku+++G1999ZWz1ag3CgoKcPbsWezfv7/X5yL52D1BROSFRx55BLNnz8bgwYNRVVWF1157DadOncLatWv7umh+dd999+E73/kOVq1ahT//+c9eH19UVIQFCxb4JDCcOXMGmzdv7nG3JPUcQwMRkReMRiN+8YtfQK/XQ61WY8KECdi6dStuu+22vi6aXwmCgL///e8oLCyEKIpuy0nL8Yc//MFnZdHpdPjzn/+M6dOn++ycJA+7J4iIiEgWTrkkIiIiWRgaiIiISBaGBiIiIpJlwIUGSZLQ2NgIDtUgIiLyrQEXGoxGI2JiYmA0Gvu6KF4TRRFVVVUQRbGvi3JdYH0HFus78FjngXU91PeACw1ERETkHwwNREREJAtDAxEREcnC0EBERESyMDQQERGRLAwNREREJAtDAxEREcnC0EBERESyMDQQERGRLAwNREREJAtDAxEREcnC0EBERESyMDQQERGRLAwNREREJAtDAxEREcnC0EBERESyMDQQERGRLAwNREREJAtDAxEREcnC0EBERESyMDQQERGRLAwNREREJAtDAxEREcnC0EBERESyMDQQERGRLAwNREREJAtDAxEREcnC0EBERESyMDQQERGRLAwNREREJAtDAxEREcnC0EBERESyMDQQERGRLAwNREREJAtDAxEREcnC0EBERESyMDQQERGRLAwNREREJAtDAxEREcnC0EBERESyMDQQERGRLAEJDevWrcPw4cMRGhqKiRMnYu/evbKO279/P1QqFcaPH+/fAhIREVG3/B4aNm/ejGXLluGZZ55BSUkJZsyYgTvuuAM6na7L4wwGAxYtWoRZs2b5u4hEREQkg99Dw+rVq/HII4/g0UcfxejRo7FmzRqkpaXh1Vdf7fK4xx57DAsXLsSUKVP8XUQiIiKSQeXPk7e2tuLIkSNYsWKF2/Y5c+bgwIEDnR63fv16nDlzBm+//TZeeOGFLt/DYrHAYrE4Xzc2NgIARFGEKIq9KH3giaIISZKCrtzBivUdWKzvwGOdB1Yw17dCIa8Nwa+hoaamBna7HUlJSW7bk5KSUFVV5fGY8vJyrFixAnv37oVK1X3xVq1aheeee67Ddr1eD7PZ3LOC9xFRFGEwGCBJkuxfIPUc6zuwWN+BxzoPrGCu7+TkZFn7+TU0OAiC4PZakqQO2wDAbrdj4cKFeO6553DDDTfIOvfTTz+N5cuXO183NjYiLS0NWq0W0dHRvSt4gImiCEEQoNVqg+4DF4xY34HF+g481nlgXQ/17dfQkJCQAKVS2aFVobq6ukPrAwAYjUZ8+eWXKCkpweOPPw7g2+YelUqFHTt24NZbb3U7RqPRQKPRdDiXQqEIyl+aIAhBW/ZgxPoOLNZ34LHOA2ug17dfryokJAQTJ07Ezp073bbv3LkTU6dO7bB/dHQ0Tpw4gWPHjjn/LFmyBJmZmTh27Bhyc3P9WVwiIiLqgt+7J5YvX46CggLcfPPNmDJlCv72t79Bp9NhyZIlANq6Fy5fvoy33noLCoUC2dnZbscnJiYiNDS0w3YiIiIKLL+HhgceeAC1tbV4/vnnceXKFWRnZ2Pr1q0YOnQoAODKlSvdrtlAREREfU+QJEnq60L4UmNjI2JiYmAwGIJyIGR1dTUSExMHbH9Yf8L6DizWd+CxzgPreqjvgXlVRERE5HMMDURERCQLQwMRERHJwtBAREREsjA0EBERkSwMDURERCQLQwMRERHJwtBAREREsjA0EBERkSwMDURERCQLQwMRERHJwtBAREREsjA0EBERkSwMDURERCQLQwMRERHJwtBAREREsjA0EBERkSwMDURERCQLQwMRERHJwtBAREREsjA0EBERkSwMDURERCQLQwMRERHJwtBAREREsjA0EBERkSwMDURERCQLQwMRERHJwtBAREREsjA0EBERkSwMDURERCQLQwMRERHJwtBAREREsjA0EBERkSwMDURERCQLQwMRERHJwtBAREREsjA0EBERkSwMDURE/ViFTo/VG4pQodP3dVGIGBqIiOTqzQ28p8cWFpXi3U+PoLCo1Ov3JPI1VV8XgIgoWDhu4ACw/KG8gBybn5ft9jdRX2JoICKSqTc38J4em5Gu9TqgEPkLQwMRkUy9uYH315t/hU6PwqJS5OdlIyNd29fFoX6OYxqIiIJYbwdKejNmgoMyiS0NRER9wNtv+BU6PdZ/WAwAWDw/F0DbDb+2wYTt+08D6HqsRGfv5023SW/GdNDAwNBARNQH2t+AuwsRhUWleG9bCQAgPjYCAPDup0cwd1oWFt41scNNv+18J5A3YQgSExM7veG37zbpqhwclEkMDUREfuLNDdhxU69tMCE+NgI5mSn4qqwSOZkp2F1cjgZjC+ZMy0JsVJjbTbv9uf9ReBj/887nGD0iCafOXkW4yoacMaOcx+RkpmD1hiLnce3L2FVrQn8dl0GBw9BAROQnhUWl2PDxYew8WIYRqfGIjQrD4vm5yEjX4kJlHXYeLMOZizV4ctFM5GSmYOfB0yg6XI6qGiMyh2uhq6xHs7kVzWYrIAE3ZiTjvx+bg8KiUmjjIrC/5CxyMlPcbv5vfHAQ1XUmVNcaMWbUYAxPHQTg2xv+6g1F2PDxYewvOYulC6Zj3aZ9qNDVAGgLCWxNoK4wNBAR+UGFTo/aBhMGxYTjxDeVOPFNJcI0agBAg7EFn+z5Gi0WK46XXUZsVBhOn7uKsnN6NJtbIYoSzl6qhdlsQ4vF6jxn2blqrPjjFrRa7QhRK3HxSgPOXKzBXd8dgwZjC3bsP43wsLb3aLXaceKbSny46wSUIZE4/s0V5GSmoLbBhIgwNYqPX0CdwYQ6Qwsy0hOcIcGb1gTOvLj+MDQQEflBYVEptu8/De2gSISoVUhOiEJMVBgajC0oLCqFpdUGAIiLDgMAlJZfgSSJEEUJAFBvaHGeSyEAggCoVUpcrTFizKhkLLxzIl56/d+oqjHirY+/AATA0mpDs7kVSoUAlVIBu03E+Uu1ePHvO3H2Yi2USgXUSiWsdjusNjtaW+3ISE/A0gXTAQDPrP0EAJytIXKukQMjry8MDUREfpCTmYL9JWeROUyLmvomZKQn4OSZqxiRGo+M9AQYTWYAgEatwoGSs2gwtkCSPJ8rPEwNpVIJs9kGURJxRd+I6jojGpssEEUJraINjkOtNhEAYBftUAiAKEk4o9OjpdUOAFCrFBAEAdq4SISEKHH67FWs27QPWcOT3AZaygkB7Mq4/jA0EBH1UvvpkBnpWuwuLsfxskqcuViDqhojGhpbYLOL2Hf0DAxGM7RxEajUG2G12bs9vygCTc1maEJUgARUVjfi5fW7YW/LB+gkawAABEGAxfrtezhCRU2DCTX1JmhClPjihA51hmaMvWEwvjmvx7HTl7DrUBm+KqvssuuBAyOvPwwNRES9UKHT46nVhSg5dRlKhYD42Ajk52XjwLGzaLFYkRAXDgAwNLW1LJhaWgEAF640yDq/IKBtICTg7NIA4AwM3XF0d3S2vcXSds5TZ68iLioMtYZmbN9/GherGmA0WZyzOThugQCuCElEJJunFRELi0pRoatB5vBE3H/7TdDGReDe//pfnDxzFVabHecu13d645ajsy4LXx8vihKMzRbnMfpaIxbeNREXq+rxx7f24I9v7en0WK4Uef1gSwMRkUyeBv659utnpGtxx5JXUV3X5DxG6u1dP4BaXboxahtbsGXPCVyqMsDSakPJqUtu6zu44oDI6wdDAxGRTJ4G/mWka5Gfl431HxbjYlU9ys5V91Xx3KiUCoiizD4MD+x2EV9XXAUACADqDCZs+PgwgI7BgAMirx8MDUREMrkO/HOsUZCTmYKX3yzC8bLLsNnFXncn+IrNLkIh+OZcEoCGxhaEh9shSaKzxQFAhwGgNLAxNBAR9YBjtUdBAPR1Tc5ZCW3THPu4cH4gSkCTqRWrN+yBUqnAv3aUoNlsRV1Ds3P2xm9/dneflpH8jwMhiYh6ID8vGymJ0aiuNcIuSlBe+1rv68CgUAiICFPDU6NBZHgIIq6tAOlvjve32kSYLTac0dWisrqxbcEolbLLY3cdKsP3lr2OXYfK/F9Q8iuGBiKiHhMgCG2rLw5JivHqSEW7vgPHa7VKAYXQ9m+lQsCcqZn44Z0TEBER0uEcg2LCoVR0/b9xpVJAmEaNwdooqJQ9769on4Wka2WOjw1H5vBEZA1PxOoNRdh1qKzDTIp1m/bhQMl5rNu0r8fvT/0DuyeIiLrh6RkLhUWlqKw24IZhiQjVqHHb5FFYvWGPs5uiO/Ex4aipN0ECoAlRITsjCRAELLxzIg4dv4D9JWdhtdoxPisV+0vOwtRshYC2m7VaJcBqk2BqbkVkhAaNJovzvBFhIdCEqGCz2aFUChAAmFutyEpMRmOTGbaWtjUfIsND0Gy29ng6qCAICA1R4vJVAy5dNcBsscJosmB/yVlcqKwH8O2ASccy1Y6/KXgxNBARdaOrqZZnLtZgx/7TANpWbgTgvLl3RRAEqNVK2GwiWq02nDpbjSd+dAsK8iehIH+SW1DJyUxBs7kV2rhIpCXH4dTZKuwvOQdDkxm3jb0BaqUCVbVGJMdH4Ykf3YLT56rRYGyBAEAhtqC60Y6Fd07A2rc/g6mlAQAQolbBbLFCutb44O0ATkmS0Gz+drGp6lojHvneFGjjIvDu1qM4c7EGFTo9MtK1mDU5E7MmZ3r3BtQvMTQQEXXDdUrhrkNlWLdpH+69dazbPtq4CIRqlDC1iJ0GBpVSgZjIUMTFhOOnD0zD6XPV+Hj3cejrTYiNDkNOZgpWbyhCTmaKcwlnAPiqrBJ/+vX3na0cz6z9BMUndBAApCXH4a1VBVi9oQjvfnoEH+0+gQuV9RiaEgfdlXpMzU6CUqmEvt6EG0cm42JVAwQALWYrwkJDYGm1QakQnCtD9lSdoRkAcOj4BZSWX0HZuasYmdb29Ew+CXPgYGggIuqG61TLp1YXYv/Rczj+TSVCVEokxkchNTkWF6saEKJWwWoTnYskCULbN3jHjAoJwI9/MNVtnYNbc0dh3aZ9WLpgOr4qq8S7nx5xa+IHgHc/PeK2nPPi+bnOnzn+7QgYjsDR9vdljBkWjaHpQ5Cfl33tsdgaDE2JQ6hGDV1lHQbFhKO5pRV2UXJb3MlbdlHCX9/bjxFpCVCrlBiaEofaBhPWf1iM7ddaYrjwU/BjaCAikqlCp0dyQjSStVEwNbei2WzFqTNVgCAgNESF0SOToK9rgu7acyUUggAoBKQmRqNS34i46HDkZKa4ndO16X5oyiAA6NDSAAC1DSa3LpL20xtdg43jfHmTRqG6uhqzZyRCoVBg8fxcZ/A4eOwc/uedz7HgjptQZ2hB0eFynLlY06t1JhpNZtTUGTE0JQ43jkzG9v2nMXdaFhbeNZELPw0QDA1ERDKt/7AYO/afxpxpWaiqacTBr85DlAAFJAxPjQcgoEpvdO4/IjUeMyeNcrYmVOhq8FVZZaf9+55u/EBbSKjQ6Z03fLnOXKzBrn3HMWv6eIwamtihxeRSlQH/PlSO8FA1DMYWDNZGo66hGYIgoMVi9bp+7HYJF640QBOiwtTxw51hgd0SAwenXBIReSk2KgwvLc/H/FnjMCYjGXOmZQEATp2pguuIhvGjUxEfG4GhKYPw0vJ8PHTPpB5/43bc8L25AW/ZU4rdxeV4es0WtymQjhaT8aNTMCI1HhW6GmSNSMIrv7wXY0Yly3pcd2eUSgVCQ1QYFBPW43NQ/8WWBiIimVyb9zPStXhy0UwUFpWitsGEfUfPwS5KmDJuKMrOVyMmsu2m6dql4E2fvqdpnt6aNzMb35y5gH8f1uGPb+1xG5i47+hZZKQn4L7bxmFkWoKzS2REagJOfHMFSkmASqmEzW6H3YtpmXa7CEOTGS+vL0KoRoW3Cg9j2k0j8OSimWxxGADY0kBEJFP7b/uOpaSPnrqIO2aMxqL87+Cln9+De24d53zSpdz+fNfHS1fo9HhqdSE2fHwYhUWlPS7vyLQEjExNgFKhwNlLtXj30yMuQSTB2V2y/KE85yDM2KgwPHTPJCQlRKHVZoN4bZCD0ssHWdhFCaYWKyqrG1FYVNqr66D+gy0NREQ9lJ+Xjf0lZ1Ghq8HsKVnOb/ENxhYAbd0YclsXXNeCAIAKXQ0y0hOc0zB72uJw25QboNREICdziHNwZUa6Fi8tz3cGCMe1OP7OSNfi1txRePnNIujrjKiqMSI+NhxXXMZryBUepsZd3x3DgZADBEMDEVEPZaRrsXTBdKzbtA85mSnOloeUxBjnOAfHAkfd8fR4aUcIab+wlDcGa2OwrGAUFAqF2+BK10GRrq8dLR7auAiEh6rxwzsn4N+HylFxQe/p9F1SCMB3J47EyLQEXKis43oNAwBDAxGRDJ2NMdhdXI7jZZXYXVyOxfNznS0P4aFqFB+/gPjYCFk3+/Y38fYrT/rrm3r763IEn2ZzK4xNFuiu1KOhsQXN5lYoFQIEBWCzyR3jIOCrssvYdegbxEWHQalse7AV12sIXgwNREQyuH7jd13l0JVrs7+ntRZ6on2Y8DXHdTkWj8rJTEFGegJKy68gbXAsFtxxE/59qBzauAjnEtb7jp6TdW5RkpxdGnWGFkwZP4zdFEGOoYGISAbXb/yuAcJ1dUZHV4Sjmf+rskqflsEXMyran8ux2JTr4lGuwWfdpn2orDZg9pRM5Odl479efB8qpQI2u7wHcwFtK2NGhKlRZ2jGhco6dk8EMc6eICKSwXXmRH5ettvCRfGxEdi+/7TbDAFHsHBsc50d4Q3X49qfszcc53LMnlg8P9ftmhwzKlwHZP7Xi+/j1JkqZA5PRPrgWDjmU8iZV2FsbsXJiiqs+OMWj4/PpuDAlgYiIi+17zLobBCj69/tuwHktha07xZp/z491f5cnrpB2reunDpTBYvVjhtHJuPspVrnctndjXCQpLanYqpVCjQ2mbFu074Oj8+m4CBIUm9WGu9/GhsbERMTA4PBgOjo6L4ujldEUUR1dTUSE9vWiSf/Yn0H1kCsb2+6Cxz71jaYsH3/aSy8a6KsG2ZvuiR8WecVOj3+68X3UXauGj+8cwJuzR2Fn636wLkehRwqpYCcrFT84j/y3KZ/DhQD8TPeHlsaiIhceHOT9mY6pOtYB2+eIeHvgZDemDA6FRNGp2Lx/FxkpGtxy3cy8P7OryDKXDHSZpdw8UodALYwBKuBGYWIiHrIm3EDrmMb5OrJMyT6g8KiUmzffxrxsRHOssdGhSE8VI0wjRoAoFR0v3JkdZ0JL79Z5Pfykn+wpYGIyEV+XjZqG0yobTB1uzBTf2oFAHw7u6L9+TyNp3DMHDl1tgqnzlxFo8kMq63rWRUCAG1cRK/LRn2DLQ1ERC46mw3RnZ7OjvAlX86uaH++9i0kjkABACe+uYKmltZuAwMACAoBaclx/aK+yHtsaSAiaqcnsxQcKynuLzmLl5bn90n3g69Xj+zqfI5AkTtuKFKTY2E0mWG12WEwmtFstnZ6To1ahazhib1eHpv6RkBaGtatW4fhw4cjNDQUEydOxN69ezvd94MPPsDs2bOh1WoRHR2NKVOmYPv27YEoJhERgJ6NO3B9cqQvvun35Ju4r8dLdHW+nMwUDE2JAwBcqmrAlRoj6gwtiI0O6/KcLRYrfvu3HdiypxS544Zyhcgg4/fQsHnzZixbtgzPPPMMSkpKMGPGDNxxxx3Q6XQe9//8888xe/ZsbN26FUeOHEFeXh7mzZuHkpISfxeViKjHHEtIP3TPpA43wp4EAF93NfRG+8d2r95QhN3F5ajQ1eDspVrMmZaFG4ZqoVQIiAgN6fZ89YYWnDxThaqaxqAbEHq983v3xOrVq/HII4/g0UcfBQCsWbMG27dvx6uvvopVq1Z12H/NmjVur1988UV8/PHH2LJlC2666SZ/F5eIqMdcB0a6DiL0pim+/fLOcr6J+3oAZHvtH9v97qdHMHdalrNlZfaUTDy5aCaeWl3ofFhXZ10USoUACUCYRoV7bx3r87KSf/k1NLS2tuLIkSNYsWKF2/Y5c+bgwIEDss4hiiKMRiMGDRrkjyISEflFZys5dneD70lfv5xjehosKnR61DaYMHdalvM6ahtMAIClC6a7LdJ0761j8du/7UCrzd7p+ezX1nSQJEBfb5JdDuof/BoaampqYLfbkZSU5LY9KSkJVVVVss7xyiuvwGQy4f777/f4c4vFAovF4nzd2NgIoC1siKL8B6r0B6IoQpKkoCt3sGJ9B9b1Vt/zZo4BIGHezDEYkRqPZQW3AADW/GMPNm49CkDCsoKZXR4nt646O8a1zguLTnT5vp0pLDqBHQdO44d3TsCI1HjsLv4G/7f3azQ0tgCQ8B/35mL9h4dgMJrx+Zfl17YDXS3XIAjAjRnJXl1jMAjmz7jcFSwDMntCENw/PZIkddjmycaNG7Fy5Up8/PHHSExM9LjPqlWr8Nxzz3XYrtfrYTabe1bgPiKKIgwGAyRJGrBLkPYnrO/AGqj1fUVvQPHxC8gdNxSDtTHO7VEaYOHcGwG0LS3skDdhCMJVNuSOG+K2vbvjutLZMa513tX7dnYNnsq75d9fID5SgYToSMSGAbv2HcPxrytgsdmRHKdGYkxst+UdFBOGp3+chyiN/GsMBsH8GU9OTpa1n19DQ0JCApRKZYdWherq6g6tD+1t3rwZjzzyCP75z3/itttu63S/p59+GsuXL3e+bmxsRFpamnP2RTARRRGCIECr1QbdBy4Ysb4Da6DW97vbT2Lj1uNotqmwrGBUt/snJiYiZ0z3+/mCa50nJyd3+r5dXUP78s677Tsor9yDEUPice/ctsWd3ttVhvLzdQgPDYGx2dLtstIpiWLA6iCQBupn3JVfQ0NISAgmTpyInTt3Yv78+c7tO3fuxD333NPpcRs3bsTDDz+MjRs34q677uryPTQaDTQaTYftCoUiKH9pgiAEbdmDEes7sAZifefnjQUgID8vu19el5w6z88bi9qGZtQ2NOPspdouxzzcNiULt03Jcr6u0OmhCVFDrVYhJSkGp89Wo6vMIAjA7dNHo+hwOdZt2oelC6Zj1uTMHl1bfzQQP+Ou/N49sXz5chQUFODmm2/GlClT8Le//Q06nQ5LliwB0NZScPnyZbz11lsA2gLDokWLsHbtWkyePNnZShEWFoaYmJhO34eIqC/4cilpf8+C6IxjFcx3Pz2C+NgIr66nsKgUuso6xMdGYPSIJJysuNrl/pIEbNx6FIVFX6Pm2oDKgRQaBjq/R6EHHngAa9aswfPPP4/x48fj888/x9atWzF06FAAwJUrV9zWbPjrX/8Km82G//zP/8TgwYOdf372s5/5u6hERB0Ecrljf6/N0NW1dPfwrc6OzclMQXhYCBqbzNiy5yQcjQxdjVoztVhR29CEtORYLF0wvYdXQ30hIAMhly5diqVLl3r82Ztvvun2es+ePf4vEBGRTIFc7tjXy0C319W1dNdi0tmxX5VVormlFaIkQa1SwGoVIEoSuntYdqgmBL97ch5bGYIMnz1BRNQFf9/IXfmiq8O1i2NEarzbz3pzLZ0dm5+Xjf0lZ1FafgWxsRFIH6zGyTNdz4iIjQxFq82O3cXlDA1BZmCO1CAi8hFfP8/B37rq4ujNtXR2bEa6FksXTEd0ZChqG0woO9d9N47FavP6/al/YEsDEdEA4o+WEdelrV1XgHTYXVyO2oZmWFrtzhUfu9JisSEmMtR57mAJZMTQQEQ0oLh2cfhqZUJH68X+krOo0NVg58EyTBidisXzc5GRrkWDsQWWVvdnTQhAl+MajCYzPtp13OvZGtS3GBqIiKhLjlaLnMwUrNu0DyWnLqPigt55w4+NCoMoAXa7CJVSgM3e/UBIUQLUaiUfjR1kGBqIiKhLrq0XQ1MGYf2HxQDawsQ/Cg/j08+/RuawBFyobECzubXb8wkCoBAETLtpBLsmggxDAxERdct1VsZvf3a38/VbhYdRWd0Io8mCFksrpO6HNAASMCItHrFRYRzTEGQ4e4KIaADw9yJU7WdlOF6PTI1HVIQGAgC5QygkAOcr6/DmR8X441t7/FJe8g+2NBARDQD+XoSq/ayM/Lxs1DaYcODYOVhabbDa7LLOIwhtS0lbbW0Jo21wJVsbggVbGoiIBoDuloHujJwWCk/PxHA8r+LcpVrY7SKSE6IQExkKZRd3FaVCcHZfRIaHIHFQJKxWu9+WzSbfY2ggogEnkM+L6Auerq+nCzfJed5FZ/vk52Vj9MhkhIWGYPqEkcgYqu3yCZeuazjcMEyLe24di3tnjeMMiiDC0EBEA46/H/zU13x5fZ21ULgGE8c+OZkpbmElI12LP/36e3h84QwAwPGyy10vzuCi7Fw13ttW4jyPJwM9/AUjhgYiGnB62lQfLHx5fZ21ULgGE8c+X5VVdggrrtMx7TLWZwDaFn4amhIHc6sNDcaWTvcb6OEvGHEgJBENOL548FN/Fojr87QcdVdLVMdGhUFeZGhTU2+C1M38zEA+LIzkYUsDERF14KkFIiNdi/y8bBQWlbp1GTj+HRGmln3+OkNbC0NsVJhXZaC+xdBARESyeeoy+ONbe7Dh48NotXa/UINSoUBIiAqCQsCooW1hgGMWggdDAxERyVKh06O2wYS507LcugxKTl1Eq9WOVmv3azVIkogQlQK5Y9MRqlHjk8++5piFIMLQQEREnXKdwVBYVIrt+08jPjbC2WVQodM7uxoAOB953Z5CEAC0PajK2NyKE+VXoKusQ0Z6AscsBBEOhCQiok65rjSZk5mC/SVnkZOZ4vbz5pa2h1QpFQJio0JhbDZ3WFJadBn0GBqigkIQkDUiCS8tz+eYhSDClgYiIuqU6/TO3cXlOF5WiQ/+fdxtDYfU5FgAgCAIuHCloctnUISolfjRvImw2cVuZ09Q/8PQQEREnfI0g+HspVps+PgwnlpdiAuVddDXNQEAbPbuB0KmD47Fpm0lMJosOPjVBY5nCDLsniAius54epaEHIvn5yI+NgI5mSlYt2kfKnQ1WLdpHxpNZtnnOHepDhAEqFUKjLthMGobTHxgVRBhaCAius709ImYrotKDU0ZhMKiUuRkpuB4WSUMTfKCw5CkGFypMUKjViI8TOMcWDmQF+MaSBgaiIiuM47HWvf0W377lgrtoAjZoWHutNG4WFWP499UYur4YZh20wjOnggiDA1ERNcZx2Ot3/30SJff8tuHg12HyrBu0z5IkoTDJ3T4145jyJs0qq3LQYbI8BBcrKrH0ZOXUF3XhH8fKsf/vbbEl5dGfsbQQER0HZLzXAdHN0ZtgwnxsRHYebAMx05dRqhGiVarHRW6Gpy7VIsQlRItrbYu30+tUsDSasP2/WVQCIBSqcCI1HifXhP5H0MDEdF1qLuHXrmu/ggA7356BDeOTEJqcgwiw0Nw8sxVSBJgFyWYrV0HBgCw2r6dWREaqsad3x2DJxfN7PV1UGBxyiUREXXguvrj4vm5WHjXRKQlx6FtaQUBapUSkeEhAABvlltQKhUIUaswMi2BMyaCEFsaiIioA9fuC0erxK5DZTh97ioiwkJwobIe4aEqWFptECVAgASbvS09KIS25aJdCQBUKiUkSUJKYgzy87J7PPWT+g5bGoiIBqgzF2vw0a7jOHOxxutjPS3q5FgRMioiFDeNHgJJAibnDMPsKTdcCw5tRAlQKRVIHBQJlVJAmEaFsFA14qJDYRdFhGrUyEjXenxiJvVvbGkgIhqgtuwpxRfHzqDZpsLyh2712Xljo8Lw5KKZWP9hMQDgwLFzENs1LURHhmLCjanYefAbSBJwc3Ya7r11LD7afQJLF0wHIG8wJvUvDA1ERAEWqGb5eTOzEa6yYdZ039yUHStCupb7vW0lmDJ+GOobm1FnaIZGrYKx2XJtemUDQlRKREVokJwQjSnjh6Mgf5LzfN0NxqT+h90TREQBFqhm+ZFpCbh31jiMTEvw+ljXR2I7eOqycIiLDsfItARMGT8MSqUCl64a8M35asTFhKHZ3Iod+0+zG2IAYEsDEZEPyWlF6O/N8hU6PZ5aXYgKXdtYiM5aAxwtDzsPluHkmSqEqFWYOn44xt2QguNllyEIAqbdNAKxUWEA+u/1knwMDUREPiTnuQ6BbJa/ojfg3e0nkZ83VnZXSGFRKSp0NchIT+jyRp+RrkV+XjaOnb4EbVwEJtyYhsXzc7F4fq5zvMPi+bmcGTGAMDQQEflQf2tFKD5+ARu3HgcgdBtUHK0kOZkpeOieSbLGXKz/sBhFhyugVikxPisVGelaVOj0zqdhckrlwMLQQETkQ/1tcF/uuKFotqlkhZiePP2ywdgCURQRHqpBTmaK23n2l5zFhcp65zLUDA/Bj6GBiGgAG6yNwbKCUVAouh/33pNWktioMCgUCjSbbfiqrBKzJmc6j8/JTMFXZZWobTA5w0h+XjZbH4IYQwMR0QDU1tVwAnkThiAxMVHWMT1pJVk8P9f5b0dYcD3PrMmZzu4KR2DwtjWD+g+GBiKiAaiwqBQbtx5FuMqGnDGjvDrWm3UkMtK1+O3P7u52H0dA6G9jPsg7DA1ERANQ201ZQu64IV4f68/WgP425oO8w9BARDQAZaRrsaxgJqqrq70+lq0B1BmGBiIicsPWAOoMl5EmIgoQT0sz96X+Vh7q/xgaiIgCpL89Crq/lYf6P3ZPEBEFiC/GCvjyCZkcu0DeYmggIgoQX4wV8OXMBo5dIG8xNBARBRG2DlBfYmggIgoibB2gvsSBkERERCQLQwMRERHJwtBAREREsjA0EBERkSwMDURERCQLQwMRERHJwtBAREREsjA0EBF1gg90InLH0EBE1Ak+0InIHVeEJCLqBJdsJnLHlgYi6tb12kzvWLK5t0+TJBooGBqIqFtspicigN0TRCQDm+mJCGBoICIZ+GRFIgLYPUFEREQyMTQQERGRLAwNREREJAtDA9F16nqdRklEPcfQQHSd4jRKIvIWZ08QXac4jZKIvMXQQHSd4jRKIvIWuyeIiIhIFoYGIiIikoWhgYiIiGRhaCAiIiJZGBqIiLzA9S3oesbQQETkBa5vQdczTrkkIvIC17eg6xlDAxGRF7i+BV3P2D1BREREsjA0EBERkSwMDURERCQLQwMRERHJEpDQsG7dOgwfPhyhoaGYOHEi9u7d2+X+n332GSZOnIjQ0FCMGDECr732WiCKSURERF3we2jYvHkzli1bhmeeeQYlJSWYMWMG7rjjDuh0Oo/7nzt3DnfeeSdmzJiBkpIS/PrXv8YTTzyB999/399FJRqwuCAREfmC30PD6tWr8cgjj+DRRx/F6NGjsWbNGqSlpeHVV1/1uP9rr72G9PR0rFmzBqNHj8ajjz6Khx9+GC+//LK/i0o0YHFBIiLyBb+u09Da2oojR45gxYoVbtvnzJmDAwcOeDzm4MGDmDNnjtu2uXPn4vXXX4fVaoVarXb7mcVigcVicb5ubGwEAIiiCFEUfXEZASOKIiRJCrpyB6vrqb7nzRwDQMK8mWP67Hqvp/ruL1jngRXM9a1QyGtD8GtoqKmpgd1uR1JSktv2pKQkVFVVeTymqqrK4/42mw01NTUYPHiw289WrVqF5557rsN59Ho9zGZzL68gsERRhMFggCRJsn+B1HPXU31HaYCFc28EIKK6urpPynA91Xd/wToPrGCu7+TkZFn7BWRFSEEQ3F5LktRhW3f7e9oOAE8//TSWL1/ufN3Y2Ii0tDRotVpER0f3ptgBJ4oiBEGAVqsNug9cMGJ9BxbrO/BY54F1PdS3X0NDQkIClEplh1aF6urqDq0JDsnJyR73V6lUiI+P77C/RqOBRqPpsF2hUATlL00QhKAtezBifQcW6zvwWOeBNdDr269XFRISgokTJ2Lnzp1u23fu3ImpU6d6PGbKlCkd9t+xYwduvvnmDuMZiIiIKHD8HoWWL1+O//3f/8Ubb7yBU6dO4cknn4ROp8OSJUsAtHUvLFq0yLn/kiVLcOHCBSxfvhynTp3CG2+8gddffx2/+MUv/F1UIiIi6oLfxzQ88MADqK2txfPPP48rV64gOzsbW7duxdChQwEAV65ccVuzYfjw4di6dSuefPJJ/OUvf0FKSgr+53/+B9/73vf8XVQiIiLqgiA5RhkOEI2NjYiJiYHBYAjKgZDV1dVITEwcsP1h/Ulf1XeFTo/ColLk52UjI10bsPfta/x8Bx7rPLCuh/oemFdF1I9xoSUiClYBmXJJRN/Kz8t2+5uIKFgwNBAFWEa6FssfyuvrYhAReY3dE0RERCQLQwMRERHJwtBAREREsjA0EBERkSwMDURERCQLQwNRgFTo9Fi9oQgVOn1fF4WIqEcYGogChIs6fYsBiig4cZ0GogDhok7fcgQoAFyzgiiIMDQQBQgXdfoWAxRRcGJoIKKAY4AiCk4c00BEQYtjI4gCi6GBiIIWB5cSBRa7J4goaHFsBFFgMTQQUdDi2AiiwGL3BBH1OxyrQNQ/MTQQkVcCcUPnWAWi/ondE0TklUAszMSxCkT9E0MDEXklEDd0jlUg6p8YGojIK7yhE12/OKaBiIiIZGFoICIiIlkYGoiIiEgWhgYiIiKShaGBiIiIZGFoICIiIlkYGoiIiEgWhgYiIiKShaGBiIiIZGFoICIiIlkYGoiIiEgWhga67jge7XzmYk1fF4WIKKjwgVV0XanQ6fHU6kJU6GoASFg498a+LhIRUdBgSwNdVwqLSlGhq0FGegLmzfTfo52JiAYitjRQv7TrUBnWbdqHpQumY9bkTJ+dNz8v2/n3iNR4VFdX++zcREQDHUMD9UvrNu3DgZLzAODT0JCRrsXyh/IAAKIo+uy8RETXA4YG6peWLpju9ndvVOj0KCwqRX5eNjLStb0+HxHR9YqhgfqlWZMzfdbCUFhUinc/PQIAzlYGIiLyHkMDDXiu4xgo8BwtPfNmjkGUpq9LQ0S9wdkT1G841k+o0Om73OYtxzgGdk30DUdLz5Y9pX1dFCLqJYYG6jccN5fColK3bRs+PnxtbYWeBwfqO/l52Vh410ROcSUaABgaqF+o0OlR22DC3GlZyM/LdrYwaOMiEKJW4vTZq9fWWHBvefBFSwT5l6OlZ2RaQl8XhYh6iaGB+oX1HxbjvW0lANpuMo4Whlc2FEFf14T0lEHIyUzBU6sLseHjw87WiPatEwwRRET+w9BA/YZdFHHg2Dk8s/YT5GSmICM9AQ2NLbCLEiaMTsXu4nKUnLqElMQYt8GNC++a6HztqYtDDoYNIqLucfYE9QuL5+ei6HA5vq6owtcVVTh1tgqR4RpEhmtww7C2AYwNxhYoFQpMGJ3qHNToulgT0BYiahtMqG0woUKnlz34kdMyiYi6x9BAfcoxHS8nMwVX9Abn9sMndACAVqsdKpUCuisNmDstCz9dML3LqZMZ6VrEx0bg3U+PID42QnYA4LRMIqLuMTRQj3mz0qJrOPiqrBLauAi88WEx9HVGSBKQnjIILWarc3+bzY6RaQk4X1kHSZJgtlhxsaoei+fndvtePQkA7VssiIioI4YG6jE5TfqOsFDbYMInn30NQQAamyzQqJWoNTQDAMI0Kugq6yBK3x4nSkBUpAahGjWqatqCxe7icozPSu325s4A4HtcipuIAA6EpF5oPwjRE9dgkZGegJp6E0wtFhiaWpz7mFtt0NebOhxbcuoyjCYLJAkQBCA2KhQ5mSm+vxDqVk8HmBLRwMKWBuoxOd/oHYGipt6IL0svwm4XIUmAze7SrCB1PE6lFJz7qNVK3DBUizpDM74qq/TpUy9JHo75ICKALQ3kZ45gsXnbMbRYrLCLEsI0aoSGfJtXPWQGAMK3PxclhGpUuPuWMbxpBYCn6adcipuIAIYGCoAKnR5R4SHO1yFqJcytti6PiY4MxWBtFGIjQxEdqYGush7xsRG8aQUAuyKIqDPsniC/KywqRaOpFQpBgChJMDSZu9xfAFB3bZCkIAAatQrfGZvOVoYAYVcEEXWGLQ3kM52tqqiNi0B0pAbhYWpZ53HtrpCkttd1hmZcqKzzXWGpU+yKIKLOMDSQz3TWrP3R7hO4WtMElcL7j5sgAFarDSfPVGHdpn0e9+ES0EREgcHuCfKZzpq1ly6YjjqDCRW6WgBt3Q+eBz+2iYsOg9FkRkSYBtpBEUhOiEaLxYalC6Z73L+79SK4xgARkW8wNJDPdDYFc9bkTLzw1x2wXBv82FVgAICmZgsS46NQU29CVY0R359zU5dTO7vrg+dzJYiIfIOhgfyuQqd3e65Ed2KjwhAXHY7a+mYIggBtXESX+3e3XgQH9hER+QbHNJDfrf+wGEaTBYIAt6mXnWkwtqC61ojQUBVMLa34aPeJXr0/B/YREfkGQwP5VYVOj6OnLkKtUiImIhTG5tZuj7GLEtJTBmHJ/VORlhyLe28d6zwXBzwSEfUdhgbyq8KiUlRWN2LimDSEhChlHSOJEkakxuPAsfNotdqdz6XgokNERH2LYxrIp9rPVHCMI9DGRWDt2/LWWQgJUeHkmSpU1xqRNSIJ+XnZqNDpUdtgwtxpWRybQETUR9jSQD7lqTXgzMUa/HrtJ9BdaXBuEwQPB19jtdrwzQU90lMG4aXl+chI12L9h8V4b1sJAHBsAhFRH2FLA/lU+5kK6z8sxvs7jkFsN88yTKNGi8UKycP8S0n6touCAaH3uE4FEfkKQwP5lMfpj4KA9umg2Wzt9BzStT+xUWHObYvn5yI+NoJdEz3AdSqIyFcYGsivFs/PxcWqeuwuLker1d7t/kqFAIVCwORxQ7F4fq5ze3drMVDnuE4FEfkKxzSQX2Wka5GWHAerjMAAtE23tNpENFvaWiI4xbL3uE4FEfkKQwP51a5DZfh49/G2B07IpFAIGJGawCmWneB6FUTUV9g9QX718pu7UV1ngiAAkeFqtFjssNvFDvsNig5Di8WKFosNyQlReHLRTOfPgrVZ3V8DEDlGgYj6CkMD+ZU2LhJKhYCYqDAYjC2wt59GcU1dYwviY8LRahMx7oYU5002mG+K/rq5c4wCEfUVhgbyC8e37KiIUGhCVFApBQwbEoczFz0v8BSiViIiPAQWqx1pyXFu5wjWqYL+urlzUCgR9RWGBvK5Cp0eT60uRIWuBjmZKbCLEqrrTKiuM3V6TKvVjstXDcgakYQGYwueWfsJGowt2LH/NGobTPjtz+4O4BX4Bm/uRDTQMDSQzxUWlaJCV4OM9ATo65vQarVBIaDDAk/t2UUJl6oacKGyHkqFgIyhWthFEUdPXUKFTh+UrQ1ERAMJQwP5nGuz/IXKOqzbtA8XKuvclpHujKHJjFHpCRg/OhVGkxlX9KE4d7EWT60udC4pTUREfYNTLsnnXNcFmDU5Ey8tz0f64DioVQoouph66XgeRUiICiPTEnDw2Hk0NpkRFalBha6GUy/b4dRLIgo0tjSQXznGNxw+oYPV1nGqpStJalujobXVhjMXazBnWhZio8Jwa+4o7C4uR22DyefdFFf0Bry7/STy88YGXSsGp14SUaCxpYH8av2HxfiiVAelQoBS2f3HTRQlVFyswaeffY2RaQlYPD8XX5VVAgC27z/t89aG4uMXsHHr0aBsxcjPy8bCuyZy6iURBQxbGsjvrDYRdrsItUoJOYtJSxIQqlEjPy/b+W167rQsv9wgc8cNRbNNFZQ33u5mZwT7lFUi6n8YGsjnHDernMwUAMDsKTfg+DeVuFrbJPscg2LC3c7hrxvfYG0MlhWMgkIx8Brd2H1BRL7m1/9T1tfXo6CgADExMYiJiUFBQQEaGho63d9qteKpp57C2LFjERERgZSUFCxatAiVlZX+LCb5mONmtW7TPmzffxppyXG467tjMHzIoG6PVSoEDB0ci4z0BLy6aR92F5fzYUs9xO4LIvI1v7Y0LFy4EJcuXcK2bdsAAD/5yU9QUFCALVu2eNy/ubkZR48exW9+8xvk5OSgvr4ey5YtQ35+Pr788kt/FpV8yHGTyslMwe7ichw9dQm6yjqEh6q7PVYCYLWLOKOrQbPFigZjS6f7urZofFVWed00w8vtduDiUkTka34LDadOncK2bdtw6NAh5ObmAgD+/ve/Y8qUKSgrK0NmZmaHY2JiYrBz5063bX/6058wadIk6HQ6pKen+6u45EOuN6uvyipRWW1AeFgILl6p7/I4hQCEadSo0jdCAqBSKhAbFdbp/o4Wjf0lZ3Ghsu3c18NNkt0ORNRX/BYaDh48iJiYGGdgAIDJkycjJiYGBw4c8BgaPDEYDBAEAbGxsX4qKflLhU6P2gYTpk8YAaPJjIbGFhiazJ3uL0qAqaXV+VoQBGQNT3Q7n+s3bNcWDUdLQ3/k6wGJfGAVEfUVv4WGqqoqJCYmdtiemJiIqqoqWecwm81YsWIFFi5ciOjoaI/7WCwWWCwW5+vGxkYAgCiKEMWu1wXob0RRhCRJQVduhzMXa7BlTynmzczGyLQEFBadwI4DpzE0JQ5Hv74Ic6sNyfER0NebIHWzpLRSIcBus2Pj/x2Bvr4J82ZmY8ueUmzcehSAhGUFMzEiNR7LCm4BAORNGgUAXtVdoOq7sOiEW7l7y/W6g+mzEuyf72DEOg+sYK5vuYPBvQ4NK1euxHPPPdflPl988QWAtm+K7UmS5HF7e1arFQsWLIAoili3bl2n+61atcpjefR6Pczmzr/V9keiKMJgMECSpKAczb9r33F8cewMwlU2NI0bCrvFhPtnZSI1KQaWZiNazFZER2iQHKdGV/9NhWnUSEmKhiAAyfHR+OJYGcJVNuRNGIpwlQ2544agurq61+UNVH3nTRji03LLdUVvQPHxC8gdNxSDtTEBe9/OBPvnOxixzgMrmOs7OTlZ1n5eh4bHH38cCxYs6HKfYcOG4fjx47h69WqHn+n1eiQlJXV5vNVqxf33349z585h9+7dnbYyAMDTTz+N5cuXO183NjYiLS0NWq22y+P6I1EUIQgCtFpt0H3gAGDW9PFotqkwa3pbq8B7u8owZ2omiveUo+JyE0RJguVyEyxWW7ctDaHn6nFzdjq+f0c2omMrkTEiBf/a/U3b+0THITExodflDVR9JyYmImfMKL+dvzPvbj+JjVuPo9mmwrKCwL9/e8H++Q5GrPPAuh7q2+vQkJCQgISE7v+HPWXKFBgMBhw+fBiTJk0CABQXF8NgMGDq1KmdHucIDOXl5SgqKkJ8fHyX76PRaKDRaDpsVygUQflLEwQhaMs+amgilj90KwAgP28sAAG1DSaUX6jBsNQEjEiNx7FTl3C+sg5KhQLmVlun57LaJZw6cxVFhysQHxuBosMVeG/bMQBAfGykzwYABnN9u46VANBh3ITjd5Cfl91vri+Y6ztYsc4Da6DXt9/GNIwePRq33347fvzjH+Ovf/0rgLYpl3fffbfbIMisrCysWrUK8+fPh81mw/e//30cPXoUn3zyCex2u3P8w6BBgxASEuKv4pKPOWZQVOj0iI+NQG2DCZ989jXsdjskCV0GBgCwWu1otdrRYGzB9v2nMXdaFu6//SYAHADo4DqLAkCHGRWccklEvubXdRreeecdPPHEE5gzZw4AID8/H3/+85/d9ikrK4PBYAAAXLp0CYWFhQCA8ePHu+1XVFSEmTNn+rO45Aeu4eH0uas48vVF2GUOEjI0mfHZFxXQDorErbmjMGuyvBk33gjmB1Z5mkXBQEVE/uTX0DBo0CC8/fbbXe4juXRuDxs2zO01DRwZ6VosXTAdL/x1B3SVdTA2t3rcT0DbAk8AoFYpoK9rQq2hGbuLyzE0ZZDPn6XQ9sCq4wCEoPtW3r4lIdjKT0TBh8+eoID5qqwSl6oa0NJF14RrZAwPVUMTokZTswUNxhY8tboQFboaAL67QQbzA6u6wodVEZE/DMyRGtTvOBZ6mjJ+GAZFh0PGrFuYLXbc8p0MjB6ZjJNnqnD67FVkpCf49Abf9sCqmQPuxuoY7xCMj/wmov6LLQ0UEOs/LMZ720qQOCgS9cYWhIao0GLpejCkBAkf7TrhHAMxekQSXlqeDwBYvaFowH6L9kUrAVeNJCJ/YGiggLHa7Dh7qQZ2EbDK2L/Vau+wLSNdi9UbivDup0dQ22BCfGzEgAsPvni2BGdOEJE/MDSQX7T/trx4fi52HDjtfLCUIACSBCgUAkRR/uDXCp3e+e25tsHkvLnm52X3iz787tZOkIOtBETUX3FMA/lF+z71jHQtfvfkPIzJSMao9ASMTEuAWqWAQsbYBqAtXFTXGlFYVOr8Fn1r7igMTYlDTmZKv+nDdy1HT8vkuL72QaNCp8fqDUWo0Ol9WWQiItnY0kB+4enb8qzJmc61FnYdKsPPVr2POkMLVErAZu+6tUGSJCTGR+HMxRo8s/YTLJ6fi6/KKnGhst75hMvaBhNqG0yo0On7rLXBn2sn8JHYRNTXGBrIL7rqU6/Q6bFu0z40NllgF8Vun0OhCVFCFCW0ttpQWFQKpULA6XNXsXTBdABwNv/Hx0bg3U+PID42wi83VTkDFP25dgK7LYior7F7ggKusKgUFboaDIoNh0qp6HL6pVIpIDsjGVabiHJdDaxWGwQBOH32qrOFoe18bWMdFt41sdubqqOZ/8zFGq/LLbe7oauuhPY/k9PtwHUXiKg/YGiggKnQ6fHM2k9w7PQlpCRG484ZoxEXHY7khKhOg4NaqUSFrtb5WpQAUQTCw0I6jGXobCxAe45jtuzxbqyB3FACtE0xfXXTPqz/sLjT93eEj87CiGuY6C9jNojo+sbuCQqYwqJSvLetBOZWG0JDVKipN0Ff39T2RDgB8DSswdxqc3u4lUopQKEALlTWY93GvVj6wxnYX3IWOZkpssvhuOnPmzkGQNsaED3penBof2yFTo+jpy7B3smskPbdDJ11OzjWtqhtMGHx/FznPmx1IKK+wtBAAeMYrNhgbEFsVBiyhifijQ+LUV1rRKPJDHtrx3UZ2pMkwNTStsrDoRM6RIQfxuETOkjSXqzbtA9LF0zv9sFWjpu/KIqorq4G0LtBhu2PLSwqRWW1ATeNHuK82Xt6f8Bz4HCdsmm12fHp518ja3hip+9HRBQoDA0UMBnpWvz2Z3e7bdPXm/CHN3Z5/FbuWMvBlShKUCoE2EUJol3E3iNn0Gq1oeTUJTSbbWg2t/boaZi9GWTYVcuBp5YA12DgKXA4Xt+aOwrvbS9BVY0R//PO58664IBIIuorDA3Up/LzsrHx0y9x4UpDh595mlUhAc6AYRclmFqsUCoEhKhVaDZbcf5SHf5ReBj6epNXzfe9WUGx/bHdncs1GHQVOAqLShGiUiItORZPPPhdt2tiCwMR9QWGBupTGelazJk2Gm9+VAyrTezROeyiBEOTGYCAOkMzVq7bBkmCcyyAp/7/tm/7J5A3YQiMFgXe/OgwAGDx/Fy/jRNwtDBo4yIQFaHBlj2lyMlMwfKH8pyDHvPzsp2BwNGdAwBTxg/n+AUi6nMMDdTnFs/PRYOxBTv2n7528/eeUiEgIS4CdruExiaz8xHbjm/17Z9TUVhUio1bjyJcZUOz7TLe21YCAF2u8dDbAYiOsgxNiUOFrgatVhvWbdqHWZMzPY5TcF17wlE2Dn4kor7E0EB9LiNdi5FpCQgPC4FSKaDO0NJhHwFAV2tAWW0iqmubcN/sHBhNZhz/phJZwxORkhiDnQfLUHS4HJXVBuw8eBp/+vX3kZOZggPHzmJ46iAoQyKRMTQBI1ITkJOZgmfWfuIcrOna8tC+W6H9AEbH9MrOWiscXQ85mSmQpL34qqwSmcO0+M8X/onDx8/DbLFCGxeBCp0ef3xrD0pOXYIA4MaRSTh66hIqqw0AOPiRiPoOQwP1C6431Bf+ugNfV1S5/VzOI63sonRtQGQrqvSNWLluG5Ljo3C+sg52uwgJwPGySqz/sBjxsRG4UFmPc5fqcLHmCiou1GDC6DR8VVbpnBaqVilx9NQlTBidisXzczuMN2g/gNG1taKrB2gNTRkEQRDQ3GLF1r2ncLW2CXZ7W9fMR7tPQF9vQmFRKSzXppq22uxobml1rk1BRNRXGBqoX3Ad3Dc0ZRDuXvo31BmavT7Phco6iKIECYDRZEFTswWCIECU2mZj2Owi9hwux/P/dSfOXqzB4VIdTpytR4vFigZjCxbPz3VOCz17qRZl566i4oLe2W3hGH9w5mINtIMinDdx1/EH+XnZHdZYWP9hMY6eugRdZR12HiyDNi4C40cPQZhGhSp9IwQB0KhVOKPT40JlHZLiI2E2W2Eyt2L2lBtQdl6PCl0Nviqr7NHsECIiX2BooH4nI12LZ34yG69sKEJVTSNEL8ZHKgQBVpcD1Col7HYRCgEQFALsdgnluhps+PgwTp+7iphQCZVXDVAolQDaFlS6WFUPfb0JN2Wl4IregHE3pLhNbywsKsWO/acBALuLy7G7uByAe7dEg7EFzRYrDhw7BwDYuPUoLFYb4qJCceKbSqiUCkwck4bahiY4ZpuaW22o1BuvXQcAoe2x4SWnK/GnX3/Pbf0GIqK+wNBA/YpjsGFtgwmNTRavAgMAWKzuC0S1ur52WXJy75EzgCQhNjUSNwzXYur4kbhYVY/3d34FUZQgAPjmfDWami2w2UT88a09uKI34OylWhTMuxn3334TTp2twqb/O+p8jwZjC0amJSA/LxuxUWEQAJRf0CMtORYhagVMzXbUNDRDABAbFYYKXY2zdQJoW+3SZpegENqWy3bMOQ3TqJwtMa6zLDggkogCjc+eoH7FdaxA5vBEqFUKKBUKhIb4Nt+aWqwwt9ogisDVmiZ8vPsEdheXQ7z2tV+hVMBqtUGSgOq6Jny06zj2l5xHZXUj1vzjc2QNT8SJ8iswmiyw2uyw2WzY+vlJvP7+QWdrBdAWWg4eO4+GxhZIaFucyi5KqK5rQpPJjLxJIxGmUUOlUjgfDy5JbU/2dGixWDvUD59BQUR9gS0N1G9U6PSobTAhd9xQNBhbMCI1HhNGpyJreCKef2272zMofEGS2roBWq02NJpa3X5mt4uwuzRSKAQBEWFtC0hZWm345SuFEOC4yUsQJaC5pRUxUaE4euoSjpdddq470dk0UmNzK3YXn4EmRIkWi3uTiihKGJWegOioUPziP251budqkETUlxgaqN8oLCrF9v2noR0UgdLyKqhVSjy+cAbe+PAQDEYzlAoBQ5JicPmqodOHQXlLlIBms7Xb/SxWOyxWOzQhKlhabc4WCaAtfGjjItBstkKlVKDigt7japaeWG12WG3fppPI8BCEqFUwW6yYOWlUh2W3uRokEfUldk9Qv+F49PSI1ASoVQpkDk+89o267bnZEoC500ZjyvhhPn1f0YsAogCgULg/x1upEGButcHU0grdlQbY7CJsdu8GYwiCAAFAi9kKs8WKiWPSPD7sioioLzE0UL/h+Bb95KKZeHzhd/GnX38PGelaPDw/F5ERIdCo2xrGfv/ze/CDueMRExka8DK2tNo6PBRDFCUYTRbna5VSgCC0P7JrCgFQqdvGMQxPHYSXludzoCMR9TsMDdTvOMKD46aprzchKjwUN2enOac1/uW/f4D7b7+pbRChsvuPcYha6Ta4sDfaN0y0b6cwNFmcuUKllJce7KIEq9UOCAKmjh8BAFi9oQgVOn0vS0tE5Dsc00D9TvtnPHT2qOnF83MRHxuBMxdrUFhUClEUr638KECpAEQREK/dvW020flvf/G01LXN3vE9B2ujEBGmwZmLNW6NFmEaFe6eme18yFb7Z1EQEfU1hgbqd9rfMDsb/Oe6dkFsVBgajN8+s+LspVqcvdj2UKjmFitESUKIWglRlLwebyCHQiEAkiRrAGRNQzPMFhviosNRZ2iGSikgIlyD/7dkLgryJwHgLAki6p8YGqjf6e6G2b4lIiNd65xlUKHT46nVhdBV1iF71GAsXTAdH/z7OM5eqsXCOyfg0PELbs91cBUXHYamZkuHR3QLwrfDGMI0arTa7M5nRTjIGUypUiralrK22dFstiI0RIWUxGhneNDXm9yuravnVxAR9QWGBup3uptW2FXTfWFRKSp0NUhPiUPW8CQMTRmEv/z3D5w/P32uGkqFgMRBEYiNCoNK+e10x5TEGPz9uQV4/tVtOHryItKSY/HNhRoYTRYoBECpUkKpFBATqkFTc6v7apO4NsfDJWCoVQooBAEWqx1KBZCTNQSSJKK0vApKQQAEwGq1Y+r44QCAMxdr8F8v/guV1Y3Oc7KLgoj6E4YGCjpdtUQ4ttU2mLB9/2nng6ZcqVVK3HPrOPzHvZPwzkefQ6e3wNBkxo0jk5GRrsVbqwoAtA1E1F05iDCNGgBgsdhgttrQ1NyKiDA1bDY7xGsLRIWFqjFYG4PqWiMsVjskqW3lx9xx6ahvbIGltW0q5TcX9LBa7QiLUiNEpUTWiCTnGIaNW4+i1WrD2HbPumAXBRH1FwwNFHS6aolwHefgeES1a5O/Y+2DBmML3vyoGGq0LfU8LjMFTy6a6TyP40mW6SmDoI2LwL6jZxEeGoJWmx1RERrnTAeFAIwekYSp44fjk8++Rk7WECxdMB0vv1mEsnNXMVgbA0EQUHLqMgDghqFt3QxpybFIS45zzgbJz8vGzoNlKDt3FRNGpzq7I9jCQET9CUMDDUiuweKZtZ84H1P925/djfjYCLy3rQQKAci7ORVKhYAJo9OcN2rHuIiSU5fbFm6yWGG1iYiO1MDSasOcaVmIjQrDgWNncaGyAVPHD3fO5HAdf7Bu075r56tB5nAtJoxOc7YqbPj4MEwtrW7l5ZMsiai/Y2ig605OZgoyhiZgZGoC5s0YjhtGDkV+3lgA3waG02evOm/0DcYWXKpqQEa6FoO1MbjvtnGYNTmzw4BM11aBr8oqcaGyHlnDk/DQPZPcwkR+Xjb2l5xFha5tqqjjOC4RTUT9HUMDDXiurQBA2w1dX2fC7CmZyMlMxewZiVAo2haIcgykzBqR5FyVsUKnx8i0BNQ2mPDJZ19j3aZ9GJoyqNObvOPBW3OnZTm7H1xlpGvx0vJ8tioQUdBhaKABr/3N3XGjnjdzDAD3qZOeFpJyHSdx+tzVDi0E7TkevLXwromdTpVkqwIRBSOGBrruOG7Yoiiiurra488Az+tByGkh4MJMRDRQ8dkTNOBV6PQ9eo6DYz2IwqJS57b2z8XwRM4+RETBiC0NNOD19DkO/moxaN+CQUQULBgaaMDr6c3fX+MO+DAqIgpWDA004Pn65t/blgKOeSCiYMUxDURe8jTWwRsc80BEwYotDUReYksBEV2vGBqIvMQ1FojoesXuCSIiIpKFoYGIiIhkYWggIiIiWRgaiIiISBaGBiIiIpKFoYGIiIhkYWggIiIiWRgaiIiISBaGBiIiIpKFoYGIiIhkYWggIiIiWRgaiIiISBaGBiIiIpKFoYGIiIhkYWggIiIiWRgaiIiISBaGBiIiIpKFoYGIiIhkYWggIiIiWRgaiIiISBaGBiIiIpKFoYGIiIhkYWggIiIiWRgaiIiISBaGBiIiIpKFoYGIiIhkYWggIiIiWRgaiIiISBaGBiIiIpKFoYGIiIhkYWggIiIiWRgaiIiISBaGBiIiIpKFoYGIiIhkYWggIiIiWRgaiIiISBaGBiIiIpLFr6Ghvr4eBQUFiImJQUxMDAoKCtDQ0CD7+MceewyCIGDNmjV+KyMRERHJ49fQsHDhQhw7dgzbtm3Dtm3bcOzYMRQUFMg69qOPPkJxcTFSUlL8WUQiIiKSSeWvE586dQrbtm3DoUOHkJubCwD4+9//jilTpqCsrAyZmZmdHnv58mU8/vjj2L59O+666y5/FZGIiIi84LfQcPDgQcTExDgDAwBMnjwZMTExOHDgQKehQRRFFBQU4Je//CXGjBnT7ftYLBZYLBbn68bGRud5RFHs5VUEliiKkCQp6ModrFjfgcX6DjzWeWAFc30rFPI6HvwWGqqqqpCYmNhhe2JiIqqqqjo97qWXXoJKpcITTzwh631WrVqF5557rsN2vV4Ps9ksv8D9gCiKMBgMkCRJ9i+Qeo71HVis78BjnQdWMNd3cnKyrP28Dg0rV670eJN29cUXXwAABEHo8DNJkjxuB4AjR45g7dq1OHr0aKf7tPf0009j+fLlzteNjY1IS0uDVqtFdHS0rHP0F6IoQhAEaLXaoPvABSPWd2CxvgOPdR5Y10N9ex0aHn/8cSxYsKDLfYYNG4bjx4/j6tWrHX6m1+uRlJTk8bi9e/eiuroa6enpzm12ux0///nPsWbNGpw/f77DMRqNBhqNpsN2hUIRlL80QRCCtuzBiPUdWKzvwGOdB9ZAr2+vQ0NCQgISEhK63W/KlCkwGAw4fPgwJk2aBAAoLi6GwWDA1KlTPR5TUFCA2267zW3b3LlzUVBQgMWLF3tbVCIiIvIhv41pGD16NG6//Xb8+Mc/xl//+lcAwE9+8hPcfffdboMgs7KysGrVKsyfPx/x8fGIj493O49arUZycnKXsy2IiIjI//zafvLOO+9g7NixmDNnDubMmYNx48bhH//4h9s+ZWVlMBgM/iwGERER+YDfWhoAYNCgQXj77be73EeSpC5/7mkcAxEREQXewBypQURERD7H0EBERESyMDQQERGRLAwNREREJAtDAxEREcnC0EBERESyMDQQERGRLAwNREREJAtDAxEREcnC0EBERESyMDQQERGRLAwNREREJAtDAxEREcnC0EBERESyMDQQERGRLAwNREREJAtDAxEREcnC0EBERESyMDQQERGRLAwNREREJAtDAxEREcnC0EBERESyMDQQERGRLAwNREREJAtDAxEREcnC0EBERESyMDQQERGRLAwNREREJAtDAxEREcnC0EBERESyMDQQERGRLAwNREREJAtDAxEREcnC0EBERESyMDQQERGRLAwNREREJAtDAxEREcnC0EBERESyqPq6AL4mSRIAoLGxsY9L4j1RFGE0GhEaGgqFgnnO31jfgcX6DjzWeWAFe31HRUVBEIQu9xlwocFoNAIA0tLS+rgkREREwcNgMCA6OrrLfQTJ8dV8gBBFEZWVlbISU3/T2NiItLQ0XLx4sdtfHPUe6zuwWN+BxzoPrGCv7+uypUGhUCA1NbWvi9Er0dHRQfmBC1as78BifQce6zywBnJ9B1+nCxEREfUJhgYiIiKShaGhH9FoNHj22Weh0Wj6uijXBdZ3YLG+A491HljXQ30PuIGQRERE5B9saSAiIiJZGBqIiIhIFoYGIiIikoWhgYiIiGRhaOhD9fX1KCgoQExMDGJiYlBQUICGhgbZxz/22GMQBAFr1qzxWxkHGm/r3Gq14qmnnsLYsWMRERGBlJQULFq0CJWVlYErdBBZt24dhg8fjtDQUEycOBF79+7tcv/PPvsMEydORGhoKEaMGIHXXnstQCUdOLyp8w8++ACzZ8+GVqtFdHQ0pkyZgu3btwewtMHP28+4w/79+6FSqTB+/Hj/FtDPGBr60MKFC3Hs2DFs27YN27Ztw7Fjx1BQUCDr2I8++gjFxcVISUnxcykHFm/rvLm5GUePHsVvfvMbHD16FB988AG++eYb5OfnB7DUwWHz5s1YtmwZnnnmGZSUlGDGjBm44447oNPpPO5/7tw53HnnnZgxYwZKSkrw61//Gk888QTef//9AJc8eHlb559//jlmz56NrVu34siRI8jLy8O8efNQUlIS4JIHJ2/r28FgMGDRokWYNWtWgErqRxL1iZMnT0oApEOHDjm3HTx4UAIgnT59ustjL126JA0ZMkQqLS2Vhg4dKv3xj3/0c2kHht7UuavDhw9LAKQLFy74o5hBa9KkSdKSJUvctmVlZUkrVqzwuP+vfvUrKSsry23bY489Jk2ePNlvZRxovK1zT2688Ubpueee83XRBqSe1vcDDzwg/fd//7f07LPPSjk5OX4sof+xpaGPHDx4EDExMcjNzXVumzx5MmJiYnDgwIFOjxNFEQUFBfjlL3+JMWPGBKKoA0ZP67w9g8EAQRAQGxvrh1IGp9bWVhw5cgRz5sxx2z5nzpxO6/bgwYMd9p87dy6+/PJLWK1Wv5V1oOhJnbfneJTzoEGD/FHEAaWn9b1+/XqcOXMGzz77rL+LGBAD7oFVwaKqqgqJiYkdticmJqKqqqrT41566SWoVCo88cQT/izegNTTOndlNpuxYsUKLFy4cMA+kKYnampqYLfbkZSU5LY9KSmp07qtqqryuL/NZkNNTQ0GDx7st/IOBD2p8/ZeeeUVmEwm3H///f4o4oDSk/ouLy/HihUrsHfvXqhUA+N2y5YGH1u5ciUEQejyz5dffgkAHh9BKklSp48mPXLkCNauXYs333wz6B777U/+rHNXVqsVCxYsgCiKWLdunc+vYyBoX4/d1a2n/T1tp855W+cOGzduxMqVK7F582aPYZo8k1vfdrsdCxcuxHPPPYcbbrghUMXzu4ERffqRxx9/HAsWLOhyn2HDhuH48eO4evVqh5/p9foOSdZh7969qK6uRnp6unOb3W7Hz3/+c6xZswbnz5/vVdmDlT/r3MFqteL+++/HuXPnsHv3brYytJOQkAClUtnhG1d1dXWndZucnOxxf5VKhfj4eL+VdaDoSZ07bN68GY888gj++c9/4rbbbvNnMQcMb+vbaDTiyy+/RElJCR5//HEAbd1BkiRBpVJhx44duPXWWwNSdl9iaPCxhIQEJCQkdLvflClTYDAYcPjwYUyaNAkAUFxcDIPBgKlTp3o8pqCgoMN/4HPnzkVBQQEWL17c+8IHKX/WOfBtYCgvL0dRURFvaB6EhIRg4sSJ2LlzJ+bPn+/cvnPnTtxzzz0ej5kyZQq2bNnitm3Hjh24+eaboVar/VregaAndQ60tTA8/PDD2LhxI+66665AFHVA8La+o6OjceLECbdt69atw+7du/Gvf/0Lw4cP93uZ/aIPB2Fe926//XZp3Lhx0sGDB6WDBw9KY8eOle6++263fTIzM6UPPvig03Nw9oR3vK1zq9Uq5efnS6mpqdKxY8ekK1euOP9YLJa+uIR+a9OmTZJarZZef/116eTJk9KyZcukiIgI6fz585IkSdKKFSukgoIC5/5nz56VwsPDpSeffFI6efKk9Prrr0tqtVr617/+1VeXEHS8rfN3331XUqlU0l/+8he3z3JDQ0NfXUJQ8ba+2xsIsycYGvpQbW2t9OCDD0pRUVFSVFSU9OCDD0r19fVu+wCQ1q9f3+k5GBq8422dnzt3TgLg8U9RUVHAy9/f/eUvf5GGDh0qhYSESBMmTJA+++wz588eeugh6ZZbbnHbf8+ePdJNN90khYSESMOGDZNeffXVAJc4+HlT57fccovHz/JDDz0U+IIHKW8/464GQmjgo7GJiIhIFs6eICIiIlkYGoiIiEgWhgYiIiKShaGBiIiIZGFoICIiIlkYGoiIiEgWhgYiIiKShaGBiIiIZGFoICIiIlkYGoiIiEgWhgYiIiKShaGBiIiIZPn/uDbJuB+wAqsAAAAASUVORK5CYII=", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "# Set the number of posterior draws you want to get\n", + "num_samples = 5000\n", + "\n", + "# Obtain samples from amortized posterior\n", + "conditions = {\"observables\": np.array([[0.0, 0.0]]).astype(\"float32\")}\n", + "samples_at_origin = approximator.sample(conditions=conditions, num_samples=num_samples)[\"parameters\"]\n", + "\n", + "# Prepare figure\n", + "f, axes = plt.subplots(1, figsize=(6, 6))\n", + "\n", + "# Plot samples\n", + "axes.scatter(samples_at_origin[0, :, 0], samples_at_origin[0, :, 1], color=\"#153c7a\", alpha=0.75, s=0.5)\n", + "sns.despine(ax=axes)\n", + "axes.set_title(r\"Posterior samples at origin $x=(0, 0)$\")\n", + "axes.grid(alpha=0.3)\n", + "axes.set_aspect(\"equal\", adjustable=\"box\")\n", + "axes.set_xlim([-0.5, 0.5])\n", + "_ = axes.set_ylim([-0.5, 0.5])" + ] + }, + { + "cell_type": "markdown", + "id": "01821d24", + "metadata": {}, + "source": [ + "The posterior looks as we have expected in this case. However, in general, we do not know how the posterior is supposed to look like for any specific dataset. As such, we need diagnostics that validate the correctness of the inferred posterior. One such diagnostic is simulation-based calibration(SBC), which we can apply for free due to amortization. For more details on SBC and diagnostic plots, see:\n", + "\n", + "1. Talts, S., Betancourt, M., Simpson, D., Vehtari, A., & Gelman, A. (2018). Validating Bayesian inference algorithms with simulation-based calibration. *arXiv preprint*.\n", + "2. Säilynoja, T., Bürkner, P. C., & Vehtari, A. (2022). Graphical test for discrete uniformity and its applications in goodness-of-fit evaluation and multiple sample comparison. *Statistics and Computing*." + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.11.10" + }, + "toc": { + "base_numbering": 1, + "nav_menu": {}, + "number_sections": true, + "sideBar": true, + "skip_h1_title": true, + "title_cell": "Table of Contents", + "title_sidebar": "Contents", + "toc_cell": true, + "toc_position": { + "height": "calc(100% - 180px)", + "left": "10px", + "top": "150px", + "width": "165px" + }, + "toc_section_display": true, + "toc_window_display": true + } + }, + "nbformat": 4, + "nbformat_minor": 5 +}