diff --git a/.gitignore b/.gitignore index 1e83934b..3249a8c8 100644 --- a/.gitignore +++ b/.gitignore @@ -78,6 +78,11 @@ target/ # Jupyter Notebook .ipynb_checkpoints +# WIP PDE tutorial notebooks -- kept locally, not tracked/pushed +/docs/tutorials/pde/pde_inference_setup_exponax_ks.ipynb +/docs/tutorials/pde/pde_inference_setup_phiflow.ipynb +/docs/tutorials/pde/grid_interpolator_examples.ipynb + # IPython profile_default/ ipython_config.py diff --git "a/docs/tutorials/pde/Kuramoto\342\200\223Sivashinsky_initial_condition.ipynb" "b/docs/tutorials/pde/Kuramoto\342\200\223Sivashinsky_initial_condition.ipynb" new file mode 100644 index 00000000..b1f0bbda --- /dev/null +++ "b/docs/tutorials/pde/Kuramoto\342\200\223Sivashinsky_initial_condition.ipynb" @@ -0,0 +1,532 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "721e3a43", + "metadata": {}, + "source": [ + "# Minimal Example: Sparse-Observation Forward Simulation + GP Initial-Condition Inference\n", + "\n", + "A trimmed-down version of `pde_inference_setup_exponax_ks.ipynb`, keeping only the two\n", + "patterns needed for this pairing:\n", + "\n", + "1. Building a `dynestyx.DynamicalModel` whose observation operator is a genuinely sparse\n", + " matrix (`jax.experimental.sparse.BCOO`), and running a forward simulation with it.\n", + "2. Treating the equation's initial condition as an unknown, Gaussian-process-distributed\n", + " latent, and recovering it from noisy observations with NUTS.\n", + "\n", + "See the full notebook for the exploratory detours this one skips: dense-vs-sparse $H$\n", + "comparisons, inferring the equation's $\\lambda$ parameter, `RegularGridInterpolator`, and\n", + "EnKF sanity checks against a naive baseline.\n" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "id": "5cfd0bd8", + "metadata": { + "execution": { + "iopub.execute_input": "2026-08-11T21:53:23.262123Z", + "iopub.status.busy": "2026-08-11T21:53:23.261975Z", + "iopub.status.idle": "2026-08-11T21:53:23.475799Z", + "shell.execute_reply": "2026-08-11T21:53:23.475314Z" + } + }, + "outputs": [], + "source": [ + "import jax\n", + "jax.config.update(\"jax_enable_x64\", True)\n" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "760d7f0d", + "metadata": { + "execution": { + "iopub.execute_input": "2026-08-11T21:53:23.477209Z", + "iopub.status.busy": "2026-08-11T21:53:23.477112Z", + "iopub.status.idle": "2026-08-11T21:53:24.477838Z", + "shell.execute_reply": "2026-08-11T21:53:24.477557Z" + } + }, + "outputs": [], + "source": [ + "import equinox as eqx\n", + "import jax.numpy as jnp\n", + "import jax.random as jr\n", + "import matplotlib.pyplot as plt\n", + "import numpyro\n", + "import numpyro.distributions as dist\n", + "from jax.experimental import sparse as jax_sparse\n", + "from numpyro.infer import MCMC, NUTS, init_to_value\n", + "\n", + "import exponax as ex\n", + "\n", + "import dynestyx as dsx\n", + "from dynestyx import DynamicalModel, Filter, LinearGaussianObservation\n", + "from dynestyx.observation import GridInterpolator\n", + "from dynestyx.inference.configs.filter import EnKFConfig\n" + ] + }, + { + "cell_type": "markdown", + "id": "b885728a", + "metadata": {}, + "source": [ + "## The equation and its (fixed, known) dynamics\n", + "\n", + "Conservative Kuramoto-Sivashinsky, as in the full notebook: a discrete-time transition of\n", + "`n_substeps` fixed ETDRK2 steps, wrapped as a bare `sample()`/`shape()` object (no\n", + "`log_prob` needed -- `EnKFConfig` only ever calls `.sample(key)` on it).\n" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "id": "c3d5036c", + "metadata": { + "execution": { + "iopub.execute_input": "2026-08-11T21:53:24.479415Z", + "iopub.status.busy": "2026-08-11T21:53:24.479251Z", + "iopub.status.idle": "2026-08-11T21:53:24.913764Z", + "shell.execute_reply": "2026-08-11T21:53:24.913422Z" + } + }, + "outputs": [], + "source": [ + "num_points = 64 # spectral resolution\n", + "domain_extent = 60.0 # periodic domain length\n", + "dt_ks = 0.5 # Exponax's own internal step size\n", + "n_substeps = 2 # ETDRK2 steps per dynestyx transition\n", + "\n", + "stepper = ex.stepper.KuramotoSivashinskyConservative(\n", + " num_spatial_dims=1, domain_extent=domain_extent, num_points=num_points, dt=dt_ks\n", + ")\n", + "\n", + "state_dim_ks = num_points\n", + "control_dim_ks = 0\n", + "x_grid = jnp.linspace(0, domain_extent, num_points, endpoint=False)\n", + "\n", + "\n", + "class KSSample(eqx.Module):\n", + " \"\"\"Bare sample()/shape() object required by EnKFConfig -- no log_prob needed.\"\"\"\n", + "\n", + " loc: jax.Array\n", + " noise_std: float = eqx.field(static=True)\n", + "\n", + " def sample(self, key):\n", + " return self.loc + self.noise_std * jr.normal(key, self.loc.shape)\n", + "\n", + " def shape(self):\n", + " return self.loc.shape\n", + "\n", + "\n", + "class KSTransition(eqx.Module):\n", + " \"\"\"Discrete-time state_evolution: n_substeps fixed ETDRK2 steps of size dt_ks.\"\"\"\n", + "\n", + " n_substeps: int = eqx.field(static=True)\n", + " process_noise_std: float = eqx.field(static=True)\n", + "\n", + " def __call__(self, x, u, t_now, t_next):\n", + " v = x.reshape(1, num_points)\n", + " for _ in range(self.n_substeps):\n", + " v = stepper(v)\n", + " return KSSample(loc=v.reshape(num_points), noise_std=self.process_noise_std)\n", + "\n", + "\n", + "state_evolution = KSTransition(n_substeps=n_substeps, process_noise_std=1e-3)\n" + ] + }, + { + "cell_type": "markdown", + "id": "fde9d38f", + "metadata": {}, + "source": [ + "## A sparse observation model\n", + "\n", + "We observe $u$ at points placed anywhere in the domain (not aligned to the grid).\n", + "`dynestyx.observation.GridInterpolator` already builds exactly this kind of observation\n", + "operator -- piecewise-linear (or nearest-neighbor) interpolation on a regular grid, with\n", + "periodic boundary handling built in -- so there's no need to hand-roll an interpolator and\n", + "recover its matrix via `jax.jacobian` (see `grid_interpolator_examples.ipynb` for a fuller\n", + "tour). Each row of $H$ has only 2 nonzero entries (`method=\"linear\"` in 1D) out of\n", + "`state_dim_ks = 64`, so we use `GridInterpolator.as_sparse()` -- a\n", + "`jax.experimental.sparse.BCOO` array -- directly as `LinearGaussianObservation.H`;\n", + "`EnKFConfig`, the filter used throughout this notebook, applies a sparse `H` directly.\n" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "id": "edaeaf54", + "metadata": { + "execution": { + "iopub.execute_input": "2026-08-11T21:53:24.915072Z", + "iopub.status.busy": "2026-08-11T21:53:24.915008Z", + "iopub.status.idle": "2026-08-11T21:53:25.474487Z", + "shell.execute_reply": "2026-08-11T21:53:25.474219Z" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "H_sparse nse (stored nonzeros): 64 out of 2048 dense entries\n" + ] + } + ], + "source": [ + "n_obs = 32\n", + "x_obs = jr.uniform(jr.PRNGKey(42), (n_obs,), maxval=domain_extent)\n", + "\n", + "interp = GridInterpolator((x_grid,), x_obs[:, None], method=\"linear\", boundary=\"periodic\")\n", + "H_sparse = interp.as_sparse()\n", + "print(\"H_sparse nse (stored nonzeros):\", H_sparse.nse, \" out of\", n_obs * state_dim_ks, \"dense entries\")\n", + "\n", + "observation_model = LinearGaussianObservation(H=H_sparse, R=0.02**2 * jnp.eye(n_obs))\n" + ] + }, + { + "cell_type": "markdown", + "id": "3e923486", + "metadata": {}, + "source": [ + "## A Gaussian process prior over the initial condition\n", + "\n", + "A GP prior over our finite grid is exactly a multivariate normal whose covariance is a\n", + "kernel matrix; `numpyro` needs no dedicated GP class for this (see the full notebook for\n", + "the derivation and the periodic-kernel formula). We build the kernel matrix $K$, its\n", + "Cholesky factor $L$, and use the non-centered parameterization $u_0 = \\mu + Lz$ with\n", + "$z \\sim \\mathcal{N}(0, I)$.\n" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "id": "49e42eac", + "metadata": { + "execution": { + "iopub.execute_input": "2026-08-11T21:53:25.475563Z", + "iopub.status.busy": "2026-08-11T21:53:25.475506Z", + "iopub.status.idle": "2026-08-11T21:53:25.724213Z", + "shell.execute_reply": "2026-08-11T21:53:25.723856Z" + } + }, + "outputs": [], + "source": [ + "def periodic_kernel(x1, x2, variance, lengthscale):\n", + " d = jnp.abs(x1[:, None] - x2[None, :])\n", + " d = jnp.minimum(d, domain_extent - d) # periodic (wraparound) distance\n", + " return variance * jnp.exp(-2.0 * jnp.sin(jnp.pi * d / domain_extent) ** 2 / lengthscale**2)\n", + "\n", + "\n", + "gp_variance = 0.3\n", + "gp_lengthscale = 1.5\n", + "gp_mean = jnp.zeros(state_dim_ks)\n", + "\n", + "K_gp = periodic_kernel(x_grid, x_grid, gp_variance, gp_lengthscale) + 1e-6 * jnp.eye(state_dim_ks)\n", + "L_gp = jnp.linalg.cholesky(K_gp)\n", + "\n", + "\n", + "def build_dynamics_ic(u0):\n", + " \"\"\"DynamicalModel for a given initial condition u0, with the sparse-H observation model.\"\"\"\n", + " initial_condition_u0 = dist.MultivariateNormal(\n", + " loc=u0, covariance_matrix=1e-4 * jnp.eye(state_dim_ks)\n", + " )\n", + " return DynamicalModel(\n", + " initial_condition=initial_condition_u0,\n", + " state_evolution=state_evolution,\n", + " observation_model=observation_model,\n", + " control_dim=control_dim_ks,\n", + " )\n" + ] + }, + { + "cell_type": "markdown", + "id": "4d7fd6d7", + "metadata": {}, + "source": [ + "## Forward simulation, with the sparse observation model\n", + "\n", + "One draw from the GP prior, forward-simulated through the KS dynamics. `dsx.simulate`\n", + "doesn't care that `observation_model.H` is sparse rather than dense -- nothing else about\n", + "the call changes.\n" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "id": "95f71f7e", + "metadata": { + "execution": { + "iopub.execute_input": "2026-08-11T21:53:25.725681Z", + "iopub.status.busy": "2026-08-11T21:53:25.725597Z", + "iopub.status.idle": "2026-08-11T21:53:27.163328Z", + "shell.execute_reply": "2026-08-11T21:53:27.163103Z" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "states shape: (1, 30, 64) observations shape: (1, 30, 32)\n", + "any nan: False\n" + ] + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAk4AAAGGCAYAAACNCg6xAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjExLjEsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvctoD+AAAAAlwSFlzAAAPYQAAD2EBqD+naQAAVgZJREFUeJzt3QeYFEXaB/B30ubAsruwAksQUQmHoGIGEcwoeoiigqIIh4rhTjlO9BTUU+A8RUWMp3eK6TNhDiiiYkIUUEREJC4578LmnenveQtmbma2qra7dzbM9P93zxxuz3RPT3V1zztV1W+5DMMwCAAAAADq5K77JQAAAACAwAkAAADAArQ4AQAAAJiEwAkAAADAJAROAAAAACYhcAIAAAAwCYETAAAAgEkInAAAAABMQuAEAAAAYBICJ5N+/vlncrlc9Pzzz1N9ffjhh2JbX375Zb231dz2oak/29q1a+mcc86hli1biv145JFHmmQ/nKZdu3Y0YsQIShRcd/7+97832vr1eb/67mtzqjONWW6xsHfvXiooKKCnnnqqyfYhXhUXF4vrdCy+U5t14PTaa6+Jiip7HH300Q23l9CsvPnmm+KYf/vtt9TcXH311bRhwwb69ddfiWcTuu666yjRLViwgC677DLq1KkTpaSkiIvRH/7wB7riiivo008/FeUQxOURft4mJSVRx44dady4cbRjx44m/RyJrKKiQpT35MmTm3pX4kpzL7cpU6ZQeno6XXnlleR0LVq0oFNPPVXb8BB+HLOzs+nmm2+miRMnUllZGcUTr52VXnrpJbr44otjvzcOceaZZ0Z8mSWSpv5s8+fPF8FBq1atyAn4QnTXXXfRmDFj6PXXX6fDDjuMAoEA/fDDDzRjxgwaOHAgLVy4sNYPm82bN4tfyvyL+f3336dRo0aJsuP1fD4fOV1967DV9RP1epDI5cbnDp9jd9xxB3m9tr5KHe9Pf/qTuIbNmjWLxo4dGzflga46SBilpaXil0tqaio5ATdx33nnnTRt2jR64okn6MgjjxS/fjMzM6l///4ikPrvf/+rDYT4tcOGDRMtTkuXLqVPPvmkUT8DQLx64YUXxDVn+PDhTb0rcSs/P5/OOOMMevzxxymeNEjgVF5eTrfeeisdfPDBoivgoIMOEr9o+Vdu0L59+0TT3T/+8Q/xi5cv+snJyfTvf/+bBgwYQCeccELENv/4xz+K1z/55JOhZdy14Ha76b777ovoNw92Q/AXRvv27emaa66hXbt2mXpvxl8gp512GqWlpYlf5bfffrv4FW8Gf3H/9a9/pc6dO4svcC4Dbg3g7iPdOKCXX35ZLFu8eLH4MuQy46ZMbgLm5mp+f47M27ZtK/aLyyP8M7F//etfYhvRXS5LliwRy/k9dIKvCz54/3v27EkPPvhg6Jcd/ze/Nzv++ONDr+UvaNVnY6tWrRKtlHyicFlzywg3c/v9/lplwPvxz3/+UxxL3odTTjmFfvnlF+2+c9lkZGSI/+byC+5XTU1NRNnec889ok5wvdmyZYt4PQcYxx13nChXDiT4/biLKxyPm+rRowdt3LiRzj33XPFe3MXFv5TY+vXrafDgwWJ9rjPhdVJn/PjxoX3lfcrJyRGtdt98802d6/Jn5s/Czd0qI0eOpCOOOKLObR1++OHi36KiIu3r9uzZI+pzbm4uZWVl0QUXXEDbtm2Tvo4/09SpU+ndd9+l3r17i2tBsJ5wGYWfpx06dBDB2+7du0Pb6Nevn3iE47LndZ555pnQsq1bt4pl06dPr/Xec+bMoV69eokuTP6MPNzAztgZq9sMX5+7jYPBfHjdDG+1l43VMVNGVvA1hMf8BQNs3v6ll15Kq1evjrh+3XLLLaLbN3jtHj16tChju2Vhts7Eqtwa4jPIvPPOO+L1bdq0qfXc008/Leo8Xw9at25NZ511Fn311VfS9+ehD927dxfvz//KrtNm6kJd51xd+xSsIw8//LC4ZnDZc/fbeeedR8uXL6eGMmDAAHHNr+va06wYFrz66qv87Wm89NJLytf4/X7j1FNPNVq0aGG8/vrrRnFxsfHNN98Yhx56qNGxY0dj586d4nV79+4V2zr33HONSy65xPj999+N3377zfj444+Ne++91/B4PMbu3bvFa2tqasT2UlNTjQsuuCD0Xi+//LLYxuLFi6X7sm/fPmPevHlGly5djDPOOCO0XPfeq1evNrKzs41+/foZv/zyi9jfhx56yBg6dKhYZ9asWdoyGj16tHHQQQcZ8+fPN8rKyoz169cbzzzzjDFx4sTQaz744AOxLX5NEJcpL7v44ouNxx9/XHz2r776ysjNzTXGjRtn3HrrrcbMmTONXbt2ifLMy8szRowYEfHe9913n9jG9u3bI5Zz+UQfN9k+RNuxY4fx5JNPGsnJycaMGTNCy2fPni3W5f2IJtvuunXrxP4eddRRYl/4sz377LPieHL5R5cBfy4uc35/Pgbdu3c3unbtKuqWTvC4Tpo0KWJ5cLsXXnih8cADD4jy4bq5detW44knnhDP/fWvfzU2bdpkrF271rjiiitE/Xv77bdD2xg0aJCow1z/+HPv2bPHuOeeewyXyyXqzVlnnSWOF9f3adOmiW1++OGHhhXV1dWiHg4fPtzIysoS+6KycuVK8R5XXXWVpffgusTrbd68OWL5hAkTxPL3339fuS6fhyeccIJRUFBgfPTRR+KzzpkzxzjvvPNEnef9DuJjzNsbPHiwOJ6rVq0yVqxYYcydO1d6nvLyzp07i3IOuuuuuwyv12uUlJSEyiczM1PUm2HDhoVe9/zzz4v3Wrp0acR7DxkyxBg1apSxZs0aY9u2bcbll18ujiuXXV14/dtuu63W5zG7zej1y8vLpXVT9XqzZWRm3SC+tqSlpYnrCNd13n8+N2666SbxPJ9f/fv3N1q2bCnOcT6+XKf5+nnwwQeHrsdWysJKnYlFuTXEZ1Dh74nLLrus1vLXXntNXBeeeuopcZ3gazZ/5vPPPz/0muD783cQb4Pfn89Jvg7x8hdeeMFyXdCdc2b2ifF6fI7xdxZ/9/G1m881/v7l78a6ymPgwIHS5/jcVB1H/q7g5/g8jhe2AifZg788GH/Z8N/8hSv78g5W8uCXXJs2bYyqqqqI1y5cuFA8xweb8RcV/83BR05OTugLlL808vPzjUAgoN1v3g6vH/wi0r33mDFjxIWZT6LogMhM4HTIIYdIT6ZwusDpuuuuq/WFlpKSUms5lwV/qZSWljZY4BTEgUSPHj1sB05/+tOfDJ/PVysQuPvuuyO2EyyDa6+9NuJ1HOTw8s8//7xegdOVV14ZsbyiokJcYAcMGBCxnOsXB2p8LIP4AsXb+Pbbb0PLuN61a9fOSE9PN7788suI5R06dBCBmh1cJ7kOTp06Vfkavhjy/kyePLlegROX2SuvvCI+Q7du3YzKykrlusHjwNcB2fklC5wKCwtFwGNG8IfQhg0bIs77N998M+ICy3WfA/HgeT9y5EjxJRz93vxlyV/cQfxlkZSUFPEjxmrgZHabsQ6cVGVkdl0O4vl1Dz/8sPI1b7zxhngNf2nKrsfBfbdSFlbqTCzKrSE+gwwHH7w+X5+j8bWav5d0gu/P14no8+Okk04S15W6vtei64LunDOzT8FrymOPPVbrOsn7U9ePtOzsbGV8EHzIjiM3XPBz/EM0XrjtDg4/EHSFHsHR9HPnzhX/DhkyJGIdbgrl7qvg80HcXBg9BoObkrlZl5tQ2ccff0yHHnqoaG7lpkke7Bpczu/LzZNBfKcXN+fz4GCPxyOeGzp0qHju999/r/O9ef+424a7lMKdf/75psqGmzi564e7mlauXGlqneh9CsdNwdxVx12H4bp27Sq6obiLKJa4u/LYY48Vzbnh3XDRZWcFl+lRRx0lmpfDBY9LdJ0YNGhQxN/cRcbCuxTs4K60cNx1x92d0XWVu8x4GX/mdevWhZZzNxqXTRCXDddLrmcnnnhixHLuijSzvzt37qQbb7xRnBvchRm80427u3VlrhsUy12I4V2u3G0Rjbsv+Dn+TNxdyHflzZs3T7w31//w9YN1n48T/83nV3S5cpnJnH322dKBs19//bU4znyeBc/TYBdM8HP36dNHdBWEXwe6desm7hbk7uhFixaJ5TwuS3Y3D3d58raD+LNyV3d96lFDbFPFTBmZxUMSGHfNqaiu3XxjAZ+70eepmbKwU2fqoyE+gwx3izG+Tsq+A7Zv3y6+r/j7iK/TKrLzg883Htrx22+/2aoLsm2a2Sfuegy/Lgfxdemkk06izz//nOoycODAWrEBP3j4iwp334aXaTyIec3lLwI+aBz4yPppo8ffcCWttVNut+j3DL9gcuDA44X4wX9zpeKgITygWLZsmRgUy33Fn332mRi4xweNx92w6urqOt+b95/7f6PJlsnwIF3+EuLAib9Uuf+b7xYIH+Okw19o4YInpmq5mcpm9s4THr/EdzlceOGFok+byyt4S3902VnBZcrHPlpwWXSdiP6ssTqxoo8371f4ftS1b9H7xYJjmmTL69pfLlu+cPMYB643PO6Dxxjwcg4YdGUeDEJl4wI4lxVvY82aNcr1ebwhv4bfg4PDxx57rM47Ebm8eNwdX0jD8Y8P/sKRkZ1jP/30kxhHxuPE+GIcPE95XAYLfm7+guDXRV8H+Lziz89/89g3HncW/cNCdby4LtWnHjXENmXMlpFZXLf4uMmuy+HHl6+dfIzNXLvNlIWdOlMfDfEZZPj8ZCUlJbWe4/G8PMaRjxuPA+XXcuAYPZ6I6b5rgvtqtS7Izjkz+xQc88nlxN/hfP7xdzE/eNxV8HoZa8EybIj6EDeBE+eQ4Wg2euAy48F5eXl5EctUd/zwhZC/ALhVgCPk008/XSznf/lCyhdNFv5L8//+7/+osrJSJCPjX6Z8AjHVF4jsvfnCEj6IMHzfzeD1+Q4BvlBx7goeuMuDDbnimwlgwlvPzCwPF7xY8G2y4fiLxYznnntO/LLgFggemB381aL7AjZbJ3RlGl0nzHxWO6KPN+9X+H7UtW/1OTYyHJx+//33dNttt4l6zMePt8UX7bou3Icccoj4EcGtLWZvXDCLz7fwX4sc2AXrNiet43MsHF+4VQOWZecYX4R5HW7drOs85esA/6LmwaPfffdd6DrAy1XXgYasRw1VN+tTRmZwSwUfN92XH58P3LrNx9jMtdtMWdipM/XREJ9Bhs9VfgSDjXAcaPA1lHsc+Acz/yjif2U3ueiuPcEg12pdkJ1zZvaJy4aDJS47/g7nG3f42hL8MSf7To+F4E1j0T0SjgqcuKmOzZ49O2L5jz/+KO6sCj5fl+AvSL47jw8atyQFl/OF/Y033hDdWIWFhbUqCHc3hAve+WQGt3Tx9qN/mbz11lumtxHcD75DggOn66+/Xlz8w+8qbAjc3cM4YAsX/GViRvQvQ97n6FvU+Y4cFn0xVOFjzgFCdOsId2kGn28K3CXMv3Ki6yrXN17G5dkYJ3N0mXMAa8akSZNEaxG3FDYGPjfCf+mGN/FbDd7MnqfB6wAnyeOL+sknnxxazr+WOajjrlxZy0Fzwp+Vv9DMnjOxuJZF3xXKXnzxReVrgudhMFAO4i5R/hFr5zytb52xWm4N8RlU+EcmX9d0uPWH0xVwiy4HP9FJgz/44IOIO4uD3zX8w5W7+xuiLqj2iVugeF+C1+XG8t1334l/+/btS44NnPgE5Sh2woQJovJyMxxnNub+WA5y/vKXv5jaDo/T4F/V3M3GYy6C3TVc8fmE49vFo5vngxcHDlY4OuYvFe4mk90uqsIXaD5RubuKb4flX0V8C6/ZX0d8YecxYPzefLJz3y6fCBzkybp0Yonfu0uXLiJ9Andl8q9L/lLlpl0zeNwBj3Ph/ed1+GLD/d3Rv+Z5fBV/ifFJz7/u6sItKtx1xdviAJp/0XAOIk4NcNFFF4nj2xQ4YOF94MCQjzv/euTgjscBcGvQ/fff36DvzxdGfnATOnczB8uF67yZBJ6XX365KFv+JXnttdeK1lkeG8XHjn9F8oUxlq0kPPaCx3jdcMMNosy4ZZP/5Qu4lcCFz1M+h/k85fOKzxUuc06tEI2vAXwt4DLhFCXBoJ3rJP8q5i55WTddc8NffHzefPHFF6a6PKyUkRncUsfXYK7njz76qKjrPOaFW+mD6Sz4tnP+8uK/ObDh48tpMS655BJxDHgsXmPXGavl1hCfQYUDDf6O2LRpU8RyHtrAQzX4HOTrI7f4cy8IBz7RaXY46L/qqqvE8eWWJh6PyIloOVVL8LyNRV0ws09cR3iYCf/Q514T/lx8LeGuQk4Hwd8rDeHTTz8VY6CjG0EcFThxRedfF5y3g4Mkbv7jL2Q+eXiAm66PPVrwghh+YeTmUR40yqK/0I855hjxpc8Vj6Nq/rXDrT6c78Is7v7gk5Q/B+e84ICHK5nZ+ZA4lxIPxOQWMu5H5s/OFZMHJTbEYMhwHMxwSxyXEedf4gefJLo8P+G4dY+/hDno5aZ9LjdOrhjM8RPEv4Y4IOOLLve7h+dxkuFWG/5VwxcuDnx525ztmi/inESuKXGOL/4cfPLysefPyq2D/EXNF+GGPl58cecLIA8uD47f41YBs8EO5yLjnFkcdHF+LW5B4x8KnCuHL+ociMXqghcMlnlcFge8/D4zZ84UF1krdZvPBz7uHKTzNvg85nNNlTlYdh3gLhluMWSqaR6aGw5k+Xzka1N0PqL6lpEZvD2uL9xNw+ciDxh+++23Q9fH4PHl8TD85cllzHWK6yZfu+2MQYlFnbFSbg3xGVS41Yavf9HXML6Oclc7f15+Px6Yzj0YHORHX0t5vzgA4/Lh6wAfDz5nw+fwi0VdMLtPzz77rPjByNdzHkvIP/b5Zgw+VrEMOoM4eP/oo4/EVFnxxMW31jX1TgAAAMQb/kHNP9ZXrFhhadoVDmI4gOGWJdldr05xzz33iCCay48TEMcLTLkCAABgAwc9nKk8PJM9mMOt5Ny6xcFjPAVNDDMTAgAA2MBddQ1900+iys7ObrA79RoaWpwAAAAATMIYJwAAAACT0OIEAAAAYBICJwAAAACTMDhcgROOcQKw4GS3AAAAjYmzBXEST87f1NB5ADlXVlVVla11k5KSQlPBOAECJwUOmuIpkykAACQmntGAEw83ZNDUqVMBbdlSe44/MwoKCsTceU4JnhA4KXBLE+uVfBF5XJGTJv5YvX+2dpm+SYOly89tGzkfUbheBfJJeLv0ipwQMlzSJfJ0+74eVyrXodnjpYsfu/NC5SpPblsvXX44dVSuM+Lg2jOGB/U9dqF0eeZ56jysgQx5tt89T6mne3n8c/m8Rz/uVh+HwjSPdPmJrSInTQ7Xq3CtdHnL1pFzHYYzDHkL5rZN6il5lm6SXzSLytT5T9I96s96cJb8GB3SNnI+wXC57WpPaMpKd6qzMf+w4n/zbYVbtqf27PVBOb4a6fITO6gnuG3XeZ10+a9LIzM1h7vrR/m0Nl9XqOcAS0tSz114gneAdPnRLdUtBd1z5F9UXVqrb3E/+MTF0uWBc/bPqCCTtHKJdPnH09Rztz20Yv/0NtHSPfKJ2dnYw+TTopw85APlOv69qdLlT/7fEOU692xQb6+qRn6L+4mp6ozjp7eWn/tHt1IfhyNOlM9T57oqclqVcIHU/ROLR9v1F/nx2VftpxM+WBL6Pmoo3NLEQdPaoocoK0t+PFRKSsqpY+GNYhsInBwu2D3HQZPHFTm5osulvhB6o14blKr5Esvwyi9EWcnq90nKkMe8vixNIrE0+Rd2ilu+z8wdFTSG3ofU66R51M9lJck/U2a6LnCSX9T8SfLlLFnxmXwu9XFIcsvLNM2jnmA0wydfJ1Ozb6rAqUyTeThVUaa6Y5eiqXNpii+/TMXn0dVHt89juS7o9jvVI3+fDE35qOpVuuLc0p2rROqued2571NsL1nTxaIqH+1nTZHvX0BxTWBJaW5L9UBXPj7FNUEXVKn2mfmr3JbriO44qIZWqI83v5fHcv1RnQ+uTPU6gVT5PlRrzjuxzUYaLpKRkSweVgQsTvCdCNDiBAAAAGQYNeJhhWHx9YkAgRMAAACQYfjFwwrD4usTAQInAAAAoIBRIx5WBNDiBAAAAE6Erjpz0OIEAAAAB7rqrI5x8juu5JA5HAAAAMAktDgBAAAAGYEa8bDCsPj6RIDACQAAAIi4m87qYG8DgRNE+Z0Wk4sik6Mle9VZjzumy5OftUlXp7JvlS/PMp1+1G7lOpWH/1W6vGL1q8p1fn38eOnytzZpEiW65Rlrj8tXrkLd2sqzjbP0TvJMvK4ydWbcwJfycvh6ydnKdbZVyBPGtUpR/1Zony4vh5yUcsvJLEuL1Z+nqlKeBG+vJgu46n0yvepjl+mrVj6X4pU/53FrxisE5D37/hp1AsxKv/y5Cr86oV9JtXydsgr1dA7JOfLs7od3X6Fc56T18uz7C7a0Vq5TXrVJ+dxWzz7p8sqA+nqhkpZWpnzO217+RVWekadcx79SXod/2CnPZM3WeuRld4RLng2etUgtlS735KuTJLrT5Fnsc5PV86a53eokkwbJy8etSWya4pHvX5bmOKR0lmdJr8jrpV5n8fPS5V+vOEa6vMzPZfADNRYMDjcHLU4AAABAxN1uAfUPLqkAWpwAAADAgfa3OHksr+M0uKsOAAAAwCR01QEAAMCBrjprLU6ErjoAAABwJAROpqDFCQAAAPj+SxvpBfyOKzkETgAAAECuQA25FGlHdOs4DQInAAAAONBVZ/GesYDzAifcVQcAAABgElqc6rCvagu5XJHxZYfU45Sv75gh7+9tnanOHJ7buUi6vPKovsp1/JXbpcs9z3yrXOfFX4dJl2/xyN+f9fYWSpf3yJFnzWU5+bvIqppl8qzLbN13R0qXr9mrzs6dpPhJkJdsKNfJV2Qq9mmyaVdWybOABwz1b5KqKnnW48oadTZkr1ue2Tjdq/61l6bIDs6SFc+5XOryqVHsd3W1JouzYnNedRJnCigyPFf5NZcrRflktN+iXKV3S/k5WbD7D8p11td8qnyu2C2v91WBFsp1PC75fqdmqDNWG63lmc1dNZXKdUpWyLOkryhW19MyQ56xP9OnXicrQ545nLIzlOu4SH7u+xTHdP86ut/88udS3Oq7xVK98rJrka2+blNn+bWRXOp6GvhKPqvCl9tOlS6vCqiPaYNAi5MpCJwAAACAXEYNuTQ/+lTrOA0CJwAAACAKBIgCfuvrOAwCJwAAADhwV52mH13C5cDB4QicAAAAYH9rk+W76vyOKzncVQcAAABgElqcAAAA4MBddda66ghddQAAAOBEroDfRuZwPzkNWpwAAACAyLAxxslA4AQAAAAO5AoELLcguZCOAGpVCpenVubwgwOKjLFEVJhWLl3eOm+Hcp3UrvIMvWXZ7ZTrpH3xlHT5nI9OUa6zZI88M3a+katcp2uuPEdHQdYe5ToeTTbrqu3yLMp7N+cp11m3tUC+Lc0voxZJ8pTVLZPU+5ZuI5t2lSJrdo1fnaW4ukbe0FutWUeZRVmbOVz9nNcjvzgahnp8Q02FPEt6lSJ7OnMpNpfiUed+SXKry1tJURc8LeTnIytsIc9+38Evz7LNNrjSlM/VUI2l7Om6jPApWaXqWQNyWsm3tWuDcp2N6+TXkrWVpZZvHWqZrF4lPV2+PSMtXbmOa698Hb+mLhqkrj8et/wYpXnU14tMn/zcz2qlngWhqu2h8id2L1Wus/HbntLlSxRlUGPIr9kNe1ed1TFOfnIa3FUHAAAAYBLGOAEAAMCBweFWE2D6HVdyCJwAAAAAXXUmIXACAAAAtDiZhMAJAAAA0OJkEgInAAAAIFfAsJxewBWwcRdsnMNddQAAAAAmocUJAAAADnTVWSyIAO6qAwAAAMdOuWJjHYeJyxanuXPn0oIFCyg5OZlOPvlkOvrooyOef+mll+iHH36IWNahQwe6/vrrLb9Xui9fZA+P2Fa6PFs0K0iXZwjPbi3PUsyM1vnS5cnrIj9DuN1vyzNwf7lVvi1WacizKHdIUWdDLkwrk++bItMuq6pQpxauVmQI37ZFng2Z7alMsZZNmzP++uT97hle9UmepMhmrc2mrcj2rZsos0qROVz3Pj63fL89LnUZJCsyoesyh/sDHsvHtapKfT6o9i9Vkzk8TXGMkjzqTOiGIveMS5OFPDtjn3R56yR1/U0JqDPc+0i+nkeTFifVJ88MnZRTolwnkNJRvs6aZcp1Vm8/Vbp8m2ercp0UypYuz0vWHLsseZkaKanKdahKvr19ivOEBQLquu3zZEqXZ3rVByInWX5tTGu7XblOTcZx0uUpK75WrrNkzQnS5es9q6XLA4a6zjcElxEgl+Y6pFrHaeJqjFNFRQX17t2bpk6dSqWlpbRmzRrq378/3XrrrRGve++99+jzzz+ngoKC0CM3Vz2tCAAAgOOJrjobDxvWr19Pa9euJb8//lqs4qrFyeVy0bPPPks9e/5vvh9ucbroootozJgx1KlTp9Dy7t270/jx45toTwEAAOIM31Fnea66gKWXP/HEE6LxgwMmwzDEvzNmzKALLriA4kVctThx11x40MS4BYpt3LgxYvmKFSvo9ttvp+nTp9OiRYsadT8BAACgNu4p4h4hbnEqKiqim266iS699FJauXIlxYu4CpxknnvuOcrMzKQjjjgiomUqPT2d3G43LVmyhI4//ni67bbbtNuprKykkpKSiAcAAICzWpxsPCzg1qb27duH/r7xxhupurqavv32W4oXcdVVF+3DDz+ke++9l5588kkRPAXdddddEd12Q4YMofPPP58GDx5Mxx57rHRbU6ZMoTvvvLNR9hsAAKC54eSXmvtNlOvUx48//ii67Dp2lN/00BzFbYvTZ599JvpEOUgaNWpUxHPhQRM777zzqEWLFvTll18qtzdx4kQqLi4OPbgJEQAAwDHq0eJUEtVjw704dSkrKxPjk/v160cnnXQSxYu4DJy++OILOuecc+iWW26pdUedCg9A0x1IHj+VlZUV8QAAAHCMegROhYWFlJ2dHXpwL44Ofx9zb1B5eTm98sorYohNvIi7rrr58+fT2WefTX/729/E4G9ZBPvzzz/TMcccE1r29NNP0969e+m0005r5L0FAACIEyIQsrEOkeilCW9w4MYIlaqqKhE0rV69WvQetW7dmuJJXAVOO3fupEGDBlFOTo7oTgtPNzBixAjq1auXiFp5ucfjocMPPzw0gn/atGnUp0+fJt1/AACA5p053OKkvcb+wMlsT00waOK76ObNm0dt2rSheBNXgZPP56M77rhD+lxqamroX+7K+/rrr+mnn36iAQMGiBantm3b2nrPXE9HcrsiMyO3TVWH5C0z9kqXJ+cWK9dxVcgrqmuBOhPwsp/Pli7fXKbJjuuR/wJoq04cTllJVZazXJftS1c+V1WVJF2+pzRDvY4im3WyJvu06rdOiiJjti7LdUDzWasVmcPdLvXFJ2BY7yH3KLaXpPk8SV511mGPIut6wO+2fOyqNRmevYrM3Zk+9b6lKzKepyaru9pVGcKNKvXncSnKtEWS+ni38BeqnwvIs/mnetV1IS1J/pm82aXKdVRHPLB2t3Kd1Xvl2bTLaL1ynVYB+WdtnSK/JrDkFvLrn+FrqVyH9srLe7cmI71Oqlee7DhLXn2F7DR5eXtby2dOYDWKxI+BJZuU6yzZJa8jpf4l0uVGAk5ncuGFF4qeo1mzZtHmzZvFg7Vr104kq44HcRU4cTRrNqnlCSecIB4AAADQPO6q27hxI3Xp0kXc2BWOp0QbOXJkXBymuAqcAAAAoPmNcTLr+++/p3iHwAkAAAAaJXBKBAicAAAAYP/AcKuBUMDiYPIEgMAJAAAADgROFgsi4LzAKS4TYAIAAAA0BbQ4AQAAwIExThYzeAec1+KEwAkAAAAQOJmEwAkAAAAwxskkBE51yPfnkdcVmXI2L0We2Zilp5ZbylLMAhvk2XZ3LTpUuc6q3fLsuLpxfbnJ8izXLZPU2Wl9bvlzVZps0aWl6lTkFYrs02VV6nmNVJIV2a91mbZVn0eso9me1SzgBukyh1ufzNKryBCuynZe13MqgYB62GONX37MaxTZ03X7oMoOrstWn5pSQVYFStX1qrpanpk6Rf1xKC+gnh4i3yXPmJ+tyZKemSbPTO3OUddTo0p+jSlfpc66vLZU/lkNzRWjFcmnz2iVuk+5TpIy47k6c3j19mzp8t2Vuqzv6udy3PJjlJOk/qzZGfLP5MrePyOFjHePPEN48bJOynVWKbKkBxTHQXd8GgRPn2L1+mSgqw4AAACciIMgq7Ga4bzACXfVAQAAAJiErjoAAADAGCeTEDgBAAAAAieTEDgBAADA/rHhFsc4Gc6bqg6BEwAAAKDFySy0OAEAAMD+O+osz1VHjoO76gAAAABMQosTAAAAoMXJJAROAAAAwKnK9z+sMJxXcAic6pDrTSFf1JQrmV75lAfMo5gao7pUnbq/epV8qobNGw5SrrOnSj6FQppmuohUxdFO96qnd1BNFVOtmXJFN21HZY18v6s167gVZ6bPxpQrXs06qj0wNFMQKKdP0ayj2p5uWh6PqgwU9c3uFDI1muPqVxwj1XJdeaf71FOuZKbIz6/k5ErlOgG/fB8CxRnKdUorUqTLfZoBDKppSFhBirzscpNV05AQZWWXSJe7cuTXBOYu3iFdvm1NO+U6m8rk9SeJ1NelghT5uZqXJp8iinmy5FPIuMrV07SUbpFPH7VLPvOO4POoj2uuP1+6vEWSeuqb9AzFMUpRT9nj3bxGunzTevVUWdsq5R8qyS0/DgFDfX43BCPgEg9r65DjIHACAAAAdNWZhMAJAAAA9reUW2xxIgd21eGuOgAAAACT0OIEAAAAGONkEgInAAAA2N9NZ7WrLuC8gkPgBAAAAPvHOGnuCJYynFdwCJwAAAAAXXUInAAAAMA0zslmuavOcFwB4646AAAAAJPQVVeHdK+LfK7I+NLnVmdz9fvlqbsrdqkzDleUpkmXby/JVr+Poh86Q57sV0jzyEfxJWkyTKuyXFf51VXHrfkFUqXITK3LPu1S/ADy2socHttfR3aygKu4Neuotud2qcvAzj7ojkONom7rJHnk2ZqTNJtKS5JnCPdqMj8HKiOz+wdV7ZWfW6y0Up45PElTR1orsmmL51Ll6+WnyrNps8z83dLlRkamch33xg3S5Ws2HaFcZ0e1vEyz3XnKddqmyT9Py0x15nBXsrw+unbIs52znVt6SZeX1qiPQ6anlfK5ViQ/5tm+CuU6SamKrPQ16jrn2rhFunzznuOU6+wz5JnDk93y4x0g9fs3CAwONwWBEwAAAIgfgroppmQM5/XUIXACAAAAjHEyCy1OAAAAICbstT7Jr+G4kkPgBAAAADbnqnM5ruQQOAEAAIDNMU4ux5Uc0hEAAAAAmIQWJwAAADiQANNie0rAeQWHwAkAAABsTrniclzJIXACAAAAjHFK1MDJ7/fTm2++SQsWLKDk5GTq378/DRw4sNbrVq9eTbNmzaLt27dT79696fLLLyefT5NWWyHZzZmEzb++vEKejXjfHnUW8L370uXrVCVbzjKd7g1Yzhzu02TgDtjJHK6ZLrs6IE8ZbZD1Xy2q7OD7n1NkMNbsm51M26rysTN4UPf+quNtJ9u4TkDTTK96Tvc+Po88y75Xk30/JUmeXdmlyZJeVZoqXV5arM7AXV4tzzae4lF/noIUdT1tmybf79ZZe5TrpB6kyKjtVWf0rl4lv8as2NNSuU6ZS541u7WRbfnzZGapM4erTi//+mrlKlt25UqX6+5ybxVoq3wuP1V+jcnwqffB41Vk6C5XZxuv2SK/Pu8qT1Wv45LX+xSSfwcESL3PDQJddYk3OLy8vJy6detGL7/8MuXn55NhGDR06FC69tprI163aNEi6tmzJy1fvpzatWtH06ZNo7PPPpsCAQd2xgIAAFjoqrP6cJq4anHyeDz0wQcf0MEHHxxaduyxx9LgwYPppptuokMOOUQsu/nmm2nAgAEiwGLDhg2jLl260GuvvUYXXXRRk+0/AAAAxLe4anFKSkqKCJrYYYcdJv7dunWr+LekpIQ+//xzuuSSS0Kv6dSpEx1//PH0zjvvNPIeAwAAxFceJ6sPp4mrFieZJ598knJycqhXr16hsU3chdexY8eI13HwtHLlSuV2KisrxSOIAzAAAADHwBinxGtxisZdb9OnT6fHH3+c0tPTQ+OgWGZm5KBQ/jv4nMyUKVMoOzs79CgsLGzgvQcAAGg+MMYpwQOn9957j0aMGEEPP/xwxLilrKws8e/u3bsjXr9r167QczITJ06k4uLi0KOoqKgB9x4AAKB5QVddAgdO77//Pl1wwQV033330bhx4yKe4wHinKbgl19+iVi+bNky6tGjh3KbvA4HVuEPAAAAxzAOZA638jDiMoyol7j7xB9++GEoaLr++uulAdCQIUPoqaeeCo1Z+uyzz2jp0qURA8YBAAAAEnpwOCez/OMf/0i5ubn0448/0ujRo0PPjR07lvr06SP++/7776dTTjmFjjjiCOratSt98sknNGHCBOrbt28T7j0AAEDzhSlXEjBwSk1NpRkzZkif42Aq6KCDDhKB1aeffko7d+6ku+++W9tNp+Nx7X+YyX7Nyirl2WRde9VpcEvK0qTLKzTZub1u+fZSFNnBWZLiOV0GblVG7xq/ugx0maT9iltXVRm4ddmxdVnA3YrN6TJt656zyk7Wbv2+KTKh29xnXXnHcr9VGcJ9qkzN/Jwiw7Oh6RKoKlNk7N+boV5HUYdVGfbreq5NWpl0eX7eTuU6ngJ5ZmpXWalynZ2/dpAuX71XPdMAkfx9DkqRZ09nBanyjOdJKf+7+ziaf7c8a3bZxnzlOttK5ccoWXUS83671RnPcxXFkOrVZOFWXE+NEnU9rdwhz3heprtuG/L6k2bIM9z7jcbNHG4Y+8c5WV3HaeIqcMrIyIhoZdLhLruzzjqrwfcJAAAgIdjJBB5wXh6nuBvjBAAAALHHrbp2Hlb99ttvYoaPU089lebPn0/xBoETAAAA7G89svOwYObMmXTuueeK4TVz584NzfoRTxA4AQAAQKO46KKL6NdffxXzy8aruBrjBAAAAA3DztxzhsXX5+erbxSIFwicAAAAoF7pCEqi5nflG7T4kYjQVQcAAAD1GhxeWFgYMd8rz/+aqNDiBAAAAPVqcSoqKoqYqixRW5sYAicAAACo1xinLAfN8YrAqQ6cvDY6gW2FJmt2aaU8g7GfJ0NUrVMlX6dGs44q27cuc3iyOxCzjNkB7ZmkPvF05aCiyhCuSSyszjYew+zguvexQ7dvqveJdSZ0qxfNuvbb65FnDk/SZA5Xbc9foz7vqirkv27LFZn8dVnxs3zqfUtRfB52UKY803b2QduV67iy5bMGGBt3KddZs+5Y6fIt5cpVKIPk15g2aepjl5Miz4SuU7EtR7p856bWynX2Vcuzl2epk5pTqlddF3KT5ccoRZM53PDLr0v+PfJyY1XFGZav2yku+X4HAoo678AJdOMBjgoAAACEWpysPqxYsGCBSHx59tlni7/vuusu8fcjjzxC8QItTgAAALA/ELI6xsmw9vrOnTvTLbfcIv771ltvDS3nweXxAoETAAAA2JpCxbA4y29eXp5oYYpnCJwAAACgXnfVOQkCJwAAAGiUzOGJAIPDAQAAAExCixMAAACgxckkBE4AAABARsD6mCVDm9QvMSFwAgAAALQ4mYTACQAAAGymI3A7ruQQONlQpZs+RTF9gC51fkWN/DAENHcreFVTcCjXUE/TopsyQ7UPuqlTdFN9qKa50FFNraJ7H+U0LdT0lNPBKPZZPBfjKWRieSeMbh88qml+FMt1qqt9yudqquXnULXi3BL75pLvQ1ZSlXKdnBT1vCZ5Obuly1MP2qlchxTdImW/qaco+XVXvnT5Pr+6TPN88utSq2TN9DI+eTlUV6rnQinZkiddvmOXfCoWsT3FtSQnyV7dVk2Zo5r+h9VUyetWdXG6cp2KUvl0OToZbkV9VBw69dFpGHy9133vqNZxGgROAAAAIAJ5y3mZAs4LnJrDD3AAAACAuIAWJwAAAMDgcJMQOAEAAAACJ5MQOAEAAAACJ5MQOAEAAAAFDLd4WBFAOgIAAABw7CS/VjOHG7irDgAAAAAU0FUHAAAAGONkEgKnOnBvb3SPb5WmKbNckalYlzG70u9RrKPmspgdXPecW/NOAcVTbpfL8jp26TKEW11Hl+XabhZuqxrrfRprv3V1TrWONru8oulflR2cVVYlWc5wn6TIJJ3sVedrzssoUT6X00qeIdyTXaFcx9ghX77lt57KddaXplqaTYC1SpGXactkdZZ01TEq02TMVmV3LylPs/w+2T51JnRdnUtXHD9Vpnhd5vCq4gzlOuWK46DruMr0eiylCK9u5Bl0RVedkZhddeXl5bR7925q0aIFpaVZz/oeDgkwAQAAIDTlitVHc7V8+XK68cYb6dBDDxXBUtu2bSk9PZ06d+5M48aNo59++snWdhE4AQAAQKjFyeqjudm8eTMNHz6cevfuTWvWrKGxY8fSO++8Q1988QW9++67dP3114vXHHvssTRkyBBat26dpe2jqw4AAAASpqvujTfeEK1KGzdupNzcXOlr/vznP9OePXvoscceo5dffpn+9re/md4+AicAAABIGOPGjTP1Oh7vNHHiRMvbR1cdAAAAJNwYJ7Z06VKqqJDfoMGDxefOnUtWIXACAAAAMgw745yoWXv//fepT58+IoAKN2/ePOrZsyfNnz/f8jYROAEAAEDCDA4Pd+2119KRRx4pgqeHHnqIKisracKECXT66aeLAeS33XYbWYUxTgAAACCCIKtdb0YzD5wyMzPp2WefpbPOOouuvvpquvPOO0VKgo8//pj69+9va5sInAAAACBh7qqTKS4upurqakpNTRU5nTIy1MlNEzJwWrVqFT355JPi33/+85908MEHRzz/8MMPi3wN4TgB1r333mv5vTzu/Y9wNZqKUuGXF2lAk09WVfF0kX8ss2nbod83apR9c2kynscyO7edbcX6+NjJhN5YdPtgpxz8ikz6OgFFhnDdvqX65Fmzk73VynVaZKszh6e0lD9n1GiuF+vzpcvXbS1QrlNSLd9elk/9PvnJ8gzUqZos6aqs6/vK1FmXKxWZw1XXReZzy/ctw6feN5/uuCqOn64uVCsyz1fsTVeuU1GRYrnOZ8iLhwIkr/PVBkbT1NfOnTvpqquuok8//VSkHrjwwgvp5ptvphNOOIEmT55Mt9xyC7nd7sQOnLiZbdasWXTmmWfS66+/Lj50dOD03XfficIKvyVRlcsBAAAAErPF6amnnqJt27bRkiVLQrHCo48+SmeffTaNGjWKAoEA/f3vf0/swGnkyJF0xx130IoVK2jmzJnK13Xo0IGGDh3aqPsGAAAQr+ykFwg088CJM4OPHz+evN7IcOecc84Rd9r98MMPlrcZd4FTx44dTb2OC+Oyyy6jli1b0imnnELnn39+g+8bAABAvErEFqdDDz1U+Vzr1q1Fy5NVCdmB6vP5xBw0AwYMEJlBr7zyShoxYoR2Hb5FsaSkJOIBAADgFImYALMhxF2LkxkPPPAA5eTkhP7mfA0nnXSSCKAGDhwoXWfKlCli/BQAAIATGeKWG4stTuS8wCkhW5zCgyZ24oknUl5eHi1cuFC5Ds9Xw7crBh9FRUWNsKcAAAAQTxolcKqqqqK7777b8nOx4vf7qbS0lDwe9W3OycnJlJWVFfEAAABwikTMHB7XgdO0adOUz3EupljZt28fvfbaaxHL/vGPf4j34VH0AAAA4OwxTg8++CC1b99e5HbixpW4GuP0888/U36+PAmczIcffkj//ve/ae/eveLvv/3tb6JrjtMUnHvuuZSUlESzZ88Wy3k0/Zo1a8QMyM8//zx17dq1AT8JAABA/ErEu+pUjj76aHHT2AsvvEA7duyg22+/nZpF4MTdY5xPKfjfPM4oHEd5PJ6IxxeZ1aVLF7r44ovFf48ZMya0/LDDDhP/cuDEBbFx40aRo4ETX3bv3l2kWLeDs2B7ouqFX1NRqhTZdu2Q59PV02XHtVO91QP/1O8T64zndrblboT3t5u5W5XxXHvs7LxPI60T6+2pLsS6C7SqznkUWamZ11MpXZ6RWq5cJyNr/w82GXeyPBN5zS711A671rWRLt+0L8ty3W4pT34tZPnkv6g9muNTpcj27a9UX+OqahTraK6LXsUxSvP4LWcbZ8ma9VSqqhQpvUn9vVGhyDauK9MMxTeuX7FKlZ0vgXrgGS4s53Gi5h04GYYhHtHZwfmGMR7/zM+5XNY+Q4MGTjxu6JFHHhG3+o8dO1b8dzgOcjiTZ69evUxvs3PnzuJRl7Zt24oHAAAAOLPFadq0abRnzx6aOnWqpeeaLHDiTJ3cOsQtS9zyZHcmYgAAAIBY4kadlBT5nINNPsaJ72ZD0AQAANB8ia46SoyuukWLFtHXX39N33zzDZWXl9fq8eLhQ//5z39s3ZzW5IPDAQAAoBmwk17AaJ6BE0+7xsHSrl27RK/X+vXrI57Pzs6mYcOG0QUXXGB52wicAAAAIKEm+R0zZox48M1ifBf+1VdfHbNtI3ACAACAhBkcXl5eTqmpqeK/hw8fXufry8rKLN15n5BTrgAAAIA1AZuP5mbmzJk0ZMgQWrBggfZ1PA0b38DGc9VagRYnAAAASJgWpxtuuIECgQCdeeaZIsF2v379REJsnkqtpKSEVq5cSZ9//jlt2bKFbr75ZpEw2woETgAAAJAwkpKSaMKECTRu3Dh68cUXxYwjTz/9tJhFhAeFc1Lsv/zlL6Ibz868tAic6joA7v0Ps6oVGXJ1MXlMs2lT4wgY+mzr6vVclstAlWk71tyK97FzfHQZsxsre3pjifU+qH7B2slI7/PWKNdJUjyXllamXMeXWqF8LlApzyRdsSdTuc6WbfLppvYqMnCzZI/8s6YolrN0r9/ysVNlDner0lyL65/H4gwE6kzbKR71sfO4DcuZyHUtI9XV8szhfr96YvhKxTq6Mk33yp+rUeybr7EzhxvWB3sHmv4SpJSenh4aJB5LCJwAAABABLi6IFfGaKZ5nBoSAicAAABo1HQEa9eupaKiIjH/bEFBQYOV/rfffkufffaZ8vnjjz+eTj75ZEvbROAEAAAAB7rqrBVEwOLrq6ur6fLLL6d33nmHDjvsMPrll1/EeKN77723QY7AihUr6LXXXotYtm/fPvr9998pJyeHJk2ahMAJAAAAmmdX3QMPPECffPIJLVu2TMxhy9OicIvPcccdR4MHD6ZYGzlypHhE4/fnlAVDhw61vE3kcQIAAIBG8d///pcuvfRSETSxE044Qcxly8sbE99ZN2jQIHr33Xctr4vACQAAAEJjnKw+GOdHCn9UVlaSLKM3d50deeSREcv57yVLllBj427Dbdu2WV4PY5wAAACADGP/wwrjwOsLCwsjlvPYocmTJ0csKy4uJsMwqGXLlhHLc3NzRY6lhsDJLrlbLhxP+vvTTz/Rv//9b/roo48sbxOBEwAAAIjxSgGbY5yKiooikkkmJydLE1MGW56i54oLPhdrPAj9jjvuiFjm9XpFV+EjjzwisopbhcAJAAAA6jXlSlZWVp1ZuLmliV+zYcOGiOX8d8eOHRvkCNx0003iEUsY4wQAAAD1GuNk1mmnnUZvvfVW6G8eC/X+++/T6aefTvECLU518Ln5Yb7T16+oRKrlTLX1xpqmJV7ppjVRPadbxw7VcbAzTYx22pkY77ed7Sk/ayPtt7bOK6bZ8CmWs6SkKvk6ydXKdQxD/VuzsjhdunzPtlzlOrvLMhTvoz77MxXzcKRoPquqHHRfelWK6UZ0v7Z11zkVj0ux35pNeTXXZFU90X1W1fQpLhtT0uhKIMUj/6zVinqleHlcmzx5skg9MHr0aDrjjDPE3XQej4duvPFGihdocQIAAADxc8/Ow4oePXrQggULRLDEQRMnweTs3nl5eRQv0OIEAAAAjTblSvfu3emJJ56I2xJH4AQAAADEPYNWewcDDiw3BE4AAABQr7vqnASBEwAAADRaV128w+BwAAAAAJPQ4gQAAAC27pIzHFhuCJwAAAAAXXUmIXACAAAA3FVnEgKnOvhchniYHRgWMGzcsqkYXJeI2cHtZNpurHJorGzjjfU+8cpOOXgVmbG93hrlOj6fPEO42+NXruOvkmeYZpVlKdLlxSWZynXKquUTm3o0ZZCmSCedotlvdeZw5SrkD8ivdIZm3+wMFNZlAbecbVxTf4KT0cpU13htZFa3/vWZpPisKe7YzUBQH7irzhwETgAAACDCNKt5mQwHlhvuqgMAAAAwCS1OAAAAILozLSfAJOflcULgBAAAAGLMm27cm0zAgX11CJwAAAAAeZxMQuAEAAAAyOOEwAkAAADMCti4qy7gwOLFXXUAAAAAJqGrDgAAAJAAM5EDp4ULF9Ljjz9Oq1atoieeeIIOO+ywWq+ZP38+Pfnkk7R9+3bq3bs3TZgwgXJyciy/l89j1Mr2aiezse4Wz4Bie43VHBjQ3E6q+qi6zMaNRVc+qmzj7hina1PVBV22czuZ0JUZ12N8HHT7ZifjuZ39Vr6PJsO0W5FJ2q3ImK17H3+1+rJYXSnP9M1KS9Oly/dVpCrXqQnIz70kTRZw1bmnyg5uN5u26jk7d1HZyYCtq4t2rj+qTOismjzS5QHDbXl7uv1OUmR9r1Z8P6i+GxoKuuoStKtu/PjxdM0111CbNm3o888/p71799Z6zZw5c2jAgAHUvn17Gj16NH3xxRfUr18/qqysbJJ9BgAAaO4Mw97DaeIucJo4cSJ9//33NHz4cO1rRowYQffccw8NHTqU3nnnHVq5ciXNmjWrUfcVAAAgXnDvg52H08Rd4JSbm6t9fufOnbRo0SIaPHhwaFnLli2pb9++oiUKAAAA1AkwrT6cJi7HOOmsW7dO/Nu2bduI5fz38uXLletxN154V15JSUkD7iUAAADEo7hrcapLVVWV+Dc1NXJAZlpaWug5mSlTplB2dnboUVhY2OD7CgAA0GzYGd9kkOMkXOAUvHOOu+zC8d+6u+p4XFRxcXHoUVRU1OD7CgAA0FxgjJNDA6fOnTtTZmamGOcUjgeUc1oCleTkZMrKyop4AAAAOAXuqnNo4OT1eumyyy6jxx57jHbt2iWWvfTSS7R27VoaOXJkU+8eAABAsxSw+XCauBscPnv2bHrooYeorKxM/D127FjRwsS5nYYNGyaWTZs2TaQf6NSpk8jlxIkyH330UerZs2cT7z0AAEDzZOcuuYADxzjFXeB0zDHH0OTJk6VddEEZGRki9QAHTzy2qWvXrmLAtx2coTY6S60+Y7V8uS4q92hmqlZRPdNYGb11GYdjPVow1tmxLWcb12W5juFnjfXnbKxys7MP2szhigzhquzg4jlF1mxdtnFVNv8aTebwqip15vCy8hT5OjXWL7M+Tfl4FVnFdee+nYz5quuP6hoX6wzhsa6/utkbagIey9nGazTPqaiOUbKi/hqN3J5jZ6y3Qc4Td4ETpxWITjWg0qVLF/EAAAAAcGTgBAAAAA3VVWctE3jAgU1OCJwAAADA1txzBgInAAAAcCI7d8kFyHnQ4gQAAABocTIJgRMAAACgxcmpCTABAAAAGgpanAAAAEB01Vm9S87A4HAAAABwIiTANActTiYy+Oqy+DZoRlsb72snQ3m80maftpOxOobHubHexy472ZpddjKrq95Hk9FblSHcTrkZAV22aPnlL6DJCF1eIc8OziqrfYr3UZ95qizcPrc8O7jYP23WfjmX9VWaNTt1QZcFXLU93bGzmu+IeRXv41ecD95Gbs7BlCvmIHACAAAA3FVnEgInAAAAwF11JuGuOgAAAACT0OIEAAAAGONkEgInAAAAwF11JiFwAgAAALQ4mYTACQAAAHBXnUkInAAAAAB31ZmEu+oAAAAATEKLUx08LoM8UVldlZm+m0lW6Hikyz7d1NG97piq9lv3eaxuS7cPdvYt1tvTZhvXZAi3um926M5VVYbwakUGcPFcjfqSWe1XZCLX7IOqTKOvOeHi8RJjpy7G+phrksgTKdbRZRs3bGRw5+8TK8vtXEfqg3P2W52rLkDOg8AJAAAAcFedSQicAAAAQAwOt9qCZMRhC2h9NXUvCAAAADQDHATZeTS0JUuW0NixY+m4446jefPmUVND4AQAAAChu+qsPhrS9OnT6YorrqBu3brRggULaOfOndTUEDgBAABAszRq1KhQi1NzgTFOAAAAsD9zuBgibl6ggbvqsrOzqblB4AQAAAD1uquupKQkYnlycrJ4JCIETgAAAHCgxYlstTgVFhZGLJ80aRJNnjy51utnzJhBL7zwgnabTz/9NHXv3r3ZHhEETgAAAEDGgf9ZYRx4fVFREWVlZYWWq1qbzjvvPOrTp492mx06dGjWRwOBEwAAANSrxSkrKysicFJp3769eMQzBE51cMluPYxxGvzGmqbFcNQ0LdanMFBOKRLjkmvO0/LEepoW9fvE9iZm1fQpAUN947C/xiNdXqWZVkX3nG56DhVV3dLVObf1mT4sv39z4LY7lY6N7ammxQnYWEc7BZFqahXF63Hbe/OE4wIAAADNMo/TV199JRJfnnzyyeLvW2+9Vfx9//33N9kRQ4sTAAAAkGHYGONkNGyrJQ8Sf/DBB2stLygooKaCwAkAAABstSAFGrjcWrRoIVqYmhMETgAAANAsW5yaIwROAAAAIEImqy1IhgPLDYETAAAAUMAwbEy5Yjiu5HBXHQAAAIBJaHECAACAemUOdxIETgAAANAs76prjhIycLrtttvovffei1jWs2dPeu655yxvy+MyxMNsRdFltG2udPus/DXhSrxMxXbEMpu2rYzDNtap67lYrhPL+qjLAm4EFOtosnnX+D2WlteVHVyVSVpHdfx02cFVY0oMzUmZaOedHXYyu2szlNsoUtUeRH/H1LW8ofD4JstjnMh5dSshA6d169aJuXDuuuuu0LL09PQm3ScAAIDmDIPDHRw4sZYtW1KvXr2aejcAAADiAsY4Ofyuunnz5tHxxx9PgwYNon/9619UVVXV1LsEAAAAcS4hW5w4Rfs111xDJ554Iq1du5YmT55M77zzjgim3G55rFhZWSkeQSUlJY24xwAAAE0LY5wcHDhNnz6dfD6f+O++ffvS0UcfTd26dRPB03nnnSddZ8qUKXTnnXc28p4CAAA0DwicHNxVFwyagrp27UqtW7empUuXKteZOHEiFRcXhx5FRUWNsKcAAADNa4yT1f85TUK2OEUrLy+nXbt2UWZmpvI1ycnJ4gEAAOBEHARZTS9gODBwSrgWJ24tmjp1KlVUVIi/+d/rrrtOtEINGTKkqXcPAACgWQq4ArYeTpNwLU6cr4lbmNq2bStSEmzevJm6dOlCc+bMocLCwqbePQAAAIhjCRc4eb1eMcj7jjvuoDVr1lBeXp64y84uzuxbK7uvLptsI2d6TRR2slw31jq67Ny652K5jh26z6raB1sZxd3W19FlAbeaHVxsT5EV2h/jzOHa7OUxnDVAX0dUmdXt1EVqtnTtGG47meftrKOdVaHhj3djXSuCuJvOapb5gAO76hIucAryeDx0yCGHNPVuAAAAxIX9I5ysdb0ZDpytLmEDJwAAADCPQyDrLU7Og8AJAAAAxEBvl8XB3gEHhk4InAAAAEAEQS6LgVDAgYFTwqUjAAAAAGgoaHECAAAAtDiZhMAJAAAAcFedSQicAAAAAIPDTULgBAAAAKLFyepgb8OBg8MRONWBc8O6rGRzjWH2YGjcTLyq/CW2smnHeJ1YZvpuDnRZwK1mB9dlCNdlAffbyDauyyRthzLDvTYLuCvusoDHWmN9VccyO7juPPY3k+8Ng/xkWLxnzCA/OQ3uqgMAAAAwCS1OAAAAcKCbDnmc6oLACQAAAA5M2Gs1cDIcV3IInAAAAODAGCdr460MB45xQuAEAAAA6KozCYETAAAAIAGmSbirDgAAAMAktDgBAAAABcR4JZeNdZwFgRMAAACgq84kBE4mMr1GZ3vVZQ9WZYa1k3FYm6HcQX3GyizgzeA2WJedDOUxPK52so3r1rOzb3aygBua80H1nC6jt50s4KrndHcV6fY7EOMs041xjbHzPvFKe+wUz8W6TBvjfKyPgGGjxclAixMAAAA4EM87Z3XuOQNz1QEAAIBzAydrLUgGAicAAABwIsPgSVcsJsA0GmvK5eYD6QgAAAAATMLgcAAAADjQ7WZ1ypWA40oOgRMAAACQYeMOOQN31QEAAIAT7R/hhBanuqDFCQAAAA4M9Mbg8LogcAIAAADLqQjsrhPvcFcdAAAAgElocWok8Tp9gcvV9OVgZzoCO+u4m3j6lMZ6n+YwzYWd6S+066imT7ExtYudfYtXzeG6pCpvu3VetT3dsdNNsxNLLostGI3dsmEYXOYBG+s4CwInAAAAsJVawEA6AgAAAHCi/akFjGaXObysrIyWLFlCgUCAevToQS1atKCmhBYnAAAAsBUEGQ0cOE2aNImeeuop6tixo+gW/Pnnn2nq1Kk0btw4aioInAAAAKBZdtXl5ubSihUrKDMzU/z9wgsv0GWXXUb9+vWjP/zhD9QUcFcdAAAANEs33HBDKGhiw4YNE/8uXry4yfYJLU4AAABQr666kpKSiOXJycniEWvz588XXXZdu3alpoLACQAAAOrVVVdYWFhrbNLkyZNrvX7ZsmW0cuVK7TZPPvlkysnJqbV8586dNHr0aDr//POpT58+TXbEEDgBAABAve6qKyoqoqysrNByVWvTwoUL6c0339Rus1u3brUCp+LiYjrzzDOpVatW9NxzzzXp0ULgBAAAAAeCJqutTob4fw6awgMnlSuuuEI8rOBuwNNPP53cbjd9+OGHEWOemkLCBk4vvvgiPfbYY7R9+3bq3bs33XvvvdSpUyfL2+HstdEZbHUj6lXZaf2arLUeRYZcw0bGX13WYzvZtJXb0uxdLLODi+1Z/AWkfR9qnOzcdrKaN4fMz7r601hZnFVZwP2K5Tq6dfxGbO+NsZOJXLUH/hjXBVV5684HO9m0VfugO96xrsPKzPMxzg5uJ8N9Yk7ya1BDCgZNbM6cOZSdnU1NzZ2oQRNHtJdffjk9//zzVFVVJW5djB68BgAAAM3XWWedRcuXL6drr72W5s2bJ7r5+FHXOKmGlJAtTnfeeSddffXVNGbMGPH3rFmzqHXr1vTMM8/Qn//856bePQAAgGZn/0Bviy1O1LAtTvn5+XTKKafQ7NmzI5Zzw0iXLl2oKSRc4LR582b67bff6IEHHggtS0tLEy1On332GQInAAAAKeuBEzVw4FTXQPKmkHCB04YNG8S/BQUFEcv5759++km5XmVlpXgEoVsPAAAcxcYYJ2rgMU7NUcKNceJJAJnXGxkT+nw+8vvVQy6nTJkiBp0FH9E5KQAAABK9q87Ow2kSLnDKy8sLJcoKt2PHDtFXqjJx4kSRJyL44JwUAAAAzhGw+XCWhAucOOUAB09ff/11xO2S33zzjTbTKCfrCuahMJuPAgAAAJwl4QInTpDFd9TNnDlT3K7IQdO//vUv0QI1atSopt49AACAZsrYP2bJyoOcN8Yp4QaHszvuuIO2bt1KPXr0EC1JnGX01VdfpYMPPthyUq9yf5XkOVeTJ8BUraPjV0zg6NFM7KhKTudxqdfxagYLqvbbq9kHr5gGwFoiSc+BsW61tuVWj3PzuAOWE/S5VesoljOXoiro3kf1WVX7XJdYJuFsvASYLsvvU12j/pw1ijqiO7+rFevs379AzBJ3Vus+q41Ejqry1p3Hdt5HdX4HNNXN7YptHVXVk1iWG6u2kajV6vdD8PunoZNMRiYXcF4gZJXLaLwj0ujKysrE3XE8tw23RFm9Ow8DxAEAoKnxmNt27do12PYrKirEMJctW7bYWr+goIDWrFlDKSkp5AQJHTjV9+68TZs2idaqvXv3iiAqehJDp+OgFOWCMkFdwTmE60rD4K9n/v5p06aN5R//doInnmXDjqSkJMcETQnbVRcLXEmDEb7rQN8KBo3LoVxQJmahrqBcUFesaay52TjwcVLwUx8JNzgcAAAAoKEgcAIAAAAwCYGTCXxn3qRJk8S/gHJBXbEO5xDKBXUFEgUGhwMAAACYhBYnAAAAAJMQOAEAAACYhMAJAAAAwCQEThrFxcU0ceJEOuWUU+j888+nt956i5ymsrKSnn/+eTrjjDNo0KBB0tf4/X6aMWOGeA0/+L95WSInR33llVdoxIgRdOqpp9K4cePot99+q/W67du301/+8hfq378/XXjhhfTJJ59QIuNjPmvWLBo2bBiddtppdN1119Hy5ctrve77778XZXfyySfT2LFjRcZhp+DryXHHHSetC3yenXPOOTRw4ED6xz/+QeXl5ZSodu3aJcoh+jFv3ryI15WWltKdd95JAwYMoHPPPZdefvnlJttngCAkwNR8CZx++univ2+77TYxYfDQoUPpP//5j7joO8WRRx5JvXr1otatW9OcOXOkr7nhhhvojTfeEJMpc7JQDhZWrFhBjzzyCCWia6+9VgTV/CXH5fLiiy9S79696ZtvvqGePXuK1/CXXr9+/eiggw6i8ePH06JFi+jMM8+kt99+m84++2xKRFwP+Ly5+OKLKT09XQQCxxxzDC1YsIC6desmXrN48WLq27evCJiGDx9OTz/9NB1//PH0448/irJMZK+99hq99957tHTpUtqxY0fEc/fdd58IEO6//37Ky8sT15xvv/2W3n33XUpEnKGa68ULL7wQMYfooYceGvG68847j7Zt2ybKZvPmzXTllVeKsuOgHKDJ8JQrUNsrr7xiuN1uY8OGDaFlN998s1FYWGgEAgHHFNmePXvEv/fdd5/RunXrWs+vXbtWlNPs2bNDy1599VWxbP369UYi2rt3b61lPXv2NK655prQ348++qiRmppqFBcXh5aNHDnS6NWrl5GoSktLay3Lzc0VdSdo8ODBxsCBA0N/V1dXG+3atTNuueUWI5GtWbPGaNOmjbF48WIxnfxLL70Ueq6srMzIyMgwpk+fHlq2cOFC8bovv/zSSESbN28Wn2/p0qXK18yZM0e8Zvny5aFld999t9GyZUujqqqqkfYUoDZ01SnMnTtXtLa0bds24tcPz1cn65Zxarp/blrn6Wm4NSWIW2K45enTTz+lRJSRkVFrGbewhM/zxPWHu6LC5zbk+rNkyZJarQ2JIi0tLeJvbkXilrkjjjgiNO8W1wnucgnyer2iBS6RuzFramrokksuodtvv50OP/zwWs9zy8u+ffsiyuXoo48W85MlcrkEWym5G+5Pf/qTODfC8Tl0yCGHRJQZn0PczcctuABNBYGTwrp168SFK1zwb34O/ldOLVu2jJjjiP87JyfHMeXEwQB3q/A4ODP1Z/369ZSofv31VzFWpUePHqJL7plnnhHjnRh/4XGAICuXRK4rf//73yk/P5+uvvpq6fPBz+60cuEu2lGjRtGECRPI4/FQnz59IromcQ2G5gpjnBSqq6spMzMzYllqamroOfhfOckyqnNZOaGceCwXD4YeM2aMaGnTlYsT6k9hYSE9+OCDtHv3bjH2i8d3HXXUUWKMU/Bzy8olUcvk448/pueee060vqkEPzvPMO+UcuFxXPPnzxcBE+MWax4XyPUleB459RyC5g8tTgrcisK/kMPt3LlT/Jubm9vwRyaOyylYVoleTr///ru4A4rvJHzssccinnNq/eEuS25xOuuss8Qddh06dKBp06aJ51q0aCG6dWXlkqhlwkETd1FyNxyXC3ffMu6244HOwbrCONh0SrlwF20waAriO1R5GETwbkKnnkPQ/CFwUuDxTdznzuMTwsci+Hw+6t69e2Mdn7goJ75lOPy2859//pnKysrEnWaJatWqVSJNBd859+yzz4qAILpc+Lb7cFx/eMxYp06dyCkKCgpCY7q4C7dr1660cOHCWuWSqHVl8uTJNHv2bNEKxw++e46NHDmSbr755lBdYeHlwkEU38mbqOUis2XLFtHqxtfYYLn88ssv4loSXld4/GRw3BxAk5AMGAfDEHeEpaSkhO504bujunfvbgwfPtyR5aO6q66mpsY49NBDjUsuuUTcbciPiy66SCzj5xLR6tWrxd2Vl156qfIz8t1CHo/HmDVrlvh769atRocOHYwbb7zRSFT33nuvuEMs6LPPPhN3Fs6YMSO07IEHHjBycnKMFStWiL/nzp0r7sD84IMPDCcoLy+vdVcdGzBggNGvXz+joqJC/D1+/Hhx91jwrtZEw3feLlmyJPT3smXLjIKCgojr6/bt242srCxj0qRJ4m+uW8ccc4wxaNCgJtlngCAEThqvv/660aJFC6NTp07iduH+/fsbu3btMpzkz3/+s3Hsscca7du3N3w+n/hvfnAagqCffvrJOOSQQ0Rg1apVK/HfvCxRnXvuueLL7+ijjw6VBz/C0xGwZ555RtQbLg8OIHg92S37iWLq1KmiDvAPjI4dOxrZ2dnG5MmTDb/fH3oNB5pXXXWVkZycLIJr/veee+4xnEIVOBUVFRlHHnmkuN5wUM7l+OmnnxqJatGiRcZxxx1ntG3b1jj88MONpKQkY/To0UZJSUnE6z766CMjPz9fXH+4PvF5xqkMAJqSi/+vadq64kNFRYW4U8hpXSxB/Nn37NlTazknxQy/k46zafNrGd8+HN115YQy4dQDwUSPQdyNyV0uPCaDB04nOu7a5s/Lg3r58wa7XaJt3bqVNm3aJM4pHvvkFHy55e6mLl26SMfpcNnxGB/u0lSVXaJ1z/G4JU6CGRz4HY3TfPA5x+kuOD0BQFND4AQAAABgUuI2CwAAAADEGAInAAAAAJMQOAEAAACYhMAJAAAAwCQETgAAAAAmIXACAAAAMAmBEwAAAIBJCJwAAAAATELgBAAAAGASAicAiJkff/yR3njjjYhlO3bsoJdffpm2bduGkgaAuIfACQBiJicnh6666iqaMWNGaNnIkSNp5syZ0rnZAADijbepdwAAEkf79u3p8ccfpyuuuIJOOeUUmjdvHn311VeiJcrj8TT17gEA1Bsm+QWAmONWpm+++YaKioro6aefpksvvRSlDAAJAYETAMTcunXrqFOnTtSnTx9asGABShgAEgbGOAFAzI0fP566dOlCixYtog8++AAlDAAJA2OcACCm/vOf/4hgafHixfT888/TlVdeSUuXLqX8/HyUNADEPXTVAUDMrFq1inr16kXTp0+n0aNHk9/vp759+4qg6a233kJJA0DcQ1cdAMTM22+/Tdddd50ImhjfScetTmlpaeLOOgCAeIcWJwAAAACT0OIEAAAAYBICJwAAAACTEDgBAAAAmITACQAAAMAkBE4AAAAAJiFwAgAAADAJgRMAAACASQicAAAAAExC4AQAAABgEgInAAAAAJMQOAEAAACYhMAJAAAAgMz5fx09vhumFO6nAAAAAElFTkSuQmCC", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "z_demo = jr.normal(jr.PRNGKey(0), (state_dim_ks,))\n", + "u0_demo = gp_mean + L_gp @ z_demo\n", + "\n", + "n_steps_demo = 30\n", + "predict_times_demo = jnp.arange(n_steps_demo, dtype=jnp.float64)\n", + "\n", + "demo_sim = dsx.simulate(\n", + " build_dynamics_ic(u0_demo), rng_key=jr.PRNGKey(1),\n", + " predict_times=predict_times_demo, n_simulations=1,\n", + ")\n", + "print(\"states shape:\", demo_sim.states.shape, \" observations shape:\", demo_sim.observations.shape)\n", + "print(\"any nan:\", bool(jnp.any(jnp.isnan(demo_sim.states))))\n", + "\n", + "fig, ax = plt.subplots(figsize=(6, 4))\n", + "mesh = ax.pcolormesh(x_grid, predict_times_demo, demo_sim.states[0], shading=\"auto\", cmap=\"inferno\")\n", + "fig.colorbar(mesh, ax=ax, label=\"u(t, x)\")\n", + "ax.set_xlabel(\"x\")\n", + "ax.set_ylabel(\"t\")\n", + "ax.set_title(\"Forward simulation from a GP-drawn initial condition (sparse H)\")\n", + "plt.tight_layout()\n", + "plt.show()\n" + ] + }, + { + "cell_type": "markdown", + "id": "be1da31d", + "metadata": {}, + "source": [ + "## Inferring the initial condition\n", + "\n", + "Pick a different GP draw as ground truth, generate a short window of noisy sparse\n", + "observations, and recover $u_0$ with NUTS running through the EnKF-marginalized likelihood.\n", + "`EnKFConfig.crn_seed` is fixed (required for `jax.grad` through the filter to be\n", + "well-defined), and NUTS is initialized at the closed-form GP-regression posterior mean given\n", + "the $t=0$ observation alone -- exact, since both the prior and the observation model are\n", + "linear-Gaussian:\n", + "\n", + "$$\n", + "\\mu_{\\text{post}} = \\mu + K H^\\top (H K H^\\top + R)^{-1}(y_0 - H\\mu).\n", + "$$\n", + "\n", + "`max_tree_depth` is capped at 6 (the default of 10 turns out to be needlessly expensive\n", + "here -- see the full notebook for that comparison).\n" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "id": "37213baf", + "metadata": { + "execution": { + "iopub.execute_input": "2026-08-11T21:53:27.164538Z", + "iopub.status.busy": "2026-08-11T21:53:27.164464Z", + "iopub.status.idle": "2026-08-11T21:53:28.154203Z", + "shell.execute_reply": "2026-08-11T21:53:28.153934Z" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "||warm_start_u0 - true_u0||: 0.09152363621851407\n" + ] + } + ], + "source": [ + "n_steps_infer = 6\n", + "predict_times_infer = jnp.arange(n_steps_infer, dtype=jnp.float64)\n", + "\n", + "true_z = jr.normal(jr.PRNGKey(42), (state_dim_ks,))\n", + "true_u0 = gp_mean + L_gp @ true_z\n", + "\n", + "synthetic = dsx.simulate(\n", + " build_dynamics_ic(true_u0), rng_key=jr.PRNGKey(7),\n", + " predict_times=predict_times_infer, n_simulations=1,\n", + ")\n", + "obs_values_infer = synthetic.observations[0]\n", + "\n", + "# closed-form GP-regression warm start from the t=0 observation alone\n", + "H_dense = interp.as_matrix() # GridInterpolator exposes the same H densely, for this linear algebra\n", + "y0 = obs_values_infer[0]\n", + "K_obs_obs = H_dense @ K_gp @ H_dense.T + observation_model.R\n", + "K_x_obs = K_gp @ H_dense.T\n", + "u0_warm_start = gp_mean + K_x_obs @ jnp.linalg.solve(K_obs_obs, y0 - H_dense @ gp_mean)\n", + "z_warm_start = jnp.linalg.solve(L_gp, u0_warm_start - gp_mean)\n", + "print(\"||warm_start_u0 - true_u0||:\", float(jnp.linalg.norm(u0_warm_start - true_u0)))\n" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "id": "a5087c51", + "metadata": { + "execution": { + "iopub.execute_input": "2026-08-11T21:53:28.155271Z", + "iopub.status.busy": "2026-08-11T21:53:28.155208Z", + "iopub.status.idle": "2026-08-11T21:54:05.791769Z", + "shell.execute_reply": "2026-08-11T21:54:05.791507Z" + } + }, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "sample: 100%|██████████| 300/300 [00:37<00:00, 8.01it/s, 63 steps of size 2.50e-02. acc. prob=0.93]" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "||posterior mean - true_u0||: 0.07786563009288042\n", + "mean posterior std across grid points: 0.006649992493998666\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "\n" + ] + } + ], + "source": [ + "def ks_ic_model(obs_times=None, obs_values=None):\n", + " z = numpyro.sample(\"z\", dist.Normal(jnp.zeros(state_dim_ks), 1.0).to_event(1))\n", + " u0 = gp_mean + L_gp @ z\n", + " dynamics_ic = build_dynamics_ic(u0)\n", + " with Filter(filter_config=EnKFConfig(n_particles=32, crn_seed=jr.PRNGKey(0))):\n", + " return dsx.sample(\"f\", dynamics_ic, obs_times=obs_times, obs_values=obs_values)\n", + "\n", + "\n", + "nuts_kernel = NUTS(\n", + " ks_ic_model,\n", + " init_strategy=init_to_value(values={\"z\": z_warm_start}),\n", + " step_size=1e-2,\n", + " max_tree_depth=6,\n", + ")\n", + "mcmc = MCMC(nuts_kernel, num_warmup=150, num_samples=150)\n", + "mcmc.run(jr.PRNGKey(1), obs_times=predict_times_infer, obs_values=obs_values_infer)\n", + "\n", + "posterior_z = mcmc.get_samples()[\"z\"]\n", + "posterior_u0 = gp_mean[None, :] + posterior_z @ L_gp.T\n", + "posterior_u0_mean = jnp.mean(posterior_u0, axis=0)\n", + "posterior_u0_std = jnp.std(posterior_u0, axis=0)\n", + "\n", + "print(\"||posterior mean - true_u0||:\", float(jnp.linalg.norm(posterior_u0_mean - true_u0)))\n", + "print(\"mean posterior std across grid points:\", float(jnp.mean(posterior_u0_std)))\n" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "id": "d8569e2e", + "metadata": { + "execution": { + "iopub.execute_input": "2026-08-11T21:54:05.792733Z", + "iopub.status.busy": "2026-08-11T21:54:05.792659Z", + "iopub.status.idle": "2026-08-11T21:54:05.899932Z", + "shell.execute_reply": "2026-08-11T21:54:05.899719Z" + } + }, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAk4AAAGGCAYAAACNCg6xAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjExLjEsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvctoD+AAAAAlwSFlzAAAPYQAAD2EBqD+naQAAjTxJREFUeJzt3QV4k1cXB/B/6u5eWrx4cXeXwXAdDBkuY2MwxgdDNhhjDNmwAWO4bTBguPtwd5cWWuru7fs9576ka0tbUmiJnd/zBOK5uUmTk3vPPVchSZIExhhjjDH2VgZvvwpjjDHGGOPAiTHGGGMsD3jEiTHGGGNMRRw4McYYY4ypiAMnxhhjjDEVceDEGGOMMaYiDpwYY4wxxlTEgRNjjDHGmIo4cGKMMcYYUxEHTqxATZo0CQqFgntZRSNHjoSVlZXGvy55vf37PJ62vod+++030e6nT5+qpd+Y5qD3AL2O9J7Iza5du8T1Tp069cHaxvKOAycNU6JECfGHozxYWFigYsWKmDdvHlJTU/P98caOHQsjI6N8v1+mf32uq8+roHG/sQ/t5s2b6d8xGzdufOPyhQsXissuXryYft66devEedu3b8/2Pr/44ov0HwqFChXK9D2W02Ho0KHitg8fPsSAAQPE95+5uTmKFCmCDh064J9//imQ7733xYGTBqJAibYQpMOTJ0/QsmVLjBkzBl999RW0zfTp08XzYLr1uuT19vw+4H5jmul///sfkpKS8vU+/f3907/D6DBlyhRx/pUrVzKdTyNwFMRVrlwZt2/fxtq1axEaGopjx46haNGi6NixI06ePAlNw4GThnN1dcVPP/0EX19f8SZLSEhQd5MYY4zpgDZt2ogRIhphUpdff/0VMTEx2LRpE2rXri1mWWjEiWZZtm3bBjMzM2gaDpy0RKlSpZCYmIigoCBxmiLyJk2awNraWrzRatSogb/++ivTbegP4pNPPoGnp6e4Trly5fDdd98hNjZWXN6lSxfMmTNHDIVmHD6lYVMl+hXQrVs3ODs7w8TEBKVLlxZv6IyjDcq8nKioKDHc6ujoiJIlS+aao6FK+3O735y87X5/+eUX0R765ZPVH3/8IS47c+ZMvj3/rFTpc0Kv0cCBA2FnZwdbW1vxOkZGRr5xf6q0LzvZvS7K56HKY2e9/dueV3aPp5wOUB5sbGzQoEED7N69G++K+qNnz55wc3MTQ/70S5Ze17S0tHd676n6OtD0BY0U04c8/a2uX78+2/blR78V1HNQxdve6/S8a9asCUtLS3G9pk2b4uzZs2/cz4YNG1C3bl3RfrqPdu3a4dq1a+mXU3/MmjVLvJ9NTU3h5OSE7t27Z/o7qV69OsqXL59tO6lvaNon49+BKm2jPqLpI/oMoPbRe2jq1KniMuozmlYtVqyY+Ftzd3fHsGHDEB4enuk+7t69KwISehwXFxeMGzcOKSkpee7rmTNnis9uakPDhg0zTZtRPxgYGIiR3KyePXsGQ0NDTJ48+a2P0bhxY9HWGTNmICIiAuoQFhYm3uPUn1l9/PHHqFWrFjQNB05a4t69e+KPlf4Q9+7di2bNmok/qhs3buDx48fiNH2BLl68OP029Afh5+eHw4cPiz/uv//+W5xPUTzZsmWLmP6jP7KMw6f0gUMuXbokPpDj4+PFBzUNof7www8i+KIPg6zoA4fmpemPesKECTk+F1XbXxD326dPH/FBTF+mWdF5ZcqUEb96Cur5v63PlT7//HNxX8+fPxe3oaTRrFO1eW2fqlR57Hd9Xlm/hJXXoy9KGrKnL0N67AsXLuS53efPnxe3Dw4Oxp49exASEoLVq1fj9OnT4u/nXd57qvQF5WF06tQJ9erVE/d36NAhHD9+HPv27Xtrm9+l3wriOeRVdu/1b7/9Fv379xftePTokbiMRsrpS//cuXPpt6Xr9+3bF61atRKvOV2P7m/+/Pnp16H7oYDl66+/xqtXr3DkyBHxQ5ACH0pfIBQM3rp1K9N9E7qc/h4osFMGnaq2TXl7ClqWL1+OBw8eoEqVKoiOjkb9+vWxY8cOLFu2THzZ03uMAi/qe+VUV0BAgLgeXU6X0fuOXsuJEyfmqX/pxw+1nYLJ69evi+CpUaNGuHPnjric7rNFixaiLVlzgGhmgt5D9HxVQQEqBYUUPKlDvXr10qfz8nvKsMBITKMUL15cqlixYvrpV69eSd988w39bJJGjhwpzitbtqxUsmRJKSUlJdNtW7ZsKdnY2EixsbFSYGCguM3SpUtzfbyvvvpKMjQ0zPay2rVrS8WKFZMSEhIynT9//nxxGz8/P3F6xIgR4rFWrFjxxn1MnDhRXJaRKu1/2/1mR9X77dGjh2Rvb5/ped29e1c81s8//5yvzz+vfa68rw0bNmQ6/8svv5SMjY0ztUXV9mUnu9clL4+d3e1ze17ZXT8nJUqUkIYOHZrn21arVk3y9vaW4uPjc7xOXt97qvRFmTJlpEqVKmW6Xlpamngsuo8nT57ka78VxHNQVU7v9fv370sGBgbS+PHj37hN9erVpaZNm4rjd+7ckRQKhTRmzJgcH+Py5cviMSZPnpzpfH9/f8nU1FTq27evOB0ZGSlZWFhIgwYNynS9SZMmif6k6+elbcTW1laytLSUQkNDM11v2rRpot3Xrl3LdL7y+axcuTJT3z5//jzT9b744gvxnJYsWSLlZufOneJ63bt3z3R+RESEaFu3bt3Sz9uxY4e47vbt29PPS0xMlFxcXDI9p+zcuHFD3Hb27Nni9IABA0TfPn36VJxesGCBuPzChQvpt1m7dq04b9u2bdne5+jRo994vytNmTJFXHblypU3LktOTpY+/fRT8RpZW1tLrVu3liZMmCAdPnxYSk1NlTQRjzhpIPqVoRyyL1y4sJi6oDwn+hUSGBgopiPat28vfqVmRMP+NIROv9ZpaNvDw0P8cqJf3fSrLS/o1zoNV9NQKY3QZES/sOhXDv2Sz4iu+zaqtr8g7/ezzz4TI3DKkTflaJOxsbEYkSrI56+qjz76KNNpmpJITk4WI4jv2r78euz8EhcXJ0YCaJSPflFnnKbKOnX5NjTKRFMZNPKTU07Eu7z33tYXNMJAowA01ZQRPY/8fD8U5HN4F1mfG42+0HRo165d37guTYlRgi+NKtBoGf3fq1evHO+bRsgJvZYZ0QgbTdsoL6epXXq8zZs3i/cSoTbQ5x2NZtH189K2jCMgDg4Oma63c+dOMQVLo1QZ0VQiPQ6NMCrbXqlSJXh5eWW6Ho3OvU//0hQrTavRyJtS27ZtxfdDxlFGGlGkdA76jMuL77//Xryf8joylh9oJS69ZjTFSKkU1M80WkuvDY2o00ippuHAScNX1dE0DA3V0tQLvcFoOoZQDkdWyvPoS5X+CA4cOIAKFSqIYXC6rGzZsmIaR/khkxtloEWJe/S4dH90oHl1ZV6Bsi2EphEpWHsbVdtfkPdLf5C0YmPFihXiNOUfrFmzRnz50VRoQT5/VVDOCn0pZKQ8rcxDyGv78vOx88unn36KBQsWiOlFWoVDwR695+mLh77Y80KZ+6f8ssyP954qfaG8T1rEkVV2572vgngOeZXde50COkJTacr3I70X6fDjjz+KKRia7sqP1ynj86PpOgoWlfld9JlHAWHGwEHVtill1za6D5p2o9tnvA8KkOm9q2wz/Z8f74Wc7iPj3zQ9/pAhQ3Dw4MH0HxpLliyBvb29WI2WF/Qj+8svvxS5Z9nlfyofj6RlyBfMSHl+1oBeVVTCgKYXaYDg8uXLIrWE2kJ5ZJqGAyctQ4mUJLsRJOV5yg81Sgan/Av6gKTRB/rlOW3aNAwfPvytj6O8D1qqSoEFfanRgf44lEFdxvuh0Zr8bn9B3S992NEfKP06pF859IuUPhgzftgW1PNXhSoFD/Pavvx87PxAq2hoxI/aSB/y9PopP5gzFoxUFSXHkxcvXuTbe0+VvlDlPvNTQTyHvMruva58TMo5Ur4f6b2Y8f1IAZsqr5NytCen55jx+dHoEI1QKHMW6X/68UOjMXlt29ueX9WqVcXtM96H8vb0Oat8ffLjvZDTfWQdCaPAkQJZymuifDcqnEnJ/++yEm38+PHieeaUH6nsx5AMgWtGFBTT+y1rG98VfS7QCB+XI2DvjX5x0dQG/aFmjfy3bt0qPgCqVauW6XyayqlTpw5mz54thntPnDiRfhmt/KD7ybrqg37dUKItfbnl9dd/fre/IO6XAif6I1+5cqUYeaJfmVQvq6Cff259nhcF2b4P8byo7+kLJ+s0ozLQzyv6sqT+oF+pOZXsKIj3Hq0EoukaSrjOiJ4bTe/kd78V1N/P+6IfZRT40rRZbmjBCr32NLKRExoRJhmn0snLly/F9LTyciX6wUOfaXQZJW/TSGbG4EfVtuVGuerv/v37uV6PVvNdvXpVjEJlRO3Ki6zvHRpVo4T3rM+dAlGaoqXPsblz54rz8jpNp0QrHClBm35Q0pRqVjRtRlPqBw4ceOMyWvFN05X0N0jv57ygdBJ6bbOivwcKFmmlo6bhESctRAEQrfagVSM0YkJvLlq2TG92Wp5Kb1wa4qQcARrGpV8CNOVHx2klFgVPSjStQx/y9MGfdXUGDfvS/DJF/nQ7muKjDwT60G7evPk7TQWp2v6Cvl8aFqZAaenSpWLEqV+/fm8MMRfU88+tz/OioNr3rvLyvOi1oC8Z6n9afUQjUMpcvqx5JKqi8gb0a5i+5Oj9T8vwaZqbfpUrVyMVxHuPphrpy3LUqFHiC4BGU2h6IbeVce/zfiiI56DcEoT+Dt4FBY+UH0OPT1+E1C76zKF8LJp6oSkl5fVopRxNMVNeDb1/KVCmv0HlKjBaxUY5ULTai3Jf6HIaTencubP44s66zJ5W6FGgROUKaNqN+uVd2pYbGoWhVAcayaKghlbNUZ7kv//+i8GDB4vcIkLlCih4pdV7tGKQ2v77779nGxjkhgIRep3pb5im4Xr06CGem7KQZEYjRowQ7Vm1apXoO5rqflfUFz4+PuL1yIoCGHr8rVu34ptvvhGrD+kzhwJKyuGiNlCb84o+u+hvgAI/+tuhHz7Ud/QeoBzC91khXGDUnZ3Ocl9Vl5NDhw5JjRo1EitAzMzMxIqijRs3pl9OqxFotUWrVq3EKgu6Hq3GmT59eqbVNHS9wYMHS87OzmJ1CL0lHjx4kH45HadVLJ6enmK1CK1a6tSpk1jxkHGlDd1/dnJaEfW29r/tft+1XzLaunWraBs970ePHmV7nfd9/tnJrc9zuq+//vrrjVUuqrYvL6vqVH3s7G6f2/PK7voBAQFi9ZCjo6NYDda+fXuxGqlmzZpSw4YNc32snNy8eVPq0qWL5OTkJJmbm0uVK1eW/vjjj0yrc97nvZfT60DvpQoVKkgmJiZiVeDq1avFCipVVtXltd8K4jkoV1nRaqbcvO29Tv1A7aIVYLTirVy5ctK4cePeWGW2Zs0aqUaNGuI1oteqXbt20tWrVzOttJoxY4bk4+Mj3tcODg7idb137162j9uxY0fR/jp16rxX2+iyIUOGZHv7mJgY6dtvvxWrKGkFGrW7fv36YoVhxs/UW7duiRWOyudGK+1oZV9eVtWdOHFC+u677yR3d3fxWPQ4586dy/F29D6n2y1atEhSRdZVdVn7iS7L7n2ufP80et2PtHqRvl86d+4sXbp0KcfHy21VXVBQkPTLL7+Iv3k3Nzdxn9Rv1Ic5reBTNwX9o+7gjTHGmPosWrRIjCJQjSPlAgmmPajGE9WjohEaTZza0jU8VccYY3qO8lpok1YOmrQPBUuUQE3lFjho+jB4xIkxxhjTQpQPRLl1lBxOuUa0kpoVPB5xYowxxrQMbUlDdbr2798vVgZz0PTh8IgTY4wxxpiKeMSJMcYYY0xFHDgxxhhjjKnISNUr6jOq0EsFzKiy6ofakoIxxhhjHwZVZqI9C2nfPuX2TzoVONFKAqrISktn3/YEs1ZjpaWbVNk1L/vpUNCUdbdrxhhjjOkW2iSadpbQmcCJRn5oB2fapoH2uKLy+7TNAu3Vo4pBgwZh7dq1GD16NObPn6/y49JIk7JDs+40zhhjjDHtRvsB0gCJ8vteZwIn2gdn/fr1uHjxolh6SXt19ezZU2x6+balmOvWrRMbNL7Lkk3l9BwFTRw4McYYY7pJlXQcrUoOX7x4sdiwkzYEpCc3fPhwFClSBMuXL8/1drQhJm0sScGTkZFWxYqMMcYY0yBaEzgFBQWJnbTr1KmT6fy6deviwoULOd6OdpSmnaVpZ2xVdytnjDHGGNPqwCk4OFj87+jomOl8Jyen9MuyQyNNNCo1YMCAPCWR03xnxgNjjDHGmNbMWylXz6WkpGQ6Pzk5GYaGhtne5siRI1i9erXYwPLp06fpI1AUCNFpCqiyM3PmTEybNi3Piet034ypi7GxcY5/C4wxxvQscFIuDwwMDMx0Pp329PTM9jb+/v6wtbVFp06dMpUWoPMpqHr06FG2XzQTJkzAmDFj3si2zwkFTE+ePBHBE2PqRLuju7m5cb0xxhgrIFq1V13VqlVRsWJF/PHHH+mjTe7u7qK8wLfffivOCw0NFVNtVMQqO5UqVUKjRo3yVI6AAicKwCIjI99YVUfdR7lX1BZVCmcxVhDofRgXFydyASl4or8Lxhhj7/89r7UjTmTy5MmiZhMFULVr18acOXPEiNGwYcPSrzN+/HicPXsWN2/e/CBtoqlD+sKioIl2qmZMXaiuGaHgiYrD8rQdY4zlP60aHmnfvj02bdqEzZs3i/pNNLJ04sQJkSCuRMdzmrojFODkpWr426Smpor/TUxM8u0+GXtXyuCdRkAZY4zp+VSdJg7h0fYvlN9UtGhRmJmZqa2NjPH7kTHGCn6qTqtGnBhjjDHG1IkDJz21e/dunDp1St3NYIwxxrSKViWHs/xD+/xRJfV69epxtzLGmLahLJuESCA+DIgKBKRUwNIZMLMBTG0AE0vaeE3drdRJHDjpIWVB0NjYWCxcuFCc17t3bzHHu2vXLgwePBgHDhwQZRY6d+6MK1euiOKKjRs3Tr+P8+fPi3pYGWtkUbI+3TfVyipevDgaNmyYa3mGZ8+eYc+ePZlWRVIb1qxZg08//ZQ3VGaMsYxSk4G4MCAuFIh6AcSHAynxgIEJoDAEwp/KwZKRuRxAWbkC5naAGR1yz9thquPASQ8FBASIoIk2PL5792560HP79m2MGjVKbIZMyXE0ItWmTRtRN8vKyipT4EQBz6FDh9IDp4cPH6JVq1ZwdnZG2bJl8csvv8DS0hIHDx6EtbV1tu3YuXMnFi1alClwOnfunCg+OmTIkALvB8YY0wpJcUDIAyDiOZAUDaSlAsYUHNkCxllqttFlFEwlRAHRr2hoCjC1BtwqAPZFaRsOdT0LncGBUwEVIlTXUnSFCkOzNLpEZR0oMMpaCJSqn1MwNHbs2Dw9Nt1n165dxXY1yjINFGjNmDEDP/74Y7a3uXr1qihImtHly5dF4EUjXIwxpveiA4GA60B0gDxyZOUGGOby+WhgCJhYyQfllB6NUD07I//vWh4w4ZqD74MDp3xGQRONzqhDTEyMGOV5X3nZEJk8fvxYjBS1aNECv/32mwge6UAVrHNLQKfAiYKtjGhaMGMwRdvZ0LQglXqgwqeqBIaMMaYT03LB94Cg2/Iokl1hOSjKK/rMtHSSc56C7sjTe+6VAGvXgmi1XuDAiWVC1abzWiDUz89P/P/ixQuEhYWln+/t7Y0aNWpkexsakbp165YYkcoaOA0dOjQ9aGrQoIEYBYuOjkbNmjWxatUqfsUYY7qN8pgCb8g5SxaO8pTc+6KpPQq+YgKBpyfkkSfHkoAhhwF5xT1WANNlNPKjDgW15QsFU8oK6UoZpyOVgRYllVNwo4r79++L4qG092DGAOzBgwfpI05//fWXGLXat2+fyMEqV64cbty4gQoVKuTTM2OMMQ1CG8WHP5GDpsQYwNYr92k5EnwXOP4TYFtIDoZcywFOPtnfjkasbDyB+AjA/4IcoLn7yjlQTGUcOOUzmkrKj+mygkbJ36rmYtHI0ZEjR9JP0wgQJX0rnycFNMWKFcO8efNE7pQSBVs0jVeyZMk37jMwMFD8nzFxfNKkSWKKTxlMXbt2TUz/EVNTU5EzRdN7HDgxxnROcgIQcA0IfQAYWwJ23pnLCcQEAQFX5etQYFS2vXw+BUxxIUBsEPDysnyegTHgXEoOorzrAI7FMz8WrbSjEajQR3JJA+9agEX+bUWm6zhw0lN16tTBd999JwIeysmi5O6c9O3bVySR9+rVC76+vqJkQXh4eHrgRCUH1q5di48++giNGjUSwQ5N2VFJA1qll13gRMEPBU3du3dH/fr1RSAWHx8PLy+v9BEsGpGigEmJjtN5jDGmc/lML68AIfcBa3c5qFGisgNnF8sBk1Js8H+BEyWBt/wBCHsCvLoJvLoFJES8Pn5Tzm3KGjgRI1PAvjAQ6Q+8uAwUqZv5cVmOOHDSU8OHDxcbIlNOEdVjoqmwwoULY8SIEW9ct3Tp0uJ6f/75p8g7olVyVG/p3r17mQIxmmajDZhplKlQoUJiqq1UqVLZPj49NiWO0/Xpsal8Ad2O6kspUXuU5RIIlUv4+OOP870vGGNMbSjxm1bNUdBE02gU0IjzU4BbfwPXNgGpSYDCQB5pcq8IeFTOfB9iiq48UKadvIou6qUcNAXdkhPBlWiKjqbllEnmdJ/0mFTmgKYHPatxuQIV8Ca/KuBNftWDAjrKd5owYQJCQ0NFkEXBGtWfYtnjTacZ0yIU5FDAQlNwVKzSOEOe6pmFwP198nH3ykDtEYC127s/FgVie8cDNPtXf1zm+0qOl6cCacrOqQT0UVQeNvnlbyCmsWjUav/+/aIAJ5UjOHr0KAdNjDHdCZqo3EDgdcDCOXPQRGgqjhK4q/QDijVCUGgEHj+8hcfPX+Kx30s8fh7w+v+XCAwJQ6UyJdCsbjU0q1sVdaqWh1mGNAeBRpVoWi45FvhnFFBrKFCsiZxHRVN0NBJF04H0P5cqyBWPOKmAR5yYtuARJ8a0RNhj4Pk5OUeJkrVfXALCHgEVuv13ndRknLvxAF9OX4Azl2+pfNdmpiaoX903PZCqVLakvP0VjSqdnCNP4ZEi9YFaIwDT17UHKbCysJfPp9woPRLFI06MMcaYhhIJ2ZcAIzM5aHp4CDj9i3yZW0WxIu7lqxB889NSrN22P33FdiE3ZxTz9kAxLw8U83YX/xcv7AkHGwucvXIbh85cwaHTl8RtD566KA6kSnkfbJg/GaWKecuJ5De3AlfXA09PyqNerX+Si2RSYnrkc3nkqVANrvGUAx5xUgGPODFtwSNOjGk4GvV59q+8ko7yjB4dBU7NlfeUK9kCCRX7Yu6qf/DDknWIjYsXN+nXuTV+GDcI7i5O8hRfSiKQHCdPu1FyuTLZOy0VksIId/3DcejCHRFIHTlzGTGx8bAwN8OCKaPRv2sbeQcGCphOzJYLYnpUAZp/J98H3Xf0SzlR3LUs9EUUjzgxxhhjGoa2O6G8JUrGpvpLT44Dp+eJoEnyaY2/w8tibOtBeOov17mrXaUcfvn2c1Sv4APEhwHhz+QAi1be0RSfdQnA0hEwpWRmSWzsq4gJQhlTS5TxsseoDjUQEBqDPpOX4vDZa/jsm1k4cOoCfpv+FeyozlPLGfJIV63h/7WR7tvcQU5ap4rltp7q6y8NxSNOKuARJ6YteMSJMQ1FIzk00hT5Qi5u+ew0cOInQEpDYuHG+PgPfxw4dUlc1dPNGT+NH4qe7ZpCkRgplxGwdJaDGJraowRuCpZy2rsuJUkubEmHmCCkRb7A7D+2YtLiLUhJSUVhTzdsmP8t6lTNZRcGKmlAuU+U72SW+yozXcAjTowxxpimoOk12mA3wg+w85JznGiaTEpDgld91Pn5Fq7cfiiSuscN7onxQ3rB0ihNzjeikaVC1QCHoqoXqDQyAayc5YNTCRjEhWH8aGc0quKDnpOW4Yl/IBr0+BzTvuiPb4Z+IrbVEmg0jGo7eVaV851or7xXtwHvmpmrmOs5LkfAGGOMFaRIPzlwokDGwEgOnir2QHzQE9ScfQs37j2Fi6M9Dq6ZC9+ShYCYV0CKEeBSWi56+b6b/NJ2KoXroqalC676FMew75dhw95/MWnO7zjy72X8s3wmLMNvA4e/k8sifDRHnkqksgS0dx5VGLdxz6/e0HoG6m4AY4wxprNouoxWqVHARHvQvebv0hSVf7wpgiYPVycc3zgfvoUs5WRt2ty3eGM5Qft9gyYl2vTXtSxsKrTGuvlTsXrqIFhamInk8UETZkNyrQC4lJETzo98DyTFvK4tJQEh9+QkdCZw4MS0xsCBA/Htt9+quxmMMaaa1BQ5yZqSwmNDgMNTxWa+T/xeokGPUbj3xB/eHq44sWE+Sjsq/sspon3jrFwKZnrMygWKYg3w6aAR2LtwPIyMDLFx5yHMX70daPQ/OZeK9sejqUQKlqzcgAh/edSMCRw4sQLTr18/TJ06Nd/uLyQkRGwuzBhjWoH2n6OilkbmwIlZonZT6Mk/RND0xC9A1GA6uWk+itulyTlFhevI02I5JX3nF1o55+6L+p0HY87XA8VZ435cgmPXngCNJwGGpnKdqcur5ZEquj5NNVKCO+PAiRUcCnQiIiLy7f5WrFiB6dOn59v9McZYgYkKkDfaNbUFTs0To04J5u6oPmk//AOCUbq4N06snwtvyyQ554n2icuvaTlV0GiWjQdGTZiO3u0aITU1Fd1GTYF/kjVQd7R8HdpkmFYC0igU1Z+iZHHGgZO+6tmzp5j2+uqrr1ClShVUrFgRCxYsyHQdSZLwyy+/oHr16ihevDg6duyI69evZ7rOxYsX0b59e5QqVQqNGzfGhg0bxPljxozB4cOHRbBDe87R4dGjR+KyP//8E40aNUKJEiXQokUL7Nu37422TZ48GWPHjkXZsmXRrZu8BcH48eMxZ86cPLUvp/vKrj+++eYbjB49GlWrVkWFChUwf/58EfiNHDkSpUuXRs2aNbFly5Y3bvu250ObFCv7gPp5yJAhCAoKeuPxp0yZIq5brVo1VK5cGT/99JN4jowxLZMUK2/cm5YG3N4htjhJNTBFw2UBeBIYDt/SxXF83Wx4mMUDDsUAr5r/bXvygSmsnLB06XJULFUYwaER6Dz8WyR61ALKd5av8PKKPAJGQV3wXSAxRi3t1CgSe6vIyEj69hL/ZxUfHy/dvn1b/J9JYkzOh6S8XDdOtevmUdOmTSUDAwNp1KhR0o0bN6R169ZJVlZW0ooVK9KvM3PmTMnBwUH666+/pOvXr0uDBg2SbGxspMDAwPTnbm9vL/3vf/+T7t27J508eVLq06ePdPnyZSksLEw8xmeffSb5+fmJQ3JysrRo0SLJy8tL2rp1q7gNPa61tbW0f//+N9o2adIk6c6dO1JQUJA4v3379tKIESNUbl9u95VTf9B90vV+++038ZoXK1ZM+vnnn8V58+bNk4yNjaXHjx+n306V5xMeHp7eBxcvXhTPo2bNmlJaWtobj//999+L9xPdn5mZmbR58+Y8va45vh8ZYx9GaookPTsrSRdXS9KRHyRpio049KtuKz5TqlUoLYWe+0uSLq6SpKdnJCk5QSNemUcXj0j2NpaijYN6tJOk+wcl6fSvkvT4hHx4dFxu84srkr59z2fFBTALqgDm1FyGXEu2AD7567/TM9zl8vnZKVwP6L/7v9M/FQPiQt+83tRI5EWzZs0QGBiImzdvpp83bdo0rF27Fg8fPkRycjIcHR3FqMfQoUPF5WlpaShTpgw6deqEmTNniuuVLFkSAQEBcHNzS78fuh5tKNm2bVsxCkMjN4Tu09XVFatWrcLHH3+cfv3//e9/uHr1Kvbs2ZPetqSkJJw4cSJTmzt06CBGbRYuXKhS+3K7r+z6w8TEJL0NhEbRypUrh7///jv9PHr8GTNmoG/fvio/n6yio6NhZ2cn+p7aq3x8Y2Nj7N27N/16Xbt2Fddbvnw5VMUFMBlTs5CHwPMz8mjT/vGiSvhv1wwxbHs4KpcriSMrf4CdYRzgUhZwryjnEGkCScK+TcvQ5pNhYqR72Q/jMKhHu8zXSYgCUuKB4k3kEgd6WgCTk8P1WN26dTOdrl+/vphOi42NxePHj8UXfIMGDdIvp2CITl+7dk2cpmCRppVoemrevHkiWFBeLzt3794Vyd3Dhw9HkSJFULhwYXh7e+O3334TQVhGNH2YG1Xap+p9KVGQlJGzs/Mb5zk5OYncrbw8H39/f3EdagddTsESfTA9fZo5X6B8+fKZTlNQpnwsxpgWoArfIq/JGpBSkWpghvOBBhj1TzjKliyCA8unwc4wHnDzBTwqa07QRBQKtOo+ENPHDRMnR06Zj/PXbsuX0arABwfkCuIpCUDwfbmop57iApgF5X8vc75MkWXFxLiHuVw3SxDyxQ3kFxphye40jdAkJsqrJ0xNTTNdh04rL6Nqs6dOnRJ5PwcPHhSjPxRs7Nq1SwQI2Y2GkHXr1omRqIyMjDK/FTON3mVDlfapel9K6dVz33KeMu9I1efTqlUrMXpFOWQeHh7iMsrJytrO3B6LMabhaIQp6K48KmNfBEFxBuiwLhWPn0WisJcnDq6YCSfTJMC9EuBann7pQeMYGOKb6XNx4eoNbD9wEp2HfYvLW+bB+dg4ub4TbRVjV1jvi2Jq4CunI0wscz4Ym+XhullK7Od0vXdw40bmIIxGamj6y97eXnyx0xd51tGbK1euwMfHJ1Og8sknn4jpqmfPaANKYNGiReJ/un3GL366HQUN9+/fT0+WVh4yTvWpQtX2FSRVns/Lly9x69Yt/Pjjj2KEj0bpaESPpvkYYzqEah+FPRajMuFRMWjR9yucue0PExtnHF4zGx4WiYBTKcC1nGYGTa8ZGJti9fo/4VO0EPwDg/HFT6sA79ryheeXy6UJCCWKU50qPaS5rx4rcMePH09fBUdTSzRiNGyYPExraWmJAQMGiBVpNNVEAdDKlStx7tw5jBgxIj1I+f777xEaKudc0bQVTZ/RqBPx9PQUQYUyeKL5Y1pRRvd58uRJcR6NuuzcuVPkLeWFKu0raKo8HwcHB5ibm6fnO1FfKfuYMaYjkhPEyjkqEiltH45fvh2Oa3cewtXJAYfWzkVh61R5tMbdt+BrNOUDGyc3bFi3DgqFAhv+OYSziqpyLSqqIP74OGDlKu+3F+UPfcSBkx6j5GP6gqcEZFpuX6dOHZHYrPTzzz+LvBsa3aFkuYkTJ2LNmjXpuTg0/UQjJ3Rbyv0pVqwYmjRpglGjRonLKUC4c+eOCB6U5QgoF4oCHkqmpvuk2/3xxx9o3rx5ntv/tvZ9CG97PjRNSCUZqBAojeRRHhTlOllY0FYGjDGdEPJAbOCbdvEPKJJj4WsRBAc7GxxcMwc+TkaAmZ2c06TqJr0aoGqdhhjUt5c4PnTGCqSV7yJfcHmVXFHcyOz1qJP+jZ7zqrqCWlWn4WgVFyV20xRSfHw8UlJSYG1tne116fKYmBgRFNAvkOyEhYWJACy7xHAaiaLpKZq+Uub+0Ao4Gn2hoCprbg8lRFO+Vda+puvTdelxVG1fTveVVXbXCw4OFlORGc+j+ksU9FhZZa65ktvzITQiRo9BbafVczSFR9dVvmeye3yqIUVF6Wj6VFXa+n5kTGvFhgKPj0K6tQ2Ku7vwIioNtdakYduK+ahWwgVISxIb7FKxSW0TEhwMn5IlEB4Zhd+mjsQQy/3yBsS+PQDf7kD0S6BYI3lDYC3Hq+pYntBUUk5Bk/Jymn7LKWgiFATktJqORlpoxCljwjRdl+4zuyCDAqDsAh0KILIGTW9rX073pcr16D6znufi4vJG0PS250OobXQ5BU2EksQzBjbZPT4917wETYyxD4xGXoJui61VpLty2ZhR+5OxYdEsVCvtBSRFA+6VtTJoIk7Ozvh+qrw/6IS5qxBVqvt/FcVppR0tdKJq4nq2iIWn6hhjjLF3EfEMCH2EuPNrYQAJm28mo2HnwahftQwQ+wpwKSdXBtdiQ0Z+Ad+ypUTC+9cbr8qBYKnWgImFXMsp8gUQFwZ9wuUI9NSmTZveKEfAGGNMRbT1yKvbSLx/BBbxLxEWL2F3XAWs/rQjEOUHOBQH3DS07EAe0EzBgoUL0bBJcyzbtBuDe/yGKhXkwr0C7WFHieKW+jM6rt2vKHtnqk5hMcYYy4KmpoLvQYoNwZ4TF5CaJmH6OWPMm/4tFDEBgJWbXK9JkwpcvocGjZuhZ5f2Ildz1NRfMteXo8T38Cfy/nx6ggMnxhhjLC8oQTr0IVYduIZOS++i4tJ4dBo2DY5maYChqbyCTk2b9haU2XMXwNLCHP9euY112w8AgTeBw9Pk+lVU9DMqAPqCAyfGGGNMVbT8Pug27jzyw8jpv4mzevbph3pVywIJ4YBrWcBKrmWnSzy9vDBp/Fhx/OuZi5F07wDgfwG4vV0uwhz2UG8KYnLgxBhjjKkq/Cnin11B6L4fUdgyCU3rVMU3Qz+Rl+bbFQEcM2+/pEu+HD8RJYp6IzAkHL+ce12/6flZsZExYkOAmEDoAw6cGGOMMVXQlNSrW7i3awHqeSTjj45WWDdvEgwTIwFjC8Ctgs7kNWWH6trNnzdfHP/f7/sQY0/FhiXg3m55X1U9KU3AgRPTCrTxMG3nos+oSjsVaWOMqQEFBCH38e+OFahkE47kVAlpNYbAzcEaoMCJ9qCj5fk67qP2HfFRiyZISUnFD8defx49PEyb3AFRL/WiNAEHTqxAUCVtqnxN/+eHZcuWoXr16tC0YC4/UEBIe9y9zY4dO8S2NowxNYh5hRfXj6N4+HFx8kS8D+o0aQtEBwD2ReXyA3pi3q+LRdHfmdtvIs6iEJCaKKqniz37qDSBjuPAqaDQGygx+sMd6PE0yOPHj0XF8OfPn+fL/WXd+qQg0dY6tNVJdgICAvDll1+Kyt9UbZ02Mp45c2bm5bl5VLt2bSxZsuQ9WswYK1CpyZBe3cblrfPgaqnA02gjNPjsByAuFDCxkkebDPWnLGLJUqXQrVN7cXzFzdfPmyqnU5K4HpQm0J9X+kOiIObJCSDxA06rmNoARRsAxqrtT0Z7x9H2ILTtBwUJ9MWfcUuUjGgfO7oOBS/ZUQYZyu1G6L6U02o0tUQjT7TViKWl5Rv7qmXdT43ui25LewPSNiXUTnrcvn37okePHm88Nj0W7VWXddNcGumix6Zgi34ZxcXFiTYotzzJDQVDp06dQqVKld64bMuWLShSpAguXbok9t47duyY2OCXiol+9dVXKo1SZSw8SnvsUVvpOVA/kYzbytDzo/5Xpd2MsQIS/gz/rP4V7QrRZ7oCBrVHwNjIAEiIBbxq6cUUXVbj//ctNm3ZhrEbb+Czn2vAonRTwMwWiH4hlyZw0t0keR5xKghpyXLQZGgCmFoX/IEehx6PHldF7du3x8CBA9GyZUsxMkRBzZAhQ0QeTcbRlXbt2qVvaktTZRcuXEi//MGDB2jYsKG4nPaqowCCRpgoYKHzSb169USgQYGPMjCaPHmyKMBJwRGN3Pz888/p93njxg3Rnvnz54vAhC4/fPjwG1N1dD/ffPONuC7dDwU7y5cvT7+c2kGX0X3TPnl0X9u2bcP7GjVqFEaPHg13d3cR2DVu3BhdunQR02g5ocBnxIgRYoSKgqLy5ctj79694jIKBu/du4fvvvtO9BMdlNN21HZ6DhT8Va5cOVPfM8Y+kMRoPL96HPF3D8NAocCttGLwrtJMXkVnXwRwKKqXL0WlypXRqnkTJKVI+OqsE+DTCjAyBYypNMEjnS5NwIFTQRJvIouCP9DjvIP169ejc+fOCAsLw/nz57Fr1y789NNP6Zf36tVLjP74+/uL0ZBq1aqhbdu26aNJFER4eXkhPDwcwcHBGDp0KA4cOCACmcuXL4vrXL9+XdyWRmrI//73P+zZswf//vuvCBAogJg3bx7++OOPTG3bvHkzLl68KKbNKLjLSnmbffv2iZGruXPnYtiwYThy5Mgbz/H48eMimOvWrVu2/UCjPtRG5YHQc8x6Xk4oSKNAMCerVq3Czp07RVBII19bt24VwSChPi9TpoyY7lM+Fo2w0fOaNGkS1q1bJ0bdZs+ejUWLFuXaDsZYPpMkSEH38NnEeei9NQbzb9mhdLep8hQdja640pYq2W/srQ8mTJws/l/59wEEBofKZ5o7ALHBOl2aQOsCpytXroiRkQ4dOuDbb78VX9q5oS8d+sL55JNPxKgHjUpkHFXRZ3Xr1sXgwYPFFF3FihUxfvx4/Prrr+n9TNNQ1HcuLi5iRIqCExo9oWBEOSJF01k04kTTT23atBGjWDmhoIHu//vvvxcjSRScFC1aFAMGDBABQkY//vijGCnKCQVONDVWq1YtMUXYvXt3dOrUSbQxo2nTpqF48dyTNrt27Zo+2kMHes/Qc8l4Xk7+/PNPHD16VIwo5YT6iUbElPdTqlSpTKNs2VmwYIEYjaJAlaYamzVrJvqJMfYBxQRh2bKlOHT+FozpM27oTBiaWsg5PBQ0mf83ra6P6jdogNo1qiIxKRmLVqwH7u8Dzi3R+dIEWhU4nT17ViTSUr4HfdnRFzt9+dMXcnYod4SmRe7evYuPPvoI9evXF7/s6csov1Z7abOsOTw0HRQUFCRGPWj6iEY+ypUrl365ubm56E+6jFDgQsEr9S0FMjR1lxt6HWh0iAICCopotMrb21sECVkDYAouckIjRC9fvkTVqlUznU9Tecq2qXI/SjTqlXF0iYLEkydPvnXE6cSJE+jXrx+mT5+Opk2b5nj/vXv3xosXL1ChQgUx4kaB1tuSyel5ZPf6MMY+kNQUPL18BGf3bICpITBz3GD4FC0kr6JzLCZP0+k5hUKBbyZMEse3/LMb0tklwIMDQFyEXJogPveBDW2lVYETfem0bt0aCxcuFCNINM1B0yQrVqzI8UWl6R76YqZpJxoNoSkgmk7ifJH/krqVaDSJ0AgUHSi4zBpg0nWUSeSffvopnj59KoJYSpamUausIz7ZocTrjEEJHWiEK6PckqFphIleW2V7s2ubKvfzPug5UMD49ddfi/dlbmhUjYLKWbNmidFO6rcmTagOSkquzzGn14cxVvDSwp7g99kTsbKdEW6PdsTnfdoD8WHyQhyXsno9RZdR248/RrkyPrgbGI9byZ7ymff3youkonVzuk5rAidadUS/8GmKTolyaWgKg/JBskNfro6OjpnOo2kn5aiFvjtz5swbpwsXLiwSwWlkib7kKfBUojwhytOhy5RcXV3FqAtNtU2dOjU9D0e5cizjl3/ZsmXFfVO+z/ugkS+afqM8qazBTMa2vStK4FauEMzO6dOnRQA/ZswY8ZxVQaN3NP1HuUqUc0WjpdSXyr7KGiTR88ju9WGMfQCJMfh9wSwMLhksTtqXbSKmzMUiHCo9oOdTdBlRv4wfL/94HL3ldfmZZ6eAlAQg4plOJolrTeDk5+cnvlxoeicjmvKhUQ9VzZkzRwRTNWvWzPE6lLRMQULGgy66du2aGC159uyZWBVGieHjxo0Tl5UuXVokjg8aNEiMzt2/f18ESM7OzujZs6e4Dk1PUY4P3Z6m4Q4ePJg+tUer2CjA2b9/v0g+p7whKj0wZcoUzJgxQ9QtotFCagNNn9KUX15MnDhRrLyjx6f7oZwheqwJEybkuR+yJoffvHlTvM+ym6o7d+6cCJo+++wzUc9JeXlu7xHqY3qOFCjRFOPGjRvFdCAFqYSKWlLQ9+rVq/THGjt2rHhNaHRVOaqqzC1jjBWsx5ePIe7CRnjbGiASVrBvMOj1XnSF5WKXLJMevXrBu5AnjtyLgj88AClNLslDVcQpUVzHaE3gpKzSTF/GGVFisqoVnGl10+LFi7Fy5Uox8pET+pKj0SzlIWuwprKURCA5ruAP9DjvgFahhYSEoHnz5iIIoKBp+PDhmfqrQYMGYiquUaNGYhTm0KFD6aNJFLj8/fff4jpUtoACAOXqOLrO0qVLRX+XLFkyvRwBBQSUoE8jVJSTRMEY5agpAzZ6DOpz8esulwKYdDsK9Ci/iFb7UakBmrqtUqWKuJxuT/eT28hRTsnh2R2UaKqX7pv6JuPlVHYhJ5R0T8EZTRfXqFFD5E/RdDGVcFAmsFPARFOdynIElMu3du1a/Pbbb+K+6bnRVB89J8ZYwUmLeoUZE7/EyOryZ4d1kzHytJORuZwQrkeFLlVlbGwsPtvJxL2vA6WHh+TvJ8p10jEK6X1KHn9AtCSeAhj6AqHcEiX65U+/5Gk5fW42bdokckvol3ufPn1yvS59cWXcAoNGE+ixaWl81urVlOz85MkTkceSXsxRCwpg0hQnBRy0eo3pjmzfj4wx1aSlYuF3Y1DV73fU9jJCjHMVWLWaDEQ8BwpVB1zLck/mgH4AF/b2QkhoGCKmFYZtWjhQbSBQtD5QsgVgkrlIsaah73n6YZrd93xWWhM605Qc1cqhJOKMgROdzrq6KiuazlGWInhb0KQc3cipSrZKKHihICYPBSnfG22wqGLQxBhj7E1Pb5zB3X3LMbKVMZJgDKuGn8sJzjaegKPuVsLODxYWFhj9+ef4dspULLuYgLEf14SCVh/SAELMK50qFKo1U3VEOWJE00uEcloocFJOAymnj6gQoxIVXqRgiSpPZ7xegaMg5kNUDVce8hg00VRl1mlPxhjTV1JSPAaP+Bz77ifgTJA5jKp9KhcXpppEbuUBo/+2SmLZGzHqc1hZWuLrna+wV9EEcPOVd7aI8NOpmk5aM+JEaFsKmpajnBk60PEffvghU34JJfZSvSdCQ26UV0LDb5SQSwclyunJriK1vti+fbu6m8AYYxpj1ZK5OHj6CkxNjOHQ6WcYFPEAIv0A90qAtZu6m6cV7O3tMXTIYPw8dx5+XLwGbZrUBczt5SriVNNJR/b006rAiVYiUVLt7du3xQokWt5Oy+EzUq50Ug4d/vPPP9neV8bCjowxxvRXwOO7+GbqTHH8uy8/Q6niheWkZmtXwNlH3c3TKl9+NRa/LlyIk5fv4MaFU6hg9hLwqCJP13HgpD4UMNHhbQERZfq3atXqA7aMMcaYNpHS0jBy+BDs7AK8SHJFu96tgaQ4eUk9raIz5pSGvPDw8EDnjh1FbnGR6z8DRiny/nU0audYUidWJWpVjhNjjDGWn/5asxzOIWdQw9MQH5cAjBQSEPsKcCohJ4WzPBsybDhSJWDttder0/3OA3HhQJycn6ztOHDKJ1pS1YHpON6DkTHVhbwKxJSJ4/FDU3lxjWHVTwEpVc7LcS5N209wd76DBg0aoFTJElh2IV4+w/88kBgNROpGTSftHzNTM5oOpK1dgoODRVVtOs6YOgJ3KgRL70Mq0KksUsoYy9mXo4biq6qJcDA3QZpdURiUaCbn4hSuI69WZu9EoVBg8JCh+GrsWNwLN0Qp+xTg1U3Axg1ILqP1059aUwBTkwtjUVVoKtDJXcnUjRZEuLu7c+DE2Fvs2fE3vh/ZDWc+s5TPaD1bLtJI26pQ4MSb+L6X0NBQeHp6YEgl4JdWZoBDcaDel3JBTPv/dmPQFDpZAFOTUU0kKo9Am+Iypi60vYyRkRGPejKmwpfksOHDsL3t6/p3NNJk4yFvX+VSloOmfODo6IguHTtg/Y6/MKeFOYzCHgHRr4BIfzk41eLZGQ6c8vFLS5V90RhjjKnX2NEjYZEUgsL2VpCMLaGo1AeICwU8qwKWjvzy5JMhw0di/aY/setBCtqXs4KCqohHBwAJEXIemZbi5HDGGGN6Y9+e3Vi+ai3uhqThToX/QdFkEpCaANi4A04l1d08nVKvXj2UKVUSI3fH4feUzkDJ5kBSPBATBG3GgRNjjDG9EB4ejoEDB4rjo/t2Qt06deTcG0r1dSknb7HC8jVJfMiQoXgRLWHRpj1yHrCJJRD+VGyorK04cGKMMaYXRo8aAS/DYIxq4Iwfxg2RA6b0mk0e6m6eTurTtx/MzExx7d5TnL96G0iOk6dFY7W3phMHTowxxnTejh07xH6ly9qZ4dfGibB4ekAuyMg1mwqUg4MDunZsD2sTwPP8VGDn53IxzOhAaCsOnBhjjOm0kJAQDB40EJ/XNEEFF0PA1AYoUk8e/aBVdFyzqUANGT4K0UnAy7BYeSubwOtAlD+Qqp0r0TlwYowxptOGDxsG44RQfN/4dfmBqv2ApBh5WbwG1hTSNXXq1kW50j5Ycfn1FizPTgNxEfKUnRbiwIkxxpjO2rx5M/7asgXzWpnDwhiAcxmgUHXAyIxrNn3QJPEh2HQzGQkpACL9gIincl0nLcSBE2OMMZ0UGBiI4cOHo3kxQ3QtawQoDICaQ4H4cMCpFNds+oD69BuAZIUZ/r7zenou4KrWTtdx4MQYY0zn0NL3QYMGITI8DMvbv953rkw7wNgCsHblmk0fmJ2dHbp37oANN14HSv4XgLgwrVxdx4ETY4wxnbN69Wrs2rULhkZGSK3SH3DzBcp1AtJS5JpNxq/zndgHM3jYSBx4lILwBEke9Qt9oJXFMDlwYowxplMCAgLwxRdfiOPThnZCsbodgRYzgIQowLE4YOOp7ibqpVp16qB06VL4cl8CdqAF4FZRK6frOHBijDGmU0aNGiV2uW9euTDGDpcrhYsRDjMbwLkUYMBffepKEh/Qvx9WX0vGzG03AAsHID5C66br+N3DGGNMpwpdbt26FZ3KGmNvh1gYPT4oj2gkRgAuZQBzO3U3Ua/17NMfhoaGOHf9Hu4/DaBsNK2bruPAiTHGmE6IiorCiBEjYGMKrOhoB8O0JHlEI+YVYFMIcCim7ibqPVdXV7Ro0gBF7RQIOjAHuLdH66brOHBijDGmEyZMmIAXL17g9472sDNKBKzcgNJt5DIEruUAQyrkxNStz6f94GljgHrmDyE9OgrEBGvVdB0HTowxxrTev//+iyVLluCjkkboWiqVMmqAuqOBhEjAsQRg5aruJrLX2nfqgmuhpvCPSoOCtr0JuqVV03UcODHGGNNqSUlJGDx4MOzNgHVdbeQzy7aXgyULRzkhXKFQdzPZaxYWFujUvi0231IWw7wuVxNPSYI24MCJMcaYVps1axZu3bqFZe2tYGecAtgWAip0kzfxdaVNfK3U3USWRZ++A8QWLER6eUUecdKSves4cGKMMaa17t69i+nTp9PEHHzKVwEMjIG6X8pfwpQMThv5Mo3TqGkzvExzxKOwNChSE4FXN7Rmuo4DJ8YYY1opLS0NgwcPElN1Lev4onyPaUCXPwBze/ngVh4wMFR3M1k2qCTBJ927/Tdd9/Ka1kzXceDEGGNMK/3+++84efIUrM1NsGTG11BQYUsjMyAlXg6azGzV3USWiz79B4rpulcxEuLNnOVEfi2YruPAiTHGmNbx9/fH11+PQ68KRngw1g1FbNIAKQ2IDgScfHiKTgtUqFgJBs4l4Tk3Gqufe8mvnxZM13HgxBhjTKskJyeje7dusEyLxpK2lnA1iACenwWiXwFWznLNJt5WRSv07tULqRKw9u+9gKmNVkzXceDEGGNMq3zzzTf498wZrOxgCRsTCXAsCZRqDUipgFsFwMRS3U1kKurV9zMYGBjgzJXbeHHvMhAVoPHTdRw4McYY0xp///035s6diyFVjdGimIG8iq7O50B8mLwXnY2nupvI8sDD0xNNG9TG7l4W8Lw+H3hxSd4iR4Nx4MQYY0wrPHjwAP3790cZJwMsaPN6VKnKp/JWKtZugHNpLnSphfr0+RTHnqWI41LAVSDqhUbvXceBE2OMMY0XFxeHLl26ICE2Cjt628HYIA3wqAwUaySXHHCrCBibqbuZ7B107NYL/zx8HY4E3wXCnwFxYdBUHDgxxhjTeCNHjsT169dRws0G3t6F5UTiWiOB+HDApSxgzXvRaSsrKytUbdASZ/xSoIAkT9fFae6mvxw4McYY02grVqzAypUrRRLxwmljYPrxXKD1LCA1AbDzlssPMK3W59P+2HRLnq5Loy1YIp4DabRZs+bhwIkxxpjGunr1KkaMGAFDBfD9sE5o3LixPDVnZA4YWwLuFQEjE3U3k72npq0+wqkgeU9BRdgjINxPHk3UQBw4McYY00gREREirykxMRGnhrliQl0TIC0ZSE4AkmMBd1/AwkHdzWT5wMjICA1bd8L5F6li30FQkriGliXgwIkxxphG7kPXr18/PHr0CBObOqCWczwUD/YDUS+B6ADAqSRgX1TdzWT5qGfvvhh/KAGN1yYj3qM2EOEHSBI0DQdOjDHGNM6sWbOwY8cOVHQ3xnf1X59ZtZ9cekBUB6cNfPkrTJdUq1UXDxIccexxPA5ffSLX5kqIgKbhdx1jjDGNcvDgQUyaNAmmhsDhQW4wkFIAz6pA0UaUASPnNXF1cJ2jUCjQoV0bcXzbobNASoJGliXgwIkxxpjGePbsGXr27Cmm6vYNLQpHRSRgZvu69ECoXHrAxkPdzWQFpGOX7vBxNEBD6TTSrv0JRPpr3HQdB06MMcY0QkJCgkgGDw0Nxax2bmjkRMnBCqDuF3IyuH0RwKW0upvJClCDxk3h4WCFT8sDEm3cTDltidHQJBw4McYY0wijRo3CxYsX4WBrhX69ewHG5nJeEwVMptbyKjrKcWI6y9jYGN5Vm+NZRBoMpWTg5WWNW13HgRNjjDG1+/3338WB8lw2zvoCLtU7AO0XAz6tgeQ4OWgyt1d3M9kH0KFTV2y/J+9VJwVcB6IDoUk4cGKMMaZWNMpEW6rYmwELRn2EFi1ayBdQoBTzSq4MbleEXyU90fKjdtj9SFRzQurLa/Kmv0lx0BQcODHGGFObkJAQdO7cGanJiTg62BXD3a4AQbcBKU1ODLbzAtx8ufSAHrGwsIC1TwMEx6bBKDUeeKlZxTA5cGKMMaYWKSkp6NmjB54/f45VXexR0T4eitRkeTsVKnJp4Qh4VgGMzfgV0jPtO3XBP/flvesQeF0eedQQWhk4XblyBfv27YOfn1+B3oYxxljBGfPllzh0+DC+qG2OT8rQhq4KoP44ORHcwFiu3USlCJjeaduhM/65l4YHoakITbWUV9elJEITaFXgFB0djYYNG6Jly5aYMWMGSpUqhalTp+b7bRhjjBWs3xYvxoKFC9GsmCHmtHi9SW+VvoBLGbnwIY00Wbvyy6CnHBwcEOtSBT4LY7HygS2QGKUx03VaFThRJdkXL17g3r17OHnyJHbt2oVp06bh2LFj+XobxhhjBefIoUMY+fnnqOZhgF29bWEACSjeBCj1kVzk0q2CXIKA6bUOnTqL/7cfOCkXwYwJhibQmsBJkiSsW7cOn332Gezt5SWpTZo0QdWqVbF27dp8uw1jjLGC8+D+fXTp2hmpqalY1KUQTBUpcqBUYwgQEwA4l5ZHnRTyqiqmv9p37i7+v3j9DsKf3ZYXC6S+zntSIyNoCX9/f4SFhaFixYqZzq9UqRKuXbuWb7chiYmJ4qAUFRX13u1njDF9FxEejnYftUJ4RBRqVSwF30FzgLtbgQpdgdggwL7w6xV0hupuKtMAXl5eqFW5PHY0fQr7a4sBRyfAqzpg5aLWdmnNiFNkZKT4XzlypOTo6IiIiIh8uw2ZOXMmbG1t0w/04jHGGHt3KcnJ6N6pHQL9nsDLzQnblv0IM0sboGp/ID4CsHQGPHgFHcusXYeOOP+CFg4AeHEJiA2BumlN4GRiIicPxsVlLoIVExMDU1PTfLsNmTBhggi6lAdeiccYY+9BkjBm+ADcuvQvrg61xoWpDeDm5CBfRlWhaWsVsYLOhruZZdKxSw9suytPz6VSFfFIPyAtDeqkNVN13t7eMDQ0fCOIodNFixbNt9sQCqpyC6wYY4ypKC0Vv/38HTZvXI8T/SxQxE4BxNySN+2latAKA6CQ+qdfmGYqXaYMbiV6IDUtDIZRfkDYUyAlHjCxVFubtGbEyczMDI0bN8bWrVvTzwsPD8fhw4fRunXr9POuXr2KI0eO5Ok2jDHGCkBKErYs/xnffjcdB/tYoJSToTwl12KGXJNHSgUKVQNsPbn7WbZo78JGrdrj5PPX03W06a+aKSRaeqZF+xnVr18fvXv3Ru3atbF06VLEx8fj/PnzIkgiAwcOxNmzZ3Hz5k2Vb/M2lBxOuU40bWdjw0PJjDH2Vokx2PTbbIyeNB27epijuqchJHN7KFrNkqfmkmIBrxqAY3HuTJarc2fOYMOXjfBLKzOkOZSAwdAT+T7ilJfvea0ZcSLVqlUTgZC5uTn279+Pdu3a4dSpU5kCoMqVK6Np06Z5ug1jjLF8FBuKdb9+h6nTp+N0fws5aDK1hqL5dMDESgRVosClQzHudvZW1WvWxL+htth8MxlXU0pA3bRqxEldeMSJMcZUFOGH1Yt/Rv9JC9C+lCG2dbeAZOkCRbOpgIWDCKrgWRlwKcu1mpjKhg8agCW/r8TATk2xfOMOHnFijDGm5WilU/A9rJg/XQRN9JvcrcpHSKs7BoqP5gCWTkBssFzs0pkLXLK86dilm/h/x9ELoniqOmnNqjrGGGMaihK9X93Gb4t/waNjG+FpDXTo0Am/ThktknuRHA/EBAIu5QC38oCBVmWJMA3QqElT2NnaiPfT06fPULxMBbW1hQMnxhhj746KVwZcxeLflsH24XbMbm6GLxo6waP/cDloSowG4kLkbVTcK3JVcPZOjI2Ncf7YfhQzj4ZhUfXmxnHgxBhjLO8oPTbSH9LLK1i4bCUqRR5AfV9jpEoKeDTsD4WRCRAXJq+ec68MuJbloIm9l5K+NeSteYwtoE4cODHGGMsb2mg16A7SAm5gweLf0Mn6CrwKGyE+zQhmLaZAQcnfVBEc0n8lB3jTXva+aIrX2g3qxoETY4wx1VEpgYCrSAq8i43L5mGYpx9MDA0QJtnAoeOPgG0heVsMGhUQxS0Lce8ynZLnwOnVq1c4ePCgqNBNVbipYFT58uXRvHlz3gyXMcZ0GY0ivbyC6CA/9Ji4FFNKPYOJoSGeGZVA4a4/AEamQMRzeQWdZzXAylndLWYs36m8tOHWrVvo2rUrChUqhM8//1wUkfT39xdVur/55hsUKVIEbdu2FRW5GWOM6ZDkBIA2WH1yAkEvn6PJiLnYc/Iy+u5Kwx3HVijcax5gaAREPJO3Tylch4Mmpt8jTlu2bMHw4cPRr18/UYXb19dXXi2RwZ07d7Bx40Z06NAB3377LYYNG1ZQbWaMMfahEsCjXgKvbgHRAQi6dwHrN23GxRshcHKwxeols1CmYll55RzVaHIsAXhUBkzUm7zLmNorh9O+bzTSZGdn99Y7jI2Nxb1791ClShXoCq4czhjTO7QaLuguEPIASI5DyLnNcAo5Ky7qutMUM36aD5+iheRVTqnJgGs5ueSAobG6W85YgX7P53nLFRpZ8vHxgaGhIfQFB06MMb2qAE7J3a9uy0FR6CMknP0dZmmx4uLlty3w0ZcL4eHiAES+AMxsAI9KgK0Xr5xjWqtAN/nds2cPOnXqhPj4+Gwvj4uLy+tdMsYY0wQ05eZ/AXh6Cgh/AunSSuDfX0TQdC8kFV9d9ETXiSvg4WAlJ4FTPlPR+oCdNwdNTG/kOXCiPKegoCA0bdoUYWFh6ec/evQIn3zyCaZOnZrfbWSMMVbQo0zhT4HHx4GQ+4C5PaSTc6B4eQWJKRImH03Akqhm+HHuItgZJgJxwYBreTkJ3NyeXxumV/IcODk6OuLIkSNwcXFB3bp1ce7cOZEIXrp0aTx+/BhdunQpmJYyxhgrmLpMLy4CT08DKQmAXWEERiVh+qlUHH6cgkrL4uHebCTmT/4cxnEBclFL79pyEjiVH2BMz7xTAUxzc3P88ccfqF69OmrVqiXqOO3YsQNt2rTJ/xYyxhgroC1T/IDAm3Ku0qPDgGdVXH4agfZD/gf/gGDY21pjy8LZaFKpMBDlJ+cxuVbgUgNMr+U5cIqOjsavv/6KOXPmiNEnqu108uRJuLu7F0wLGWOM5a+kOLFlipiWo1VzV9cCMUGIf3gCzWaHITwmEaWLe2Pn4qko4WQo5y951QYcivKqOab38hw4LV26FIsXL8YPP/yAgQMHwsjICNOnT0fDhg1FvacWLVrofacyxpjGigoQW6aI5O77+4BHR8TZ4akW6LY2COExqWhZvzo2zRgMO2szwKkE4FRKXj3HGMt7OYKHDx/C09NTTNdltHLlSowYMUIUwWzfvr1OdS2XI2CMaT36qA9/Ary4Ary4DFzfCMSFiot2vXREj1VPEJsMfNGnHWaP7AAje2/AtSxg7c4r5pjOi8pDOYI8jziVKFEi2/P79+8PV1dXseWKrgVOBb6VQWoS/5pjjBXsqrngu0DANSDmFXB2oTg7ycwZn+2Mx7p/n8DUxBirp/bDp13ayIUs7YsCRib8qjD2LoETDUpl3WIlO5QcTpv9qnp9RkljL+UtDQrX5V91jLH8R1W9aZ+5oNty6QDbQkDxZngSHIu6M04jICwGHs722DbnC9Ro0FwOmszfvksEY/pKpXIECxcuxBdffCE29c1NYGCgqOPEtZzyOHweFwYkRuXlVowxploSuN954PY2wNgcMLMVP2xn33ZBiXEHRdBU27cELm75FTXa9gW8anLQxFh+jDj17t1bbNxbvHhxNGjQQBxo2xWaD6RVdg8ePMCJEydEfaeePXti1qxZqtwtU0qIkIMnM1vuE8ZY/kiIBF5cAq5uAG5uFQUrY+tOwKCJ87Bx5yFxlc86Nsai2T/A1KsCYGLJPc9YfieH+/n5YdmyZdi7d6/Y+DcxMRHGxsYoU6YMWrZsicGDB+eYA6XNCio5/PqRv2C4fwKczdPgMmAj4F0z3+6bMabHYoIBv7PAuaXA05PirDCHKmi04BFuPPSDkaEh5k8chuFjJkBhw8nfjEUV5Ca/Wfels7Cw0PkeL6jAafFPUzE8bp58ot0vgG8PwNgs3+6fMaZn6OM84hnw5CTw7wIg5B4kKHAi2RfNZ59BcnKKyGfasHQuGrbryZW/GfsQm/xmpA9BU0Fq2bkPzvmniuOx944D8eHqbhJjTFulpshVwG/+DRz+TgRNaYam+PaKKxr9cFIETR2b1cL1c8fRsGM/DpoYe0fvFTiRGzduoG/fvqhSpYrYeqVdu3bYs2fP+96tXqCcsfMRcl5T2P1zQGyIupvEGNNGyfGA/0Xg5WXg0mogNgixBjaovyoBM/65DwszE/z+w1hs3bEbjkUrqLu1jOlv4BQQECAqhVerVg3r1q3Dtm3bRD0nKoS5b9++/GulDjMq0VT875L6Egh7DKTJI1CMMaYSWljy7F8g9D5g44G46iNxLdoeRWa9wL+Po1GtbFFc2b8en439HgoLB+5UxtSxya/SmTNnxCa/o0aNSj+vZMmSeP78uUggb9Wq1fu2T+fV6TgYzzdtg7etARIenYJZ0foAf7gxxlQR6S8HTbSFSokWCAyNQKtBP+HanWeilt6Ezzpg2nffw9itLGDw3hMMjLH3HXGiFXQXLlzAo0ePMiVY0VQdlStgb+dbvQ6Ov5Sr8768fkz+9cgYY2+rBB50F7i5DTg0BTj7G17eOoV63Ubi2p2HcHW0xdFVM/HD3MUw9ijPQRNjmjLi5Ovri/Hjx6N27dpi7zoqTRAUFCRqOQ0ZMiT/WqnDFAYGiHOrhYOPjuGRQTyGRr2QN9VkjLGcVs5RFfBrG4GLfwDJcUg2tsHgSXPx6Fk4ing44eDaX1CiTlveyomxAvBe5QiUUlNTRdXwpKQkeHh4wNTUFLqkoDf5PX14L+o1awNbawsEHVsBk3JtAVOrfH8cxpiWo49r2nPu+Gzg9t/idLS5F2r88gx3A2JQvkQh7N/4OzwqNeZ95hjTlE1+jx49KlbSVapUSYw42dnZwdDQEJ6ennm9K/Za7YbN4Opkj1ch4Th25jJaFKvDgRNj7M2gKeQ+cPQH4M4/4ix/q0ooP/UcImMTUbuiD3Zv2QD7YpV5ao6xApTnwCkiIgJLlizB/fv3kZaWhiJFiqBixYoikKL/adrOzc2tYFqrowyMjNG+VVP8s30rnl46BLTvBth5q7tZjDFNQqtuL69ND5qumdVBtQkHkZKSilb1q2LLn3/C0q2YulvJmM7Lc3J4x44dcefOHTGs9e+//2LcuHHi/OnTp6NHjx74+eefC6KdOq9buxbw/9IKg70fIS3gJpCSpO4mMcY0KWiiOk2eVQGf1jinqILK/9svgqae7Zpix659HDQxpumr6iwtLcXo0vDhw7F9+3YcOnQI9evXx+TJk/O3hXqiftueuBwoH39y9h+uIs4Yk4U/BfwuAIYmkCwcMeWsGWpNPQZKTx3RtxvWbdkFExsn7i3GPpB8K+zRqFEjFC1aFEeOHMmvu9QrJhZWeKKQp+cSn18G4kLV3STGmLpF+AGn5gEXliPVxBrDJs3BdwtWi4umfDUcC1ZsgIEJ72/JmEYHTteuXcOBAwcQHBz8xmWU23Tp0qX8apvesa/WSfxfzDQcUshDuVYLY0w/Rb0ETs2Xt1Dxv4Al33+BpRv/EYUtF/84CVNnL4TC0FDdrWRM7+Q5OfzixYtiS5XExERRekC5uo42/KWk8VmzZhVMS/VA7U7D4Td9HrxsDfD0/C4UKVYfMLdXd7MYYx8aFcI9swi4uIKW02HbMxuMWncLJsZGWL94FroM+AJQKPh1YUwbAqfPPvtMbOpLCeJXr14VBwqmaASKCl/26dOnYFqqB6wc3XAizhletqEIvXMSRejDkwMnxvRLUixwfjlwdrEImv5+bIYua/1hY2mOHasXolHHflxugDFtqxxuZGSEChUqiAMHSvnLuFRzIGoTPKQAICoAcCyez4/AGNNYtJr21jbg1FxASsO2h4bosj4ILo622LdxOSo17cxBE2Nqxrs+apgqnUaj744ElF8UhSf3bwFJcepuEmPsQ6CcxpdXgENTgZQEnPIHum8MR7FCLvh31wYOmhjTEBw4aRhH79LwMyuDsHgJ2/afAOJ501/G9GYrldAHeObdGRcDgY/XR6GCTxGc3rURxWq04pEmxjQEB06axsAAnT5uK47+ffQiEBui7hYxxgpaxDMg4BoeBMah9ucrUH1pFIoWKYzDf6+Ba4WGHDQxpkE4cNJAHTp2xtg6JphZ2R/Bd04DqSnqbhJjrKBEvwKO/wT/x/fQ5LNJCAgKRYUShXBg4xLY+dQCDLjkAGNanxzOClahkhXwaVVrVHBIxonjW+Fcsytg5cLdzpiuiY8ATv8KXF0P2yQgNToapYt44NC6eXAs2xAwNFZ3CxljWfCIkyYyNEKYTTlx1DT0tvyLlDGmW1ISgasbIZ1bIk4uvpAIS3tXHF49Ey6+zWg7AXW3kDGWDQ6cNFThJgPE/752cYh4epWn6xjTtWTwB4eRduR7KKRUbLuTjKV3bHBkxVR4VG4BmNupu4WMsRxw4KShitTrCr9oA5gbK3D94CZeXceYLgm+i5TdX8EgORaXA1Ix/l8LHF46EV5VmwPWbupuHWNMlwKnDRs2oGbNmihSpAjatWuHmzdv5nr927dvY9CgQShXrhwqVqwotosJDAyExjOxwDPjkuKoYeAVIPbNvQEZY1ooNhQpe76BUcxLBMemYeABE+z69WsUrdIEsC+q7tYxxnQpcNqyZQv69esnAqHdu3fDyckJjRo1QlBQULbXT01NRffu3VGrVi389ddfWLFihQi0mjZtivj4eGg6t4b9xP+V7GIQ9fQakJaq7iYxxt5HcjzSnp3F/WvnxMlh+4G1s8bCp1ojwLUs7z/HmBZQSBJNtmsH2lC4evXqWL58uTidkpIiNhqmUaQpU6Zkext6erSbuNLDhw9RsmRJHD16VARdqoiKioKtrS0iIyNhY2ODD0WKi8DDCUVxIzARBhW7o8PnPwFWzh/s8Rlj+V0Z/BLGTpiM+ev3oZWPKcZ9NQYNm7YAvGoBxmbc3YypSV6+5w206Uldu3YNzZs3z7RnHo0enThxIsfbZQyaSGJiovjfxMQEmk5hbou1qe3Q+c94rDl0g6frGNNmoY+waOFCzFm3D6kS0KvfEDSsVwfwqMxBE2NaRGsCpxcvXoj/XV1dM51Pp1++fKnSfdDo0zfffAMfHx8xcpUTCq4oUMt4UAuFAp27dhNH9/57AzEv7/F0HWPaKCYIj1f0R/LVP2FsAMz4oi96taoFuPsCZrbqbh1jTFsKYI4ePRqbN2/O9Tpnz54VieDKGUUaZcqITlMukyq+/vprMTp1/PhxGBvnXFhu5syZmDZtGjSBb416KOHtJj54T+7fjtYl6/F0HWPaJCkOD/+cjBIJN/BFLRMku1fB2N6NAJeygK2XulvHGNOmwOm7777DhAkTcr2Os7Nzpv9DQjLv3RYcHAwXl7dX1Z40aRKWLl2K/fv3i1yp3FCbxowZk36aRpy8vNTzAacws8UvXQqjjXUcTvsdlqfrOM+JMe2Qlgr/42vh/GAjYAr89cwBXwwdBIV9EcClDCeDM6aF1Bo4USIWHVRBgRONPJ06dQrt27dPP//kyZPo0KFDrredPHkyfv31V+zbtw+1a9d+62OZmpqKg0ZQKFC0Xjfg2jRUtI1C3ItbsHAuxftXMaYFwh9eRNiOb1DICbgaYoSWg6bC2MoWcPMFjDQ/z5IxpsU5TmTkyJH4/fffcenSJTE9N3fuXPj7+2Pw4MHp1/nqq68yrZajKbf58+eLoKlOnTrQRqVb9seLaMDKRIFrBzYC8eHqbhJj7C2So4JweEZn+DqlISIBcG03BTbmRoB7RcDCgfuPMS2lVZv80vTZq1ev0KBBA6SlpYlRKKrtVLp06fTr0FJC5XReWFgYpk6dCjMzM3Tq1CnTfVHQ1atXL2gDhbkdHqIIPPEUaf4X5Ok6Syd1N4sxlpPUFCz7ujtGFI0WJ8PLD0BRJ1vAuTRA03SMMa2lVXWclKh+U3R0NOzs7N4oN0D5SMnJyXB0dBQJ5RRoZYemCM3NzTW6jlNGt7b+hHI3ZiAqUYJJ91Uw8+0AGGjVgCFjemPVrz9g529TsOJjc4Q5VEGxFoMBSxfaS4lLDzCmgfLyPa9VI04ZV9LZ29tne1nGJ0xBlZubbuz7VLb1QLw8MwMeVgqcO7gRNUs2BCwd1d0sxlgW54/tx9BxU5GYlIK6LVpjTMtOgIER4FGJgybGdAAPWWgJhbk9Hkje4nji0/NcDJMxDRTo/wydenyCxKRkfNysLr4YNgBITZLzmnh6nTGdwIGTtlAoYFN3IDpujkPXDcFIDH4ob+HAGNMISYmJGNC9HXZ8nIBB9d2x9ueJMIh5BTgWBxx4817GdAUHTlqkYss+OB9mjaCIeBw6cpJX1zGmQUYPH4RP3R+iqochFrS1gY0iBjC3BVzLcfkQxnQIB05axMDCHp1byaUWthw8A8RlLgbKGFOPZYsXIvTsRvQob4w0KGDaYDSQkgi4luctVRjTMVqZHK63FAp069Qejs/3oLXbVSQF3IWJY0leXceYGv17+jS+n/glrgwyE6cNfLsBxpaAYzEuPcCYDuIRJy1Tu0kbDKlmihoeCtw4shlIiFB3kxjTW4GBgejcuQMWtjSGk4UBJPtiQIlmPEXHmA7jwEnLGFo54W6ypzie+PgMEMvTdYypAxXh7fdpH7Rwi0T70saQFIZQ1B4BpCTxFB1jOowDJ21cXVe1qzha1iIUKa/u8+o6xtRgwa+/Yv/BQ6jjbSz/aVb6BDAy5Sk6xnQcB05ayPfj4XgVK8HOTIEbhzfydB1jH9j169fx9fjx4nhalQFA8++Bog0AMxteRceYjuPASQsZ2bjgzuvpupgHp7kYJmMfUHx8PHr17ImkpCS0bVAFQz/tCriU4VV0jOkJDpy0kUIBx1qfiKMVrMIQ9+I2T9cx9oGM//prvHp6Bxu62mDl9JEQu2VGB/IUHWN6ggMnLVW+7VA8CFdg251k7N+3j4thMvYB7NmzBwsWLsSvrczQsyzgdPN3IC6Up+gY0yMcOGkphaUj1qS2xYB/EvDHHt67jrGC9urVK/Tv3x/tfIzQs4IxoDAAqnwKJMfyKjrG9AgHTtpKocAnvT8VR/eduYHgJ9d5uo6xAiJJEgYMGICEyCAsb28pn1muI2BiKRe5pANjTC9w4KTFSlesgapli8LXWcKhHX/ydB1jBWTRokVimm5uSwu4WkiAtTvg0wYwtgBcyvJedIzpEQ6ctJmZLRZ29sSlwVawCzzBq+sYKwC3bt3C2LFj0biIIT6r/HqXqlpU6DJWDposHLjfGdMjHDhpM4UCJZv1E0fruibg8dWTPF3HWD6XHujevTsSExMxt52jfGapNnKwZOsNOBbn/mZMz3DgpOUcq3dCcIIhbEwVOL9nPU/XMZaPxowZI0acXB1t4dFxOlChu5zbZGQiJ4QbylXDGWP6gwMnbWdmhxDr8uKoeeh1SDGv1N0ixnTC1q1b8dtvv4nja6cPg0uR0kClXkByPOBUGrByVncTGWNqwIGTtlMoULjFUHG0SaEUXDy+l6frGHtPz549w8CBA2FqCGweUQ3NmzWTL4gJBGw8AGcf7mPG9BQHTjrAonxrhCQaw9pUgWsHNwPxYepuEmNaKzk5GT179kRERASWdXVDN6f7wKn5QGIM/VKRp+hoM1/GmF7iwEkXmNkhyrGyOGoafhfJES/V3SLGtNbUqVNx5swZNCphgT6l4uUzC9eVV606lwKs3dTdRMaYGnHgpAsUCni1HYfOf0vo/3c0Du35B0hLVXerGNM6R44cwcyZM8UU3Y7ejlBAAoo2BByKyTlNTj7i740xpr84cNIRxkVqwrNCPaRKwLptvHcdY3kVFBSETz75RFQJ3zbEBzZSpBjNReU+gJQqT9GZWHDHMqbnOHDSFaY2+KTzx+Lo9qMXEPPqibpbxJjWSEtLQ79+/RAYGIjuNT3Ryvn16tRaw4EUWkXnA9gWUnczGWMagAMnXaFQoEbtuvijsz0ufWaMnVs28HQdYyqaP38+9u7dC3MzU6xqbyZP0RVrBNgXBSydAZfSPEXHGBM4cNIhCjtvdCxjjNJOhnh4ZicXw2RMBWfPnsX48ePF8blf9YFZg9GAe0Wg0qcZpuheb+zLGNN7HDjpEnNbSIXri6OljF4i8PEtdbeIMY0WFhYmtlRJSUlB19YNMKRzY8CjEtD8e3kvOqeSPEXHGMuEAycdY1+3r/j/o5JG2LJhNZCarO4mMaaxeU19+/bF8+fPUbZYIfwxsgkUVNySxAS9nqIrw1N0jLFMOHDSNUUaIFKyhKWJAo/O7YYUHajuFjGmkebMmYNdu3bB1NQER0aXgdXFX4D7e4GkOCAthafoGGPZ4sBJ15jbwrRUU3G0kUs0/j28R90tYkzjnD59GhMmTBDH/xrfBq5h5+XFFBaO8rYqPEXHGMsBB046yKzG6+k6HyNsWrcaSIhUd5MY0xghISEiryk1NRWDO9ZHW7OLAK2iK9lSXkVn5cJTdIyxHHHgpIu8aiLcsgQWnE/CzuMXEfb8nrpbxJjG5DX16dMHL168gE8xbyxqnAhFYhRg6wVU7v16iq4cr6JjjOWIAyddZGoNuy7zsfqJG56FJ2PNqhVAaoq6W8WY2s2aNQv79u2DmZkpjo6tBqPQu4CBMdBgHBAfATiVkIMoxhjLAQdOOkphWwhDOjUUx5du3AWJVgkxpsdOnDiBSZMmieNrJ3SDR9BR+YLqAwEjU3kvOpeyvIqOMZYrDpx0lbUbPunQAm1KmaOhYzBOHtyp7hYxpjYBAQHo0aOHPFXXoQU6t2kCVOkLFKkPFK5HPzXkopdc6JIx9hZGb7sC01LG5rAxM8TuHsaITzbC52vWoEG7XmIajzF9EhcXh48//lgET2VKFsPiMZ1ElX04FgdSkoAof8CjCqCs4cQYY7ngESddVqwx4kycYG6sgHnQZYQ8vaPuFjH2QdEI06effoqLFy/Cwd4OB6Z2gJW9M2BoDEgSEP0SsC8i70XHGGMq4MBJl1k6waJCO3G0X0VDrF65gjf+ZXqFcpq2bt0KY2Nj7P15EAo9Xgcc+BaIDQHiQuQRWLcKciDFGGMq4MBJlxkYApV6I1UyQBV3Q5w6sB1STLC6W8XYB7Fq1SrMnDlTHP9jxheoEbUPSEkAzGzlZPDkOMDdF7Bw4FeEMaYyDpx0nUsppHlUFkebuUfh2L4d6m4RYwXu+PHjGDx4sDg+cdQA9Ha+DUQ8A0ysgHpfyqNNzqUAuyL8ajDG8oQDJ11nag1j327i6CcVjPHHqjXyXlyM6agHDx6gU6dOSE5ORtePW+O7JmbA05PyhfW/AqRUwMoVcCkHGPBHIGMsb/hTQx+U+QiJxrZ4FZOG61fOI+jJTXW3iLECERYWhrZt24r/a1StjDWjGsDgxmb5woo9AScfQGHwuvSABb8KjLE848BJH1i7wbTNTPQ57obrgSlY9Qcliaepu1WM5aukpCR06dIF9+/fh7eXF3bMGwOz238CqcmAZ1WgfGcgPhRwLS/+Jhhj7F1w4KQPaMWQZxUM6ShXEl+2cQfSaFURYzoiJSVF7EF39OhRWFlZYefS7+BmkQY0+RYo1gioNwaICgAcisujTowx9o44cNIX1m7o8VEjONtawDIpGId3b1N3ixjLt6Cpd+/e+PPPP0XZgT+X/AhfFwVgW0iU5ED9sUBijLylCk3RGXLdX8bYu+PASV+Y2cHSGHg6yhT7PrHA7ytXcZI405mRps2bN4ugaeuqRWjtHACEPPyvNlN8OGBoCNDqUq6czxh7T1r30ys6Ohrbt2/Hq1evUKFCBbRs2VLl2966dQvbtm1DzZo10bx5c+gVhQLwrg1jE1O4G6Yi6fklPL5+BsWqNVV3yxh756Cpb9++2LRpkwia/lq5CO2KJgFHVsg1mij527MKkBAFeNXgvCbGmP6NOPn7+4tgacGCBXj06BH69++Prl27QqKtE94iNjZWXPenn37C7t27oZdsC8G4RBNx9LPKRpg9Z668VxdjWiY1NVUETRs2bICRkRH+WrUE7ctaAGeXyEGTUyk5ITw6UNQyg2MJdTeZMaYjtCpwGj9+PJydnXH69GksWbJEJILSCBJtqfA2w4cPx0cffYRixYpBbxmbAZV6iqNtShrh5LFDCHxwVd2tYizPQVO/fv3Sg6Y/V/2G9uWtgXNLgUg/uTJ4w2+AmCCANvN1rcD1mhhj+hc40YclTdHRhp00LE9KlSqF+vXrY8uWLbnedv369bh+/TpmzJjxgVqrwbxqQnLzhYFCgVHVDDBv3hwgNUXdrWJM5c8BGmlet26dCJo2r/oNHX3tgHO/AYHXAEMToOEE2t4XMLORk8HpBwNjjOlb4PT8+XPExcXBxyfzUmI6fffu3Rxv9/DhQ4wZM0YETyYmJio9VmJiIqKiojIddIaFIxQVu4ujfSsaY+v2fxDhn3P/MaZJdZroh9PatWthaGiITauWolMlR+D8MuDFJcDACGg8UR5loh8DHpV4HzrGmG4lh+/cuRPXrl3L9TojRoyAvb09YmJixGlbW9tMl9vZ2aVflt0Hbffu3TF58mSULVtW5XbRxqDTpk2DTqItJkq2gnR+OUzCnqK2eyoW/TofE2cv4+kMprEiIiLENio0PS8HTcvQuYozkBQLOJUE/M8DDcfLxS2jXwIeVQBbL3U3mzGmg9Q64kR7SSUkJOR6UCZ+W1paiv+zjv5ERkamX5YVDedTEnl4eDimT58uDrQa7/z58+J4TknlEyZMEPerPPj5+UGn2HhAUbU/dpl3xrrryZj/x5+IC3qi7lYxluNoc7169dKLW+7euhFdqrkCidGAtQfg2x1ovwTwrCYHTZQI7lxaXknKGGO6NOJEvyDpoApvb2+Ym5uLqbcWLVpk2tCTcp2yU65cOYwcOVIEYEoULFGehDIoU2Tz4WpqaioOOouWaRdvijbGtijqtQ9P/ALw+5IF+HzqPP6yYRrlypUrYlFHQEAAPDw8sHvTClRyTAbu7gbKfPzf+9XaFYjwA+wLy/WauMglY6yAKCRV1vJrCJp2o1+fJ0+eFImhNJpUunRpkfPQo0cPcR0qNfDixQsMHjw42/uoVKkSGjVqhPnz56v8uDTKRVOENPpkY2MDnRAfATw8jCVbjmDWL0tgYGmPuzeuw8TBU90tY0zYt2+fKCFCU/Hly5XDnpVz4GUcBtzbC9zZATiXAlrNoo8xIOI5YOclFj/w5r2MsbzKy/e81iSHk1mzZonAqXHjxiLhm/5v1aoVunXrln4dKk/w66+/qrWdWsHcDrD3xqAiz/Hwc2vUcYzGhlXL1N0qxoTff/8dbdu2FUFT00YNcGr9j/AyCBDBvgiaSOG6ctBEJQhsPYFC1TloYowVOK0KnIoUKYKbN2+KLRYoYZwCpB07dsCAEp5fow/bIUOG5HgfgwYNyjTVp9fsisDI0hFGBsDXdU0wa+HvSIsNVXermB6jFa2UY0h/pzSl/mn3TtgzbyRsE14CV9YDN/+Sr1ipN1C2AxDlL0/TUdBkaqXu5jPG9IBWTdWpi05O1RF66e/vg/RnXyhSE9FqXSwGf/MTOg0co+6WMT1DH0P0I2js2LFiCp5M/mIQpvZpAEVskFynKeaVXHKgxmCgZEsg0h+wchJbCYmil4wx9o50dqqO5TNKrHWvCEWReuLk13VNMXP+Yki0txdjHwiVJGnatCk6duwogiY3V2dsnDsB03rXhcLCCbiyTg6arFyBNj8DPq3lkSZLRzmniYMmxtgHxIGTvqMvI9/ukBQGaFLUCAh/gsM7X0+HMFaAgoKCxCKOypUri1IDpqYmmDi0Ox5s/QE9mvgCdoUBU0ug3higSH2g7S+AQzE5aDK3l4Mm+p8xxj4gDpz0nYEh4F0TCqqBA2BcHVPMnPMrkBSn7pYxHRUdHS022y5RogSWL18upum6t6qHe1tmYvrwrrAyNwWC78rvTUIlBqi4pZEpEP7sv6DJwkHdT4UxpofUWseJaQgqIlihi6i+3KyYEQbtvI5ju/5Eo0791N0yll1eWkoikBIPJCfI/9Pp+Ej5uMJAzgPKeKCaRgpDwND49cFEPojLlMcL/jfUjRs3sGTxYqxbv14ET6Rq2WKYP6Yn6tWuIT+3m1vkcgNSKmDlBriWk29MFcJjX8kjTrT/HO1DxxhjasCBEwOMTIDiTYCq/TD9r0eIStyNUeOn4HKzj2Bs48w9pG4UNMSFyXk+scGvA6dEIC1ZDjaIAQVERvJpKe31+Wn/naZ8tvR1IAp5NIeuT7ejAMrEUh7Jof+NzQEj89f/m2UfVIn7pfujg0K+/2yKySbERGHrnxuwZNnvOH3uUvr5PoXdMKH/x/i0W3sYUDtvUMC0B0hNlK/gXklUuBdiQ4DkOMCtIuBSVn6/MsaYmvCqOn1eVZcRTc09PIiwiBj4tB2B0PBIzJ86BqMn/8zVxD80CkgSo4C4UCAqAIgJApJi5NEkYwt5hIiCB/qfznuX+09LyXxITZKDMQqyCAVTNDVGgRMFUDQCJIKlVCAtQ0CmRO1QBk4Ghrj35AVWbNmPlVv3ISRcXmxgZGSIDk1qYFivDmhcryYUKQnAjT+Bu7sAOk6cSgGVe8uBEwVlUS/l50kb9tJoE2+jwhhT8/c8B0753KFaLeAG8PIKlu27jqWL5uFhtCnuXTkHtxIV1N0y/ZAYA0QHyAUdaYSJRlloRMjURq5R9C5B0rtKTZZHf8TIVsp/o0rZ/S9GnYAXgcHYvPcU1u8+gct3HqfflZe7Cwb3bIfPun0EdxenzMH61gFyUOhYEqj0CeBZVb5fCuQiXwBWzvJ5Vi4f7rkzxvROFAdO6utQrd+G5e4eSKfmIi3kAXyXxKB6/eZY9dfu/xJ1Wf6ikZv4cHmftfAnQEKkPKpEwRKN9ChHWGi6ii6nICM5PvOhRFN5NIZEBwJ+5+Ql+mZ2r/9/fSiA1zA8Mhp/7zuO9f8cxLGzV9M3zqbRpZb1a4iAqU2D6jAKeyC3K/g20Oqn/9ry6DBgbAV41fjvuVJ/0MGxuJzPZGqd7+1mjLF3/Z7nHCeWeRsW55JQmFjAUAHMbm6Kjzbsx+AD21GnVWfuqfxE011U2DHsiTzCRInelGNkX/T1CI/0XyBBU1lUADIHE1cdx7rriTAwUKBjKQPMrRfz5sNJQHSaGXZHl8WtZC9xXRuDBBQ2DkdEqjmCE40QEq9AfGIyEpKSkJD43yElJRVpUhpSU9OQliYhNS1V/E/n37j/GElJyemPU6+aL/q1a4BO9UrBXhEtB0t/LwIS5WRwgVbMKZO+izfNPAJFeVw0ukarPJ195GR2xhjTIBw4sczsigDluwCBN9GmJNCsWBJGjPkaFxu1gqGZJffW+0pJkqfjwh7L/1OeEBV5pGm4+3vFVClCHgB1P8cDRQnsOXYWT68cwezKEm4FpSE4TkJ0ooToJDpAHN97yx/PA+V8o/MKQ2y0NYGLpSL94GShgKGBAraGCfhr/ylsv0uBGdClrBHGdbVIb1piigS/KAnPI9PwPCUNS+8k46x/qrisuocBupUzhr25AvZm8sHOTAHbJqawt7DAQak+anzUB0UKucujSCcmZn7eNGpEwRCVEXAo+ua0YEygXB3FpTTg5CMH8YwxpoE4x0kFejNVR2iq5flZ4Pgs4PFR3AqW4LskGgtmfovh479Td+u0F42mUKJz6AN52o1GUowtAf8LcqDx6mamq/9+0xiDtsr7BtLon6kR4OjkgmLeHnBzdoCrkz1cHV//7+QAFyc78dJFx8QhOjYOMXHx/x2PjYVRcjRsDeIRnGyO2FRjpEkSypgGoaX9E9gbJsLWKBEGWRbFHTRqgQBzHxgZGqJ48h3UjD2Q8/OrPw4o1lA+7nceOLNAnh50rywHSy5l3pwqTEsF4kLkxHAbT/k6VJCVE8AZYx8Y5zipsUN1Am30e+cf4MBEsRT+s3/i8fcjE9y/eQ3OXsXV3TrtQjlLlL9EI0wJ4XI+j4UjokMDYbF3JAyl5PSptCNPUrDxZjJOPEvFw7A0GBsboUH1imjTqBbaNK6FUsW8oSiooIKmBykhnaYPY4LlsgcUCFEgQ6j9j48BJlZyyQL63/T1cSpdQEnclJul0mOlAgkRct9YOsujTLbecnkExhhTAw6c1NihOsP/MnDmV+DWNgTHG6DovAj07NIBy9f/zSMCb0NDP1RKgKpchz+VE7pTEvHo9hXM2PMcpy/dwP0nfjg30BI2psDqa8lYdz0Z/lGSWIHWqmFNESw1rVMV1lYqBiPa0CfUDxQwUX4XTcVRPpdjMTkJnjHG1IiTw9n7cyoBlG4LPDsNczNDuFtHYcWmHRg0eD9qNGrFPZwTWg0WfF8OmBKjkfrqLsKu74Fz8gt4pkjYtisaEa9LFg06bIliJX1QpXopLO3vg8rlSsLd2VFeii/KAcQDUZH/FYXM1n/lAMQUlyiEafL+tZ7yC03D0WpNWv1Ho1O0+s+2EGDpwoUsGWNaicfGWfZoSwuq0lxzOKzcK6LOhZ/w8O/9GPH5Fzh3+QYMjHi1Uya0Ko6ms4LuAIHXkfTsIuB3FiaKFFDt9dQ0CUefpuGz9g3QrFU7VC1fCs6OrxOgKUiigpc0IhMRkyHwMZUTx+m1oGkwYzN56xRRGZwKUqbJ017KwpR0PzT9RfcltmGJlZPRlYFVxi1WlIEVnfe+qA0i2FMGfEn/VTWnQM7SCfCsItdi4tICjDEtx8nhKtDLqTpCowSPjoig4FW8EXya9kJUTBzmzfgWX/yPE8XTg4ZIfyDotsgNenVlL1z996R34f3QVGy+ZwjTMi3Rp1dPuQAkBRTUtyLASZCDFwqOrNzlXKH0bU8oUHqHnCa6fwpe6DHo/ulAgR2VBBABWtx/gY4offAaJW+LLVgMM+91R6NaFAhRUJT+f2qGgEy5D97rYM/EWg6QKNijBHELxw+yFx5jjL0rznHKZ3obOJHQR2K6jqZWzm2ahSFLTuBmCLD/n61o2qYD9BolUQffRdzLe/jr6BUs334C9+7cwo1hlth1PwXHQ53QqG1P9OrQAuYmRvKIEgUuFKxQUEF1m2jaytxBzvn5UDWLUlP+C6iUQZUItOLk4/S/2H4l9b+RI7ERsDFgYvFfgrjYksUUMHz9PwV6vI8cY0wLceCkxg7VOTS6QKupziwU/9+KtkH5uf6wt7XGuTNnULLM60KG+oRGcoLv4dqJ3fA7sR6pEf7osFEu8GhoaIiOTWtiSO8uaFrLF4qkaDlgosEZWoVm7SavVLNwkKuDa+rSe0rgTp9yS5ODIxE88cgRY0z3cHI4yz80CuJcGijeBHh6GuWsozC+dVHM2vsEH3/cFmcvXIGtnZ4UK5QkpEX4YfOKhbhzeD0Gl4lCRRcq2qhAt+quqNSwPfp1bgl3K4UcLNHSfpqCcy0vL7unYElbVpBRgGRgRslu6m4JY4xpFE4OZ29HxQmpkGGJJsD9/fihbgKO3nfE+YdP0bNrB+zcd1iMtOi0xBgc/2c9Fs75AV+UDUHPmvSnY4DgRBMEerfHxg29YUCjS7SqztBe3piWkqJpGo7rEzHGmM7gcXemwrvEAHAuBZT+GHAoDoPkWBwZ7A4bCxPsPXQc33w1Wnd7MS0N984fQee2zXBgwedY1ywcdb2NkCgZIaZ0Nzh/tgkVmnSGQeRzeVrLo7I8OudaVl5FxkETY4zpFA6cmGooCKA9xKr0EcnBlnHPceHbWuKin39ZhNUrlulcTwb7PcbI/t1Rrk5z7Dt+DgMqm8DUSIFE54ow7bocVlW6yfvN0Wo1CpQoYHKvICdOM8YY00kcODHVUBIzFcWkrTFqDhVn+aQ9xI8ju4rjg4ePxJnTp3SiNxNiojBr8jiUq1ABi9dsQWpqGpo2qAOj+l8A9b+Cacup8qoz2meNql8XawwUqsYb0zLGmB7gHCemOkpudiop9q9DzWEiWBhn4Yxz915i28HT6NixIy5cugwvLy+t7FUpNQXbN6zA2InT4KEIwsneZtju54ka3cehce0q8pWowGSUH2D9elNaa3deacYYY3qER5xY3lCuk53X65ViLjAwMMCaud/C16cIXgWHoFmTRrh//7529aok4cbZo2hWvyYGDhmGbyqF42R/S5RyMsTXDW3QuGZFuSxDxHO57pFHVaBofcDWk4MmxhjTMxw4sbyhQocelQAzayA2WJxlFX4LJ8dVhpebA+4/fIyaNarj4IEDWtGzIc8fYPinXVCpblMUiruBuyOtMKiKiXxhyRZQfDRHLlpJo0wUKBVrBLiVl/uBMcaY3uHAieUdVbx2ryQXSAx7ChyZDhv/w7ixoDdq+ZZERGQUWrdpgwW/zIdEVac1UFJ0GH79fjxKlq+Mg3u243Afc6zuYA5nCwVg5w20miXncsWFypW+vWoBhevKJQYYY4zpLQ6c2Luh4II2AVakARV7irNsb6/HsSVj8Wn7pkhNTcXnX3yJoYMGICmJNprVDCkJsVi1YCZKlSmD0ZN/QkR0LCr5eKNhEdr41hSo0g9o+4tc3TvqJWBfRB5loinKD7UlCmOMMY3Fm/yqQK+3XMlNShLw/AwQ/hS4shbwvyA2d5XqjcWc3Xfx9U/LxIhTw7q1sGX7Tjg5qW+0Ji05CZtX/4apM37C/acvUM7ZAKGww9TR/TGwe1sYPjkGuJaTN6SNeiHvx+ZWAXAoKm9iyxhjTGfxXnVq7FC9Ex8BPD0FxIbI+9mFPqB4XIxC7QrxRq8vv0d0bDyKenvin392oXzFSh+0eVJaGrZvXInJ383AzftPUNbZAAvbWqGhlwKJLX+GuXupDM8lHEiIAOyKyHlMtIqQMcaYzovKw/c8T9Wx92NuB7hXBIzMgCaTAJ9WFK4A1zagbWlTnNn6G4p5ueHJ8xeoUr06PuneBf+ePl3guU/xsTHY9McSVK9UFp16D0SA31Msa2+FG8Os0dibBpEMYR7zXL4y5TDRijlaOVeoBlC4DgdNjDHGssVTdSrgEae3oCAo8Cbw8jJg6wU8PQGEPEgvlBkSFoHeX36H/Scvpt+kkm95DB8xCr0++QSWlvlTaZuCsdPHj4gq5n9u34WomDgYGQBj6lpiaiMzmBsky1f0rgNU7Q/YuMt1meJD5cKeNDXHyd+MMaZ3ovIw4sSBUz53qN6i0Rq/c0DoY7nOk0GG2qoJUYD/eVyKK4RFa/7Cxt3HkJAoBzG2tjbo368/Bg8ZAh8fn3faLPjp06dYs2IZ1qxbh0dP/dLP9/Zwxqn+FvAyCpfPoCrf1QcB7r5ye2MC5XY6l5GTv41elyFgjDGmV6I4cFJfh+o12rPN/yIQ6SevSjO2oCQj4PA04MUloEh9oHIfhKZaYuWmHViy4R88fhGUfnMqpuns7Aw3N7f0g6urK1xcXJCcnIyIiAhxoNchIiwUEeFhCAsLw/1HT9Lvw8bCDJ1bN8SnndugQY2KMLi7E7jxp3hclGgubx1DJQaSYuTRMVoZaO2qpg5jjDGmCThwUmOH6j3awy3wBhB8FzCxlms+3fobuLxaDqIUBnI9pPJdkGZfFPuPnMSiNVux78x1sSfcu1AoFGhZyxdTPi6CGoa3YFCtvxykKfOXqE208S4FdrS/HLWJAia7woAh7zrEGGP6LopHnNTXoYyClTQg7BEQcA1ITZFziSjn6dpG4MV/eU7wqCLXgHIpg5T4WIQEv0JgUDACg0MQGBSGV2GRCAyNxKvQSJiaGMPO2gJ2Npaws7WDnZ0dbO3s4WBtBl+Dh7DzOyCPJBE3X6DlD/89DhXqjA4ADE0AJx/AsQRgasUvFWOMsTx/z/PPbZb/DAzkzYBNrYGXV4DwZ4BDMaDZVCDsMXBzK/D0pJxM7llVBE5G5pZwc3GCm6uzXEOJptSUxAo8SR6tUqK6US+OAY+eyqNKxNwBKN8ZKNlSPp2WKo8wJSe8LthZBrBy5lecMcbYO+PAiRUcazd5Wo5GnsIfi02BRQDVYJycc0T5R8ogh1zfBNz5R07YNrOTp9So3AEFYLT6renU/wIqCsZCH75+HHc5YCreVK7uTVOCVFeK9pizcpGDMwqcuJAlY4yx98SBEytYZjaAV005+Am6IxeZtHCSgypa4ZZRcpz8P40g0UgRHTKiTYUpECJl2slbodDIFgVkFFDRyFRcmFzEkiqAUz0mCph4Q17GGGP5hMsRqIBznPIBBTW091v4EyDyhZx3RCNKpjaZp+WoTICygjf9T5XJabTJ0hnwqiEneWd336IeU5g8QkV5TJT4bWKRHy1njDGm46I4x4lpHAqObD0BGw85iZum2iKeyfvc0WgUBVE0lUZTbTSqpBxZygkFWDRClRwPpMTLeVGUbE57y3HiN2OMsQLCU3XswwdQVJ2bDjTNRhvqhj4CImn7E4Wc3yQOhpmP06gSBUkULFEOE51vbC7fD9WMoqDMzJZfTcYYYwWKAyem3vwnOlBFbyoXkBT73wgSrYSj6Tw6JNOqOQVgailPwVlQwriNPMpkbMavIGOMsQ+GAyemfhT80BRbdvWg0pLlaTlC+U0Z86EYY4yxD4wDJ6bZ9aAMTHlVHGOMMY2RoaIgY4wxxhjLDQdOjDHGGGO6HDilpaUhJiYmz7dLSkpCSsrr7TkYY4wxxnQ5cJIkCRMnThQbvDo4OKBo0aLYvXv3W2936dIlNGjQANbW1nBzc8PIkSMRF/e6SjVjjDHGmC4GTr/88gsWLVqEAwcOiMBn2LBh6NSpE+7fv5/jbW7fvo2GDRuiVq1aCA8PR2BgICpUqIBbt2590LYzxhhjTPtp1ZYrxYsXR8eOHfHzzz+nn0ejThQ8zZkzJ9vb0GX+/v44f/78Oz8ub7nCGGOM6a68fM9rzYhTcHAwHj9+jHr16mU6n6bgzp07l+1tKJ9p37596Nq1qzgdGxv7QdrKGGOMMd2k1jpOlOCdkJCQ63Uol8nAwEAETsTZ2TnT5XT67Nmz2d6WbhMfHy8CptKlS+P58+cwNjZG37598dNPP8HMLPuq04mJieKQMRJljDHGGFNr4DRhwgRs3Lgx1+tcuHBBTMcpXleMTk1NfWNUiQKr7Chvs2DBAhw8eBBVqlTBjRs30KxZM5iammL27NnZ3m7mzJmYNm3aOz4rxhhjjOkqtU7VUUATEhKS64GCJuLh4SH+f/XqVab7oNPu7u7Z3r+TkxNMTEzwySefiKCJUGI4jTjt3Lkz14CO5jmVBz8/v3x81owxxhjTVlqz5Qolbfn6+uLQoUPpOUs0+nTkyBEMHTo00/QfjUJRyQIjIyPUrVv3jdIDdNrc3DzHx6LRKDooKfPnecqOMcYY0z3K73eV1stJWmTz5s2SsbGxtGrVKunevXvSwIEDJQcHBykgICD9Op999plUrly59NOHDx+WLC0tpQ0bNkjPnj2T/vzzT3F6wYIFKj+un58f9SQfuA/4PcDvAX4P8HuA3wPQ3T6g7/u30apyBGTNmjWYP3++mKKjabdZs2ahYsWK6Zd/+eWXuHjxIk6ePJl+3q5du0Te0rNnz+Dt7Y2BAwdiwIABeapU/vLlS1FAU5k3lZ9RrpeXl5gOfNsSSMb9xu+3D4//Rrnf+P2m+3+nkiQhOjpapAXllDetpHWBk67hGlHcb/x+02z8N8r9xu83zReVhzpM70tr6jgxxhhjjKkbB06MMcYYYyriwEnNaPXelClTMq3iY9xv/H7THPw3yv3G7zfNZ/oBv0s5x4kxxhhjTEU84sQYY4wxpiIOnBhjjDHGVMSBE2OMMcaYijhwUpMnT55gyJAhaNiwIXr37i2KdrLMkpKSxCbQrVq1QvPmzXMsTrpkyZL068ybN09suaPPzp8/j2HDhqFp06bivXXgwIE3rkN9RH1FfUZ9R31IfanPHj58iDFjxohNwLt06YIVK1a88V5KTEzEjz/+KK7z0UcfYdWqVWprr6aJj49Hu3btUKtWLVFLJ+s2V9999x2aNGkirrNhwwboM39/f9FPWQ9nz559ozbRxIkT0bhxY3To0AHbtm1TW5s1yYEDB9CrVy/xGTdt2jTx3ssoPDwcX3/9tei3Tp06Yc+ePfnbgDzve8LeG20R4+rqKnXt2lXas2ePNGrUKMnMzEy6fPky924GlStXlrp37y4NGDBAsrW1zbZvxowZIzk5OUlr1qyRNm7cKHl4eIhtd/TVpk2bpFq1aklLliyRDh06JP3444+SiYmJOJ0R9WmhQoVEn61evVpydHSUvvrqK0lfPXjwQKpSpYr066+/SkeOHJGWL18u/kYHDRqU6XqdOnWSSpQoIf3111/iOtbW1tL333+vtnZrksGDB0vly5cX21YEBwdnuqxly5ZS2bJlpa1bt0qLFy+WzM3NpXnz5kn6/H6jfvr777+lM2fOpB8iIiLSr5OamirVrl1bqlq1qrR9+3bRX7Tl2MqVKyV9NmvWLLFt2s8//ywdO3ZMmjlzpjRixIj0y5OSkqRKlSpJderUkf755x/xGWhkZCS2W8svHDipwddffy0VLlxYSklJST+vYcOGUocOHdTRHI2l/BChfQWzC5woADU0NBRf/ko7d+6UFAqF+GDSR9HR0W+c9+WXX0olS5ZMP33//n3xob13797089auXSs+XAIDAyV9lJCQIL6oMpo7d26m992FCxdEv9EXnBK9NykIiIqKkvQZBZIVKlSQtmzZ8kbgRIEonXfjxo308+jLjvqW+l2fA6cnT57keB0KquizjPZYVRo/frz4cZj1vaov7ty5Iz7z161bl+n8+Pj49OP0I5o+y4KCgtLPo8CKfvDkF56qU4PDhw+jdevWMDQ0TD+vffv2OHTokGo7M+sJKp+fm2PHjiE1NRVt27ZNP69ly5YwMTERfayPrKyssj2Ppj2VqG/MzMzEdFPG9x9NS1Gf6iOq/ZJxfyrqi3///TfTPpj09+nk5CSmVDL2G00T0HX11dOnTzFq1Cgx/ZZdDR16vxUuXBjly5fP1G80nafvKQqDBg0S0000tX7r1q03+o3ef7S/asZ+o31T79y5A320ceNG8b3Qo0ePTOfT51nGfqO/UWdn50z9RlPxtF9tfuDASQ3oxaONBDOi0zExMWJulqnej7QnUcZgwdjYWPzB5NcfiLYLCQnBsmXLRH6EEvWNi4sLjIyM0s+jDaypH/W932iT8Bo1asDd3R0RERGZckqy+7ul69HG3/rabxRgUq7JN998kykwyoj6xtPTM9N5yn7U134jjRo1Qv/+/TF27Fjxw6Zy5co4evToW78nlJfpo9u3b6NSpUrYvn072rRpI/IMf/jhB5FD9yH77b9PTvbBJCcnv/HLzNzcPP0y9u79qOxL7kc5IZcCJldXV8yYMYP7TQW0YKNz5864efOmSDqlhOb58+fn+H6j4JMO+vp+mzRpkvjx8vnnn+d4Hf68exONwB05ckQE3YRmIGgEbty4cemjcNRvFhYWmW6n798TCQkJuHTpEhYuXCj6ikZ7J0+ejP3794ugk0aNP8T7jQMnNXBwcEBYWFim80JDQ8WLbmdnp44maW0/0ggdTW8qP4CUfeno6Ah9Rh8oNDxN/UMfKJaWlrm+/6gP6Tx977fSpUuL/+vVqydGkyjwpJV2NF2SXb9FR0eLD2N97bfVq1eLqZPatWuL0zRKR2il5oABAzB8+HDRb3fv3s10O/obJfrabzQynhVNnY8ePTr984z6jablMtL3fnNwcBB/c3///Tfs7e3FeYUKFRJTc1evXkWVKlVy/H7Nz37jqTo1oBf3woULmc47d+4cypYty3vW5bEfaaqA/mCUHj16JP5oaNhbX9GvMgqaXrx4IX7V0rRc1n6jDxIqiaFEv+IoX0yf+y0rCpwyfuhSv1E+D01/Zvy7Jfrab7t37xYlGWhUjg6Us0NohJPKDij7jQInSkXI2m807cJkgYGBYmRE+SOQ+u3atWuZRkmo32iEM6dpUV1XrVo10UfKoCnj36kyzYX6jT7PMuYLU7/R7UqVKpU/Dcm3NHOmMipBQCsDjh49Kk7fvXtXrDDR5+W5uclpVV1aWppY/kxLxJWrTPr06SMVLVpULEnVR7RKSbn0O6cVcomJiVKRIkWkfv36idPUd+3btxeroqhP9dGOHTukixcvZlqdSO8rKtmgfC/ReS4uLtLo0aPFaTq/SZMmUt26ddXWbk1Dq1qzrqoLCwuT7OzspIkTJ6avgKJl9vQ+1Vfr16+Xbt++nX76ypUrkoODgzRkyJD08/z9/cWKzdmzZ4vTtHKT/kapRIu+CgoKkmxsbKQ//vgj00pDe3t7KTw8XJx++PChKNugLMFC7z8fH598LVPDgZOaTJ8+XTI1NRUvKNXZoRc1Y3kCJpdtqFmzpviSp0CTjtOBltNnXJ5aqlQpydnZWXJzcxPX1ed6WL/88ov44qL3lbK/lIeM769Lly6JvqI+o76jPqS+1Fc3b96UGjVqJLm7u4svJ6oTQyVCrl+/nul6J0+eFMvBPT09xYe1r69vrkvK9U12gRM5ePCgCDq9vb3Fj6Dq1atLL168kPTV2bNnRX0mLy8v8bdK3wW0ZD4uLi7T9ah+E73P6Mcg1QyrX7++FBISIumzgwcPis8tKrFC/UelfahmXdZ6dhRgFStWTLKwsJCaN28uRUZG5lsbFPRP/oxdsbyifACaLqGMf0rgZZndv3//jblq4uvrmylpkt7CNBVAla/LlCmTaVm5vgkICMhx5UjGZfSE+ouWNVN/UW5PxjwxfUXTcJRXQnkTlCuRHZoepn6jBFQfH58P3kZN/0yjv0WaUsm4apPQlBP1G/3tlihRQm1t1CT0XqOk8KJFi2ZaUp+1Wj31KSXh0/UY0v8G6b1UpEiRTKV9MuZ53rt3T0zrUTJ+fuLAiTHGGGNMRfr705wxxhhjLI84cGKMMcYYUxEHTowxxhhjKuLAiTHGGGNMRRw4McYYY4ypiAMnxhhjjDEVceDEGGOMMaYiDpwYY4wxxlTEgRNjjDHGmIo4cGKMMQAPHz7E5s2bkZSUlGnbBjqPtkZijDHCW64wxhiAqKgoVKxYEZ06dcKcOXNEnwwfPhyHDh3ClStXYGlpyf3EGOPAiTHGlE6fPo1GjRph9+7dYnPVzp07i/OqV6/OncQYEzJvX80YY3qsbt26+N///oe+ffuKHdinTp3KQRNjLBOeqmOMsQzi4uLg6uoKa2tr+Pn5wdDQkPuHMZaOk8MZYyyDKVOmwN7eHtHR0Vi+fDn3DWMsEx5xYoyx144cOYKWLVvi6NGjePDgAUaOHInLly+jVKlS3EeMMYEDJ8YYAxAWFgZfX1/069cP06dPF33SrVs3PH78GGfOnIGxsTH3E2OMp+oYY4zs2bMHbdu2FVN1SkuXLkWZMmVw/Phx7iTGmMAjTowxxhhjKuLkcMYYY4wxFXHgxBhjjDGmIg6cGGOMMcZUxIETY4wxxpiKOHBijDHGGFMRB06MMcYYYyriwIkxxhhjTEUcODHGGGOMqYgDJ8YYY4wxFXHgxBhjjDGmIg6cGGOMMcZUxIETY4wxxhhU83/wAtZEmKwTewAAAABJRU5ErkJggg==", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "fig, ax = plt.subplots(figsize=(6, 4))\n", + "ax.plot(x_grid, true_u0, color=\"black\", label=\"true $u_0$\")\n", + "ax.plot(x_grid, posterior_u0_mean, color=\"tab:orange\", linestyle=\"--\", label=\"posterior mean\")\n", + "ax.fill_between(\n", + " x_grid, posterior_u0_mean - 2 * posterior_u0_std, posterior_u0_mean + 2 * posterior_u0_std,\n", + " color=\"tab:orange\", alpha=0.25, label=r\"posterior $\\pm 2$ std\",\n", + ")\n", + "ax.set_xlabel(\"x\")\n", + "ax.set_ylabel(\"$u_0(x)$\")\n", + "ax.set_title(\"Posterior over the initial condition, recovered by NUTS\")\n", + "ax.legend()\n", + "plt.tight_layout()\n", + "plt.show()\n" + ] + }, + { + "cell_type": "markdown", + "id": "0b16516c", + "metadata": {}, + "source": [ + "## What this notebook set up\n", + "\n", + "- A `dynestyx.DynamicalModel` whose observation operator is a real\n", + " `jax.experimental.sparse.BCOO` matrix -- built with `GridInterpolator.as_sparse()`,\n", + " not hand-rolled -- used for both forward simulation (`dsx.simulate`) and filtering\n", + " (`EnKFConfig`, via `dsx.condition`/`dsx.sample`) identically to how a dense `H` would\n", + " be used.\n", + "- The initial condition treated as a 64-dimensional Gaussian-process-distributed latent,\n", + " recovered with NUTS running through the EnKF-marginalized likelihood -- initialized at a\n", + " closed-form GP-regression warm start rather than a hand-picked value.\n", + "\n", + "See `pde_inference_setup_exponax_ks.ipynb` for the full derivations, the dense-vs-sparse\n", + "comparison, and the additional inference problem ($\\lambda$, the equation's nonlinearity\n", + "strength) this notebook leaves out.\n" + ] + }, + { + "cell_type": "markdown", + "id": "fb675450", + "metadata": {}, + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "dynestyx", + "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.13.14" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/dynestyx/observation.py b/dynestyx/observation.py new file mode 100644 index 00000000..b610db84 --- /dev/null +++ b/dynestyx/observation.py @@ -0,0 +1,310 @@ +"""Grid-interpolation utilities for building observation operators. + +`GridInterpolator` observes a state that lives on a regular N-D grid at arbitrary query +points -- not necessarily grid-aligned -- via nearest-neighbor or (multilinear) +piecewise-linear interpolation. Both interpolation kernels are linear in the grid values, +so the whole operation reduces to a fixed matrix applied to the flattened state; this module +computes that matrix (and, for the ``"constant"`` boundary mode, an additive bias) directly +from the interpolation weights, rather than recovering it after the fact via autodiff. + +The intended pairing is `dynestyx.models.observations.LinearGaussianObservation`: + +```python +interp = GridInterpolator((x_grid,), query_points, method="linear", boundary="periodic") +observation_model = LinearGaussianObservation(H=interp.as_matrix(), R=obs_noise_cov) +``` + +or, more directly, `interpolation_observation_model(...)` below builds that +`LinearGaussianObservation` in one call. +""" + +from __future__ import annotations + +import itertools +from typing import Literal + +import equinox as eqx +import jax.numpy as jnp +from jax import Array +from jax.experimental import sparse as jax_sparse +from jaxtyping import Bool, Float, Int, Real + +from dynestyx.models.observations import LinearGaussianObservation + +BoundaryMode = Literal["periodic", "constant", "edge"] +InterpolationMethod = Literal["nearest", "linear"] + +_VALID_METHODS = ("nearest", "linear") +_VALID_BOUNDARIES = ("periodic", "constant", "edge") + + +def _row_major_strides(shape: tuple[int, ...]) -> tuple[int, ...]: + """C-order (row-major) strides for flattening a `shape`-shaped array.""" + strides = [] + acc = 1 + for size in reversed(shape): + strides.append(acc) + acc *= size + return tuple(reversed(strides)) + + +def _axis_candidates( + q: Float[Array, " n_query"], + x_axis: Real[Array, " _"], + method: InterpolationMethod, + boundary: BoundaryMode, +) -> tuple[ + Int[Array, "n_query k"], Float[Array, "n_query k"], Bool[Array, "n_query k"] +]: + """Per-axis candidate grid indices, weights, and validity for one query coordinate. + + `k` is `1` for `method="nearest"` and `2` for `method="linear"` (the two grid points + bracketing the query coordinate), kept uniform across axes -- including degenerate + single-point axes, where the second `"linear"` candidate is a zero-weight duplicate -- + so that axes can be combined by a plain outer product later. + """ + n = x_axis.shape[0] + k = 1 if method == "nearest" else 2 + + if n == 1: + idx0 = jnp.zeros_like(q, dtype=jnp.int32) + idxs = jnp.stack([idx0] * k, axis=-1) + if method == "nearest": + ws = jnp.ones_like(q)[:, None] + else: + ws = jnp.stack([jnp.ones_like(q), jnp.zeros_like(q)], axis=-1) + valid = jnp.ones((q.shape[0], k), dtype=bool) + return idxs, ws, valid + + x0 = x_axis[0] + dx = x_axis[1] - x_axis[0] + extent = n * dx # assumes x_axis is uniform and, for boundary="periodic", does not + # repeat the periodic image of the first point (e.g. jnp.linspace(..., endpoint=False)). + + if boundary == "periodic": + q_adj = (q - x0) % extent + x0 + elif boundary == "edge": + q_adj = jnp.clip(q, x0, x_axis[-1]) + else: # "constant" + q_adj = q + + idx_float = (q_adj - x0) / dx + + if method == "nearest": + idxs = jnp.round(idx_float).astype(jnp.int32)[:, None] + ws = jnp.ones_like(idx_float)[:, None] + else: + idx0 = jnp.floor(idx_float).astype(jnp.int32) + frac = idx_float - idx0 + idxs = jnp.stack([idx0, idx0 + 1], axis=-1) + ws = jnp.stack([1.0 - frac, frac], axis=-1) + + if boundary == "periodic": + idxs_final = idxs % n + valid = jnp.ones_like(idxs, dtype=bool) + elif boundary == "edge": + idxs_final = jnp.clip(idxs, 0, n - 1) + valid = jnp.ones_like(idxs, dtype=bool) + else: # "constant" + valid = (idxs >= 0) & (idxs < n) + idxs_final = jnp.clip(idxs, 0, n - 1) # safe for gather; masked out via `valid` + + return idxs_final, ws, valid + + +class GridInterpolator(eqx.Module): + r"""Linear observation operator for a regular N-D grid, observed at arbitrary points. + + Precomputes, once at construction, the sparse "corner" structure of either + nearest-neighbor or (multilinear) piecewise-linear interpolation: for each query point, + the up-to-$2^d$ bracketing grid points and their weights. Both kernels are linear in the + grid values, so `as_matrix()`/`as_sparse()` expose that structure directly as an + observation matrix $H$, and `__call__` evaluates the same structure directly rather than + via a matrix-vector product. + + Attributes: + grid_shape: Shape of the regular grid, e.g. `(64,)` for 1D or `(64, 64)` for 2D. + The flattened state this operates on has length `prod(grid_shape)`, in C + (row-major) order. + method: `"nearest"` or `"linear"` (multilinear in N-D). + boundary: How out-of-range query points are handled -- a single, global setting + for now (not per-axis): + + - `"periodic"`: wrap around (query points are taken modulo the domain extent + per axis). No bias term; every row of the resulting $H$ sums to 1. + - `"edge"`: clamp query points to the grid's extent (flat/Neumann-like + extrapolation). Also bias-free. + - `"constant"`: query points (or interpolation corners) outside the grid + contribute `fill_value` instead of a grid value. This introduces an additive + **bias** (see `bias()`) alongside $H$ -- pair with + `LinearGaussianObservation(H=..., bias=...)`, not `H` alone. + fill_value: Constant used for out-of-range contributions when `boundary="constant"`. + Unused otherwise. + + Note: + Grid axes are assumed uniformly spaced. For `boundary="periodic"`, an axis is + assumed *not* to repeat its periodic image as an explicit grid point (matching + `jnp.linspace(..., endpoint=False)`), consistent with how periodic PDE grids are + built elsewhere in this codebase. + """ + + grid_shape: tuple[int, ...] = eqx.field(static=True) + n_query: int = eqx.field(static=True) + method: InterpolationMethod = eqx.field(static=True) + boundary: BoundaryMode = eqx.field(static=True) + fill_value: float = eqx.field(static=True) + corner_indices: Int[Array, "n_query n_corners"] + corner_weights: Float[Array, "n_query n_corners"] + corner_valid: Bool[Array, "n_query n_corners"] + + def __init__( + self, + grid_axes: tuple[Real[Array, " _"], ...], + query_points: Real[Array, "n_query ndim"], + *, + method: InterpolationMethod = "linear", + boundary: BoundaryMode = "periodic", + fill_value: float = 0.0, + ): + """ + Args: + grid_axes: Per-axis 1D grid coordinates, e.g. `(x_grid,)` for a 1D grid or + `(x_grid, y_grid)` for a 2D grid. Any number of axes is supported. + query_points: Observation locations, shape `(n_query, len(grid_axes))` -- + including for 1D grids, where this is `(n_query, 1)`, not `(n_query,)`. + method: `"nearest"` or `"linear"`. + boundary: `"periodic"`, `"constant"`, or `"edge"`. See class docstring. + fill_value: Constant used for out-of-range contributions when + `boundary="constant"`. + """ + if method not in _VALID_METHODS: + raise ValueError(f"method must be one of {_VALID_METHODS}, got {method!r}.") + if boundary not in _VALID_BOUNDARIES: + raise ValueError( + f"boundary must be one of {_VALID_BOUNDARIES}, got {boundary!r}." + ) + query_points = jnp.asarray(query_points) + if query_points.ndim != 2 or query_points.shape[-1] != len(grid_axes): + raise ValueError( + "query_points must have shape (n_query, len(grid_axes)); got " + f"{query_points.shape} for {len(grid_axes)} grid axes." + ) + + self.grid_shape = tuple(int(axis.shape[0]) for axis in grid_axes) + self.n_query = int(query_points.shape[0]) + self.method = method + self.boundary = boundary + self.fill_value = float(fill_value) + + ndim = len(grid_axes) + strides = _row_major_strides(self.grid_shape) + + axis_results = [ + _axis_candidates(query_points[:, i], grid_axes[i], method, boundary) + for i in range(ndim) + ] + k = axis_results[0][0].shape[-1] # uniform across axes by construction + + flat_indices = [] + weights = [] + valid = [] + for combo in itertools.product(range(k), repeat=ndim): + corner_idx = jnp.zeros((self.n_query,), dtype=jnp.int32) + corner_w = jnp.ones((self.n_query,)) + corner_valid = jnp.ones((self.n_query,), dtype=bool) + for axis, sel in enumerate(combo): + idxs_i, ws_i, valid_i = axis_results[axis] + corner_idx = corner_idx + idxs_i[:, sel] * strides[axis] + corner_w = corner_w * ws_i[:, sel] + corner_valid = corner_valid & valid_i[:, sel] + flat_indices.append(corner_idx) + weights.append(corner_w) + valid.append(corner_valid) + + self.corner_indices = jnp.stack(flat_indices, axis=-1) + self.corner_weights = jnp.stack(weights, axis=-1) + self.corner_valid = jnp.stack(valid, axis=-1) + + @property + def state_dim(self) -> int: + state_dim = 1 + for size in self.grid_shape: + state_dim *= size + return state_dim + + def __call__(self, values: Real[Array, " state_dim"]) -> Float[Array, " n_query"]: + """Interpolate a flattened grid state at the query points.""" + gathered = values.reshape(-1)[self.corner_indices] + grid_contribution = ( + jnp.where(self.corner_valid, gathered, 0.0) * self.corner_weights + ) + fill_contribution = ( + jnp.where(self.corner_valid, 0.0, self.fill_value) * self.corner_weights + ) + return jnp.sum(grid_contribution + fill_contribution, axis=-1) + + def as_matrix(self) -> Float[Array, "n_query state_dim"]: + """Dense observation matrix $H$ such that `H @ values == self(values) - self.bias()`.""" + n_corners = self.corner_indices.shape[-1] + query_idx = jnp.repeat(jnp.arange(self.n_query), n_corners) + col_idx = self.corner_indices.reshape(-1) + valid_weights = jnp.where(self.corner_valid, self.corner_weights, 0.0).reshape( + -1 + ) + H = jnp.zeros((self.n_query, self.state_dim)) + return H.at[query_idx, col_idx].add(valid_weights) + + def as_sparse(self) -> jax_sparse.BCOO: + """Sparse (`jax.experimental.sparse.BCOO`) observation matrix $H$. + + `LinearGaussianObservation.__call__` computes `H @ x`, which accepts a `BCOO` + operand directly -- pass `as_sparse()`'s result straight to + `LinearGaussianObservation(H=..., R=...)` when the query points are sparse enough + (few corners per row relative to `state_dim`) for that to be worthwhile; use + `as_matrix()` otherwise. + """ + n_corners = self.corner_indices.shape[-1] + query_idx = jnp.repeat(jnp.arange(self.n_query), n_corners) + col_idx = self.corner_indices.reshape(-1) + valid_weights = jnp.where(self.corner_valid, self.corner_weights, 0.0).reshape( + -1 + ) + indices = jnp.stack([query_idx, col_idx], axis=-1) + return jax_sparse.BCOO( + (valid_weights, indices), shape=(self.n_query, self.state_dim) + ) + + def bias(self) -> Float[Array, " n_query"]: + """Additive constant from `boundary="constant"` fill-value contributions. + + Zero for `boundary in ("periodic", "edge")`, where every corner is always valid. + """ + invalid_weights = jnp.where(self.corner_valid, 0.0, self.corner_weights) + return jnp.sum(invalid_weights, axis=-1) * self.fill_value + + +def interpolation_observation_model( + grid_axes: tuple[Real[Array, " _"], ...], + query_points: Real[Array, "n_query ndim"], + R: Float[Array, "n_query n_query"], + *, + method: InterpolationMethod = "linear", + boundary: BoundaryMode = "periodic", + fill_value: float = 0.0, +) -> LinearGaussianObservation: + """Build a `LinearGaussianObservation` that observes a regular-grid state at arbitrary, + non-grid-aligned points via nearest-neighbor or piecewise-linear interpolation. + + A thin convenience wrapper around `GridInterpolator`: constructs the interpolator, then + returns `LinearGaussianObservation(H=interp.as_matrix(), R=R, bias=interp.bias())` + (the `bias` argument is only nonzero for `boundary="constant"`). + """ + interp = GridInterpolator( + grid_axes, query_points, method=method, boundary=boundary, fill_value=fill_value + ) + bias = interp.bias() + return LinearGaussianObservation( + H=interp.as_matrix(), + R=R, + bias=bias if boundary == "constant" else None, + ) diff --git a/pyproject.toml b/pyproject.toml index 80edc6a1..abd38492 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -75,7 +75,11 @@ dev = [ "nbconvert>=7", "seaborn>=0.13.2", "imageio>=2.37.2", - "optax>=0.2.0" + "optax>=0.2.0", + "exponax>=0.2.0", +] +pdes = [ + "exponax>=0.2.0", ] [tool.ruff.lint] diff --git a/tests/test_observation.py b/tests/test_observation.py new file mode 100644 index 00000000..186261b0 --- /dev/null +++ b/tests/test_observation.py @@ -0,0 +1,234 @@ +"""Tests for GridInterpolator and interpolation_observation_model.""" + +import jax.numpy as jnp +import numpyro.distributions as dist +import pytest + +from dynestyx.models import DynamicalModel +from dynestyx.observation import GridInterpolator, interpolation_observation_model + +# --- 1D fixtures ----------------------------------------------------------- + +NUM_POINTS = 8 +DOMAIN_EXTENT = 8.0 # unit spacing, dx = 1.0, for easy hand-checked arithmetic +X_GRID = jnp.linspace(0, DOMAIN_EXTENT, NUM_POINTS, endpoint=False) + + +def _values_1d(): + # a smooth-ish, non-constant field so interpolation is a nontrivial check + return jnp.sin(2 * jnp.pi * X_GRID / DOMAIN_EXTENT) + 0.5 * X_GRID + + +def _periodic_linear_interp_reference(x_query, u_values): + """Independent reference implementation (mirrors the notebook's hand-rolled version).""" + dx = X_GRID[1] - X_GRID[0] + idx_float = (x_query % DOMAIN_EXTENT) / dx + idx0 = jnp.floor(idx_float).astype(jnp.int32) % NUM_POINTS + idx1 = (idx0 + 1) % NUM_POINTS + frac = idx_float - jnp.floor(idx_float) + return (1 - frac) * u_values[idx0] + frac * u_values[idx1] + + +# --- basic correctness ------------------------------------------------- + + +@pytest.mark.parametrize("boundary", ["periodic", "edge", "constant"]) +def test_linear_exact_at_grid_points_1d(boundary): + values = _values_1d() + query = X_GRID[:, None] # exactly at every grid point + interp = GridInterpolator((X_GRID,), query, method="linear", boundary=boundary) + assert jnp.allclose(interp(values), values, atol=1e-6) + + +def test_nearest_exact_at_grid_points_1d(): + values = _values_1d() + query = X_GRID[:, None] + interp = GridInterpolator((X_GRID,), query, method="nearest", boundary="periodic") + assert jnp.allclose(interp(values), values, atol=1e-6) + + +def test_nearest_picks_closest_point(): + values = jnp.arange(NUM_POINTS, dtype=jnp.float32) + # query points offset by 0.3*dx from grid point 3 -> should round down to 3 + query = jnp.array([[X_GRID[3] + 0.3]]) + interp = GridInterpolator((X_GRID,), query, method="nearest", boundary="periodic") + assert jnp.allclose(interp(values), jnp.array([3.0])) + # offset by 0.7*dx -> should round up to 4 + query = jnp.array([[X_GRID[3] + 0.7]]) + interp = GridInterpolator((X_GRID,), query, method="nearest", boundary="periodic") + assert jnp.allclose(interp(values), jnp.array([4.0])) + + +def test_linear_matches_manual_periodic_reference(): + values = _values_1d() + key_points = jnp.array([0.0, 1.75, 3.4, 6.999, 7.999, -0.5, 8.5])[:, None] + interp = GridInterpolator( + (X_GRID,), key_points, method="linear", boundary="periodic" + ) + expected = _periodic_linear_interp_reference(key_points[:, 0], values) + assert jnp.allclose(interp(values), expected, atol=1e-6) + + +def test_linear_partition_of_unity_for_bias_free_boundaries(): + query = jnp.array([0.3, 2.6, 5.1, 7.9, -1.2, 9.4])[:, None] + for boundary in ("periodic", "edge"): + interp = GridInterpolator((X_GRID,), query, method="linear", boundary=boundary) + H = interp.as_matrix() + assert jnp.allclose(jnp.sum(H, axis=1), 1.0, atol=1e-6) + assert jnp.allclose(interp.bias(), 0.0) + + +# --- boundary modes ------------------------------------------------------ + + +def test_edge_boundary_clamps_beyond_domain(): + values = _values_1d() + # far outside the domain in both directions + query = jnp.array([[-5.0], [50.0]]) + interp = GridInterpolator((X_GRID,), query, method="linear", boundary="edge") + expected = jnp.array([values[0], values[-1]]) + assert jnp.allclose(interp(values), expected, atol=1e-6) + + +def test_constant_boundary_far_outside_returns_fill_value(): + values = _values_1d() + query = jnp.array([[-50.0], [500.0]]) + interp = GridInterpolator( + (X_GRID,), query, method="linear", boundary="constant", fill_value=3.25 + ) + assert jnp.allclose(interp(values), jnp.array([3.25, 3.25])) + # every corner is invalid this far out, so as_matrix()'s row is all zero and bias + # alone carries the value + assert jnp.allclose(interp.as_matrix(), 0.0) + assert jnp.allclose(interp.bias(), jnp.array([3.25, 3.25])) + + +def test_constant_boundary_partial_straddle_mixes_grid_value_and_fill(): + values = _values_1d() + dx = X_GRID[1] - X_GRID[0] + # query point half a grid cell past the last grid point: one bracketing corner + # (the last grid point) is valid, the other (one past it) is not. + query = jnp.array([[X_GRID[-1] + 0.5 * dx]]) + fill_value = 10.0 + interp = GridInterpolator( + (X_GRID,), query, method="linear", boundary="constant", fill_value=fill_value + ) + expected = 0.5 * values[-1] + 0.5 * fill_value + assert jnp.allclose(interp(values), expected, atol=1e-6) + + +def test_call_matches_matrix_and_bias_for_every_boundary(): + values = _values_1d() + query = jnp.array([0.3, 2.6, 5.1, 7.9, -1.2, 9.4])[:, None] + for boundary in ("periodic", "edge", "constant"): + interp = GridInterpolator( + (X_GRID,), query, method="linear", boundary=boundary, fill_value=-1.0 + ) + reconstructed = interp.as_matrix() @ values + interp.bias() + assert jnp.allclose(interp(values), reconstructed, atol=1e-6) + + +# --- sparse export --------------------------------------------------------- + + +def test_as_sparse_matches_as_matrix(): + query = jnp.array([0.3, 2.6, 5.1, 7.9])[:, None] + interp = GridInterpolator((X_GRID,), query, method="linear", boundary="periodic") + assert jnp.allclose(interp.as_sparse().todense(), interp.as_matrix(), atol=1e-6) + + +# --- N-D generality --------------------------------------------------------- + + +def test_2d_bilinear_exact_on_affine_function(): + nx, ny = 6, 5 + x_axis = jnp.linspace(0.0, 5.0, nx) + y_axis = jnp.linspace(0.0, 4.0, ny) + + def f(x, y): + return 2.0 * x - 3.0 * y + 1.5 + + xs, ys = jnp.meshgrid(x_axis, y_axis, indexing="ij") + values = f(xs, ys).reshape(-1) + + query = jnp.array([[0.5, 0.5], [1.3, 2.7], [4.9, 3.9], [0.0, 0.0]]) + # non-periodic domain (no wraparound for an affine test function): use "edge" so + # in-range queries are unaffected by boundary handling. + interp = GridInterpolator((x_axis, y_axis), query, method="linear", boundary="edge") + expected = f(query[:, 0], query[:, 1]) + assert jnp.allclose(interp(values), expected, atol=1e-5) + + +def test_2d_state_dim_and_matrix_shape(): + nx, ny = 4, 3 + x_axis = jnp.linspace(0.0, 3.0, nx) + y_axis = jnp.linspace(0.0, 2.0, ny) + query = jnp.array([[0.5, 0.5], [1.5, 1.0]]) + interp = GridInterpolator((x_axis, y_axis), query, method="linear", boundary="edge") + assert interp.state_dim == nx * ny + assert interp.as_matrix().shape == (2, nx * ny) + + +# --- validation -------------------------------------------------------- + + +def test_invalid_method_raises(): + # dynestyx runs jaxtyping+typeguard checks package-wide (see pyproject.toml), so an + # out-of-Literal value is rejected as a jaxtyping.TypeCheckError before __init__'s own + # `if method not in _VALID_METHODS` check would run; either way it must raise. + with pytest.raises(Exception): + GridInterpolator((X_GRID,), X_GRID[:1, None], method="cubic") # ty: ignore[invalid-argument-type] + + +def test_invalid_boundary_raises(): + with pytest.raises(Exception): + GridInterpolator((X_GRID,), X_GRID[:1, None], boundary="reflect") # ty: ignore[invalid-argument-type] + + +def test_wrong_query_points_shape_raises(): + with pytest.raises(Exception): + GridInterpolator((X_GRID,), X_GRID) # missing the trailing ndim axis + + +# --- integration with LinearGaussianObservation / DynamicalModel ----------- + + +def test_interpolation_observation_model_builds_valid_dynamical_model(): + values = _values_1d() + query = jnp.array([0.3, 2.6, 5.1])[:, None] + obs_model = interpolation_observation_model( + (X_GRID,), query, R=0.01 * jnp.eye(3), method="linear", boundary="periodic" + ) + + dynamics = DynamicalModel( + initial_condition=dist.MultivariateNormal( + loc=jnp.zeros(NUM_POINTS), covariance_matrix=jnp.eye(NUM_POINTS) + ), + state_evolution=lambda x, u, t_now, t_next: dist.MultivariateNormal( + loc=x, covariance_matrix=jnp.eye(NUM_POINTS) + ), + observation_model=obs_model, + control_dim=0, + ) + + assert dynamics.state_dim == NUM_POINTS + assert dynamics.observation_dim == 3 + + obs_dist = obs_model(values, None, jnp.array(0.0)) + interp = GridInterpolator((X_GRID,), query, method="linear", boundary="periodic") + assert jnp.allclose(obs_dist.mean, interp(values), atol=1e-6) + + +def test_interpolation_observation_model_with_constant_boundary_carries_bias(): + values = _values_1d() + query = jnp.array([[-50.0]]) + obs_model = interpolation_observation_model( + (X_GRID,), + query, + R=jnp.array([[0.01]]), + method="linear", + boundary="constant", + fill_value=7.0, + ) + obs_dist = obs_model(values, None, jnp.array(0.0)) + assert jnp.allclose(obs_dist.mean, jnp.array([7.0]))