diff --git a/README.md b/README.md index e2bb7b3..87241fe 100644 --- a/README.md +++ b/README.md @@ -23,4 +23,7 @@ git clone git@github.com:LSSTDESC/BPD.git cd BPD pip install -e . pip install -e ".[dev]" + +# Might be necessary +conda install -c nvidia cuda-nvcc ``` diff --git a/experiments/exp4/README.md b/experiments/exp4/README.md new file mode 100644 index 0000000..2f85725 --- /dev/null +++ b/experiments/exp4/README.md @@ -0,0 +1,8 @@ +# Experiment 4 + +Same as experiment 3 but we let deviations from true centroid `(dx, dy)` be free parameters in +the fit. + + +For now we assume that the true and interim prior on the deviation is the same, so it does not get +accounted for in the shear inference. diff --git a/experiments/exp4/figs/contours.pdf b/experiments/exp4/figs/contours.pdf new file mode 100644 index 0000000..851bca1 Binary files /dev/null and b/experiments/exp4/figs/contours.pdf differ diff --git a/experiments/exp4/figs/hists.pdf b/experiments/exp4/figs/hists.pdf new file mode 100644 index 0000000..97726ff Binary files /dev/null and b/experiments/exp4/figs/hists.pdf differ diff --git a/experiments/exp4/figs/scatter_dxdy.pdf b/experiments/exp4/figs/scatter_dxdy.pdf new file mode 100644 index 0000000..811e67f Binary files /dev/null and b/experiments/exp4/figs/scatter_dxdy.pdf differ diff --git a/experiments/exp4/figs/scatter_shapes.pdf b/experiments/exp4/figs/scatter_shapes.pdf new file mode 100644 index 0000000..1698cb3 Binary files /dev/null and b/experiments/exp4/figs/scatter_shapes.pdf differ diff --git a/experiments/exp4/figs/traces.pdf b/experiments/exp4/figs/traces.pdf new file mode 100644 index 0000000..1a2c2bf Binary files /dev/null and b/experiments/exp4/figs/traces.pdf differ diff --git a/experiments/exp4/get_interim_samples.py b/experiments/exp4/get_interim_samples.py old mode 100644 new mode 100755 index 7780deb..ed98dea --- a/experiments/exp4/get_interim_samples.py +++ b/experiments/exp4/get_interim_samples.py @@ -1,46 +1,122 @@ #!/usr/bin/env python3 -"""Check chains ran on a variety of galaxies with different SNR, initialization from the prior.""" - -import time from functools import partial import jax.numpy as jnp import typer from jax import jit, random, vmap -from jax._src.prng import PRNGKeyArray from bpd import DATA_DIR -from bpd.chains import run_sampling_nuts, run_warmup_nuts from bpd.draw import draw_gaussian -from bpd.initialization import init_with_prior +from bpd.initialization import init_with_truth +from bpd.io import save_dataset from bpd.pipelines.image_samples import ( get_target_images, get_true_params_from_galaxy_params, - loglikelihood, logprior, - logtarget, + pipeline_interim_samples_one_galaxy, sample_target_galaxy_params_simple, ) -def sample_prior( - rng_key: PRNGKeyArray, - *, - mean_logflux: float = 6.0, - sigma_logflux: float = 0.1, - mean_loghlr: float = 1.0, - sigma_loghlr: float = 0.05, - shape_noise: float = 0.3, +def main( + seed: int, + n_gals: int = 1000, # technically, in this file it means 'noise realizations' + n_samples_per_gal: int = 100, g1: float = 0.02, g2: float = 0.0, -) -> dict[str, float]: - k1, k2, k3 = random.split(rng_key, 3) + lf: float = 6.0, # ~ SNR = 1000 + hlr: float = 1.0, + shape_noise: float = 1e-4, + sigma_e_int: float = 4e-2, + slen: int = 53, + fft_size: int = 256, + background: float = 1.0, + initial_step_size: float = 1e-3, +): + rng_key = random.key(seed) + pkey, nkey, gkey = random.split(rng_key, 3) + + # directory structure + dirpath = DATA_DIR / "cache_chains" / f"exp4_{seed}" + if not dirpath.exists(): + dirpath.mkdir(exist_ok=True) + + # galaxy parameters from prior + pkeys = random.split(pkey, n_gals) + _get_galaxy_params = partial( + sample_target_galaxy_params_simple, g1=g1, g2=g2, shape_noise=shape_noise + ) + galaxy_params = vmap(_get_galaxy_params)(pkeys) + assert galaxy_params["x"].shape == (n_gals,) + assert galaxy_params["e1"].shape == (n_gals,) + assert "lf" not in galaxy_params + extra_params = {"f": 10 ** jnp.full((n_gals,), lf), "hlr": jnp.full((n_gals,), hlr)} - lf = random.normal(k1) * sigma_logflux + mean_logflux - hlr = random.normal(k2) * sigma_loghlr + mean_loghlr + # now get corresponding target images + draw_params = {**galaxy_params, **extra_params} + target_images = get_target_images( + nkey, draw_params, background=background, slen=slen + ) + assert target_images.shape == (n_gals, slen, slen) + + # interim samples are on 'sheared ellipticity' + true_params = vmap(get_true_params_from_galaxy_params)(galaxy_params) - other_params = sample_target_galaxy_params_simple( - k3, shape_noise=shape_noise, g1=g1, g2=g2 + # we pass in x,y as fixed parameters for drawing + # and initialize the function with deviations (dx, dy) = (0, 0) + x = true_params.pop("x") + y = true_params.pop("y") + true_params["dx"] = jnp.zeros_like(x) + true_params["dy"] = jnp.zeros_like(y) + extra_params = {"x": x, "y": y, **extra_params} + + # more setup + _logprior = partial( + logprior, sigma_e=sigma_e_int, free_flux_hlr=False, free_dxdy=True ) - return {"lf": lf, "hlr": hlr, **other_params} + # prepare pipelines + gkeys = random.split(gkey, n_gals) + _draw_fnc = partial(draw_gaussian, slen=slen, fft_size=fft_size) + pipe = partial( + pipeline_interim_samples_one_galaxy, + initialization_fnc=init_with_truth, + draw_fnc=_draw_fnc, + logprior=_logprior, + n_samples=n_samples_per_gal, + initial_step_size=initial_step_size, + background=background, + free_flux=False, + ) + vpipe = vmap(jit(pipe)) + + # compilation on single target image + _ = vpipe( + gkeys[0, None], + {k: v[0, None] for k, v in true_params.items()}, + target_images[0, None], + {k: v[0, None] for k, v in extra_params.items()}, + ) + + samples = vpipe(gkeys, true_params, target_images, extra_params) + e_post = jnp.stack([samples["e1"], samples["e2"]], axis=-1) + fpath = dirpath / f"e_post_{seed}.npz" + + save_dataset( + { + "e_post": e_post, + "true_g": jnp.array([g1, g2]), + "dx": samples["dx"], + "dy": samples["dy"], + "sigma_e": shape_noise, + "sigma_e_int": sigma_e_int, + "e1": draw_params["e1"], + "e2": draw_params["e2"], + }, + fpath, + overwrite=True, + ) + + +if __name__ == "__main__": + typer.run(main) diff --git a/experiments/exp4/get_posteriors.sh b/experiments/exp4/get_posteriors.sh new file mode 100755 index 0000000..fcfebf9 --- /dev/null +++ b/experiments/exp4/get_posteriors.sh @@ -0,0 +1,7 @@ +#!/bin/bash +export CUDA_VISIBLE_DEVICES="0" +export JAX_ENABLE_X64="True" +SEED="43" + +./get_interim_samples.py $SEED +../../scripts/get_shear_from_interim_samples.py $SEED exp4_$SEED "e_post_${SEED}.npz" --overwrite diff --git a/experiments/exp4/make_figures.py b/experiments/exp4/make_figures.py new file mode 100755 index 0000000..8aed6b7 --- /dev/null +++ b/experiments/exp4/make_figures.py @@ -0,0 +1,154 @@ +#!/usr/bin/env python3 + +import os + +os.environ["CUDA_VISIBLE_DEVICES"] = "" +os.environ["JAX_PLATFORMS"] = "cpu" +os.environ["JAX_ENABLE_X64"] = "True" + + +import jax.numpy as jnp +import matplotlib.pyplot as plt +import numpy as np +import typer +from jax import Array +from matplotlib.backends.backend_pdf import PdfPages + +from bpd import DATA_DIR +from bpd.diagnostics import get_contour_plot +from bpd.io import load_dataset + + +def make_trace_plots(g_samples: Array) -> None: + """Make trace plots of g1, g2.""" + fname = "figs/traces.pdf" + with PdfPages(fname) as pdf: + fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(10, 5)) + g1 = g_samples[:, 0] + g2 = g_samples[:, 1] + + ax1.plot(g1) + ax2.plot(g2) + + pdf.savefig(fig) + plt.close(fig) + + +def make_contour_plots(g_samples: Array, n_examples=10) -> None: + """Make figure of contour plot on g1, g2.""" + fname = "figs/contours.pdf" + with PdfPages(fname) as pdf: + truth = {"g1": 0.02, "g2": 0.0} + g_dict = {"g1": g_samples[:, 0], "g2": g_samples[:, 1]} + fig = get_contour_plot([g_dict], ["post"], truth) + pdf.savefig(fig) + plt.close(fig) + + +def make_scatter_shape_plots(e_post: Array, n_examples: int = 10) -> None: + """Show example scatter plots of interim posterior ellipticitites.""" + # make two types, assuming gaussianity and one not assuming gaussianity. + fname = "figs/scatter_shapes.pdf" + + n_gals, _, _ = e_post.shape + + with PdfPages(fname) as pdf: + # individual + for _ in range(n_examples): + idx = np.random.choice(np.arange(0, n_gals)) + e1, e2 = e_post[idx, :, 0], e_post[idx, :, 1] + fig, ax = plt.subplots(1, 1, figsize=(7, 7)) + ax.scatter(e1, e2, marker="x") + ax.set_title(f"Samples ellipticity index: {idx}") + ax.set_xlabel("e1", fontsize=14) + ax.set_ylabel("e2", fontsize=14) + pdf.savefig(fig) + plt.close(fig) + + # clusters + n_clusters = 50 + fig, ax = plt.subplots(1, 1, figsize=(7, 7)) + ax.set_xlabel("e1", fontsize=14) + ax.set_ylabel("e2", fontsize=14) + fig.suptitle(f"{n_clusters} galaxies plotted") + for _ in range(n_clusters): + idx = np.random.choice(np.arange(0, n_gals)) + e1, e2 = e_post[idx, :, 0], e_post[idx, :, 1] + ax.scatter(e1, e2, marker="x") + pdf.savefig(fig) + plt.close(fig) + + +def make_scatter_dxdy_plots(dx: Array, dy: Array, n_examples: int = 10) -> None: + """Show example scatter plots of interim posterior ellipticitites.""" + # make two types, assuming gaussianity and one not assuming gaussianity. + fname = "figs/scatter_dxdy.pdf" + + n_gals, _ = dx.shape + + with PdfPages(fname) as pdf: + # individual + for _ in range(n_examples): + idx = np.random.choice(np.arange(0, n_gals)) + dx1, dy1 = dx[idx, :], dy[idx, :] + fig, ax = plt.subplots(1, 1, figsize=(7, 7)) + ax.scatter(dx1, dy1, marker="x") + ax.set_title(f"Samples ellipticity index: {idx}") + ax.set_xlabel("dx", fontsize=14) + ax.set_ylabel("dy", fontsize=14) + pdf.savefig(fig) + plt.close(fig) + + # clusters + n_clusters = 50 + fig, ax = plt.subplots(1, 1, figsize=(7, 7)) + ax.set_xlabel("dx", fontsize=14) + ax.set_ylabel("dy", fontsize=14) + fig.suptitle(f"{n_clusters} galaxies plotted") + for _ in range(n_clusters): + idx = np.random.choice(np.arange(0, n_gals)) + dx1, dy1 = dx[idx, :], dy[idx, :] + ax.scatter(dx1, dy1, marker="x") + pdf.savefig(fig) + plt.close(fig) + + +def make_hists(g_samples: Array, e1_samples: Array) -> None: + """Make histograms of g1 along with std and expected std.""" + fname = "figs/hists.pdf" + with PdfPages(fname) as pdf: + fig, ax = plt.subplots(1, 1, figsize=(7, 7)) + + g1 = g_samples[:, 0] + e1_std = e1_samples.std() + g1_exp_std = e1_std / jnp.sqrt(len(e1_samples)) + + ax.hist(g1, bins=25, histtype="step") + ax.axvline(g1.mean(), linestyle="--", color="k") + ax.set_title(f"Std g1: {g1.std():.4g}; Expected g1 std: {g1_exp_std:.4g}") + + pdf.savefig(fig) + plt.close(fig) + + +def main(seed: int = 43): + # load data + pdir = DATA_DIR / "cache_chains" / f"exp4_{seed}" + e_post_dict = load_dataset(pdir / f"e_post_{seed}.npz") + e_post_samples = e_post_dict["e_post"] + g_samples = jnp.load(pdir / f"g_samples_{seed}_{seed}.npy") + + e1_samples = e_post_dict["e1"] + dx = e_post_dict["dx"] + dy = e_post_dict["dy"] + + # make plots + make_scatter_shape_plots(e_post_samples) + make_scatter_dxdy_plots(dx, dy) + make_trace_plots(g_samples) + make_contour_plots(g_samples) + make_hists(g_samples, e1_samples) + + +if __name__ == "__main__": + typer.run(main) diff --git a/notebooks/exp4_check1.ipynb b/notebooks/exp4_check1.ipynb new file mode 100644 index 0000000..08096b2 --- /dev/null +++ b/notebooks/exp4_check1.ipynb @@ -0,0 +1,175 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 1, + "id": "4a213c34-b49c-41ce-86e2-442af71fff46", + "metadata": { + "tags": [] + }, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/pscratch/sd/i/imendoza/miniconda3/envs/bpd_gpu3/lib/python3.12/site-packages/tqdm/auto.py:21: TqdmWarning: IProgress not found. Please update jupyter and ipywidgets. See https://ipywidgets.readthedocs.io/en/stable/user_install.html\n", + " from .autonotebook import tqdm as notebook_tqdm\n" + ] + } + ], + "source": [ + "#!/usr/bin/env python3\n", + "\n", + "import os\n", + "\n", + "os.environ[\"CUDA_VISIBLE_DEVICES\"] = \"\"\n", + "os.environ[\"JAX_PLATFORMS\"] = \"cpu\"\n", + "os.environ[\"JAX_ENABLE_X64\"] = \"True\"\n", + "\n", + "import jax.numpy as jnp\n", + "import matplotlib.pyplot as plt\n", + "import numpy as np\n", + "import typer\n", + "from jax import Array\n", + "from matplotlib.backends.backend_pdf import PdfPages\n", + "\n", + "from bpd import DATA_DIR\n", + "from bpd.diagnostics import get_contour_plot\n", + "from bpd.io import load_dataset\n" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "edea5496-4787-4e10-bb0e-f0c1369b589c", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "data_dict = load_dataset(\"/pscratch/sd/i/imendoza/data/cache_chains/exp4_42/e_post_42.npz\")" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "id": "889c5a4c-4625-4919-80a0-6b231c0d42e5", + "metadata": { + "tags": [] + }, + "outputs": [ + { + "data": { + "text/plain": [ + "dict_keys(['e_post', 'true_g', 'dx', 'dy', 'sigma_e', 'sigma_e_int', 'e1', 'e2'])" + ] + }, + "execution_count": 8, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "data_dict.keys()" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "id": "2a50e56d-2151-469c-8cf9-158c20824e9f", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "dx = data_dict['dx']\n", + "dy = data_dict['dy']" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "id": "09d935ac-795d-4fad-935b-7f6021341eed", + "metadata": { + "tags": [] + }, + "outputs": [ + { + "data": { + "text/plain": [ + "((1000, 100), (1000, 100))" + ] + }, + "execution_count": 5, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "dx.shape, dy.shape" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "id": "c6a91f9e-a7b9-4d68-b36b-69c820456cd8", + "metadata": { + "tags": [] + }, + "outputs": [ + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 7, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkkAAAGdCAYAAAAGx+eQAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjAsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvlHJYcgAAAAlwSFlzAAAPYQAAD2EBqD+naQAANrJJREFUeJzt3X9wFGWe+PFPhpCQGCYBIRnRoFCwsPw4QZDsAK7nmSLBlCtK1bK5rAssyqK4dygHC3cm5HQ9OHX39vRYf9x9JVbdKpqqU1eEeDnCjzOJASJRIZHSFQSECSfZZAIL5Md8vn/Mps1AJ0wymR89835VdQ3T/XTneTpD+jNPf56n41RVBQAAAD5s4a4AAABAJCJIAgAAMEGQBAAAYIIgCQAAwARBEgAAgAmCJAAAABMESQAAACYIkgAAAEzEh7sCVubxeOTUqVMydOhQiYuLC3d1AACAH1RVWltbZdSoUWKz9dxfRJAUgFOnTklmZma4qwEAAPrhxIkTcsMNN/S4nSApAEOHDhUR70m22+1hrg0AAPCH2+2WzMxM4zreE4KkAHTdYrPb7QRJAABYzNVSZUjcBgAAMEGQBAAAYIIgCQAAwARBEgAAgAmCJAAAABMESQAAACYIkgAAAEwQJAEAAJggSAIAADBBkAQAAGCCIAkAAMAEQRKsz+MRKSoSqaryvno84a4RACAK8IBbWF9xsciTT3qXLk88EbbqAACiAz1JsL7c3N7fAwDQDwRJsL6yst7fAwDQD9xug/UVF3tfc3O9AVLX+4Hg8XiP1/3YNr5bAEAsiFNVDXclrMrtdktqaqq0tLSI3W4Pd3UQDEVFvrlOhYXkOwGAxfl7/eYrMdAb8p0AIGYRJAG9Id8JAGIWOUlAb4KZ7wQAiGjkJAWAnCQAAKyHnCQAAIAAECQBAACYIEhCbOJ5bwCAqyBxG7GJ570BAK6CniTEJuY/QiygxxQICD1JiE1m8x/Nnh2eugDBQo8pEBCCJMQm5j9CLMjN9Q2Q6DEF+oQgCbHJZvv2GzU9SIhW9JgCASFIAoBoRY8pEBBm3A4AM24DAGA9zLgNAAAQAIIkAAAAEwRJAAAAJgiSAAAATBAkAQAAmCBIAgAAMEGQBAAAYIIgCQAAwARBEgAAgAmCJAAAABMESQAAACYIkgAAAEwQJAEAAJggSAIAADBBkAQAAGCCIAkAAMBE0IOkzZs3y0033SRDhgyRrKws2bdvX6/lS0tLZeLEiTJkyBCZOnWqbN++3We7qkpRUZFcd911kpSUJNnZ2fL555/7lGlqapKCggKx2+2SlpYmy5Ytk3PnzhnbL168KEuWLJGpU6dKfHy8LFiwYMDaCwAAokNQg6Q33nhDHnvsMdmwYYN89NFHcvPNN0tOTo6cOXPGtHxVVZXk5+fLsmXL5ODBg7JgwQJZsGCBHDp0yCjz9NNPy3PPPScvvvii1NTUyDXXXCM5OTly8eJFo0xBQYEcPnxYysvLZdu2bbJ3715Zvny5sb2zs1OSkpLkb/7mbyQ7Ozt4JwAAAFhWnKpqsA6elZUlt956q/zbv/2biIh4PB7JzMyUn//857Ju3boryi9atEjOnz8v27ZtM9Z973vfk2nTpsmLL74oqiqjRo2S1atXy9/93d+JiEhLS4tkZGRISUmJ/OhHP5KGhgaZNGmS7N+/X2bOnCkiImVlZXLXXXfJyZMnZdSoUT4/c8mSJdLc3Cxvv/12n9vndrslNTVVWlpaxG6393l/AAAQev5ev4PWk9TW1ia1tbU+PTU2m02ys7OlurradJ/q6uorenZycnKM8kePHhWXy+VTJjU1VbKysowy1dXVkpaWZgRIIiLZ2dlis9mkpqYmoDZdunRJ3G63zwIAAKJT0IKkb775Rjo7OyUjI8NnfUZGhrhcLtN9XC5Xr+W7Xq9WJj093Wd7fHy8DB8+vMef66+NGzdKamqqsWRmZgZ0PAAAELkY3dYH69evl5aWFmM5ceJEuKsEAACCJGhB0ogRI2TQoEHS2Njos76xsVEcDofpPg6Ho9fyXa9XK3N5YnhHR4c0NTX1+HP9lZiYKHa73WcBEIU8HpGiIpGqKu+rxxPuGgEIg6AFSQkJCTJjxgzZuXOnsc7j8cjOnTvF6XSa7uN0On3Ki4iUl5cb5ceMGSMOh8OnjNvtlpqaGqOM0+mU5uZmqa2tNcpUVFSIx+ORrKysAWsfgChWXCzy5JMic+Z4X4uLw10jAGEQH8yDP/bYY7J48WKZOXOmzJo1S37zm9/I+fPnZenSpSIi8pOf/ESuv/562bhxo4iI/O3f/q3cfvvt8qtf/Ury8vJk69atcuDAAXn55ZdFRCQuLk5WrVolv/zlL2X8+PEyZswYKSwslFGjRhlzHX33u9+V3NxcefDBB+XFF1+U9vZ2eeSRR+RHP/qRz8i2+vp6aWtrk6amJmltbZW6ujoREZk2bVowTwkAK8jN9QZH3d8DiDlBDZIWLVok//d//ydFRUXicrlk2rRpUlZWZiReHz9+XGy2bzuzZs+eLa+99po8/vjj8vd///cyfvx4efvtt2XKlClGmbVr18r58+dl+fLl0tzcLHPnzpWysjIZMmSIUeZ3v/udPPLII3LnnXeKzWaThQsXynPPPedTt7vuuku++uor4/306dNFxDtZJYAYV1Z25fvZs8NTFwBhE9R5kqId8yQBUcrj8d5iy831BkjFxSI2xrkA0cLf6zdBUgAIkgAAsJ6wTyYJAABgZQRJAAAAJgiSAAAATBAkAQAAmCBIAgAAMEGQBAAAYIIgCQAAwARBEgAAgAmCJAAAABMESQAAACYIkgAAAEwQJAEAAJggSAIAADBBkAQAAGCCIAmAtXg8IkVFIlVV3lePJ9w1AhCl4sNdAQDok+JikSef9C5dnngibNUBEL3oSQJgLbm5vb8HgAFCkATAWsrKen8PAAOE220ArKW42Puam+sNkLreA8AAi1NVDXclrMrtdktqaqq0tLSI3W4Pd3UAAIAf/L1+c7sNAADABEESAAQTUxYAlkVOEgAEE1MWAJZFTxIABBNTFgCWRZAEAMHElAWAZXG7DegPj8d7G6X7MHQb3zlggikLAMtiCoAAMAVADCsq8s0xKSwkzwQALIIpAIBgIs8EAKIeQRLQH+SZAEDUIycJ6A/yTAAg6pGTFABykgAAsB5ykgAAAAJAkAQAAGCCIAkAAMAEQRIAAIAJgiQAGGgej3fC0aoq76vHE+4aAegHpgAAgIFWXOydkb37rOzMyA5YDj1JADDQmJEdiAoESQAw0JiRHYgK3G4DAH95PN5bad1nWreZfNdkRnYgKjDjdgBiesZtfy8WQDQpKvLNMyosJNcIsCBm3EZwdSWmzpnjfeWbMmIBuUZATCFIQv9wsUAsItcIiCnkJKF/zC4Ws2eHpy5AqJBrBMQUcpICQE5SMTlJAADL8ff6TZAUgJgOkgAAsCgStwEAAAJAkAQAAGCCIAkAAMAEQRIAAIAJgiQAAAATBEkAAAAmQhIkbd68WW666SYZMmSIZGVlyb59+3otX1paKhMnTpQhQ4bI1KlTZfv27T7bVVWKiorkuuuuk6SkJMnOzpbPP//cp0xTU5MUFBSI3W6XtLQ0WbZsmZw7d86nzCeffCK33XabDBkyRDIzM+Xpp58emAYDAADLC3qQ9MYbb8hjjz0mGzZskI8++khuvvlmycnJkTNnzpiWr6qqkvz8fFm2bJkcPHhQFixYIAsWLJBDhw4ZZZ5++ml57rnn5MUXX5Samhq55pprJCcnRy5evGiUKSgokMOHD0t5ebls27ZN9u7dK8uXLze2u91umTdvntx4441SW1srzzzzjBQXF8vLL78cvJMBAACsQ4Ns1qxZunLlSuN9Z2enjho1Sjdu3Gha/oc//KHm5eX5rMvKytKf/exnqqrq8XjU4XDoM888Y2xvbm7WxMREff3111VVtb6+XkVE9+/fb5TZsWOHxsXF6ddff62qqr/97W912LBheunSJaPML37xC50wYYLfbWtpaVER0ZaWFr/3AQAA4eXv9TuoPUltbW1SW1sr2dnZxjqbzSbZ2dlSXV1tuk91dbVPeRGRnJwco/zRo0fF5XL5lElNTZWsrCyjTHV1taSlpcnMmTONMtnZ2WKz2aSmpsYo8/3vf18SEhJ8fs6RI0fkj3/8o2ndLl26JG6322cBAADRKahB0jfffCOdnZ2SkZHhsz4jI0NcLpfpPi6Xq9fyXa9XK5Oenu6zPT4+XoYPH+5TxuwY3X/G5TZu3CipqanGkpmZad5wAABgeYxu64P169dLS0uLsZw4cSLcVQIAAEES1CBpxIgRMmjQIGlsbPRZ39jYKA6Hw3Qfh8PRa/mu16uVuTwxvKOjQ5qamnzKmB2j+8+4XGJiotjtdp8FAABEp6AGSQkJCTJjxgzZuXOnsc7j8cjOnTvF6XSa7uN0On3Ki4iUl5cb5ceMGSMOh8OnjNvtlpqaGqOM0+mU5uZmqa2tNcpUVFSIx+ORrKwso8zevXulvb3d5+dMmDBBhg0bFmDLAQCA5QU7g3zr1q2amJioJSUlWl9fr8uXL9e0tDR1uVyqqnr//ffrunXrjPKVlZUaHx+vzz77rDY0NOiGDRt08ODB+umnnxplNm3apGlpafrOO+/oJ598ovfcc4+OGTNGL1y4YJTJzc3V6dOna01NjX7wwQc6fvx4zc/PN7Y3NzdrRkaG3n///Xro0CHdunWrJicn60svveR32xjdBsvq7FQtLFStrPS+dnaGu0YAEDL+Xr+DHiSpqj7//PM6evRoTUhI0FmzZumHH35obLv99tt18eLFPuXffPNN/c53vqMJCQk6efJkfe+993y2ezweLSws1IyMDE1MTNQ777xTjxw54lPm7Nmzmp+frykpKWq323Xp0qXa2trqU+bjjz/WuXPnamJiol5//fW6adOmPrWLIAmWVVioKvLtUlgY7hoBQMj4e/2OU1UNb1+WdbndbklNTZWWlhbyk2AtVVUic+Z8+76yUmT27PDVBwBCyN/rN6PbgFhUVtb7ewCAxIe7AgDCoLjY+5qb6w2Qut4DAAzcbgsAt9sAC/F4vMFg98DQRmc6EIu43QbECo9HpKjIm2dUVOR9jysVF4s8+aQ3F+vJJ+k9A3BVBEmA1XHx909ubu/vAeAyBEmA1XHx9w/J6gD6iMRtwOrMLv4M578SyeoA+ojE7QCQuI2IQEIyAPSJv9dvgqQAECQBAGA9jG4DAAAIAEESAACACYIkAAAAEwRJAAAAJgiSAAAATBAkAQAAmCBIAjAweIYcgChDkASEQiwEEDxDDkCUIUgCQmHDBt8AYsOGcNdo4PEMOQBRhiAJCIWTJ3t/Hw14gCyAKMMDbgH0rC/PheMBsgCiDEESEAo33ND7+0jVlWf05JPfrnviCfOyNtu322bPDnrVACDYuN0WaWIhwTcW/eM/ihQWilRWel//8R/DXSP/kGcEIIbRkxRp+vLNHdZh1V4WszwjK9UfAAJAT1Kk4Zs7IklxsW8PmNXzjOipBdAH9CRFGr65R5a+JC5HI6v2gPWEnloAfUCQFGkYIRRZuKhGl9xc398lPbUAehFDX4ktouub++zZ3ler9lpEy22NUN7+7O85G+hzHS2/OzPM5QSgD+hJQnBESw9MKG9/9vecDfS5jpbfnRl6agH0gUW7KRDxoiUBPZSJy/09ZwN9rqPld2cmWnpqAYQEfyEQHNFyWyOUF9X+nrOBPtfR8rsDgABxuw3BwW2NvuvvORvoc83vDgBERCROVTXclbAqt9stqamp0tLSIna7feB/QKwPP480/D4AICr4e/2mJymSRXMCrRXx+wCAmMLX4EgWzQm0VsTvAwBiCkFSJCOBNrLw+4hu0Tw/FIB+4XZbJIuVBFqr5PrEyu8jVnE7FcBlSNwOQNATt2NFUZHvhamwkIsTQq+qSmTOnG/fL1ki8v/+X2QG7AAC4u/1m//9CD9yfRAJLr99WlJCbyEQ4wiSEH7k+iASFBd7e4+6I2AHYho5SQg/cn0QCWw2kcxM33XBfFYfgIhHTlIAyEkCooxVBhEA0S7I/xf9vX4TJAWAIAmIMQRRQGgEeUAPidsAMNC6pgmYM8f7yq1hIDgiZEAPQRIA+CtC/nADUS9CBvSQuA0A/jL7w01iNzDwImRADzlJASAn6c/I00Cs4LMORAUSt0OAIOnPmDEbQCgQpGKAkLiN0CFPA0AokDiPECNIQuAiJMEOQJTjCxlCjMRtBC5CEuwARDkS560him6LkpMUAHKSACCEoujiG9UskKdK4nYIECQBAHCZqipv3liXysqI6/ELe+J2U1OTFBQUiN1ul7S0NFm2bJmcO3eu130uXrwoK1eulGuvvVZSUlJk4cKF0tjY6FPm+PHjkpeXJ8nJyZKeni5r1qyRjo4OnzK7d++WW265RRITE2XcuHFSUlLis33v3r1y9913y6hRoyQuLk7efvvtgWgyAACIojzVoAVJBQUFcvjwYSkvL5dt27bJ3r17Zfny5b3u8+ijj8q7774rpaWlsmfPHjl16pTcd999xvbOzk7Jy8uTtrY2qaqqkldffVVKSkqkqKjIKHP06FHJy8uTO+64Q+rq6mTVqlXywAMPyPvvv2+UOX/+vNx8882yefPmgW84AACxrLjYe4utstL7auU8VQ2C+vp6FRHdv3+/sW7Hjh0aFxenX3/9tek+zc3NOnjwYC0tLTXWNTQ0qIhodXW1qqpu375dbTabulwuo8wLL7ygdrtdL126pKqqa9eu1cmTJ/sce9GiRZqTk2P6c0VE33rrrX61s6WlRUVEW1pa+rV/1OrsVC0sVK2s9L52dkbGsQAAUP+v30HpSaqurpa0tDSZOXOmsS47O1tsNpvU1NSY7lNbWyvt7e2SnZ1trJs4caKMHj1aqqurjeNOnTpVMjIyjDI5OTnidrvl8OHDRpnux+gq03WMQFy6dEncbrfPAhMDOZcJ86IAAMIkKEGSy+WS9PR0n3Xx8fEyfPhwcblcPe6TkJAgaWlpPuszMjKMfVwul0+A1LW9a1tvZdxut1y4cKHfbRIR2bhxo6SmphpLZmZmQMeLWgM5lwnzogAAwqRPQdK6deskLi6u1+Wzzz4LVl3Dbv369dLS0mIsJ06cCHeVItNAJu1FUQIgAMBa+jSZ5OrVq2XJkiW9lhk7dqw4HA45c+aMz/qOjg5pamoSh8Nhup/D4ZC2tjZpbm726U1qbGw09nE4HLJv3z6f/bpGv3Uvc/mIuMbGRrHb7ZKUlHTVNvYmMTFREhMTAzpGTBjIySWZqBIAECZ9CpJGjhwpI0eOvGo5p9Mpzc3NUltbKzNmzBARkYqKCvF4PJKVlWW6z4wZM2Tw4MGyc+dOWbhwoYiIHDlyRI4fPy5Op9M47lNPPSVnzpwxbueVl5eL3W6XSZMmGWW2b9/uc+zy8nLjGAgBm+3bicMCnRtjII+FgcOkfgBiQbAyx3Nzc3X69OlaU1OjH3zwgY4fP17z8/ON7SdPntQJEyZoTU2NsW7FihU6evRoraio0AMHDqjT6VSn02ls7+jo0ClTpui8efO0rq5Oy8rKdOTIkbp+/XqjzJdffqnJycm6Zs0abWho0M2bN+ugQYO0rKzMKNPa2qoHDx7UgwcPqojor3/9az148KB+9dVXfWojo9sQswoLVUW+XQoLw10jAPCbv9fvoAVJZ8+e1fz8fE1JSVG73a5Lly7V1tZWY/vRo0dVRHTXrl3GugsXLujDDz+sw4YN0+TkZL333nv19OnTPsc9duyYzp8/X5OSknTEiBG6evVqbW9v9ymza9cunTZtmiYkJOjYsWN1y5YtV2wXkSuWxYsX96mNBEmIWZWVvkFSZWW4awQAfvP3+s1jSQLAY0kQVfpyC80Cz2YCgJ74e/3uU04SgCjWNSdV9+Cnp8CHhHoAMYCepADQk4SoYoGHUgLAQAj7A24BWAxzUgGAD4IkeHk83jyTqirvq8cT7hoh1ML9UEo+gwAiDLfbAhBVt9tIxEW48RkEECLcbkPf8Iw0XC7UPTt8BgFEGIIkeJGPgst1jXabM8f7Guzbb3wGAUQYpgCAF0O6cbncXN/bX8Hu2eEzGLl4DA1iFEESvHhGGi5n1rMTzM9GJH4GrRIcBLuefZlDC4giBEkAzNGzY53gINj1DHWvIhAhIvArEYCI0NWzM3u29zUSe1CCzSrJ5MGuJ/liiFEx+FcPAPwUjuCgP6MKg13PcM+hBYQJ8yQFIKrmSQJwpXDkJPVnviir5E4BEcLf6zdBUgAIkgAMOJ6hBwQdk0kCgBWR/wNEDEa3AUAkiYZRhdz+Q5TgdlsAuN0GACZ4Dh8iHLfbAADhYZWpE4CrIEgCAAws8qoQJchJAhB65KxEt2jIqwKEnKSAkJME9BM5KwDCiJwkAJGLnBUAFkCQBCD0yFkBYAHkJAEIPXJWAFgAOUkBICcJAADrIScJAAAgAARJAAAAJgiSAAAATBAkwVo8Hu8cO1VV3lePJ9w1AgBEKUa3wVqKi72TEHafiJBJCAEAQUBPEqyFSQgBL3pVEUx8vkSEniRYjdkkhLNnh6cuQDjRq4pg4vMlIvQkwWqKi73P+aqs9L4yCSFiFb2qCCY+XyJCTxKsxmb79tsMPUiIZfSqIpj4fIkIQRIAWBOPdkEw8fkSER5LEhAeSwIAgPXwWBIAAIAAECQBAACYIEgCAAAwQZAEAABggiAJQN8wE685zgsQdZgCAEDfMBOvuWg+Lx6Pt33dh4Pb+I6N6MenHNbFN/fwYCZecz2dl2j4nHYFgHPmeF9jdM4cxB6CJFgXf7jDw2wm3kgVygClp/MSDZ9TAmPEKG63wbpyc31vbfCHOzSsNBNvKG+B9XReouFzumPHle9j8BEViD0ESbAuni0UHlZ6fl4oA5SezgufU8CyCJJgXVbq0cC3QpkEHAkBSnGxiKrIyZPe96rec2ClxOf580V++Uvf90AMIEiCdVmpRwPfioRbYP3VnwDPZhOJixMpKfl2XVyctUa+RUKwCYQBD7gNAA+4BfqhqsqbxNylstI6F9yiIt/grrDQv2AnUtrc3148pgBAlOEBt4g80TAUGoGz0ui4y/V3lFektLm/I+26em1nz/a+EiCFHn8/w4LbbQidaJ5sD/6zci5Zf287RUqbo2GkXazi72dYECQhdPgDDRFr55L1N9iJlDaTW2Rd/P0MC4IkhA5/oGF1kRLs9Fek9Gih7/j7GRZBvbHc1NQkBQUFYrfbJS0tTZYtWybnzp3rdZ+LFy/KypUr5dprr5WUlBRZuHChNDY2+pQ5fvy45OXlSXJysqSnp8uaNWuko6PDp8zu3bvllltukcTERBk3bpyUdB9ZIiIbN26UW2+9VYYOHSrp6emyYMECOXLkyIC0Gz0oLvYmulZWel/5Aw2EFrlF1sXfz7AI6v+QgoICOXz4sJSXl8u2bdtk7969snz58l73efTRR+Xdd9+V0tJS2bNnj5w6dUruu+8+Y3tnZ6fk5eVJW1ubVFVVyauvviolJSVSVFRklDl69Kjk5eXJHXfcIXV1dbJq1Sp54IEH5P333zfK7NmzR1auXCkffvihlJeXS3t7u8ybN0/Onz8/8CcCXvyBhlWQJItIw9/P8NAgqa+vVxHR/fv3G+t27NihcXFx+vXXX5vu09zcrIMHD9bS0lJjXUNDg4qIVldXq6rq9u3b1WazqcvlMsq88MILarfb9dKlS6qqunbtWp08ebLPsRctWqQ5OTk91vfMmTMqIrpnzx6/29jS0qIioi0tLX7vg37q7FQtLFStrPS+dnaGu0aIZoWFqt5pH71LYWG4awRgAPl7/Q5aKFpdXS1paWkyc+ZMY112drbYbDapqakx3ae2tlba29slOzvbWDdx4kQZPXq0VFdXG8edOnWqZGRkGGVycnLE7XbL4cOHjTLdj9FVpusYZlpaWkREZPjw4T2WuXTpkrjdbp8FIRINDwmFdfBAVwASxNttLpdL0tPTfdbFx8fL8OHDxeVy9bhPQkKCpKWl+azPyMgw9nG5XD4BUtf2rm29lXG73XLhwoUrfq7H45FVq1bJnDlzZMqUKT22aePGjZKammosmZmZPZbFAOOihVCKlHmNAIRVn4OkdevWSVxcXK/LZ599Foy6Bs3KlSvl0KFDsnXr1l7LrV+/XlpaWozlxIkTIaohuGghpEiSBSD9mAJg9erVsmTJkl7LjB07VhwOh5w5c8ZnfUdHhzQ1NYnD4TDdz+FwSFtbmzQ3N/v0JjU2Nhr7OBwO2bdvn89+XaPfupe5fERcY2Oj2O12SUpK8ln/yCOPGEnlN9xwQ6/tSkxMlMTExF7LIEgYuoxQsvpQfwADos9B0siRI2XkyJFXLed0OqW5uVlqa2tlxowZIiJSUVEhHo9HsrKyTPeZMWOGDB48WHbu3CkLFy4UEZEjR47I8ePHxel0Gsd96qmn5MyZM8btvPLycrHb7TJp0iSjzPbt232OXV5ebhxDRERV5ec//7m89dZbsnv3bhkzZkwfzwRCiotWbOAZYQAiSTCzx3Nzc3X69OlaU1OjH3zwgY4fP17z8/ON7SdPntQJEyZoTU2NsW7FihU6evRoraio0AMHDqjT6VSn02ls7+jo0ClTpui8efO0rq5Oy8rKdOTIkbp+/XqjzJdffqnJycm6Zs0abWho0M2bN+ugQYO0rKzMKPPQQw9pamqq7t69W0+fPm0sf/rTn/xuH6PbgAHGqDIAIeDv9TuoQdLZs2c1Pz9fU1JS1G6369KlS7W1tdXYfvToURUR3bVrl7HuwoUL+vDDD+uwYcM0OTlZ7733Xj19+rTPcY8dO6bz58/XpKQkHTFihK5evVrb29t9yuzatUunTZumCQkJOnbsWN2yZYvPdhExXS4v1xuCJGCAVVb6BkmVleGuke/0E48/7l2YiqL/mM4DEcDf63ecqmqYOrEsz+12S2pqqrS0tIjdbg93dQDrKyryfT5VYWH4H+J5eZ26C2f9rHprMhJ/x4g5/l6/LfA/CjGLWY9jTySOKuttuolwTkVh1bnDmM4DFkKQhMhl1YsA+m+gH70wEIF2b9NNhHMqCqsGG0znAQvp8+g2IGRyc3275a1yEUDk6Aq0u3+O+nprp/v0Ezt2eP89f374p6Kw6lPhmc4DFkKQhMhl1YsAIsdABNo9TT8R7s+iVYMNpvOAhRAkIXJZ9SKA/glGInJ/A20rJEUTbABBx+i2ADC6LYZZ4SIaSfw5X8EY9dTf3xMjsICoxug2IJhIKu8bf85XMBKR+5sIbtWk6EjEKFVYGEES0B+xehHt7wXPn/MVSaOeIqkuVscXClgYOUlAf8RqUnl/R4v5c74iKQetL3Xh1mvvGKUKCyNIAvojki7oodTfC54/5yuSEpH7UpeBmGYgmsXqFwpEBYIkoD8i6YIeSv294EXz+aKnpHex+oUCUYEgCYD/uOBdiZ6S3kVzgIyoxxQAAWAKAADkJAHW4+/1myApAARJAICIRyB/BX+v39xuAwAgmjG4oN9iO5QEACDaxeq8bgOAIAkAgGjG5Kj9xu02AACiGaNS+43E7QCQuA0AgPXwgFsAAIAAECQBAACYIEgCetLfJ94DsCb+z+MyJG4DPWFukcjGBHkYaPyfx2X4iwL0hLlFIlvXBW3OHO8rI3YQKP7P4zIESUBPmFsksnFBw0Dj/zwuw+02oCfMLRLZzC5oPGUefdX9tq2qyOOPi8yfz/95iAjzJAWEeZKAMCInCQOhqMg3B6mwkDykGODv9ZsgKQAESQBgcVVV3ry2LpWV9EjGACaTBADgashDQi/ISQIAxC5yD9ELbrcFgNttAABYD7fbAAAAAkCQBAAAYIIgCQAAwARBEgAAgAmCJAAAABMESQAQTB6Pd1bnqirvq8cT7hoB8BPzJAFAMBUXex970f3RFzz2ArAEepIAIJhyc3t/DyBiESQBQDDx2AvAsrjdBgDBxGMvAMvisSQB4LEkAABYD48lAQAACABBEgAAgAmCJCDWMG8PAPiFxG0g1jBvDwD4hZ4koLtQ97KEo1eHeXsAwC/0JAHdhbqXJRy9Ombz9syeHdyfCQAWRE8S0F2oe1nC0atTXCxSWChSWel9Zd4eADBFTxLQXah7WcLRq2OzfdtbRQ8SAPSIIAnoLtSzIzMbMwBELGbcDgAzbgMAYD3MuA0AABCAoAZJTU1NUlBQIHa7XdLS0mTZsmVy7ty5Xve5ePGirFy5Uq699lpJSUmRhQsXSmNjo0+Z48ePS15eniQnJ0t6erqsWbNGOjo6fMrs3r1bbrnlFklMTJRx48ZJSUmJz/YXXnhB/uIv/kLsdrvY7XZxOp2yY8eOAWk3AACwvqAGSQUFBXL48GEpLy+Xbdu2yd69e2X58uW97vPoo4/Ku+++K6WlpbJnzx45deqU3Hfffcb2zs5OycvLk7a2NqmqqpJXX31VSkpKpKioyChz9OhRycvLkzvuuEPq6upk1apV8sADD8j7779vlLnhhhtk06ZNUltbKwcOHJC/+qu/knvuuUcOHz488CcCAABYjwZJfX29ioju37/fWLdjxw6Ni4vTr7/+2nSf5uZmHTx4sJaWlhrrGhoaVES0urpaVVW3b9+uNptNXS6XUeaFF15Qu92uly5dUlXVtWvX6uTJk32OvWjRIs3Jyem1zsOGDdP/+I//8LuNLS0tKiLa0tLi9z4AACC8/L1+B60nqbq6WtLS0mTmzJnGuuzsbLHZbFJTU2O6T21trbS3t0t2draxbuLEiTJ69Giprq42jjt16lTJyMgwyuTk5Ijb7TZ6gaqrq32O0VWm6xiX6+zslK1bt8r58+fF6XT22KZLly6J2+32WQAAQHQKWpDkcrkkPT3dZ118fLwMHz5cXC5Xj/skJCRIWlqaz/qMjAxjH5fL5RMgdW3v2tZbGbfbLRcuXDDWffrpp5KSkiKJiYmyYsUKeeutt2TSpEk9tmnjxo2SmppqLJmZmb2cAQw4HswKAAihPgdJ69atk7i4uF6Xzz77LBh1HXATJkyQuro6qampkYceekgWL14s9fX1PZZfv369tLS0GMuJEydCWFsYj/CYM8f7ypxC0YuAGEAE6PNkkqtXr5YlS5b0Wmbs2LHicDjkzJkzPus7OjqkqalJHA6H6X4Oh0Pa2tqkubnZpzepsbHR2MfhcMi+fft89usa/da9zOUj4hobG8Vut0tSUpKxLiEhQcaNGyciIjNmzJD9+/fLv/7rv8pLL71kWr/ExERJTEzste0Iotxc32ec8WDW6BWOZ9oBwGX63JM0cuRImThxYq9LQkKCOJ1OaW5ultraWmPfiooK8Xg8kpWVZXrsGTNmyODBg2Xnzp3GuiNHjsjx48eNXCGn0ymffvqpTwBWXl4udrvduFXmdDp9jtFVprd8IxERj8cjly5d6tsJQeiYPcIjUtETEphwPNMOAC4XzOzx3NxcnT59utbU1OgHH3yg48eP1/z8fGP7yZMndcKECVpTU2OsW7FihY4ePVorKir0wIED6nQ61el0Gts7Ojp0ypQpOm/ePK2rq9OysjIdOXKkrl+/3ijz5ZdfanJysq5Zs0YbGhp08+bNOmjQIC0rKzPKrFu3Tvfs2aNHjx7VTz75RNetW6dxcXH63//93363j9FtIdbZqVpYqFpZ6X3t7Ax3jXpWWKgq8u1SWBjuGlkL5w9AEPl7/Q5qkHT27FnNz8/XlJQUtdvtunTpUm1tbTW2Hz16VEVEd+3aZay7cOGCPvzwwzps2DBNTk7We++9V0+fPu1z3GPHjun8+fM1KSlJR4wYoatXr9b29nafMrt27dJp06ZpQkKCjh07Vrds2eKz/ac//aneeOONmpCQoCNHjtQ777yzTwGSaoiDJCsFCPD+nrpf5Csrw10ja+HzDiCI/L1+8+y2AIT02W1FRb75GYWF5GhEMn5fABCxeHZbtCFHw1qKi72BUWWl97W3kXjkLwFAROrz6DaEiVnS8uzZ4akLrs5m+7bn6Gq/J0ZyAUBEIkiyiq6eiNxcb4DEHEHRg6kNACAiESRZRV96JmAt9BICQEQiSALCjV5CRCOPx/tZ7v65tpEGC2thdFsAQjq6DQCshBGeiGCMbgMAhA8jchEFCJIAAAPPSo8RAnpAThIAYOCRa4coQE5SAMhJAgDAeshJAgAACABBEgAAgAmCJCCceG4bAEQsEreBcOK5bQAQsehJAsKJuWQAIGIRJAHhxFwyABCxuN0GhBNzyQBAxGKepAAwTxIAANbDPEkAcDlGEwLoA263AYgdjCYE0Af0JAGIHYwmBNAHBEkAYgejCQH0AbfbAMQORhMC6ANGtwWA0W0AAFgPo9sAAAACQJAEAABggiAJAADABEESAACACYIkAAAAEwRJAAAAJgiSACDS8cw5ICwIkgAg0nU9c27OHO9rXybBJMAC+o0gCQAiXSDPnAskwAJiHEESAES6QJ45x0N9gX7j2W0AEOmKi0VURU6e9L5X9d42s/nxPdcswJo9e8CrCEQjgiQAiHQ2m0hcnEhJybfr4uJEnnji6vvyUF+g33jAbQB4wC2AkKmq8uYVdamspEcI6CcecAsA0SSQvCQA/cLtNgCwAm6bASHH7bYAcLsNAADr4XYbAABAAAiSAAAATBAkAQAAmCBIAgAAMEGQBAAAYIIgCQAAwARBEgAAgAmCJAAAABMESQAAACYIkgAAAEwQJAEAAJggSAIAADARH+4KWFnXs4HdbneYawIAAPzVdd3uuo73hCApAK2trSIikpmZGeaaAACAvmptbZXU1NQet8fp1cIo9Mjj8cipU6dk6NChEhcXF+7q9MrtdktmZqacOHFC7HZ7uKsTVLQ1OsVKW2OlnSK0NVpZoa2qKq2trTJq1Cix2XrOPKInKQA2m01uuOGGcFejT+x2e8R+aAcabY1OsdLWWGmnCG2NVpHe1t56kLqQuA0AAGCCIAkAAMAEQVKMSExMlA0bNkhiYmK4qxJ0tDU6xUpbY6WdIrQ1WkVTW0ncBgAAMEFPEgAAgAmCJAAAABMESQAAACYIkgAAAEwQJFlEU1OTFBQUiN1ul7S0NFm2bJmcO3eu130uXrwoK1eulGuvvVZSUlJk4cKF0tjY6FPm+PHjkpeXJ8nJyZKeni5r1qyRjo4OnzK7d++WW265RRITE2XcuHFSUlLis33jxo1y6623ytChQyU9PV0WLFggR44cicq27t27V+6++24ZNWqUxMXFydtvv92ntm3evFluuukmGTJkiGRlZcm+fft6LV9aWioTJ06UIUOGyNSpU2X79u0+21VVioqK5LrrrpOkpCTJzs6Wzz//3KeMP+fzk08+kdtuu02GDBkimZmZ8vTTT/epXVZo58WLF2XJkiUydepUiY+PlwULFgTUxkhu6+7du+Wee+6R6667Tq655hqZNm2a/O53v4vKth45ckTuuOMOycjIkCFDhsjYsWPl8ccfl/b29qhra3dffPGFDB06VNLS0gJqp0hktvXYsWMSFxd3xfLhhx8G3N4+UVhCbm6u3nzzzfrhhx/q//7v/+q4ceM0Pz+/131WrFihmZmZunPnTj1w4IB+73vf09mzZxvbOzo6dMqUKZqdna0HDx7U7du364gRI3T9+vVGmS+//FKTk5P1scce0/r6en3++ed10KBBWlZWZpTJycnRLVu26KFDh7Surk7vuusuHT16tJ47dy7q2rp9+3b9h3/4B/2v//ovFRF96623/G7X1q1bNSEhQV955RU9fPiwPvjgg5qWlqaNjY2m5SsrK3XQoEH69NNPa319vT7++OM6ePBg/fTTT40ymzZt0tTUVH377bf1448/1h/84Ac6ZswYvXDhgt/ns6WlRTMyMrSgoEAPHTqkr7/+uiYlJelLL73kd9us0M5z587pihUr9OWXX9acnBy95557+tU+K7T1qaee0scff1wrKyv1iy++0N/85jdqs9n03Xffjbq2/uEPf9BXXnlF6+rq9NixY/rOO+9oenq6z//taGlrl7a2Np05c6bOnz9fU1NT+93OSG7r0aNHVUT0f/7nf/T06dPG0tbWFlB7+4ogyQLq6+tVRHT//v3Guh07dmhcXJx+/fXXpvs0Nzfr4MGDtbS01FjX0NCgIqLV1dWq6r3g22w2dblcRpkXXnhB7Xa7Xrp0SVVV165dq5MnT/Y59qJFizQnJ6fH+p45c0ZFRPfs2RPVbe1rkDRr1ixduXKl8b6zs1NHjRqlGzduNC3/wx/+UPPy8nzWZWVl6c9+9jNVVfV4POpwOPSZZ54xtjc3N2tiYqK+/vrrqurf+fztb3+rw4YNM86DquovfvELnTBhgt9ts0I7u1u8ePGABElWaGuXu+66S5cuXdr3Rv6Zldr66KOP6ty5c/veyD+L9LauXbtWf/zjH+uWLVsCDpIita1dQdLBgwcDal+guN1mAdXV1ZKWliYzZ8401mVnZ4vNZpOamhrTfWpra6W9vV2ys7ONdRMnTpTRo0dLdXW1cdypU6dKRkaGUSYnJ0fcbrccPnzYKNP9GF1luo5hpqWlRUREhg8f3seWWq+t/mpra5Pa2lqf49tsNsnOzu7x+Ferz9GjR8XlcvmUSU1NlaysLJ92X+18VldXy/e//31JSEjw+TlHjhyRP/7xj1HTzoFmtba2tLT06/+kiLXa+sUXX0hZWZncfvvtUdnWiooKKS0tlc2bN/erfVZqq4jID37wA0lPT5e5c+fK73//+8Aa3A8ESRbgcrkkPT3dZ118fLwMHz5cXC5Xj/skJCRccb86IyPD2MflcvkEDV3bu7b1VsbtdsuFCxeu+Lkej0dWrVolc+bMkSlTpvjfyG71tkpb++Kbb76Rzs5O0+P31q7eyne9Xq3M1c6nP+fGX5HczoFmpba++eabsn//flm6dKmfrfNlhbbOnj1bhgwZIuPHj5fbbrtNnnjiiT620iuS23r27FlZsmSJlJSUDMiDYyO5rSkpKfKrX/1KSktL5b333pO5c+fKggULQh4oESSF0bp160wT07ovn332Wbir2ScrV66UQ4cOydatW33WR2NbASvYtWuXLF26VP793/9dJk+eHO7qBM0bb7whH330kbz22mvy3nvvybPPPhvuKg24Bx98UP76r/9avv/974e7KkE3YsQIeeyxxyQrK0tuvfVW2bRpk/z4xz+WZ555JqT1iA/pT4OP1atXy5IlS3otM3bsWHE4HHLmzBmf9R0dHdLU1CQOh8N0P4fDIW1tbdLc3OzTw9LY2Gjs43A4rhjF0DUirHuZy0eJNTY2it1ul6SkJJ/1jzzyiGzbtk327t0rN9xwQ1S3ta9GjBghgwYNMj1+b+3qrXzXa2Njo1x33XU+ZaZNm2aUudr57OnndP8Z/orkdg40K7R1z549cvfdd8u//Mu/yE9+8pO+N/LPrNDWzMxMERGZNGmSdHZ2yvLly2X16tUyaNCgqGlrRUWF/P73vzcCQFUVj8cj8fHx8vLLL8tPf/rTqGmrmaysLCkvL/evcQOEnqQwGjlypEycOLHXJSEhQZxOpzQ3N0ttba2xb0VFhXg8HsnKyjI99owZM2Tw4MGyc+dOY92RI0fk+PHj4nQ6RUTE6XTKp59+6vNhLS8vF7vdLpMmTTLKdD9GV5muY4h4/6M+8sgj8tZbb0lFRYWMGTMmatvaXwkJCTJjxgyf43s8Htm5c2ePx79afcaMGSMOh8OnjNvtlpqaGp92X+18Op1O2bt3r8+Q6fLycpkwYYIMGzYsato50CK9rbt375a8vDz553/+Z1m+fHlUt/VyHo9H2tvbxePxRFVbq6urpa6uzlieeOIJGTp0qNTV1cm9994bVW01U1dX5xN4hURY08bht9zcXJ0+fbrW1NToBx98oOPHj/cZLnny5EmdMGGC1tTUGOtWrFiho0eP1oqKCj1w4IA6nU51Op3G9q5h8fPmzdO6ujotKyvTkSNHmg6LX7NmjTY0NOjmzZuvGBb/0EMPaWpqqu7evdtnqOaf/vSnqGtra2urHjx4UA8ePKgior/+9a/14MGD+tVXX121XVu3btXExEQtKSnR+vp6Xb58uaalpRkj7u6//35dt26dUb6yslLj4+P12Wef1YaGBt2wYYPpUNu0tDR955139JNPPtF77rnHdKhtb+ezublZMzIy9P7779dDhw7p1q1bNTk5OaApACKxnaqqhw8f1oMHD+rdd9+tf/mXf2n8LvsrUttaUVGhycnJun79ep//k2fPno26tv7nf/6nvvHGG1pfX69/+MMf9I033tBRo0ZpQUFB1LX1cgMxui1S21pSUqKvvfaaNjQ0aENDgz711FNqs9n0lVdeCai9fUWQZBFnz57V/Px8TUlJUbvdrkuXLtXW1lZje9dwyV27dhnrLly4oA8//LAOGzZMk5OT9d5779XTp0/7HPfYsWM6f/58TUpK0hEjRujq1au1vb3dp8yuXbt02rRpmpCQoGPHjtUtW7b4bBcR0+XyctHQ1l27dpm2dfHixX617fnnn9fRo0drQkKCzpo1Sz/88ENj2+23337Fcd588039zne+owkJCTp58mR97733fLZ7PB4tLCzUjIwMTUxM1DvvvFOPHDnSp/Opqvrxxx/r3LlzNTExUa+//nrdtGmTX+2xWjtvvPFG099ftLV18eLFpu28/fbbo66tW7du1VtuuUVTUlL0mmuu0UmTJuk//dM/+VyQo6WtlxuIIClS21pSUqLf/e53NTk5We12u86aNctnmpdQiVNVDV2/FQAAgDWQkwQAAGCCIAkAAMAEQRIAAIAJgiQAAAATBEkAAAAmCJIAAABMECQBAACYIEgCAAAwQZAEAABggiAJAADABEESAACACYIkAAAAE/8f/5g+jFW2FlkAAAAASUVORK5CYII=", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "plt.scatter(dx[3, :], dy[3, :], marker='x', color='r', s=5)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "21fe7d36-3f21-419c-a609-3c387fa70dd4", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "bpd_gpu3", + "language": "python", + "name": "bpd_gpu3" + }, + "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.12.8" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +}