From f47a33945e90249c55b3cbc5da73d9b50c5854bf Mon Sep 17 00:00:00 2001 From: jrob93 Date: Tue, 28 May 2024 17:23:39 +0100 Subject: [PATCH 01/14] add some missing docs --- src/adler/science/PhaseCurve.py | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/src/adler/science/PhaseCurve.py b/src/adler/science/PhaseCurve.py index dc53468..715c0cc 100644 --- a/src/adler/science/PhaseCurve.py +++ b/src/adler/science/PhaseCurve.py @@ -49,6 +49,16 @@ def __init__(self, abs_mag=18, phase_param=0.2, phase_param2=None, model_name="H print("no model selected") def SetModelBounds(self, param, bound_vals=(None, None)): + """By default the sbpy model uses "physical" boundaries for the phase parameter(s). + Use this function to change the boundaries of a parameter when fitting, or remove them by setting to None. + + Parameters + ----------- + param : str + Adler PhaseCurve parameter to adjust boundaries of: abs_mag, phase_param, phase_param2 + bound_vals: tuple + The (lower, upper) boundaries to use when fitting - set (None, None) for no bounds + """ model_sbpy = self.model_function param_names = model_sbpy.param_names x = getattr(model_sbpy, param) From ec47db16d88201fa2aec763319d992b878f02b00 Mon Sep 17 00:00:00 2001 From: jrob93 Date: Tue, 28 May 2024 17:34:15 +0100 Subject: [PATCH 02/14] add an example nb --- docs/notebooks/adler_phasecurve_example.ipynb | 350 ++++++++++++++++++ docs/requirements.txt | 1 + 2 files changed, 351 insertions(+) create mode 100644 docs/notebooks/adler_phasecurve_example.ipynb diff --git a/docs/notebooks/adler_phasecurve_example.ipynb b/docs/notebooks/adler_phasecurve_example.ipynb new file mode 100644 index 0000000..9098084 --- /dev/null +++ b/docs/notebooks/adler_phasecurve_example.ipynb @@ -0,0 +1,350 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": null, + "id": "d591f5d8-9148-46ff-a62b-0f2a29eb806c", + "metadata": {}, + "outputs": [], + "source": [ + "from adler.dataclasses.AdlerPlanetoid import AdlerPlanetoid\n", + "from adler.science.PhaseCurve import PhaseCurve\n", + "\n", + "import numpy as np\n", + "import matplotlib.pyplot as plt\n", + "import matplotlib.gridspec as gridspec\n", + "import astropy.units as u" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "645efb98-567d-481e-a79c-b1cfdc828726", + "metadata": {}, + "outputs": [], + "source": [ + "# ssObjectId of object to analyse\n", + "ssoid = \"8268570668335894776\"" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "10b36aab-b322-49b8-8ff3-49bef68d7416", + "metadata": {}, + "outputs": [], + "source": [ + "# retrieve the object data via adler\n", + "\n", + "# # here we use an offline SQL database which contains the observations of the sso\n", + "fname = \"../../tests/data/testing_database.db\"\n", + "planetoid = AdlerPlanetoid.construct_from_SQL(ssoid, sql_filename=fname)\n", + "\n", + "# alternatively we can retrieve the object data directly from the RSP\n", + "# planetoid = AdlerPlanetoid.construct_from_RSP(ssoid)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d9a0623d-0dc7-49c1-99dd-a76ef970a3ff", + "metadata": {}, + "outputs": [], + "source": [ + "# inspect the whole object\n", + "planetoid.__dict__" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "1d360360-025b-4a77-acf5-325b2f2d1873", + "metadata": {}, + "outputs": [], + "source": [ + "# inspect just the ssObject table\n", + "planetoid.SSObject.__dict__" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "be8f8d63", + "metadata": {}, + "outputs": [], + "source": [ + "# retrieve all observations in the r filter\n", + "obs_r = planetoid.observations_in_filter(\"r\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "da1e483d", + "metadata": {}, + "outputs": [], + "source": [ + "# inspect the fields available in the observations table\n", + "obs_r.__dict__" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "9d7dc125-06c1-49ad-8854-17d8c8b6954f", + "metadata": {}, + "outputs": [], + "source": [ + "# plot the observations as a phasecurve\n", + "x_plot = \"phaseAngle\"\n", + "y_plot = \"reduced_mag\"\n", + "\n", + "x = getattr(obs_r, x_plot)\n", + "y = getattr(obs_r, y_plot)\n", + "xerr = obs_r.magErr\n", + "\n", + "fig = plt.figure()\n", + "gs = gridspec.GridSpec(1, 1)\n", + "ax1 = plt.subplot(gs[0, 0])\n", + "\n", + "ax1.errorbar(x, y, xerr, fmt=\"o\")\n", + "\n", + "ax1.invert_yaxis()\n", + "ax1.set_xlabel(x_plot)\n", + "ax1.set_ylabel(y_plot)\n", + "\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6efe3b5a-09dd-4d5e-9f41-20ea6e1b43df", + "metadata": {}, + "outputs": [], + "source": [ + "# retrieve the phase curve model parameters provided in the ssObject table\n", + "\n", + "sso_r = planetoid.SSObject_in_filter(\"r\")\n", + "\n", + "r_H = sso_r.H\n", + "r_G12 = sso_r.G12\n", + "\n", + "pc = PhaseCurve(abs_mag=r_H * u.mag, phase_param=r_G12, model_name=\"HG12_Pen16\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "80f552f1-8907-4cc9-b57c-2e667eab459c", + "metadata": {}, + "outputs": [], + "source": [ + "# what sbpy model is being used?\n", + "pc.model_function" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "24c1955e-95cd-4d77-ad05-aa5b8d18620a", + "metadata": {}, + "outputs": [], + "source": [ + "# set up an array of phase angles to plot the model\n", + "alpha = np.linspace(0, np.amax(obs_r.phaseAngle)) * u.deg\n", + "alpha" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c3f30fe0-0d89-4ffa-8237-9c71181d44ee", + "metadata": {}, + "outputs": [], + "source": [ + "# calculate the model reduced magnitude over these phase angles\n", + "red_mag = pc.ReducedMag(alpha)\n", + "red_mag" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "04be98a1-e4dc-4216-bcd9-ef777f6053fb", + "metadata": {}, + "outputs": [], + "source": [ + "# plot the observations with the model phase curve\n", + "x_plot = \"phaseAngle\"\n", + "y_plot = \"reduced_mag\"\n", + "\n", + "x = getattr(obs_r, x_plot)\n", + "y = getattr(obs_r, y_plot)\n", + "xerr = obs_r.magErr\n", + "\n", + "fig = plt.figure()\n", + "gs = gridspec.GridSpec(1, 1)\n", + "ax1 = plt.subplot(gs[0, 0])\n", + "\n", + "ax1.errorbar(x, y, xerr, fmt=\"o\")\n", + "\n", + "ax1.plot(alpha.value, red_mag.value)\n", + "\n", + "ax1.invert_yaxis()\n", + "ax1.set_xlabel(x_plot)\n", + "ax1.set_ylabel(y_plot)\n", + "\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "9815543d-6140-4bdb-8bad-8296994723f4", + "metadata": {}, + "outputs": [], + "source": [ + "# plot the observations as a lightcurve\n", + "x_plot = \"midPointMjdTai\"\n", + "y_plot = \"reduced_mag\"\n", + "\n", + "x = getattr(obs_r, x_plot)\n", + "y = getattr(obs_r, y_plot)\n", + "xerr = obs_r.magErr\n", + "\n", + "fig = plt.figure()\n", + "gs = gridspec.GridSpec(1, 1)\n", + "ax1 = plt.subplot(gs[0, 0])\n", + "\n", + "ax1.errorbar(x, y, xerr, fmt=\"o\")\n", + "\n", + "ax1.invert_yaxis()\n", + "ax1.set_xlabel(x_plot)\n", + "ax1.set_ylabel(y_plot)\n", + "\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "de462b92-3914-4091-b0af-bddd9e9c1ef1", + "metadata": {}, + "outputs": [], + "source": [ + "# do a different phase curve fit to the data\n", + "# adler is able to fit different models, and perform more sophisticated fits" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f92891c9-6ccf-4dac-8887-9545f633ba90", + "metadata": {}, + "outputs": [], + "source": [ + "# create a new PhaseCurve object with a different sbpy model\n", + "pc_fit = PhaseCurve(abs_mag=pc.abs_mag, model_name=\"HG\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "db24432b-6d05-4ff2-9d98-e52d8c2e4342", + "metadata": {}, + "outputs": [], + "source": [ + "pc_fit.model_function" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "9039e2e2-27d9-4d21-b2f6-9504a5b85ce4", + "metadata": {}, + "outputs": [], + "source": [ + "# use adler to fit this new phase curve model to the data\n", + "pc_fit.FitModel(\n", + " phase_angle=obs_r.phaseAngle * u.deg,\n", + " reduced_mag=obs_r.reduced_mag * u.mag,\n", + " mag_err=obs_r.magErr * u.mag,\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5d4b7144-ee72-45e0-9606-c40f83c443c6", + "metadata": {}, + "outputs": [], + "source": [ + "# plot the observations with both\n", + "x_plot = \"phaseAngle\"\n", + "y_plot = \"reduced_mag\"\n", + "\n", + "x = getattr(obs_r, x_plot)\n", + "y = getattr(obs_r, y_plot)\n", + "xerr = obs_r.magErr\n", + "\n", + "fig = plt.figure()\n", + "gs = gridspec.GridSpec(1, 1)\n", + "ax1 = plt.subplot(gs[0, 0])\n", + "\n", + "ax1.errorbar(x, y, xerr, fmt=\"o\")\n", + "\n", + "ax1.plot(alpha.value, pc.ReducedMag(alpha).value, label=pc.model_name)\n", + "ax1.plot(alpha.value, pc_fit.ReducedMag(alpha).value, label=pc_fit.model_name)\n", + "\n", + "ax1.invert_yaxis()\n", + "ax1.set_xlabel(x_plot)\n", + "ax1.set_ylabel(y_plot)\n", + "ax1.legend()\n", + "\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b7f39ed4-8334-4e10-a97c-a9471105225b", + "metadata": {}, + "outputs": [], + "source": [ + "# # now we would add our calculated values back into planetoid\n", + "# planetoid.AdlerSchema.r_H = pc_fit.abs_mag\n", + "# planetoid.AdlerSchema.r_G = pc_fit.phase_param" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "1ca4bbfd-1954-469f-8608-40c52838d300", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "adler-dev", + "language": "python", + "name": "adler-dev" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.10.14" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/docs/requirements.txt b/docs/requirements.txt index 3979f83..5f633ff 100644 --- a/docs/requirements.txt +++ b/docs/requirements.txt @@ -8,3 +8,4 @@ jupyter astropy sbpy matplotlib +numpy \ No newline at end of file From a17addcc9436ad6154727b39d539d931229c785b Mon Sep 17 00:00:00 2001 From: jrob93 Date: Tue, 6 Aug 2024 12:35:29 +0100 Subject: [PATCH 03/14] update docs int -> float --- src/adler/dataclasses/AdlerPlanetoid.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/adler/dataclasses/AdlerPlanetoid.py b/src/adler/dataclasses/AdlerPlanetoid.py index d6ceb40..b9d31b8 100644 --- a/src/adler/dataclasses/AdlerPlanetoid.py +++ b/src/adler/dataclasses/AdlerPlanetoid.py @@ -29,7 +29,7 @@ def __init__( filter_list : list of str A comma-separated list of the filters of interest. - date_range : list of int + date_range : list of float The minimum and maximum dates of the desired observations. observations_by_filter : list of Observations objects @@ -76,7 +76,7 @@ def construct_from_SQL( filter_list : list of str A comma-separated list of the filters of interest. - date_range : list of int + date_range : list of float The minimum and maximum dates of the desired observations. schema : str or None @@ -157,7 +157,7 @@ def construct_from_RSP( filter_list : list of str A comma-separated list of the filters of interest. - date_range : list of int + date_range : list of float The minimum and maximum dates of the desired observations. """ @@ -215,7 +215,7 @@ def populate_observations( filter_list : list of str A comma-separated list of the filters of interest. - date_range : list of int + date_range : list of float The minimum and maximum dates of the desired observations. service : pyvo.dal.tap.TAPService object or None From 9f1077c896b9e283a63f4f7c30542c47e1127f8f Mon Sep 17 00:00:00 2001 From: James Robinson Date: Tue, 6 Aug 2024 15:02:42 +0100 Subject: [PATCH 04/14] 142 colour functions (#154) * add x,y,z pos of observation * include topocentric coords * add tperi to mpcorb query * add ecliptic coordinates to query * include MPCORB fullDesignation * update fullDesignation unit test * fix fit covariance matrix * add more tests * work with quantity objects * no units test * update PhaseCurve parameter names * fix InitModelDict * initial apparition finder * get obs as df func * initial commit for colour code * initial colour function * update colour func * add colour demo nb * colour test data * add unit colour unit tests * [pre-commit.ci lite] apply automatic fixes * demo nb * add description to doc string --------- Co-authored-by: pre-commit-ci-lite[bot] <117423508+pre-commit-ci-lite[bot]@users.noreply.github.com> --- notebooks/colour_functions_demo.ipynb | 391 ++++++++++++++++++ .../adler_demo_testing_database.db | Bin 512000 -> 512000 bytes src/adler/dataclasses/AdlerData.py | 2 +- src/adler/dataclasses/AdlerPlanetoid.py | 6 +- src/adler/dataclasses/MPCORB.py | 10 + src/adler/dataclasses/Observations.py | 40 ++ src/adler/science/Colour.py | 126 ++++++ src/adler/science/PhaseCurve.py | 69 +++- src/adler/utilities/science_utilities.py | 115 ++++++ tests/adler/dataclasses/test_MPCORB.py | 5 +- tests/adler/dataclasses/test_Observations.py | 4 +- tests/adler/science/test_Colour.py | 112 +++++ tests/adler/science/test_PhaseCurve.py | 3 + .../adler/utilities/test_science_utilities.py | 8 + .../df_6098332225018_g_r_app_3_N_ref_1.csv | 16 + .../df_6098332225018_g_r_app_3_N_ref_3.csv | 16 + .../df_6098332225018_g_r_app_3_N_ref_5.csv | 16 + tests/data/testing_database.db | Bin 57344 -> 512000 bytes 18 files changed, 931 insertions(+), 8 deletions(-) create mode 100644 notebooks/colour_functions_demo.ipynb create mode 100644 src/adler/science/Colour.py create mode 100644 tests/adler/science/test_Colour.py create mode 100644 tests/data/df_6098332225018_g_r_app_3_N_ref_1.csv create mode 100644 tests/data/df_6098332225018_g_r_app_3_N_ref_3.csv create mode 100644 tests/data/df_6098332225018_g_r_app_3_N_ref_5.csv diff --git a/notebooks/colour_functions_demo.ipynb b/notebooks/colour_functions_demo.ipynb new file mode 100644 index 0000000..902f481 --- /dev/null +++ b/notebooks/colour_functions_demo.ipynb @@ -0,0 +1,391 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": null, + "id": "d1ef57c2", + "metadata": {}, + "outputs": [], + "source": [ + "%matplotlib notebook" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "61f50e0d", + "metadata": {}, + "outputs": [], + "source": [ + "%matplotlib inline" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "82537a61", + "metadata": {}, + "outputs": [], + "source": [ + "from adler.dataclasses.AdlerPlanetoid import AdlerPlanetoid\n", + "from adler.dataclasses.AdlerData import AdlerData\n", + "from adler.science.PhaseCurve import PhaseCurve\n", + "from adler.science.Colour import col_obs_ref\n", + "from adler.utilities.plotting_utilities import plot_errorbar\n", + "from adler.utilities.science_utilities import apparition_gap_finder, get_df_obs_filt\n", + "\n", + "import numpy as np\n", + "import pandas as pd\n", + "import matplotlib.pyplot as plt\n", + "import matplotlib.gridspec as gridspec\n", + "import astropy.units as u\n", + "from scipy.interpolate import splrep, BSpline, CubicSpline\n", + "from astropy.modeling.fitting import SLSQPLSQFitter\n", + "from astropy.time import Time\n", + "from astroquery.jplhorizons import Horizons" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b9f6dcda", + "metadata": {}, + "outputs": [], + "source": [ + "# ssObjectId of object to analyse\n", + "ssoid = \"6098332225018\" # good MBA test object" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c6f705c1", + "metadata": {}, + "outputs": [], + "source": [ + "# retrieve the object data\n", + "# fname = \"../tests/data/testing_database.db\"\n", + "fname = \"/Users/jrobinson/lsst-adler/notebooks/gen_test_data/adler_demo_testing_database.db\"\n", + "planetoid = AdlerPlanetoid.construct_from_SQL(ssoid, sql_filename=fname)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "21e0b93e", + "metadata": {}, + "outputs": [], + "source": [ + "# check orbit parameters\n", + "e = planetoid.MPCORB.e\n", + "incl = planetoid.MPCORB.incl\n", + "q = planetoid.MPCORB.q\n", + "a = q / (1.0 - e)\n", + "Q = a * (1.0 + e)\n", + "print(a, e, incl)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ee4064e7", + "metadata": {}, + "outputs": [], + "source": [ + "# plot and fit the phase curves in g and r filters\n", + "\n", + "fig = plt.figure()\n", + "gs = gridspec.GridSpec(1, 1)\n", + "ax1 = plt.subplot(gs[0, 0])\n", + "\n", + "adler_cols = AdlerData(ssoid, planetoid.filter_list)\n", + "\n", + "for filt in [\"g\", \"r\"]:\n", + " # get observations and define a phase curve model to fit\n", + "\n", + " sso = planetoid.SSObject_in_filter(filt)\n", + " obs = planetoid.observations_in_filter(filt)\n", + "\n", + " H = sso.H\n", + " G12 = sso.G12\n", + " pc = PhaseCurve(H=H * u.mag, phase_parameter_1=G12, model_name=\"HG12_Pen16\")\n", + "\n", + " alpha = np.linspace(0, np.amax(obs.phaseAngle)) * u.deg\n", + "\n", + " pc_fit = pc.FitModel(\n", + " np.array(getattr(obs, \"phaseAngle\")) * u.deg,\n", + " np.array(getattr(obs, \"reduced_mag\")) * u.mag,\n", + " )\n", + " pc = pc.InitModelSbpy(pc_fit)\n", + " red_mag = pc.ReducedMag(alpha)\n", + "\n", + " adler_cols.populate_phase_parameters(filt, **pc.ReturnModelDict())\n", + "\n", + " # add this phase curve to the figure\n", + " fig = plot_errorbar(planetoid, filt_list=[filt], fig=fig)\n", + " ax1 = fig.axes[0]\n", + " ax1.plot(\n", + " alpha.value,\n", + " pc.ReducedMag(alpha).value,\n", + " label=\"{}: H={:.2f}, G12={:.2f}\".format(filt, pc.H, pc.phase_parameter_1),\n", + " )\n", + " ax1.legend()\n", + "\n", + "ax1.invert_yaxis()\n", + "\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d0df3c51", + "metadata": {}, + "outputs": [], + "source": [ + "# inspect the r filter phase curve model\n", + "adler_cols.get_phase_parameters_in_filter(\"r\", \"HG12_Pen16\").__dict__" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c46fb9eb", + "metadata": {}, + "outputs": [], + "source": [ + "# inspect the g filter phase curve model\n", + "adler_cols.get_phase_parameters_in_filter(\"g\", \"HG12_Pen16\").__dict__" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "73871e56", + "metadata": {}, + "outputs": [], + "source": [ + "# get the boundary times for each apparation of the object in the survey\n", + "df_obs_all = pd.DataFrame()\n", + "for filt in [\"r\", \"g\"]:\n", + " obs = planetoid.observations_in_filter(filt)\n", + " _df_obs = pd.DataFrame(obs.__dict__)\n", + " df_obs_all = pd.concat([df_obs_all, _df_obs])\n", + "df_obs_all = df_obs_all.sort_values(\"midPointMjdTai\")\n", + "\n", + "t_app = apparition_gap_finder(np.array(df_obs_all[\"midPointMjdTai\"]))\n", + "print(t_app)" + ] + }, + { + "cell_type": "markdown", + "id": "63980010", + "metadata": {}, + "source": [ + "Now we can inpsect how the colour of the object varies (or not) as a function of time. The adler function `col_obs_ref` will compare the latest observation in a given filter with observations in another filter. By setting parameter `N_ref` one can set how many past obsevrations to use when calculating the latest colour." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "bd0e0185", + "metadata": {}, + "outputs": [], + "source": [ + "# Here we simulate observations coming night-by-night and the calculation of a g-r colour for the object\n", + "\n", + "# define colour function parameters\n", + "\n", + "# set number of reference observations to use for colour estimate\n", + "# N_ref = 5\n", + "N_ref = 3\n", + "# N_ref = 1\n", + "# N_ref = None\n", + "\n", + "# observation and filter field names\n", + "x_plot = \"midPointMjdTai\"\n", + "y_plot = \"AbsMag\"\n", + "y_plot = \"reduced_mag\"\n", + "yerr_plot = \"magErr\"\n", + "filt_obs = \"g\" # observation filter\n", + "filt_ref = \"r\" # reference filter (we are calculating a filt_obs - filt_ref colour)\n", + "\n", + "# define colour field names\n", + "colour = \"{}-{}\".format(filt_obs, filt_ref)\n", + "colErr = \"{}-{}Err\".format(filt_obs, filt_ref)\n", + "delta_t_col = \"delta_t_{}\".format(colour)\n", + "y_ref_col = \"{}_{}\".format(y_plot, filt_ref)\n", + "x1_ref_col = \"{}1_{}\".format(x_plot, filt_ref)\n", + "x2_ref_col = \"{}2_{}\".format(x_plot, filt_ref)\n", + "\n", + "fig = plt.figure()\n", + "gs = gridspec.GridSpec(1, 1)\n", + "ax1 = plt.subplot(gs[0, 0])\n", + "\n", + "col_dict_list = []\n", + "for app_i in range(len(t_app) - 1):\n", + " # consider only one apparition\n", + " if app_i != 3:\n", + " continue\n", + "\n", + " time_min = t_app[app_i]\n", + " time_max = t_app[app_i + 1]\n", + " print(time_min, time_max)\n", + "\n", + " _df_obs_all = df_obs_all[\n", + " (df_obs_all[\"midPointMjdTai\"] >= time_min) & (df_obs_all[\"midPointMjdTai\"] < time_max)\n", + " ]\n", + " print(app_i, len(_df_obs_all))\n", + " _time_max = np.amax(_df_obs_all[\"midPointMjdTai\"])\n", + "\n", + " # get the phase curve model and observations for each filter\n", + "\n", + " # get the stored AdlerData parameters for the observation filter\n", + " ad_g = adler_cols.get_phase_parameters_in_filter(filt_obs, \"HG12_Pen16\")\n", + " pc_g = PhaseCurve().InitModelDict(ad_g.__dict__) # make the PhaseCurve object from AdlerData\n", + " # get the phase curve model for the reference filter\n", + " ad_r = adler_cols.get_phase_parameters_in_filter(filt_ref, \"HG12_Pen16\")\n", + " pc_r = PhaseCurve().InitModelDict(ad_r.__dict__)\n", + " # get the observations in both filters\n", + " df_obs = get_df_obs_filt(\n", + " planetoid, filt_obs, x1=time_min, x2=_time_max, col_list=[y_plot, yerr_plot], pc_model=pc_g\n", + " )\n", + " df_obs_ref = get_df_obs_filt(\n", + " planetoid, filt_ref, x1=time_min, x2=_time_max, col_list=[y_plot, yerr_plot], pc_model=pc_r\n", + " )\n", + "\n", + " ax1.errorbar(df_obs[x_plot], df_obs[y_plot], df_obs[yerr_plot], fmt=\"o\", label=filt_obs)\n", + " ax1.errorbar(df_obs_ref[x_plot], df_obs_ref[y_plot], df_obs_ref[yerr_plot], fmt=\"o\", label=filt_ref)\n", + "\n", + " # simulate stepping through each filt_obs observation\n", + " x1 = time_min\n", + " for xi in range(len(df_obs)):\n", + " x2 = df_obs.iloc[xi][x_plot]\n", + " print(xi, x1, x2)\n", + "\n", + " # do the colour finding function here\n", + " col_dict = col_obs_ref(\n", + " planetoid,\n", + " adler_cols,\n", + " filt_obs=filt_obs,\n", + " filt_ref=filt_ref,\n", + " N_ref=N_ref,\n", + " x_col=x_plot,\n", + " y_col=y_plot,\n", + " yerr_col=yerr_plot,\n", + " x1=x1,\n", + " x2=x2,\n", + " )\n", + " col_dict_list.append(col_dict)\n", + "\n", + " # plot some lines to show the colour and mean reference\n", + " ax1.vlines(df_obs.iloc[xi][x_plot], df_obs.iloc[xi][y_plot], col_dict[y_ref_col], color=\"k\", ls=\":\")\n", + " ax1.hlines(col_dict[y_ref_col], col_dict[x1_ref_col], col_dict[x2_ref_col], color=\"k\", ls=\"--\")\n", + "\n", + "# store running colour parameters as a dataframe\n", + "df_col = pd.DataFrame(col_dict_list)\n", + "df_col = pd.concat([df_obs, df_col], axis=1)\n", + "\n", + "ax1.set_xlabel(x_plot)\n", + "ax1.set_ylabel(y_plot)\n", + "ax1.legend()\n", + "ax1.invert_yaxis()\n", + "\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a7734122", + "metadata": {}, + "outputs": [], + "source": [ + "# recorded colour parameters\n", + "df_col" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "0cd867d9", + "metadata": {}, + "outputs": [], + "source": [ + "# Plot how the colour changes as a function of time\n", + "# find filt_obs - filt_ref of newest filt_obs observation to the mean of the previous N_ref filt_ref observations\n", + "# colour code by time diff between obs and most recent obs_ref\n", + "\n", + "x_plot = \"midPointMjdTai\"\n", + "y_plot = colour\n", + "y_plot_err = colErr\n", + "c_plot = delta_t_col\n", + "df_plot = df_col\n", + "\n", + "fig = plt.figure()\n", + "gs = gridspec.GridSpec(1, 1)\n", + "ax1 = plt.subplot(gs[0, 0])\n", + "\n", + "\n", + "s1 = ax1.scatter(df_plot[x_plot], df_plot[y_plot], c=df_plot[c_plot], zorder=3)\n", + "cbar1 = plt.colorbar(s1)\n", + "ax1.errorbar(df_plot[x_plot], df_plot[y_plot], df_plot[yerr_plot], fmt=\".\", zorder=1)\n", + "\n", + "# TODO: leave the most recent obs out of this mean?\n", + "# loop through new obs and calculate current mean and std?\n", + "obs_ref_mean = np.mean(df_plot[y_plot])\n", + "obs_ref_median = np.median(df_plot[y_plot])\n", + "obs_ref_std = np.std(df_plot[y_plot])\n", + "print(obs_ref_median, obs_ref_std)\n", + "\n", + "ax1.axhline(obs_ref_mean, c=\"k\")\n", + "ax1.axhspan(obs_ref_mean - obs_ref_std, obs_ref_mean + obs_ref_std, zorder=0, color=\"k\", alpha=0.2)\n", + "\n", + "ax1.set_xlabel(x_plot)\n", + "ax1.set_ylabel(y_plot)\n", + "cbar1.set_label(c_plot)\n", + "\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "05abfee3", + "metadata": {}, + "outputs": [], + "source": [ + "# These colours can then be run through the previously written outlier detection functions" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "605a043a", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "adler-dev", + "language": "python", + "name": "adler-dev" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.10.14" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/notebooks/gen_test_data/adler_demo_testing_database.db b/notebooks/gen_test_data/adler_demo_testing_database.db index 1c3cdb5eca3624872b5c8052c309a211a5f99258..455b28b12e546bf0c177958a4e41c6884dac604b 100644 GIT binary patch delta 818 zcmW;I&rj1}7zgmZ?c3JGxxTPO`MvT(4+ibR3_muMqA}5E!c-6)17@O$2QXDp5Ye=~ zwA`jk^g0ft!3)zydfmx`C;kND)eByih53DZ_QmpC9K7sT10#O=xLbNJ?}_Yq$6 zp1~>aE;PIhoc6llb&tRs?sGWe&cHEu6kc$9;fzb+R%I-;zBJ^lqh4{A;IeZQ=AB`9 z&q>1j_9lE_KY_FMEjVZE@S%Mc&htw=dI}!j0q;kE(ka1g;v0-5gSfVnnA_auH+|== z!40nj*Q1HwdgG|yctdc_WAL@R1z$yTzjS9&pLY$|;qsX#+zz-HO`oj@`p*1x~#pu=z2X`8UiF2e3;`X2iV>RwxeeRdqCHh{+O(4PkblVCh(=xqimRN|&C z74RZivCND#foid0n7K+|u1qiBX%*qk+mS(}6XCzV8|gv#&+bD~h=%kd0|@UVhD0OP z#1D&>eR@n;NNbYH!epBKoAPN%WIkn5b(cyg>KrCwk0CVGkd%bh{@0dHlu%~MiRx2L zG){90%|{?bFu8xjjgtyojsvqA?^}fB6(s{f@%*Vk`h!B^fRODV&VDBzZB;Mh&Ik7s+;pGA z4R;#8cZcCSw+pVj1ip1%z&FkeeC>?Fey0a^I24Z8$5I<>PwWlU7wt9JYd?UU_8{!C z+hMo01+Q37VUP6?UbWINWnF?QAL80GaQ!j3Jp`0a31$=DVXPFy)xAU^Tji7H+;uqN z>hMmhPfc`YA%CQ&JMiS>ifR4fO^Q$;F!aE8gM$`Oso5yx{%&m>a+QoimkzN ztM9zc&&=9=mBP&P|Sw#jAF(kU+C4QQ; z7Se{?O!|0cUzjL|a)z8uhV+;=_YaOG%6#fcETnm%nfzvosC1#ZzniC@O`O$)v!qW| zQS+%mu-8xo7RsK&S8ECuOl6dXSwkVr&*=)tb2Mz_GRnWCpH)$28TNnW@Em3{^l;y& zk%aFuT%xJ~OyY%9B}o>xXN13a_YvQ3a89Hsq>;=&YltIw*O7df_Ngvnk8xCrHW(Gy fR~f^|Cs>n@IY*ob6Uhyc_v`V~H9ldaG?V`UX>+!M diff --git a/src/adler/dataclasses/AdlerData.py b/src/adler/dataclasses/AdlerData.py index 86c2e2a..e3231f6 100644 --- a/src/adler/dataclasses/AdlerData.py +++ b/src/adler/dataclasses/AdlerData.py @@ -101,7 +101,7 @@ def populate_phase_parameters(self, filter_name, **kwargs): # update the value if it's in **kwargs for model_key in MODEL_DEPENDENT_KEYS: - if kwargs.get(model_key): + if model_key in kwargs: setattr( self.filter_dependent_values[filter_index].model_dependent_values[model_index], model_key, diff --git a/src/adler/dataclasses/AdlerPlanetoid.py b/src/adler/dataclasses/AdlerPlanetoid.py index b9d31b8..4592e18 100644 --- a/src/adler/dataclasses/AdlerPlanetoid.py +++ b/src/adler/dataclasses/AdlerPlanetoid.py @@ -240,7 +240,9 @@ def populate_observations( observations_sql_query = f""" SELECT ssObject.ssObjectId, mag, magErr, band, midPointMjdTai, ra, dec, phaseAngle, - topocentricDist, heliocentricDist + topocentricDist, heliocentricDist, heliocentricX, heliocentricY, heliocentricZ, + topocentricX, topocentricY, topocentricZ, + eclipticLambda, eclipticBeta FROM {schema}ssObject JOIN {schema}diaSource ON {schema}ssObject.ssObjectId = {schema}diaSource.ssObjectId @@ -291,7 +293,7 @@ def populate_MPCORB(self, ssObjectId, service=None, sql_filename=None, schema="d MPCORB_sql_query = f""" SELECT - ssObjectId, mpcDesignation, mpcNumber, mpcH, mpcG, epoch, peri, node, incl, e, n, q, + ssObjectId, mpcDesignation, fullDesignation, mpcNumber, mpcH, mpcG, epoch, tperi, peri, node, incl, e, n, q, uncertaintyParameter, flags FROM {schema}MPCORB diff --git a/src/adler/dataclasses/MPCORB.py b/src/adler/dataclasses/MPCORB.py index 289084d..844cef9 100644 --- a/src/adler/dataclasses/MPCORB.py +++ b/src/adler/dataclasses/MPCORB.py @@ -4,10 +4,12 @@ MPCORB_KEYS = { "mpcDesignation": str, + "fullDesignation": str, "mpcNumber": int, "mpcH": float, "mpcG": float, "epoch": float, + "tperi": float, "peri": float, "node": float, "incl": float, @@ -31,6 +33,9 @@ class MPCORB: mpcDesignation: str Number or provisional designation (in packed form) + fullDesignation: str + Number or provisional designation (in readable form) + mpcNumber: int MPC number (if the asteroid has been numbered; NULL otherwise). Provided for convenience. @@ -43,6 +48,9 @@ class MPCORB: epoch: float Epoch (in MJD, .0 TT) + tperi: float + MJD of pericentric passage + peri: float Argument of perihelion, J2000.0 (degrees) @@ -71,10 +79,12 @@ class MPCORB: ssObjectId: str = "" mpcDesignation: str = "" + fullDesignation: str = "" mpcNumber: int = 0 mpcH: float = 0.0 mpcG: float = 0.0 epoch: float = 0.0 + tperi: float = 0.0 peri: float = 0.0 node: float = 0.0 incl: float = 0.0 diff --git a/src/adler/dataclasses/Observations.py b/src/adler/dataclasses/Observations.py index e99a758..eb92af2 100644 --- a/src/adler/dataclasses/Observations.py +++ b/src/adler/dataclasses/Observations.py @@ -12,6 +12,14 @@ "phaseAngle": np.ndarray, "topocentricDist": np.ndarray, "heliocentricDist": np.ndarray, + "heliocentricX": np.ndarray, + "heliocentricY": np.ndarray, + "heliocentricZ": np.ndarray, + "topocentricX": np.ndarray, + "topocentricY": np.ndarray, + "topocentricZ": np.ndarray, + "eclipticLambda": np.ndarray, + "eclipticBeta": np.ndarray, } @@ -53,6 +61,30 @@ class Observations: heliocentricDist: array_like of floats Heliocentric distance. + heliocentricX: array_like of floats + x-axis component of the heliocentric distance. + + heliocentricY: array_like of floats + y-axis component of the heliocentric distance. + + heliocentricZ: array_like of floats + z-axis component of the heliocentric distance. + + topocentricX: array_like of floats + x-axis component of the topocentric distance. + + topocentricY: array_like of floats + y-axis component of the topocentric distance. + + topocentricZ: array_like of floats + z-axis component of the topocentric distance. + + eclipticLambda: array_like of floats + The ecliptic longitude. + + eclipticBeta: array_like of floats + The ecliptic latitude. + reduced_mag: array_like of floats The reduced magnitude. @@ -71,6 +103,14 @@ class Observations: phaseAngle: np.ndarray = field(default_factory=lambda: np.zeros(0)) topocentricDist: np.ndarray = field(default_factory=lambda: np.zeros(0)) heliocentricDist: np.ndarray = field(default_factory=lambda: np.zeros(0)) + heliocentricX: np.ndarray = field(default_factory=lambda: np.zeros(0)) + heliocentricY: np.ndarray = field(default_factory=lambda: np.zeros(0)) + heliocentricZ: np.ndarray = field(default_factory=lambda: np.zeros(0)) + topocentricX: np.ndarray = field(default_factory=lambda: np.zeros(0)) + topocentricY: np.ndarray = field(default_factory=lambda: np.zeros(0)) + topocentricZ: np.ndarray = field(default_factory=lambda: np.zeros(0)) + eclipticLambda: np.ndarray = field(default_factory=lambda: np.zeros(0)) + eclipticBeta: np.ndarray = field(default_factory=lambda: np.zeros(0)) reduced_mag: np.ndarray = field(default_factory=lambda: np.zeros(0)) num_obs: int = 0 diff --git a/src/adler/science/Colour.py b/src/adler/science/Colour.py new file mode 100644 index 0000000..0144684 --- /dev/null +++ b/src/adler/science/Colour.py @@ -0,0 +1,126 @@ +import numpy as np +from adler.science.PhaseCurve import PhaseCurve +from adler.utilities.science_utilities import get_df_obs_filt + +# TODO: add colour estimate function + + +def col_obs_ref( + planetoid, + adler_cols, + filt_obs="g", + filt_ref="r", + N_ref=3, + x_col="midPointMjdTai", + y_col="AbsMag", + yerr_col="magErr", + x1=None, + x2=None, +): + """A function to calculate the colour of an Adler planetoid object. + An observation in a given filter (filt_obs) is compared to previous observation in a reference filter (filt_ref). + By setting N_ref one can control how many reference observations to include in the colour calculation. + Note that the observations are considered to be in reduced magnitudes, hence why an adlerData object with phase curve models is passed. + + Parameters + ----------- + planetoid : object + Adler planetoid object + + adler_cols : object + AdlerData object + + filt_obs: str + Filter name of the new observation for which to calculate a filt_obs - filt_ref colour + + filt_ref:str + Filter name of the reference observations when calculating the filt_obs - filt_ref colour + + N_ref: int + Number of reference observations to use when calculated the reference absolute magnitude. + Set to 1 to use only the most recent of ref observations. + Set to None to use all past ref obs. + + x_col: str + Time (or phase angle) column name + + y_col:str + Magnitude column name (e.g. AbsMag or reduced_mag) + + yerr_col: str + Magnitude uncertainty column name + + x1: float + Lower limit on x_col values (>=) + + x2: float + Upper limit on x_col values (<=) + + Returns + ---------- + + col_dict : dict + Dictionary containing the colour information for the most recent obs in filt_obs + + """ + # define fields to be recorded as a colour dictionary + colour = "{}-{}".format(filt_obs, filt_ref) + colErr = "{}-{}Err".format(filt_obs, filt_ref) + delta_t_col = "delta_t_{}".format(colour) + y_ref_col = "{}_{}".format(y_col, filt_ref) + x1_ref_col = "{}1_{}".format(x_col, filt_ref) + x2_ref_col = "{}2_{}".format(x_col, filt_ref) + + # get the stored AdlerData parameters for this filter # TODO: do this bit better + ad_obs = adler_cols.get_phase_parameters_in_filter(filt_obs, "HG12_Pen16") + # make the PhaseCurve object from AdlerData + pc_obs = PhaseCurve().InitModelDict(ad_obs.__dict__) + ad_ref = adler_cols.get_phase_parameters_in_filter(filt_ref, "HG12_Pen16") + pc_ref = PhaseCurve().InitModelDict(ad_ref.__dict__) + df_obs = get_df_obs_filt(planetoid, filt_obs, x1=x1, x2=x2, col_list=[y_col, yerr_col], pc_model=pc_obs) + df_obs_ref = get_df_obs_filt( + planetoid, filt_ref, x1=x1, x2=x2, col_list=[y_col, yerr_col], pc_model=pc_ref + ) + + # select the values of the new observation in the selected filter + i = -1 + x_obs = df_obs.iloc[i][x_col] + y_obs = df_obs.iloc[i][y_col] + yerr_obs = df_obs.iloc[i][yerr_col] + + # select observations in the reference filter from before the new obs + ref_mask = df_obs_ref[x_col] < x_obs + + # set the number of ref obs to use + if N_ref is None: + _N_ref = len(df_obs_ref[ref_mask]) # use all available ref obs + else: + _N_ref = N_ref + + # select only the N_ref ref obs for comparison + _df_obs_ref = df_obs_ref[ref_mask].iloc[-_N_ref:] + if len(_df_obs_ref) == 0: + print("no reference observations") # TODO: add proper error handling and logging here + return df_obs + + # determine reference observation values + y_ref = np.mean(_df_obs_ref[y_col]) # TODO: add option to choose statistic, e.g. mean or median? + yerr_ref = np.std(_df_obs_ref[y_col]) # TODO: propagate ref uncertainty properly + # determine the ref obs time range + x1_ref = np.array(_df_obs_ref[x_col])[0] + x2_ref = np.array(_df_obs_ref[x_col])[-1] + + # Create the colour dict + col_dict = {} + col_dict[colour] = y_obs - y_ref + col_dict[delta_t_col] = x_obs - x2_ref + col_dict[colErr] = np.sqrt((yerr_obs**2.0) + (yerr_ref**2.0)) + col_dict[y_ref_col] = y_ref + col_dict[x1_ref_col] = x1_ref + col_dict[x2_ref_col] = x2_ref + + # TODO: + # could also record phase angle diff and phase curve residual? + # need to test error case where there are no r filter obs yet + + return col_dict diff --git a/src/adler/science/PhaseCurve.py b/src/adler/science/PhaseCurve.py index a40d1cb..f0b0463 100644 --- a/src/adler/science/PhaseCurve.py +++ b/src/adler/science/PhaseCurve.py @@ -158,9 +158,18 @@ def InitModelDict(self, model_dict): """ - model = PhaseCurve() + # clean the input dictionary + del_keys = [] for key, value in model_dict.items(): - setattr(model, key, value) + if not hasattr(self, key): + del_keys.append(key) + model_dict = model_dict.copy() # make a copy to avoid changing the original dict + for key in del_keys: + model_dict.pop(key, None) + + # initialise a new model + model = PhaseCurve(**model_dict) + return model def InitModelSbpy(self, model_sbpy): @@ -221,7 +230,8 @@ def ReducedMag(self, phase_angle): Parameters ----------- phase_angle : float or array - value or array of phase angles at which to evaluate the phasecurve model, must have astropy units of degrees. + value or array of phase angles at which to evaluate the phasecurve model. + Must have astropy units of degrees if sbpy model uses units, otherwise pass radian values. Returns ---------- @@ -233,6 +243,59 @@ def ReducedMag(self, phase_angle): return self.model_function(phase_angle) + def ModelResiduals(self, phase_angle, reduced_mag): + """For a set of phase curve observations, return the residuals to the PhaseCurve model. + NB that units must match the sbpy model. E.g. phase_angle should be passed with units of degrees, or be in values of radians + + Parameters + ----------- + phase_angle : float or array + value or array of phase angles at which to evaluate the phasecurve model. + + reduced_mag : float or array + value or array of reduced magnitudes at which to evaluate the phasecurve model. + + Returns + ---------- + + residuals : float or array + The residuals of the observations minus PhaseCurve model values + + """ + + residuals = reduced_mag - self.ReducedMag(phase_angle) + + return residuals + + def AbsMag(self, phase_angle, reduced_mag): + """For a set of phase curve observations, return the absolute magnitude from the fitted phase curve model. + I.e. this is the model residuals, shifted by the fitted absolute magnitude + NB that units for phase_angle and reduced_mag must match the sbpy model. E.g. phase_angle should be passed with units of degrees, or be in values of radians + + Parameters + ----------- + phase_angle : float or array + value or array of phase angles at which to evaluate the phasecurve model. + + reduced_mag : float or array + value or array of reduced magnitudes at which to evaluate the phasecurve model. + + Returns + ---------- + + abs_mag : float or array + The residuals of the observations minus PhaseCurve model values shifted by the model absolute magnitude + + """ + + # TODO: add option to pass model & filt instead of reduced_mag & phase angle - will calculate all absmags automatically + # probably not possible as observations are contained in a separate object + # if (phase_angle is None) and (reduced_mag is None) and (filt is not None)... + + abs_mag = reduced_mag - self.ReducedMag(phase_angle) + self.H + + return abs_mag + def FitModel(self, phase_angle, reduced_mag, mag_err=None, fitter=None): """Fit the phasecurve model parameters to observations. starts with a phase curve model as an initial guess for parameters. diff --git a/src/adler/utilities/science_utilities.py b/src/adler/utilities/science_utilities.py index 68beaf7..90ff1df 100644 --- a/src/adler/utilities/science_utilities.py +++ b/src/adler/utilities/science_utilities.py @@ -1,10 +1,15 @@ import numpy as np +import pandas as pd from astropy.stats import sigma_clip as astropy_sigma_clip +import astropy.units as u def outlier_diff(new_res, diff_cut=1.0): """Test whether new data point(s) is an outlier compared to the model by considering the residual difference. + Parameters + ----------- + new_res: array The residuals of the new data points compared to the model diff_cut: float @@ -29,6 +34,9 @@ def outlier_diff(new_res, diff_cut=1.0): def outlier_std(new_res, data_res, std_cut=3.0): """Test whether new data point(s) is an outlier compared to the model by considering the standard deviation of the residuals. + Parameters + ----------- + new_res: array The residuals of the new data point(s) compared to the model data_res: array @@ -59,6 +67,9 @@ def zero_func(x, axis=None): """Dummy function to return a zero. Can be used as the centre function in astropy.stats.sigma_clip to get std relative to zero rather than median/mean value. + Parameters + ----------- + x: Dummy variable axis: @@ -70,6 +81,9 @@ def zero_func(x, axis=None): def sigma_clip(data_res, kwargs={"maxiters": 1, "cenfunc": zero_func}): """Wrapper function for astropy.stats.sigma_clip, here we define the default centre of the data (the data - model residuals) to be zero + Parameters + ----------- + data_res: array The residuals of the data compared to the model. kwargs: dict @@ -90,6 +104,9 @@ def sigma_clip(data_res, kwargs={"maxiters": 1, "cenfunc": zero_func}): def outlier_sigma_diff(data_res, data_sigma, std_sigma=1): """Function to identify outliers by comparing the uncertainty of measurements to their residuals + Parameters + ----------- + data_res: array The residuals of the data compared to the model. data_sigma: array @@ -106,3 +123,101 @@ def outlier_sigma_diff(data_res, data_sigma, std_sigma=1): outlier_flag = np.abs(data_res) > (std_sigma * data_sigma) return outlier_flag + + +def apparition_gap_finder(x, dx=100.0): + """Function to find gaps in a data series. E.g. given an array of observation times, find the different apparitions of an asteroid from gaps between observations larger than a given value. + + Parameters + ----------- + + x: array + The SORTED data array to search for gaps + dx: float + The size of gap to identify in data series + + Returns + ----------- + x_gaps: array + Values of x which define the groups in the data, where each group is x_gaps[i] <= x < x_gaps[i+1] + """ + + ## TODO: check that x is sorted? + + # find the difference between each data point + x_diff = np.diff(x) + + # select only differences greater than value dx + x_diff_mask = x_diff > dx + + # find the values in x that correspond to the start of each group + x_gaps = x[1:][x_diff_mask] + + # add the very first observation to the array + x_gaps = np.insert(x_gaps, 0, x[0]) + + # add the very last observation to the array, including a shift of dx to ensure the boundary inequalities work + x_gaps = np.insert(x_gaps, len(x_gaps), x[-1] + dx) + + return x_gaps + + +def get_df_obs_filt(planetoid, filt, x_col="midPointMjdTai", x1=None, x2=None, col_list=None, pc_model=None): + """Retrieve a dataframe of observations in a given filter. Has the option to limit the observations to a range of values, e.g. times/phase angles, if required. + + Parameters + ----------- + + planetoid: object + Adler planetoid object containging the observations + filt: str + The filter to query + x_col: str + Column name to use for ordering values and limiting observations to a range: x1 <= df_obs[x_col] <= x2 + x1: float + Lower limit value for x_col + x2: float + Upper limit value for x_col + col_list: list + List of column names to retrieve in addition to x_col, otherwise all columns are retrieved. + N.B. if AbsMag is included in col_list then the absolute magnitude is calculated for the phase curve model pc_model + pc_model: object + Adler PhaseCurve model used to calculate AbsMag if required + + Returns + ----------- + df_obs: DataFrame + DataFrame of observations in the requested filter, ordered by x_col and with any x_col limits applied + """ + + # get observations in filter as a dataframe + obs = planetoid.observations_in_filter(filt) + # TODO: add option to get all filters? calculating AbsMag gets complicated... + # TODO: split the dataframe functions into a separate function? + df_obs = pd.DataFrame(obs.__dict__) + + # if no list of columns is provided, retrieve all columns + if col_list is None: + col_list = list(df_obs) + else: + col_list = [x_col] + col_list + + # calculate the phase curve absolute magnitudes (i.e. remove phase curve variation from reduced_mag) + if "AbsMag" in col_list: + # calculate the model absolute magnitude + # TODO: add robustness to the units, phaseAngle and reduced_mag must match pc_model + df_obs["AbsMag"] = pc_model.AbsMag(obs.phaseAngle * u.deg, obs.reduced_mag * u.mag).value + + # select only the required columns + df_obs = df_obs[col_list] + + # sort into x_col order + df_obs = df_obs.sort_values(x_col) + # limit to range in x_col if limits x1, x2 are supplied + if (x1 is not None) & (x2 is not None): # TODO: add functionality to set only upper/lower limit + x_mask = (df_obs[x_col] >= x1) & (df_obs[x_col] <= x2) + df_obs = df_obs[x_mask] + # reset the dataframe index + df_obs = df_obs.reset_index(drop=True) + + return df_obs diff --git a/tests/adler/dataclasses/test_MPCORB.py b/tests/adler/dataclasses/test_MPCORB.py index d139d00..ee43b3d 100644 --- a/tests/adler/dataclasses/test_MPCORB.py +++ b/tests/adler/dataclasses/test_MPCORB.py @@ -12,7 +12,7 @@ def test_construct_MPCORB_from_data_table(): test_query = f""" SELECT - ssObjectId, mpcDesignation, mpcNumber, mpcH, mpcG, epoch, peri, node, incl, e, n, q, + ssObjectId, mpcDesignation, fullDesignation, mpcNumber, mpcH, mpcG, epoch, tperi, peri, node, incl, e, n, q, uncertaintyParameter, flags FROM MPCORB @@ -25,6 +25,9 @@ def test_construct_MPCORB_from_data_table(): assert test_MPCORB.ssObjectId == 8268570668335894776 assert test_MPCORB.mpcDesignation == "2014 QL4" + assert ( + test_MPCORB.fullDesignation == "2011 2014 QL433" + ) # N.B. that there are known DP0.3 issues with mpcDesignation and fullDesignation, https://dp0-3.lsst.io/data-products-dp0-3/data-simulation-dp0-3.html#known-issues assert test_MPCORB.mpcNumber == 0 assert_almost_equal(test_MPCORB.mpcH, 19.8799991607666, decimal=6) assert_almost_equal(test_MPCORB.mpcG, 0.15000000596046448, decimal=6) diff --git a/tests/adler/dataclasses/test_Observations.py b/tests/adler/dataclasses/test_Observations.py index 44a3717..39e1a57 100644 --- a/tests/adler/dataclasses/test_Observations.py +++ b/tests/adler/dataclasses/test_Observations.py @@ -15,7 +15,9 @@ def test_construct_observations_from_data_table(): test_query = f""" SELECT ssObject.ssObjectId, mag, magErr, band, midPointMjdTai, ra, dec, phaseAngle, - topocentricDist, heliocentricDist + topocentricDist, heliocentricDist, heliocentricX, heliocentricY, heliocentricZ, + topocentricX, topocentricY, topocentricZ, + eclipticLambda, eclipticBeta FROM ssObject JOIN diaSource ON ssObject.ssObjectId = diaSource.ssObjectId diff --git a/tests/adler/science/test_Colour.py b/tests/adler/science/test_Colour.py new file mode 100644 index 0000000..08e64c6 --- /dev/null +++ b/tests/adler/science/test_Colour.py @@ -0,0 +1,112 @@ +import numpy as np +import pandas as pd +from numpy.testing import assert_almost_equal, assert_array_equal, assert_array_almost_equal +import pytest +import os +import astropy.units as u + +from adler.science.Colour import col_obs_ref +from adler.utilities.tests_utilities import get_test_data_filepath +from adler.dataclasses.AdlerPlanetoid import AdlerPlanetoid +from adler.dataclasses.AdlerData import AdlerData +from adler.science.PhaseCurve import PhaseCurve +from adler.utilities.science_utilities import get_df_obs_filt + +# set up a good test object (main belt asteroid) +ssoid = 6098332225018 +test_db_path = get_test_data_filepath("testing_database.db") +planetoid = AdlerPlanetoid.construct_from_SQL(ssoid, sql_filename=test_db_path) +adler_cols = AdlerData(ssoid, planetoid.filter_list) + +# Set time boundaries for a given apparition (apparition index 3 for this object) +test_app_i = 3 +t_app1 = 61469 +t_app2 = 61678 + +# define columns +column_dict = { + "x_col": "midPointMjdTai", + "y_col": "AbsMag", + "y_col": "reduced_mag", + "yerr_col": "magErr", + "filt_obs": "g", # observation filter + "filt_ref": "r", # reference filter (we are calculating a filt_obs - filt_ref colour) +} + +# Load the stored colour results for this apparition from a df csv in tests/data, see also lsst-adler/notebooks/colour_functions_testing.ipynb +# test_data_path = "/Users/jrobinson/lsst-adler/tests/data" # pytest does not seem to like relative file paths, use absolute +test_data_path = "/".join(os.path.abspath(__file__).split("/")[:-3]) + "/data" +df_col_fname = "{}/df_{}_{}_{}_app_{}".format( + test_data_path, ssoid, column_dict["filt_obs"], column_dict["filt_ref"], test_app_i +) +df_N_ref_1 = pd.read_csv("{}_N_ref_{}.csv".format(df_col_fname, 1), index_col=0) +df_N_ref_3 = pd.read_csv("{}_N_ref_{}.csv".format(df_col_fname, 3), index_col=0) +df_N_ref_5 = pd.read_csv("{}_N_ref_{}.csv".format(df_col_fname, 5), index_col=0) + +# fit adler phase curve models +# TODO: replace this with the stored alder_cols when they are implemented in the database +adler_cols = AdlerData(ssoid, planetoid.filter_list) +for filt in [column_dict["filt_obs"], column_dict["filt_ref"]]: + sso = planetoid.SSObject_in_filter(filt) + obs = planetoid.observations_in_filter(filt) + + H = sso.H + G12 = sso.G12 + pc = PhaseCurve(H=H * u.mag, phase_parameter_1=G12, model_name="HG12_Pen16") + + pc_fit = pc.FitModel( + np.array(getattr(obs, "phaseAngle")) * u.deg, + np.array(getattr(obs, "reduced_mag")) * u.mag, + ) + pc = pc.InitModelSbpy(pc_fit) + + adler_cols.populate_phase_parameters(filt, **pc.ReturnModelDict()) + + +def test_col_obs_ref( + planetoid=planetoid, + adler_cols=adler_cols, + column_dict=column_dict, + N_ref_list=[1, 3, 5], + df_ref_list=[df_N_ref_1, df_N_ref_3, df_N_ref_5], +): + # gather the observations for the apparition + ad_g = adler_cols.get_phase_parameters_in_filter(column_dict["filt_obs"], "HG12_Pen16") + pc_g = PhaseCurve().InitModelDict(ad_g.__dict__) + df_obs = get_df_obs_filt( + planetoid, + column_dict["filt_obs"], + x1=t_app1, + x2=t_app2, + col_list=[column_dict["y_col"], column_dict["yerr_col"]], + pc_model=pc_g, + ) + + # test multiple N_ref + for N_ref, df_ref in zip(N_ref_list, df_ref_list): + col_dict_list = [] + + # simulate stepping through each filt_obs observation + x1 = t_app1 + for xi in range(len(df_obs)): + x2 = df_obs.iloc[xi][column_dict["x_col"]] + + # do the colour finding function here + col_dict = col_obs_ref( + planetoid, + adler_cols, + N_ref=N_ref, + x1=x1, + x2=x2, + **column_dict, + ) + col_dict_list.append(col_dict) + + # store results as a dataframe + df_col = pd.DataFrame(col_dict_list) + df_col = pd.concat([df_obs, df_col], axis=1) + + # compare results df to stored df + for x in list(df_col): + print(x) + assert_array_almost_equal(np.array(df_col[x]), np.array(df_ref[x])) diff --git a/tests/adler/science/test_PhaseCurve.py b/tests/adler/science/test_PhaseCurve.py index 77e0385..d67fc3c 100644 --- a/tests/adler/science/test_PhaseCurve.py +++ b/tests/adler/science/test_PhaseCurve.py @@ -143,3 +143,6 @@ def test_PhaseCurve_FitModel_HG_bounds(): assert pc_fit.G.bounds == (0.0, 0.1) assert pc2.phase_parameter_1 == 0.1 assert pc2.phase_parameter_1_err is not None + + +# TODO: test absmag diff --git a/tests/adler/utilities/test_science_utilities.py b/tests/adler/utilities/test_science_utilities.py index 3a2d753..b8a0588 100644 --- a/tests/adler/utilities/test_science_utilities.py +++ b/tests/adler/utilities/test_science_utilities.py @@ -220,3 +220,11 @@ def test_sigma_clip(): def test_outlier_sigma_diff(): outliers = sci_utils.outlier_sigma_diff(y_res, y_err, std_sigma=5.0) assert_array_equal(outliers, outliers_y) + + +# TODO: test apparition_gap_finder, correct number of groups identified? correct boundaries? +# check boundaires are in ascending order + +# TODO: test get_df_obs_filt, check filter and sorting/time range constraints, +# check the column list behaviour +# check AbsMag calculation diff --git a/tests/data/df_6098332225018_g_r_app_3_N_ref_1.csv b/tests/data/df_6098332225018_g_r_app_3_N_ref_1.csv new file mode 100644 index 0000000..9d0ab62 --- /dev/null +++ b/tests/data/df_6098332225018_g_r_app_3_N_ref_1.csv @@ -0,0 +1,16 @@ +,midPointMjdTai,reduced_mag,magErr,g-r,delta_t_g-r,g-rErr,reduced_mag_r,midPointMjdTai1_r,midPointMjdTai2_r +0,61485.3588,17.82546893645328,0.13199999928474426,0.6294446491054693,7.961480000005395,0.13199999928474426,17.19602428734781,61477.39732,61477.39732 +1,61500.36601,17.548929910716545,0.052000001072883606,0.4263135941635241,3.041850000001432,0.052000001072883606,17.12261631655302,61497.32416,61497.32416 +2,61504.33885,17.545916408017177,0.05700000002980232,0.42122628677132923,0.024510000002919696,0.05700000002980232,17.124690121245848,61504.31434,61504.31434 +3,61524.2819,17.60109262744087,0.06199999898672104,0.4764025061950221,19.967560000004596,0.06199999898672104,17.124690121245848,61504.31434,61504.31434 +4,61525.31814,17.530157869823057,0.04699999839067459,0.4042173365957282,0.024610000000393484,0.04699999839067459,17.12594053322733,61525.29353,61525.29353 +5,61536.27685,17.474481419144603,0.039000000804662704,0.3952621411742143,4.973070000000007,0.039000000804662704,17.079219277970388,61531.30378,61531.30378 +6,61556.37692,17.262436542715754,0.03200000151991844,0.39311468820022455,0.024060000003373716,0.03200000151991844,16.86932185451553,61556.35286,61556.35286 +7,61558.20527,17.330018817677058,0.03200000151991844,0.46069696316152786,1.8524099999995087,0.03200000151991844,16.86932185451553,61556.35286,61556.35286 +8,61558.20674,17.302025752576462,0.029999999329447746,0.43270389806093235,1.8538800000023912,0.029999999329447746,16.86932185451553,61556.35286,61556.35286 +9,61589.25295,17.086183244116974,0.028999999165534973,0.5580614970225177,1.968010000004142,0.028999999165534973,16.528121747094456,61587.28494,61587.28494 +10,61589.25519,17.035176049202192,0.032999999821186066,0.5070543021077363,1.9702500000057626,0.032999999821186066,16.528121747094456,61587.28494,61587.28494 +11,61590.002,17.013050114521022,0.04399999976158142,0.4279341627624582,0.023699999997916166,0.04399999976158142,16.585115951758564,61589.9783,61589.9783 +12,61590.24357,17.091348238945375,0.02199999988079071,0.5062322871868119,0.26526999999623513,0.02199999988079071,16.585115951758564,61589.9783,61589.9783 +13,61616.20243,17.44818378972448,0.032999999821186066,0.5004836109665014,0.19327999999950407,0.032999999821186066,16.94770017875798,61616.00915,61616.00915 +14,61648.07405,17.627385699986323,0.054999999701976776,0.4917527634000258,0.024980000001960434,0.054999999701976776,17.135632936586298,61648.04907,61648.04907 diff --git a/tests/data/df_6098332225018_g_r_app_3_N_ref_3.csv b/tests/data/df_6098332225018_g_r_app_3_N_ref_3.csv new file mode 100644 index 0000000..d55b95f --- /dev/null +++ b/tests/data/df_6098332225018_g_r_app_3_N_ref_3.csv @@ -0,0 +1,16 @@ +,midPointMjdTai,reduced_mag,magErr,g-r,delta_t_g-r,g-rErr,reduced_mag_r,midPointMjdTai1_r,midPointMjdTai2_r +0,61485.3588,17.82546893645328,0.13199999928474426,0.6646031133749091,7.961480000005395,0.1387787184136571,17.16086582307837,61469.38552,61477.39732 +1,61500.36601,17.548929910716545,0.052000001072883606,0.3998204721333032,3.041850000001432,0.061730424176583504,17.149109438583242,61477.39732,61497.32416 +2,61504.33885,17.545916408017177,0.05700000002980232,0.41306100220004893,0.024510000002919696,0.06719475757019573,17.13285540581713,61503.35033,61504.31434 +3,61524.2819,17.60109262744087,0.06199999898672104,0.4682372216237418,19.967560000004596,0.07148521046953224,17.13285540581713,61503.35033,61504.31434 +4,61525.31814,17.530157869823057,0.04699999839067459,0.41017516413599964,0.024610000000393484,0.047603929903510755,17.119982705687057,61504.31434,61525.29353 +5,61536.27685,17.474481419144603,0.039000000804662704,0.3696556612160329,4.973070000000007,0.0435304535193636,17.10482575792857,61524.3074,61531.30378 +6,61556.37692,17.262436542715754,0.03200000151991844,0.23760932081133745,0.024060000003373716,0.11609811521344783,17.024827221904417,61525.29353,61556.35286 +7,61558.20527,17.330018817677058,0.03200000151991844,0.30519159577264077,1.8524099999995087,0.11609811521344783,17.024827221904417,61525.29353,61556.35286 +8,61558.20674,17.302025752576462,0.029999999329447746,0.27719853067204525,1.8538800000023912,0.11556284964731135,17.024827221904417,61525.29353,61556.35286 +9,61589.25295,17.086183244116974,0.028999999165534973,0.483080272053364,1.968010000004142,0.09278157251016264,16.60310297206361,61562.31723,61587.28494 +10,61589.25519,17.035176049202192,0.032999999821186066,0.43207307713858256,1.9702500000057626,0.09410855558372905,16.60310297206361,61562.31723,61587.28494 +11,61590.002,17.013050114521022,0.04399999976158142,0.4132712746425149,0.023699999997916166,0.05221424730174361,16.599778839878507,61589.277,61589.9783 +12,61590.24357,17.091348238945375,0.02199999988079071,0.4915693990668686,0.26526999999623513,0.03569772593630112,16.599778839878507,61589.277,61589.9783 +13,61616.20243,17.44818378972448,0.032999999821186066,0.5299196193560611,0.19327999999950407,0.041117958024899405,16.91826417036842,61610.20336,61616.00915 +14,61648.07405,17.627385699986323,0.054999999701976776,0.5006338385771478,0.024980000001960434,0.08026750494410353,17.126751861409176,61620.16575,61648.04907 diff --git a/tests/data/df_6098332225018_g_r_app_3_N_ref_5.csv b/tests/data/df_6098332225018_g_r_app_3_N_ref_5.csv new file mode 100644 index 0000000..f083e76 --- /dev/null +++ b/tests/data/df_6098332225018_g_r_app_3_N_ref_5.csv @@ -0,0 +1,16 @@ +,midPointMjdTai,reduced_mag,magErr,g-r,delta_t_g-r,g-rErr,reduced_mag_r,midPointMjdTai1_r,midPointMjdTai2_r +0,61485.3588,17.82546893645328,0.13199999928474426,0.6646031133749091,7.961480000005395,0.1387787184136571,17.16086582307837,61469.38552,61477.39732 +1,61500.36601,17.548929910716545,0.052000001072883606,0.40214961118913806,3.041850000001432,0.06408282905927162,17.146780299527407,61469.38552,61497.32416 +2,61504.33885,17.545916408017177,0.05700000002980232,0.411249516593152,0.024510000002919696,0.06403779113904624,17.134666891424025,61497.32416,61504.31434 +3,61524.2819,17.60109262744087,0.06199999898672104,0.4664257360168449,19.967560000004596,0.0685261888982895,17.134666891424025,61497.32416,61504.31434 +4,61525.31814,17.530157869823057,0.04699999839067459,0.4033930271697166,0.024610000000393484,0.05524451285749118,17.12676484265334,61503.35033,61525.29353 +5,61536.27685,17.474481419144603,0.039000000804662704,0.3506598671390968,4.973070000000007,0.05091030826493725,17.123821552005506,61503.35078,61531.30378 +6,61556.37692,17.262436542715754,0.03200000151991844,0.20073869280633616,0.024060000003373716,0.10276020693938347,17.061697849909418,61504.31434,61556.35286 +7,61558.20527,17.330018817677058,0.03200000151991844,0.2683209677676395,1.8524099999995087,0.10276020693938347,17.061697849909418,61504.31434,61556.35286 +8,61558.20674,17.302025752576462,0.029999999329447746,0.24032790266704396,1.8538800000023912,0.10215507815432866,17.061697849909418,61504.31434,61556.35286 +9,61589.25295,17.086183244116974,0.028999999165534973,0.40181661190363016,1.968010000004142,0.12432445391335133,16.684366632213344,61558.23271,61587.28494 +10,61589.25519,17.035176049202192,0.032999999821186066,0.3508094169888487,1.9702500000057626,0.12531787533089708,16.684366632213344,61558.23271,61587.28494 +11,61590.002,17.013050114521022,0.04399999976158142,0.4366844446633422,0.023699999997916166,0.057457452706149006,16.57636566985768,61586.24199,61589.9783 +12,61590.24357,17.091348238945375,0.02199999988079071,0.5149825690876959,0.26526999999623513,0.04300417290467258,16.57636566985768,61586.24199,61589.9783 +13,61616.20243,17.44818378972448,0.032999999821186066,0.6651807942177115,0.19327999999950407,0.17001002587811542,16.78300299550677,61589.27924,61616.00915 +14,61648.07405,17.627385699986323,0.054999999701976776,0.5616368620560408,0.024980000001960434,0.10406682109366068,17.065748837930283,61616.22684,61648.04907 diff --git a/tests/data/testing_database.db b/tests/data/testing_database.db index 7a220e3bcffdbc65f92490b11c691334da75d150..455b28b12e546bf0c177958a4e41c6884dac604b 100644 GIT binary patch literal 512000 zcmeF42UHZ<_vf3OX#*&Nq0FKTlEFk>6%=#MsQ^hz4h9shCMf1Sh8eZYSx0Bgp`JNo z&Z*521qIBy_ZGD2=ivU&?tk~}9ysT{Pk&!^-S@utRdt09bm`PKA~A#ujfsy6Oyt~E zrBr%)D&RPkO0^103Kml=R#^ICF~p+#)6+dDhO1Q0@r!9?L&M!F>&_NEj342HlE0D+ zlw_bJ1OKZTNU!5!Qoce3z07)vfkBZW?K(E^(79P*ZQQ)Gf77o1T-T<}+WK?V3LC3& zb;#s^p9aA8{~bB-hBA3d-d*VVsgS6xGw#K8E(Z<>N*qGFR0L*l<_3i&2PNNh}S z_%}KpF{re436J?MqjBG?6dw{A5+4#Bj6Hpsg!m{N$(Mm8 zByyJA~(vQPR`e9g#{=S2x?}mkHK+=B> z3!8sj>nFn+{$uyUf9PJ(-v=80UH7O44FAvWvH8cfe$svNkKHH#(7mF+50v~}_oxOW z|7Z8u{Nq~RcW>^PG0a(cxtY*hIWi?Abcsod55}##kOn`#+!Qt|FE_Z!7pfm;Go)y9 zZ5Bd1@Rt_?U31&OsGx!0E(i;Z3=GE2yV*BS*0r_$rY$@qG9o58BswuZBG@k?0grrN zVt>1+d(VGc(Cgn8^!}Fx{}$=r0(~0?k89ncBjS=meoRkncwhqVu3?ei?#8k4_!1JF z7&6c=B>4L|`?926U|95zd;fnf{FjwG2V&swUvs~#+d1(2mE!RPU|>=po(Oz<0-|f} z{O#2J+YX%=gXbDQovCk@d^b_wwtt)6Z`;33?>Ft=&eAu{{}QDrh`D20z28jASE-^m zrlzv#f7zVBJkbAs!(M)V z67i+`AGfEdKs=54vL~SWz6)%iF1i;E`0Mi`26l|WGmv(J26hen_S`5M&qoT!^z|&~ z+bNF^9Dt+#;rJc@<0cgNgRO+UR|KhRX?XTw54b<-H-_*6$Rn-n^TeVtkrZ%wtX#Ltcx3G)<)f4`b^-D5P zl7W&8lw_bJ10@+K$v{a4N-|KAfszcAWS}Gi|1%kAJgz)Gtx4u6tbH>>u;w${Vhx$@ zSU1hI$GTZ2y|LYVY&h0#V>@8&J=P2BT4O6=-C-=f1wL}jV5|p^>5TQDF%7Ve7{g&5 zKE?#=u#9M|Lo>Qz9g@))>wy{7u@25K$2urI4(q`5o>&X%Jl6fwYhmr5PVbYaq$Of) zlhy}o+q7m_mrHZRx_nwGtZSqW#X2&zKh^_M{ju(t>V$Qt)UsH2)+A%yMH7g1S4{xc z-8A*F?yf0^bq{(kzvn=#dbPo7gd0{S73liEq4jF4Rb{`q{=cKXuD+-~tv;qcpx&e2 zrrxApt6r{_)brFc)sxkk>Qwbe^$>NedXRdcy1%-oy0f~ix`i6l4b>iMXLTKQHFZUG zIdv(ux!TbBll2?xm)6g$A6RRxuUenC&aysYz0Z1=^%m>D3OlCGOa4kSP?CX?43uP` zBm*TGD9J!c21+tel7W&8lw_bJ179+bJdQq#K0K2yk%wgpSPxY$--l%S;Nzss8dxW0 z(uH!u*pXPrD;M!`W4q$xnqvWLyRme+++eI7*8P-A`G#X+@v+|+`gFVd7`kYV9^-*^ zjB-IAJH`?pn`I2c+B}0Uo-HyuVr{8h)>~!N!N=AaHdw1Ouu_#ukHp$Txw!XC_ru4f z(_OIkO0S4@i*y=)OXU(jAWgu>tj9~&FpY+CnVhLuNknXF>9fDwr!(jd1x5- znKgHa1N{d#BwUYAJTO=!ruL$HB`QAZ)Io)xj~&qB139;agG9XlD1$WAXnTREm$xN9 zt_5F^W(ug}CS#SaG|0ntKM_fXen{ewe9kxvM*7j+IMVX%ZQ~ItjbGW~ZJ0+?ot~kd zb{Q(vej3&2O+JJ@ThXe?%$Ly7>{h;HJr3e?=QsMRLvPTZ;tCj$k9G29NIA{gMIsQ5mEnTYg#ox2B+a)Dp6!t;qf z5^7%RzgR?S{L1{ohE@1YY4pUiP^H0<0KFG8A;K+qzTt1N&|Y0(L(TpS5h=`DC8lD_@CuzgkO`G6X#*49YKA{%+H-AVCnpdk>4 z*|!jh!aBNRrIKHDd?FFGEpKUwe>8l}kQ%?zzp!EE`e00#bq5*OYWn1P&ch+$-u$gL zCVs5N+UU+>;*+7_j)DF=O-6zK6Auvs@=5k{6s9&!mq-OYU8g^->}7Nl3)q{UOkFNw@J!|1yE3 zlOA3^6PkiSl|Y68`IMdE0#Y8CUJi&c*Ga16cU5J2Airyv!H9S3F(Q<8T_qxQv(aG? zDLKY}9f#zTCiWoBcVit9D*0*FBY~Ln`*uZNkXA_eK;viq3mf+Hj+|a?ejF^kwiTH4 z9V55;00>FA`tJ{J{}F zOjyr|_fOLuT@^SiUd|zPHvV}5sN@v?y$BWGz3xE?DSK&@jYod5VyKq%^UL!E;`v+3 zrm69>D&B_4wigkLmWH}Rvg`z#EHi1E1#CXdd6cDBbtXN4L-V%(^$D`#bazO^#_ ze87?fi8{!rpTLL-hio_`)7!UEOZsYs(h`!N9Xyl>&V&8{q|Qc#0u(AfZnK|2OuDw1 zBMzuKQAEnG@)NYgI}Zwzwf3$PwQKw=i??CyUF-hVeskQwzMShSVfStb`!u$Vd7rWP zg4<}RRzDCDOg0^9aAiK|k2y*K`REPJv`9G(Y74}KI));WsrIzukn(f+t(JJoNhA$* z=RG6pZ1fw1LM4y>=Q2a`!OFV}D&F(&8XUFf@771)EIgM4B){Q~|Ic>X}bg~TljIzn9F?)(S7-I>ES?(*>uOQ9wg zc-pA32nGuuVLTO|-uDp&mP2eC2*kjMKPBS*)Okdv*X;?A&c4!(H24^}B8qH)-a=mv z6+g0{zepT^$O56_<3I9>K6hskhnlI3X9H4xIz!pzHGXEr+c0Q)+v}`nN_o&bR$z(ezg7zMjecBtH}|g?<@x5+p^0g z;(M2!jCc?JAPqhh?g_}k*#LuOJ*ZV4(yxb>I8EJv5yc&=L?qKYHU>z)SRKZQ=TD87 zkQzVJU)ZSfX=_mB>b1dm&ZL44>nB3^V}pt_&3i!G%@K{QO5cNoMYC$vNDBc2^$QLI z@*~5#2uL~Ql#*9DHe{QY^m|9;cq_kkX~2=bV^|7-#|@YM60*pK-U%O4vgX4IAWo=n ziom3Fb`_9JbvTkJKgdxjJTfAWq1O1B6mLWSO6zWW#=Y6VjGf!IijX!#TuYOqAL8pU zheN}L#P+BL6<2nwGh_E@Fxd4-#DIK;o2!=CIdUdZ_6_U`q>r?JNo0EG${_HSe|=wu z1pR4iNf+%KubV070P=`%2Q9Yp&Da$IDf@X1QuKk#or!YZ&WjxJUUV-6+A#iwjf;z( zIxfE6fmNDfDM#lCkU0EMe5LF-*4p9C*taM4LbD2E_dUz}0}SfT5HTP>VwtBv6fMFK zDn8tQC?mdqe4Hp39J8?#IMnDrlb1}d!^a_~;)NC;2qd+y29X3$|mtLHK8;vr#X zsdo2#u7dsmXNCd!P_{~ol+#O@h{S|39~qMAHJ&Vye)Mv^!u@-DiKxNGF^O1YBTu_S zKqbe;RRtuUuxUCY_AIwiM9Qxo?H7nwxJ-uRC!{9=QsZa%3mZ1uBUUA?uLh-SG`>`G z#9oNf^A0$aHwT*bE1lAR84pPp>zG<}z5@pKe_=d48R?Ktlt;cW5{a?4)fw^KkY`%r z+lNaPedkm#N)^~`4bUQsY!pP)p;mct(hX|mhb$Vv5qmD#Nm_n7Y6>IfoXJLw`D`J+bxDtY)tS0W#k zznT%-CwD>MKC#C^AYL3-xarPwn$Mxu`04+`hNr{fcgKd=@h`nRmTLpHz~BpRDf2GP z1V1;yIXviojL|18hbGMvh&k%RB9f`1a}@owzA_N- zHfhIUdy$P#?Ufx~PM+qcJS#ogP+`)AHWEHdE(mC)CBCxhi@<$j$V@qp@Ic- z>bWo7$hfXOCi$%tAmQnx%sQ9GvKH4$`#0{G3yrK6?r%|P1n6gPWf+i8N~)_Rrq-#g z=%}uaK>BlXFj3AuQ(hoG=w(j?E`I|cbsJT{6;g7i>pLalFpDUKp{_kd)O^m2OMsLg z)sN+fw}TRiR#=S9CabLPTOGH&V{zDg4T_S#k_?n&pdNd`(XP?CZFLm5c5 zH})|xFwrY~A4CV!!ax0@R^!6Km=>X^71nk-$EfM4FJy;5`0cn&XfrS@?I>Uqf3NaL zTEwq+ZVpIrIh8Mx_6vHVaE#_YN>x4n)sz}a!(E#y3g?^l!FF7sQ7u674xh%5<`eoi z#xZKT?{Zln&e_q72rhm4h@>mDoqg1x5L+5{*GnIeT@s+Yb z8yjd2!<>2ZOi+oV=7`uF7Ma`i9CSzWO*hQQ=BgzPX;3T>R! zjL$m~#98>}VX z?|g~Pgg+o13EW+=5Mj8r&PFi%5RlEo+4C;8~yU);d+|G zNQEyRhdNlLw`<0~-n5Hh6aTPlCxwSKtwhw|S-J2%03mXngqmM8U?&2N-fuY>sNwER ze~x&_OW9ArqiKNB9uWOmqV}0tuSBFxEB(_t_PM%Oy2Z7{%=p?fY~t^3>Matt*X@bG`T66iM7rmQ-lX~IEmKHW;5-sWamP1MNqtCrZs6zE>@VP_E(L?@5%gG{fu~SPe65O5} z3rN29uriD|A~RQu3@*A2kmJycB!BzJUXi$U@pC|e+mU?;72k1jV}(=w#&Gz+SGoST zvQnwj)W+81tSzl3qA2+*$v{a4N-|KAfszcAWS}GiB^fBmKuHElGVuR<22zYGs;sQ6 zObkpNO;nXkOiWd#CZ;CXV5KrPvEuh?-@Z4poU3}W$j!yom2-1(@#cDUbZr2o_t!Y$ z6=I@m;P<<`1Lg``TCV4av3o}C$kq3P2p5wRZrx6S(IMfX?3AIq69(Cm*Y$pVPmMCt zbzkpi-9O0Lkg;xZqp#h&dT^cEdSG9rH=lp>UPYr)TE6t}e!Hba{UlQ$ zhxQ-3&EL6oyN%V^NxHjDklnjN|EHio!cf=$f7mzrZ?Io@KSyVNgh8>LZ*5cOcD&vV zo%t@VE-py{o(X|a$A2nYXK-29@YQ>Hb0F_M%Dl+jJ~RX3%g!F~Zp;BN+HZN(_T|08 zTPBdbb)X#EbifGxVmm)}o_~DaJowwr!##TR3WO@vsy2&x=t7-S!R{_b){a&aEnir6w%B4(!F;6IEwg5(%S=s81{-G?yBf_jd}G+xKsKnR zpP~0iuMMv9f8kGKU%^IgfH$p-mG2K!!M)SO#Hf0FxP*7C-J`u#bO8Z%uhx^N12<=Yh%4wa=jH9=?CRCP0|#gGv%y7Ju4>P` zf17z}`rAm~s~>;ka9z@t!{ETODZIOyYhPRUB^&U}U5$%X{MxmqfCR5990CVtpL#2# z+rh7cOLUpmCfm2KE+C-pFa zb^gBZ=A_=+LHtMvtun0?czW()J> zUk4Xn>3Z0$yLEK|0rdd8%`2hSkju29Zg8$%lsTN$3T1vaws?~#?d6yA^4p)VeeQek zWrBdKqJBN=X1kpg)P08RES=QWT=#t#&kOrRT+Yv(J6A%2=b2GJ%)T6!-cGO4*Re(P zX@1^w$W>iHKz(L2cOIGz-pFu;FJp6acXo64z}=zrPsUa~%yL;T=3RDwa+ANA@2eTg zZCcd(nx@=p-n?$bzXr_FUN$fKk_&ic%@#2@VzL1pq+2G<$nqh2yap67-Ugk?l1Y1-Cj|Ll11 zxi>|?6}}A41EU3L^!QB%NoqzmkbfeHO zUkUMiF`v=$b(?-e;Na?)))KRCq`Kkz!Cwa#?#e|tl}*%aNXB8z9{J^LpsyW>QKfaO8HPcC`Ux~NY&+Nm7p3~cJ>`#Nl>f$pmo-pRKCm-7dT4qo0% z9Y{m=11aXA-O8_n3%_>sx7ZOgbO8bNK)bUH%6)(Ea&z}~c5!w0^3qf=`N`m_M%YF^ zXAP~|SUa@Y;+tdeO2E!gw{1S`Kzn{=%k@Xk{AQ{90tFuSYiV&gzgc6);$=IsLTe`kx z?hS4I?*s|H=j%R@a8=YcFDq=w3qAZa*kiqZfBm9wMSxew!uLhj*2!o1Jn-Vn)AQlE zZCHkx{O8xPB_3w&QUbf_0s`tD&I7H%WzQ#CK{u;juCB_F6%PCVM-en#Bx`}(@?I(Q9YA}#;3=zPg@b+i^WWM8XQYLR!| z*TIGV)qnN+8P3;I+Dmh|Bp^<;HSI@GPi}prcpBmdnw5zYZ?)a+kLc zw^`@{0_vP5v(|&#ab+_o${ZfOa8_O1X<7|_Hn<2=uikua>{YAU^>ciWZ3Mv1P@fNd zJDB9W{OFwFADgZF z4gsE%kJMspeyQ#>*mMR%ji+;F=F8ive(FEnhY|40`QP&ZjHW*2yU|4QQJ;p6P) zt66^jd&`dFN?(6Q193#Ftk7Dp^~u+T|m+G6YK! zmP9NGSmLq7VX2A54od?p{jfB|(g;gqEKRWZV&Sm>mZn&mVQG%V4U0FHT39+@F~;JD z#T`pDmKZFtSj@1PW3j+Oeyp%qV^L!%g~bC4`7Mpb3rhNd`(X zP?CX?43uP`Bm*TGD9J!c21+vUe?0?qk!+*tsZ#q`$6Bql)LPnEbTrR2W2SFRT};A_ zC8Ntmr3?cMMExy#&-AL{6aOdwG%YfQ{rvY4bie#(w4utf3k+{^8a<6z^YooPB&>CR zY?r)&HGSF6&*t|T&~#y~564b@;NL$#AYrfk*cXTP`2r*U zvvw2_;7vCXsr!@*pI7Q3hnk;tB?FLh^6L0utO7e-M#hoD^(Z&QSAW^HKs*_UXG; zBzD^I7=d>#XJm53+;>d{;=MKHMI<^IQoIv`xjnn|vF*|ktgI`SN^N8ViPj;vy7cP9 znm@9bb1!QS_^mwAwO5kNzqi$A6gl(36Nxy&Em|bTCO**;-v!qfkn*MUvIv#DeNz(- zi9hMLPeOii^0e~{lq!Dmu4Y6zb=X=UhSpa}NdCU#2n7D-)X=4*<%0%#q|r%I@lNzR z&419Nl-FS}w6}_zQ{EZkU#hNFIX{y%ynO1}%Te^;_h;Nk;BPqE_R}IYei?~BccS|-%+Fpkq_%&I=*28tSh2GD!+m8{uKLcnoS zsdLv3$3pn-=zEhwpMp{DFh-Hnp3&Vjl^h>o3&harcNON(2x6!~@8)bRl7Cue9YO_; zodZQ=(bwxfwYnss;@c!=Bk(uH9Lgv<%r=!nE#JH_Vwz?Ql#g7^^O&~@F z#)`zi9%Z$}w|NT?DrWNJm4uXE_A8u`noVX4sEa;UW^7-Xl-w^qoDtL8Y?Fvl!y1W5 z{?VRXA~VpM08-B0He8EDAMwTe_!wONYFWn-&zRl7Wepb!^&m(LwYGl#C-@JTxqsr6 zrx1BMb>|(eDH#9NT!}oaTt$Uq`gDdp$gz!gOGsvLcP5d~w;90^FCIFIP(dm0Y67xo zK0eA_sg-YgysE;H6;FtygSQ5gmfx03VWeO9)e(uuR`1p#(MQ}bd^9{YCg;+cqs-QA zNBYfsbs!?fZt*Jj84ytGQ*d0DUXU=QQZL@u7|d2tVEnPXcBLhxEFLMNB@Sxch{%jp zOBv~cE!Bbeo3HH{hZ z_JkE85`Dz}!iQZ--~8O3DSXcDnFhO7q(JP?jJG)94ksr0^tw0>ir6KSKQm;vheygWsl9^VkOIT{m5-oN zs#?$H%D1}&q+}gA2ys1K=ekWA!5pkxo=<9b8mhPJ`XVoO0a*6zfB|tn?6*irS@0Sx z5UX8Nwq*IX$`?>Gv-xx<9k;LDYl+}eZI^^B@?myZIYRIeq3I&F^5Y+0RG93#i;=#X z|4EBvCcSPT@K>o02a-k~QN{as_j~O4hG$>L@!Qglqzr9>zcj~JmK}O@XEwc^y)CES zfoc~Q^%->X09a1#%V9)b(;Ey(69HX-tG^r!1SYL0nM(`_s%2XxB!97#ox&?) zF9K@xF}Qdi`Pq-h9?L6Z&B_VqLVHA5L*T2{7fM@2fM2yu@ip6MA=1O!GirM*nCKV! zh@Y;d$g*bFM}@Ix?ntPa(eJkyl7G=45`n)2m{6TGxXd|1)cF`zO+qE7O%4^2eDa1X zT4L~~F{I`9A03H&UMn*#@qE_+35h-i{ldqX5p@Heb?L#%R}E?|gztvHHxHM8tTq?? zs^nYs;5R~~?~5H5>;y1=p$`~QPB?W-B5K5Zff&2p9+1ok&!&)k&ei7<@uc%|((qda z8KLu`S~Z)LoU&jTkx%J)oFj%BEYc$7ck)0W{i?$hj(GOU9s!9yB7Wh+vD_=YCa0b; zTcbY*3Nt!Eki|u9pX2?&FMax6pXM7NvRS0Jc^KXxJ|k6&s7NIw47?+E;{>5)?OJNgJO-bY?;X!rcG&;9s2H+jFX2XPR4t71P*>3#gAYKum(fT`8xrX zEIupE5yLlF0I`L$t-^t+7bMioz_dJvlyl4tNu!Uj;(fflRx0k@6ua(FCem^74XdWm zvZvSh#r^6r$of4x?wt_?&#$Y#Wc4?gjFK5e78AN6@QnA+IgTjI_SYgob<~H*f3(>s z5TCERCnCXk{$vSR=!5^<;IV*;AN;8bQ6ARqxz(Hp>W`Xd1&@|v!d5h}UWg=ix0vuz`g-q6UELo)M5$pYzu#hXY&h5Sb% zQs?7kH-<_UE7TQ`d|cH`go+>d`#3G}w4INLWCpdoh;q)T+An*dO3P^ z@QMG$Kg`p>U*st|PPJ#o_C#jgs`=IQ8lNEQa#K~G!6x7{zq`=BNmWQ%Y`swoS_bCP z${xWFsdGWZQOF~+J4(cmO7j?!8Mu!Vi7Im*fiq^t015cndI?#SC!;;Rw5a&xV_g|B zwO^*9<1%*&sO17pG?A&||KUjI4s0bN(MNFcK62Zs%d~eqRF#!pr#~!4;{ZWd5U#ME8E0x`~1038#)&&n0irJY)~s5qF8Hkb^Ge&4Tu559 zqsHB(!@)dM3mB14dfJ{Lt_arOKdD4NycBsEbAEEPMYnI}mXQF_>*O8+KIAq~`@E>OO6HxK7@7^dJy0h@Sp?BmCY{pm7SeQ?Pl6$_HKd(9WWoG0 z&+BocwnNKiH4gJR24E3ahr@_`Sa>cYj_zGaVdM)FBA6aClaO-3wHpZho!l|891?0Q z3uVY6ABHPuim3P~&3geUi|0}qanMZ-hXmEZu?YMJ0-@am;-w+VY@?5U#rw#xyg{EjY?Pq+IBEw^r&E&2I_#i(^0M&`rU`%(}5Xsc9( zkw-RisKHd)O_X16^3f76m&qdn_i+kC>U=bv4ygG4YnmbOY*xKSA`W`w07y{HPZLQO zq$Vo7XrSyz=%a7(KAycj7{7DsOB4By$IXWhS(PEKOVl%N@JVKG<9>BQ|LIUIFDP+I z&lh0Xb%7Rp;5D}T0#X)s-I9p4e%CVM?X0pQl9^3CtLTEA+fm|m0lX5CMLx_cyJ}JK zv97=oCv@1zh{-931>)7l@3cr}vc`q9oL}xQ(&(d4@jjm9O%9rIeu)b2^Uc&W%e{j4 z`JSHBJ06Dj`KFZ(3hjjV`4XzmObf^Rd_#dEr#3B$?&Agyp?t0b_Q|hz-~P?f&Tbvi;4+U}@flBJv5GR$+3% zVuAFaLDwZD7}t3&BAF^^1w+dB)+TC^;51r!g({j4<6Ao=TAKM}wut0oX9oi@sJ0hF zg5mQs5>kF?*hEV_HMS#yeDwN-kBjcMen%>_V%D3&eOfJ7;qk+-$D+Sig3sbLmMp#k zBuRgXNux%9xuz4ti2RV_bt$quGH9j3f#Y)oe8BXtR#d1O*_0^f_Fce`;MAazh%EA9 zWW0hDA9poRBBstdhCsc~kRj##6OSe0=Yxf_oO5X^wWE)o#rt^NYwT6$zw0*w+eW2p zw1}+?iH+wDXnnd9^J};A$%=v3pvB?R<&)bwfknS+7*NHBj9ad7^ubyj>8NTVLxP!p z8Pan8qRLw0`Q|+t5*!=87m+$2K{Gg1{P68n1*AOs+8r%1Vni|#^qear@ZZj?Y{U_- zq>&^aJ&O17C_X>t_^BOEc-`y-w*tfPy7{iPRlqX5Zq8=6#-!tQ^W!X`VPm{*Hl>Jo z-8>qI8jYzyjNM*MMB;UGHzHm)PX^-2-@Fib-R$7RkVQUJtIi0hL$mXPS9`!9Uhhd3`i-nTuJE;aS=lZqaYWcumK`9K2}@cSj9 z;YkVFjcQr?Mvn?$X>(Yl$O)@93B=?nXA!u6IR2qUg3fqNcZletTk$>~%yOM|DC|-r zFdLF;Fu2w&Xg;mIDtXs5=2yDi&N@N0An1Lp*_B!Ez{FO+~3R|~K66L-9jQ1jhy z`-;SOVRHZps)ySUI4=k6h*1CU9>gLadQn!S@G7Yi5Qk_#Xo-Dd>;&Swj~5Vly&0G) z67L%u3rO_QwRj)*mz*u16Bc`xfBEdgl0$#)h1huq%NCe)g3|HNwY^;Sk2hoEpaJ}*jFpKn|h z3q|6SyA6SO*J&kb^wFt!A9pV8-}<2g7b+jNb2-;waTEOaH~dV#o&&4c#_yp=<4I6M zOlvb>Pff5)ao{i_uW7kJK+0XSUnsh!GA|J*H zRe_c!G|zn-M!+*5C~NZP6Xo3H;n{vuB^2IqO}88+pg z_qOiGtO6l+{>&kjA9+Ek0lzy`&Z-MG8RL&u+C2{}cbAuNw7jO?Yes~G_8hT&=UfE- zQlg2XnR&Z947Hs9Ih{j7ty9lLq%KbbrjU{o%nP^vbnCW2I&z4kKs@>=n?o|g1vMDy z+%n411n0ZMFMK#}u5zP^V9spfxaiZae?oBNhEruay@Qs|%2ju|a|ogi-ss=30Un>i zdI=OcWzrf@IMYxhMn{cgNT#7EfXIUU>*M24qb$6{NU9|swmW{Ni_*4gbgi3C7 zUXKym8kJM{K|Y~H%|8kIh*0qdtG$+x=%an{K5idO3G_d-EQFO`?r{G4vcDkcNyEl- ze9A+M_%1KDTy=#ge!$I;Ls!6fMq%cXo!c;sERVjzi^Q0$&j`G-s&HLE@~eu%G{x_5@&5tMbiwCoDNm`0L_~HbH zlr^@`Ny~$ch7*}_kEcL7zshihkGKCLq6W*w3y4MAmEN2y5-K^l-*IZ?6He@Ax0nBGE_N;(gpaQuW{*t9{*B`F(w--_sm~AZxRu1yV5hwQ}#6 z<4_tRd)yr=xs3wgOrIm!5pLJj;5q+|aqi@+`Vz8i<+UzWP0C7yn~PeMYS3LS~M z?dp{vppp~oA4cHSo=qlcKFoT#gp}W2t*Isb(tn9aJmY^vK%$Q}zwqG_r7v_{lf}yK zsI*~!nKlq~cd=2xc_(NgTnY<%-5#QRA_sSKtp&!j$}o&5CtDr{;^^!|ju>-wBWY$( zn-h_I!G|ZL`Ky!dNP}hP9SEI|kND4%Q1S7N&vC@WS^FhoZ@n%=`R%sCS${Tk2~hi? zy}<$!eY7s#$F#35qdyP>qdQ4R(EIHG0?*$rw%~}Toi>O_^wH`UJ`PX%lw}^UgITRfd${I(Ux+fD zAs(9k4>Y-H7&vU}JxJEggyPa3lHHBxS4PgDRhtsQtCYd}^29^;X6&KXVGuqDCJ9 z#rwF@f7HZTQ))%9a^trpnq58ugVHB`aJh2@{L|m$N$Ots>qp-I+jL#Qs=|E8Yw>bOa=mHb;BVoad4B@ z0#g1IHH9PQRE-gd_qS>Qi9T8u@8fc-n0`Aagw~XII5o~XYHtH^y^U+!t(3s*r$7AL z?ov9`aK3+JS^Y|2nP!Ip@i^GF4@1hGFEuCfPLIMk(#Ptx2P8AEkcGhWg{{LyB-A{w zjHJki(VN1(YiMM1dKNF19wZQB4_sozOI@BxNM=y*p0u1BvQ$K(j~2!Ixb%GSypWku zSDCtM5_`9IHbhNWQEzPC6lj>T!R%gicSst0`0kAd55RoaEQ-hvTR$0r`}?2g6bAjR z(&7WAx8W>PcwrMD<(FM-5h|$ihlPOD?H}DeNy))OVnrlB%3!3z=x&!K)N=lS&OrPS zRe->M={O*Sh(7#(;bZrf7vmz&HUZ=5>2uC+wS~9~JKCQ)eTCJQEPd7$JOvlY^s3PV z53pPp&M+dM60;MCF^8KfI%si%qCaGm6S19{v`~D>xpR~=NpLDCFJXI;57P==IaGX1 zQ!6b}9{bXrv>dN`3`qX*=*t3<869XRA?2J_2N@E5`2E7iqPdeQ-WdN1%7u=4epP-9 zi3vH+&D`{ue;2>J-iL=nz|x+B7B`sz7Bil57?B_Jm%Bt9HR%sVj2zQciv-hN*NF04 zV{d_YbwU~uZ1vtUq;5X4I%rYJolk@*6dgZHq$AoK0wkzroD@jsJNPowErEyd^GMC>@aL6Jb7Pmc7V3f%MQJ&PkG7yLK zESwjX*o4}dSym`f&R2O!Nc7RHcpn#3*+bX#cJhXD^>&usiodClxa!V{(2C)#MR4)6rb&(`+=i_K+4jL1h_8w|vhO^rEXq;+d85=^djlaTWJs(~W$rs#se|N5zC z0!ZC_EVtC6l7~1gR5(1bw8Eggg%Y-dp7|I>zp8YVw0wG=^7@NDnilWle280G%cgAx zvvONkpUdaQ!l0;`o!6KKg8%jxVO$wcNURcEX8)u8VCGeqBJ&|m%624cjudVokw+E? z*v?FByJ(Se&Xchs@j;*0fCR2j8-^_Mp?7W~hl=kUPzizA*xQ8>2RCumBIQpl>I=jd z)siIQeV2&{nh*Ge58vE6bLW98*w$G&>*_`)NNhA|=&rnm%&*n9aqp&7g%-QAm!8&s z0E^z&wHT2P>5|5XqfejVh*7m)0TRqitw_uHDT}ql^ZKnA5*#<|;*do?jD{Z&Q1K%s zxfA8lH-{^XXgmQW&PS6%A8(EuDeYHwE+nE4{ue&HYu>ZC5LT7X-m+oTpH<32tmLxp zx_bnu$C(b9+q?zX9Jb6H<`DyyyGmm~JbpNTpvb&$l`RU}?eS80EjLfZc4pq8CQ&Xh zDbJ8lGj}!VA|J*cN42Qr#6M$!nAW(SL>yc#OiMhnxsHHjhC_`JIJ^080Es?)i}!JE z%90Omk-~2IiTT6wHnn3R&Y&>Hxym*X5;Qb43 zsYLW$Q20v!M877rGt+kgM7dywGD4_dwv2?j$cNQzak z^7|^rFx1RE|FeLU-?Zt9ARkSN_i^@4f@!;qR?9 z)fz$DTOEyVzZd})_eyaXk&p9Mwg7%)n|ss)ZrQc&Xt_tql~E+2V>E>C7dj%!iLLwXP3h?BWP2t2PI#WN(+ z^FF6gb#0u8TFx0-7mzsLjf(ei#=q&wlid!q#{1K%u#WrNL9q4o1;*KJ@V?cnleNyT zf~b()*G?P^0poH00!2>F9>Nhv2VK+>W5!KoNTz>ch=}B0Uyu}D2P*+J*wu3+7S1>S z@o*av6+hIly~0?}+6p^e_>*C~{Lb%`7Rf&y^H?Arobi!EqK}5f`#3Fx-S6>x&(dJ< zcl*(Q1k8l^&K{oI^=>h@xxxtV>;>R4?d;}%;_rgxfteT(_dnwX6qyf6`%7WJZs`)~ zye;jB%(UP&hm>D!4g@6VamseA+io|cN!TiDR`da4zy9MzVyOqVL_Vjj7^46xE zwX)tD%C{|g<(7H61mb#4PB}ASCv)i5(|JH557on5_lKOv&uRCH7bvo(Vl0Q01xt-c zthzB%OVl!RKr*vM%Kj}E?C~TGF69SG$RZzRYyYHHK1>tIh~x8=dEk?6UvQ+a)eb@6 z*$3?9kaFJo>jDyec>Ti1DG!x4;J1D9_2?;GEIJN{xZq^=z~eQm)ZX@NRO&LkPFuDi zDCs&_4y-O>M7&?Qiy`H{xqp+E&Fx!q#OtdT3rJ?RKS@ivVDCaAIKQZe(D^W3f0Y#9 zc-9G#IL>_p5QqKaDiE(Xu#%9>xI0gjUmfKHB>M3Dg%87)b4by-wuA; zYV=yAjj^D`wR|C;pz;>+GA5%`@Y_Jl*C z50Bz~oH!rq@W*9OXZ$?I+U|${$;8iNc!GK9oA`N*1-ou~cY&5o4|Kb1mxiCmsHj8^ z3p>g10e&7MP|=Yu>S*x+ejejD0g0c-sD!}JV`Nw{B-B_rNkA5DSB5L!N~rj#{T~qc zc?_*BBM!Rh&mr;i82<>QUx!s^#7je`Ymw-~{TDv6CQS_8f3gC8uSDFv=fVQ~Udg0) znQ!&+dnNA=_~gCA@0F~1m>V}8zgOZ-5#?d!9!o?qy{|xwZ7cy2zgN;mMDqEz%_ZW+ zBV7OqrJQ3#WYK(ll#NO8?aQYy;>a?d2)wVcxxR$N@0HBtNWW-yP9z>%c?*!}!|fM7 z{yF)#<(*V57|sr_r=GPC;u|^|)?Hx1T>YQC7`dY_xSzE+(Q3ODSnjXRQRLBw3`An+ z>DEB(-=~K_%-=klG&4Q6L_o@~hW{lZ`RuC7mZ6)EB{!9&qEizsE|sJHnvB5v;HLFB zB>y~ou@=en^6C-gr*@tUi9TF^;iJlqH$%b(UI&BUUGJ4QiiY^E?{{As@|d}?Jq`5n zzJSMyhkXNQ?!)8YHVGs0qPemK@F9JZInw>SCX2+pIn^X2GyR~CKp(sOh+q&}ARu)< z{uU%uvSxKOwaUG>-5|}YYP}=P=kyQPBAMRdS`3ND!HpaeeYpI>hxMZ+f0rM#k?Fni ziAW7-0sXAcrbPDL0iGQ;Pru`E5+eNllLqWC0HeF@MT(p;-&i7s!We-l+^L}G$N4ct zFbK)hBKa3}>WQTB7omx|?W(jIQ1L14XA4MqplghlD7f8c#K&u|OGy4rpBg~Cvm}{_ zKI#|mevFM4$&B`I6_9v8rmKWRAI`tdJLmNr@#=%AL}uRcJ|N`+>!utM zYL;1p(D^X7awH`u9aCE6)H$^T(h+A`X^Dr2TteV=TIcf;>7466fJ7g^74PFHv)|aV z`l6Lkro-rK1^{!VP;H>#AFHZSKW*3D_c>Ihq>Tv#Z z9oJq=YWv;e@A~K`$`8%R%11w+j+M=~yQ80IxPkNf($r5hIlWpJ^b`F`Eux?30|}}# zs>V}BKhdEwqd70SQV-|wrmkIN-w)@W}p1`F#jGixRzi8GQTX{i`xEed34@1vB} z?08=VqMv9v68%JF4Wf?vf3D-o(rzgOj=yOws-Y2;qSuaMp{)*`J9uIPw;M1u^|bvx z=Con0?J@hI=x_dedmGs|S*G96bwN!9iY*;aN>Dv=Cz|9|ZVzP2**n@&W)4ob6UYbPu2CP-cd#eou$c-DBlWN+0cMefUy<+MJH*QFJ zirM-5W%7CKdI2&UpH`r-4lzIG94!^YhMi$#c*FsT9Q?gKhqz9Zq)`4eX`k*DUWX&2 zj=F!YA`~h@7mmQ?Z(X8UXsuv z_FR(z`pNWOtLIRIXBsj$n^H}H;=3&xz);wu)-lMMd=4Uu)<>@u0+Jf*=c+()wC`og zV&^kQDRcBM*C=PRMT~sX#FDZt{`{A}|JNzzb!zvaUoWlwT2|)3@=AX#4rQIco?3e2 zeh6m+LgLnsNM=5#S4CB^aAD?)^?FsgW^WnKT8!{A*T|98`*X+)>+6q;KPT(&Tk(r$ zCw;lyuksQe|NT0Nzm`^}Bubg~0#v7)`N(8&pZ5|Jvnu8T?6pVrSIBprLS=~cUHi{< zT)uWIdfYyjGR!R4q+Ji!yDT(q_0|OUi@aXstp_ioKQdQ2x6ZiDE1Aui?~GcWdBRm7 z13OLRr1Y_aMCMfy7~~cXR{^ejtlKC+R_)1K2^Owz(~YYYBsK7Il0ZhaQe@IUJ5GV( zhi-bE+}OJcW%X0m3Jy_6tzzmJ8$Kj?&c!-n`geml%AK3IpLYBCE6#DO%wF54^3y}C ztci0(FQX&O=1NN$EsBw5ZW7t!Mmi_^)W{XcB*#`Hw`$l*hHBQH?i{iP3u84{ScgSG zl!T;KoSve)U~nUW@`RI@0nSHTHs+9<%{Gytnlbge1W`xLKi6^T#QFTk`2*9r`PFu{ z+7>is-cKV}9hp6XxjcP!*=R^B=5J{5EY4>cGfi=&7S-sFr8pVfw~9vk9C{^_={1ri zh%;z?D2tzmae(hTr&2Xo_`dq-V9RLf$;DxDdPtcin*~RIhk0c zXnR{nqfu9KJD_*aME{>H4+!tN)l9G%asyY5t| zWuN=60QR$`et$m;?zWvHGfS79462!tU1jn?)2{+#&Wo>0uyB2t)%BB+#7J%j@cDSY zL3hlDSdtl@H5Q>1hGD<|G5q|*D0yBK`<9jZ@k;0YO~w8pSI+udR2z; zvIDL(4Q%4aoF2ax0*Bzd?$VJ8TEzXdfz-09IPh%JnI$dpg~Q(1wypY|lDoW1LJ z9hKc?OUOl^kFwD+k~%YfG{E*f^9v`(BtMkN`{CCVi2G@q0qmbP&8LhyoQkRA?i7ax z>y7rYk`LG4S$ofsx$n>()HpPmcQhFzHgK-We5%zrQDC=`nYSs>TaLUrTE-9J-6AVyYn;Kpr75- zMh7mHMnAh!`fu=`7eJz)-AhIW?pwtn`q_P!DSz9wUm@QwiQteK7zZ#|xIPTWl~It? zke)vUC{FoRhq5>};F$*1f};l+`7yf3D;F&&r>sJzJ8e{<`z#e8Qqj zENJ+TJE7yovi1w+FHAoFnN4i@x#66}m6-MNK?+(FXKtFrp&GV#DP^_)=~RKt^@@=o zHyeCfp`7J@pF!r5FiM7nby)nIrXZD1$X@n2i8PIJqbj65)&x1aQN{Lpgcn9t;zZSB1}GV?%1 zp_b*R-lRCS%r4y#$MyP{f%_Xq%kvZNY2??2&t>wbL4g2^-dC3YP!hvi94FPVuMFU5 z_Sj2O%zAlzm(BQ$c%G(Q2VdCfnQVp{HAM&8+Cs6o_GxtKaaPfo9A`|BVt z(`n73vM1)UQ77$pdQU&kI?Z@jdWval7NoT{`(zx+%!9(IW!0~`zCP5bZTkeup;mqz za%--z*9di1WWfgYxdeF0UIaV*}Lz5@Z=khK|d2{;B7|Mn%zu6*h zjvTu^cklwH4c-^7naG+@>y2W41ThhcSwfY*O+ z3k|9f6HCaXf7K2oGg#iALGekeJc)enu#7V5aQJf_r|*6ej7ns7!uTq*#|&Z+l z-~Bh2XJ%yM{>mnCp09i0{>sMX-uv{%{gt)i{~lcl_g8w#9oEdn{gqcaS`@>D3XF_s z@>(MOzic3Je`S;cxnapx45}#x-AUGHp$LlVu;}WjA*nI9a%rhJBVn3A2H!MRATx8R zLs|XS=6AghT(@h;sH5Va>p1ns!;HTl<-<$51O<)ldyob1Y85}q^%QfvU%{Yx)C@Lp z>F|(GWtKCmzt%CdC{Ag2Qi19;>s=b@Q?m(&+%V6F>CQ0)SRa)dN|06e3T3cxeHib! zOiR`7gS!G8-Hdpxks)nvNKnmc6Qq#&bqhH8)p07IIw}-X$437RBQx49;N@?hNbc|` zij8`Cy;O3nE9=zb{L+;73t6yh^Y8)BnlrQ2?c=+JvML+eD&+W-brR&pO^L}%tdEonwmtx|J^x(jrX&A?#VWr&xfnt2(mFs28Y-0}NK9@TbxDewQ*vZs# zWOqXW;<`$x0!5Fx7c}zvEEfS{9D{`<#xa;qGRv>-8Z2BNCc7RfNUHpD6)jc$qHF=i zF{sp1chvh^BsVrs*Pxna`He%=QSQ%moGb`$l|B{UaqUv;_glJk0}Idj98-IwJNGmw zNU+RwX5A}l&lCTt%&b2SqZY-`d4m`^y}Y4J238(NnVD%10FEp+S<7VNiVy*^8aqyD zu&@r(5;Yi-I_JVpfah=5MuVc$*)kGj#^tJUsHO*&(8ya0Hj}8szL+}tmW^1SQ}4X^ zX4eb1vQ==N@114E=&Vla$sVU1>|4HIj-C5Y3d(h1HhW!Vv?xaQIL)D2w|$FrI(`u`p-JYAt$l^-3%Op3- zHB_LQW@1CqO8YN=|1bOJ*Xj7>hbxW+?W~Ub|1UlMkkuLY|941O(R3zjxv9N{Ww!-% z|G)3(MoGB;zm|q)D^9BGEkSkqCykT7m1_XJ)0(tbNCT^l9IEM7CM0Y4GzApBPNwJ2 zNl0SE)&yEAM)kd@kfHNik!rSlkwM)5-)WoFI;Jh%Vzbii1Fz5~ePEZD zDQwj4TVqe<4`d#f7e_qsN@GEJ8|$V&>%+`;){&{@DR0U!Qtm!p&mmQ!IWjk1cMsqQ z?0H9xeA?}A%FJkLCk~6YBlCO6NTSan{XB*^b;VPO^2iBx8dUSzM#_}GSe=r{Yaf;Y zs-sLXbu^f8aonX!?&`BkHIA+x|BZ$E{ye(*a|W;c_tpk|-ECM6|2s|x``~_Bzmqbx z9O-(VLDhSazQw)pVZ(ZX^7H5o0M|F~FD2FNOVt#}8rEt;x$yI6(LvuzsvYK(7O+%| z>6EIGlT*VO`Lg2{0dmuC=QOB(^D_ceN9jM;acp$ZrqcJGnBe+NYvxnADqX*^s#ESC zT)#P9^V|;^*KZOHVwYgNx&K_haTm}suHVe&M%JrLP1g^_da0ckaPOCInkg}0|^<|ZzK+J{ie18|9{tSO8vQxqnVfX?zOOGyiDJS zVN+bfna>P^s%hMgd3sh|yZc!q7JTc*SMU1em|34(z2#Xg2g~#urgu!$4PJSjmUDwU z>6FDEsoOO&vFAPkGMBHt7%W^L2A^LjNaE~xKMktmTLj3YkM~T4EV$d16jR4#00Z%A zfrkRIK1vo-$MT5Csp%a*iK!;mfxFXwu#ki8_l|4)i&@67ANIEfuyT%z4m&tuJgrmL z1hgndx}A{8CJj1p(kmL^i69YcD4H{D&cpZe=w1&@`3mf-v5m+#+q zG{X1I{Z7#_q4>V}uIz)zzwv#uH$#i~zS)ICeBXR4k^SRh6f(CP{*Py1Hh84~@qP0g zWqjXksllSp$4?^xNu4^`fs>2&jFHGmu9+mhZ`K7^yJMCL5Z^Zka)>%?i>bqU#^2*A zyTvop{Nt`?%h<3m^X2_#mG|Hc3@>f}5@W#{x9&1Mpienwb77W(7R6bP@dD|;-(4al zuNr#JNj4_A#qfI^su`b8NRU-n`kZpn`Y4s6A)y~?u$c;`0 z2vGfGW1v9PVe{uY4(IGM`0?e_d{(N@2hr#KLpIT5$cPT*zi`hrt@qp4l34fF2B%JU zb7D4C>q}@+3~wAxEsGH~hYF;Ba$CyGbjxQBxuL^8QcbDwg=CG6FOgu;`tW@(BdIZK z7tm61=E;(b9JfQ?dNZ@mL6p^RgZ~l82Z3cJh&rr`sbg7O+b_{AzHvLdjNu8DkFil< z(^d?q?7%$K@8{PK!+DL3Nr}q%XlCB;txPQkeDG4ps2Ob}GStFDf!yr(D6?XcI$uB* zx6jnKx2!4il#q+o$3$xlNnO@^3NR2`mot?p*Vy9BAvdy#lcAU${!$>H_;>)S!|Km< z9Qtcz+kfnu+TwnUaf^;Pti}BpdtctVwwifFpRtQwI1Kk=bPkI8SsC|Z>{h7d(IsqU zi2E_#3uLg(a|z;pj4YXQnp-VSCf}+hKvq3`l?DszFe(v7NgQo^MR&%nZ!+b;wks&( zevI1^ndVfNk;$j^cLD0KET#^x{x{FRw2NV7+Q)us`nCq^Y7*p&WX(WS%&@_LQenQ6HnjVgnokoA=r7Jj|O{CMjd9w)DdY}viUE^O{~;h8+E?tI(+~CY5Tb5 zOML&o;EU&w;e8%5gpbmqz&0Sv_bw_+s7nmXr@g>{%3d2u9l>Fdf2iv7!6mdJYJ!YOmZxtRdl ztAKj~`DCR2sY4wmf3D*|ub(y_tizw8A49K#ak(eakD+GYG^f+(#}M1uuYWi6V~DTq zeu|?XLl>D^_71M5k@DtWfch~c$q@Y*&I?dYTX>0)R~)BF5d9cV39zsZ%ONK?k{G-3 zy#&=^1wRzBW2t>4`Y{~TDF4WmbQ2!u$jGR}_|J9h+jhn1!;@ww#ndw<9@T0+V~SMb4w+cS6Y#%FH zDlSNz!N~~w)&TvR>nzY9`Zrf#5dE9=R|$2@{Bs?bkC04y1x7F{nr3iC9vI(OEt?fo91VzMN!uD zS56&$zL|uqx^0Z3IN3tGWwY%DNGq z-%0pEyc6)0L2l*rScYoOjOrRh9n=3@$IcPe-jD2L_L!AAHExPUa!cmXrrELv*Y@)c z(aNm@rxr1v)@r|fi$*i^?iPB>ky8Q`v@A|e7^^$1xgDciaBG}I{&Eb|$gi2vlv#zG zd!Vom@z>Bg0+Kj=@J?E)N?xxda&U_W0u*z~Sn6g>9?vN!70d%v$Fx7!v8{_=uW~1R zt#LekcVCkoyI5x{w;uA|d%V5P<2z+iHZh-A+vL#v<~W{SBBMndPy1_7oElw?ky3B{ zoz4u_TIxA}q`5-Nza%68)c0f|hric{_^Ya!LQ6w}CX(u;3C{tJr{i)s6tkkAF!Eb; zCxuMu)1O1s5m8JXbtA)5st-NMN_kX`TT|v3bAMZM*4~=lyu;mqm|hp}F`xQ*alS9WX70812a}@GdL`mJWa$G}ZRUWmbTvW$~A+%JSS&ji* z4?NM0lYQ+XWGLp$uBVV0(o=x-@!&XR)G_tXb!?4noz=Y7+orfLcv8aK3UzQ_@N04R zky_lNPwIn)OZqU+?U&t0#@EDs!LMa%Ip9zQjhwkSNhX8l`bm&k3^CE5nrAv4;Jn|0 zni6F8mWvcv`1vrtGlL_keOwhr%1zd4q<{CJB<>4V8RaZTQ;B?dv8Mo0$G^-YQ#3z@ z68!b0fz)BF`1$nT5s}-yao=O@cdcgp%}Oq=XzT9zg_XV1FvNR91KjsmmRiJpkKZJ+ zU7flz*86VI5}8yaXgQe73Vp&hCCtAp>h|*7LQv zN)j^edz{XonwIS)L)7svv&j6Xjy7Ja#y&qW0_T+;oSk?(4(FBL#rR+}{`LzlM0Xq& zj`K=yqw9J!$9bjh)H2R1xlqe$*uK&l**{J{a>03}jsh~yE6q}9d6xTO$~do7O@f7W zSp2wvfb&XIN-N|dB~N#fYfFZV^GdeB0Oyr@OAzOkY6}o`O#1WtYKvRdUF$1H&ro07 zd{81f)QE+Sc72qyt{bl$qBwVcGM?3Fx7B-E=Okw1%cw;$vOy~aszW!-qpX?^I>X7A zH5tHm_sVAh;{KT13SjS_Rr%DA*j=B@8kw0JV79h@N)S#@UPw!2cHQ1BGVbS{V7bPSy3^z|h zHOBE3Wi_mUg#yJ7A0uSQ&8##Hs%f*%D-d-QGk?tH>F=Y$D^ zxDgC6W5?lmiTv9AtOl8}aW26A!}_J3P3o?bkbl=9{uEl%QgK?PWsDrSM*s$D#|swwz2X^&U_iVSoO7Z2D04=~-DQHohnrg~P+2FpO$EWm@ka4})mqRi8Wk*In@0O=PX1S{s<-$7h zXTPB&hPFN|lM_3i*Yy^XI5Ms`$CKhG<+?^*$vFt9j?h2Xv2jDc#((vmT#?zgI-EIj z)+mfCx%|OCtv$w-Tr%J9eP`C`%EL2(pVnerNqd=E_KlfFEsL}M8NS|Lf;xwpB{`Mb{tCG%))M*B?7<8QFgjU4IN1vcZ?GKPqSu z*B=K1y8f8RN#Dw~DC7EL2aPna+8{t&e>9dMu0MX^u&9pns(>U$Y+J06QGLI0GIV}R z3F7+WWd?EmQPgv~&od1fbrdsy%!YM;Cr@~?=8F2ZOuMNA%8g*5Jq-DZv(tEmte$lo zH>WbEWzM&a+iYMqy`BhYQH(5+$DumX{g6i18+KJ96F1$IA-9~Jr0487Uxv(F`oEB2 zVI5YJycHyI!s|B@x$yXCPDY&UtdY+K4(E`Y7yRB|W{uq#=g>Svr<)g z1xKqc>cFjR)Y-0IvlFbCM~>x}QXdRia0#=LJ2Aej*_PMTqB^z14S|$hU&~}j$qEc| zoS55k--cz(?vBfoJ}xC$c;2NM4R0^xp4KTH+tlvOy4M}q=ac^j zX8oy!LM?}Gj;EHz>E^8&8EE$m;0#@Sss!;3eh()TkLddiR_KTBDahkIo%Y zLD$W{G@^m%vTYk>X527Ef@-?|L`L3Pmaafs{Q2*G|Bnv(^LAuysg^S`22{JkEBFuG z7IfR0jhgW0deWvA%tMTyP|A2d3;H|a>f{3$-|W~=4bN7c{QC%3l~;Rk$|0SY1iA4K ze^ShRuu&u9J@n%?W>;p4j9j$7p9e8|)NOn}PENnPLC>T6uac3~ywbZl6u%zXqICNA8#f0ye$7V@Ui ztA;so;@uWk>PJjAWDdJZ-X9TEo7wDJLoK5pk^cIs)wi`mL_eb5KNRxX`2!51AJG5} zs@Xo?5@Z$)BFUmUw)Gc~)X2*Eae*4()SFZ8cv@pn4D>!HLvCywMp;d>s7e`i6f-}~ z8mEE}FIz@@<#z4nT{hpWu~B2DzH2#Z3iH@h>6pP`Ulz1%Y6@?DhMD*BK!Z4fpHxL7 zqbAzvhMKn^x!I9m4#ni&i4wW{L=0utv}K3_3+pg#by`MJm)YuDBe7de3q*r$jm@c`$ieFsxU9mUK~vua+{7|*D-?ySO@S#piGms$Nit1NmmVHq#9_oX?P%g=rrdiYYbpJ*N0(dx${G&9>$% zu&@q8;rGAtSviYisXDFN1bS5UKUG?W%phi%0L7;U2Puosr{5t_M=|r$tYQJ1YhS+D zftMTUUDxt>w1HfbKoT-Hp~U{M{D=P@L8_N*QX6eCu=*Y!C*K|)r4{A?gYG1GH6 z!1bHHT_uP*ikY8gWkR>rbNZ*RR8!jdPW*f4dluTH;qfLj=5f1$M{DID7{!EBg0J=S z7tE%81gDne%n2G)eFitAtTx?rT_%%n=-)`VW%o7;|stdCfS8MSXPbA$ejS`ItZms(aMPdyaKkh{-i$W3DhQ5L_3UEt)k z*NGZrF5^2gShPMyACQsM+54+WP@Kl@0bEbZn!%y^<5fq6a%NL~?<_uQs_zw1M=|rq ztUU7IS+2?U{y3g){G(CJi#VSCmN))pCmc_^82XLwhU00Uo-aFv((yE--!wJVR#!H- zD^VVwqHl3A{nliU4qv$z1?k7iPwX0QSuyB2tHMuGv zsT0S4=TMAv-l379H?HVD?ehx}=P~M7NDy~PCu~-_YcTteJ zdF2n3apdAhvWE8l3M{O{LNTVL;>3+1wG@>qY?_4*gK3w9nzma=gXIEY;L}7%tld-=~!mD&uiQltdux$xFWaf zlk)1o_Z!S<{pDe=h6UojU`v5oj_msm;J)BFw-m}Pg~^Qak&9*k?<=>LoO1S#rj(h3 z!*l@_t`C#^I0;Fe_4paU`3!?tPWt*T5y*SX-!aII8{W|%?h7^oR7WxM$1DqQem=F1 z6~>3%*R9#9N4r?)GiR$KJNNS1KWt0XO`XVG)JN-P4VuPmj{SZo&-igjLd#*^|y zn;gj{)h`)a0o75={4ud#npPj^GWtHvV=#2u#DN&kCo%JHlX#lPAY%V`KbptjWL(}M zn#X`q%Mi+gZT;y7yb( zi(?)Ui_3*SpUuJ>-re?YbO&Zx^P)m6``qgTP{+GYj0|?@p+Rom)I))4nws4GBWBI`sdF_?Cgt2 z$G6)wj$3|yb>8c64D)XAt8?O&-pu*ksO?ocMlruRWU*IxJp8%WBBt zD1r2?y;32+Nv;esqd_qm6hC!5Dp9Tyd`p18>kxn5_^Ba@V^B?; zFj1zQr!1u`CLUmvQAaWJ$1J&WxWttBjC!n`f9bt;KTELiNrPRvu*}r0Ek{PVJ3h@2IAe^%LaOEcnqK;zbk6HA2N9On^rps7aqfKL) z%?xMZ*@G{1%g*GTb9nEePcO0VHbGrhE{kQ>?;lc&V(5f>y3=deFv@{72awF{eF4Bb zBkGE7VnjKHA6N|+BMlbTVe+#kM^fjUP!+PiN-qD0X z=1{3XfJN&&|G0#t&WhK!Mq+6GK?NE6rxpTvZ~0L``=^_f)wGy=8Df1GGd_(^gFO#d zE$g-x^Bxae_qrF(zYTH?>)ksL^B%JyQ+@hl-s9f;TZedK-s2t|Ez-Ql)H3Eh4pB&_ z7d;sH#-Nh~G4F9pnezY4dwfAaE?VDLXVRmh|5#Ii>;E5{N@Rq-Il<9umVW$$d5<-Q z<(T)_7f>C=j87w73;S_$Gz2f77+o?KU>#o>5sV&euG~x^E1)R;|=4 zW@B-jS`;VPG~?tn)7JtS5Uw9DFtf(DXgN3NcT9q6(z2=~b8fVW!NNMsF6-}GHRj@A zj-}$bjQRAaII`_^%FN_(Z^~+V#tnhI+fJa2I*J*eMmkgRZBBapVlnelJ+;*FHW)A1 z&AZN|Z_GHabZl_j3#{VJua9S4^Tc?;6BV?G@q$M&a&Sz2nHkxlx{~+LO)XsODj2Qk_BS##$sz5by zzzt40?`$i|V&b{Jlu<`9x0w`&8q!N+Oe4sRdCg1W6YOI7Sy`#9T9IiDu8 ziH%;1Ml*3g>D*rewLEjcDS+djI(dv7u>TDw^Nf25kednm_EF8M?v04`5xY!-h3mt@ zeUF5s&e&-rlMAZZ>V{X+{}wlhOVeoiuc$aqer;8jGU_O1e43#vcU&u7^Shs@hCIJ~ zY1QA@UL8K(u=_NOvvc_B*QiL0v$JdR7Jin-*TciUPUSU_TkE+1&~>$n4*D;jpNVL%rhA*84FivgR z`x5!r!u>Mj#_r7&%4u~C1&BI|89!#o&q;L#1U0$M%6VRpqRufE?(;$GSa%ioFizdk zW5Z15S=lns&o+=*#{_ePL>a$KPtcv;_% zF!Lp~Wn}edmo|*l2E_`Hm8mDnu&@r}#5o*Eo#n7jBPWz!C6PnkPb8V)k~acWGdvn{ z^3Kj5lu<{mKOavIUK46F@9C#x^{!2gl+?|qaQ~R=U)8U^;gvm3&zaWP3in@5{<6|% zIJ23YrJzMIvco$Kssr}@exFsXaEw!a@M{LiE%SX9sAjA75@hvS3NkFL!z}SWM^eYd z<^uyUx|X{_hI?<*^Yh3rBsV#CmqRtfbvTJSiWxs<@PgAzw>5auj+JR!dZESNW|*h# z@j!3S9PT-(?b3$3ELr#U_tN)_P?(KLc{GSm=Dq0}In8y3LI%XX1UL@ukwtO?=aC9j zlRD_nkkxPfKte96V@hd_9*tSKkQDtZm6FI|o{we7Ocv>%J~jPHfP$>vks~-n9mR|v zGid3>_4l9jF%mDXRCdS#n0K|MSZ_%3TVCd>T5V5{zc3GV_HgYi`kOCL)X*YkSaYP7 z)h^Gj%9L9kPb1Ys0~!N-&Rd5t%Guk_a>y#?_mE)G=Obei0zMaZHXMqfUI7w0dcy`r z-g6kjAx7s45}=wkGlfJQ#f%>_@JuRyIN)KSjgWm`c{j}{BSXhVU zBx?lp8ymy~a*53dMoymIOM@6EdWo*V^LvP>W70YaqK;z5j~N)odp9m|q7N%m-)KqX zv*9d!$ihws2cX|o^Ny=F9)8F=4{f~sQL_OXnEuU9+@n0?HAN*C7*@P~T(V`gHhHI#K&_QhTR$GbJm z1z1>zdANgwr22gH(x5oEfe6s6z@{rBUu>)>LvGqST7zoF;r9|m9mR|v(@%+7{rHU8 z7K~%)ReMm7O7mT}tXMma=DQy5wE8uTV^}Za;wzf(T0x5#$FMhp>I~s`J-&T!%gC7T zx)G`Vx-e5B(?0ehS)+Y<92VALW_DLW!hF|tG%~7uIX#E=wv&-D-*r_6F^=I}J*Thv zBp{=XV#bf@^K{(q{a$~+!TcCu9~XFFKe#u$dg;{Fm>*+Trekb>)-|+dsKvr$%#Tq; zL5pJW%jp1TgpPdCb5IHWH)DQ`+8UOtKb*=-WJ2UT1&RfwS8-Tahw-ke0+KqZzlZKr z-o^{%Hr_uW{Cd|6VW`-6;sYii8%43%3Ie29?37J_&W@}K*Xqh08x11`G z7=QdNS@gcjw2_g-1@rV*Lk;URRzX(1pDtvO8TWG5^N#_}8nXH-VFdVp@cloknDIBh zW+e7@JyyRwt8ld4ulO^=+1TW8+e(8svR)Q@E}nT%f`z>gDE};aEi+otPQf!4LzefK zpc<7mSSBYJT2sb6R5umMX`Zecd4G(Z09m8{=?XN=Jed{RNz8HBfq;K)(n^40|D-yk zT5ZcofcNeE_X_!_%UKDEO0UM0Bi*Bl8GqyJolm>B{p=Vg{@na+#h3P9S>U0=^R3k! z=KjL$@$nl~*~A7z=Z<)`gBgcR(_4=63E*g1o#5A+k$(3hbu(g(Wy-$}R00?=cjHr% zHOZeW!@@d#)Y(l-ReAakTB=Uv({#t3sjg7Y>_0^#)8gGZ`Am)jR7Wx6Z+x9#cT1kS zejKY>i!Y2@_kj8L*pj((Wg6>#Iep0cI=5MPm%C+M*f3^P+dxB$;sC$u)Uq0#`bs7z z{8L_m+;G38L|NPV7#OIV!>eeJ)p*mC!@@eUTrCiAoHOQ~?&P2Pc3YIx{Tddi-^TjL zkG4GycZs!;@B&U2tt7a~SKlqSfdl%zx`BUh3y9)@_+(zn5!bS$JsM zWv#5UnBkL2ddroCfft_~!{OWmr&QGfftISHA5Z2`oK)$$Mh>%DB|-JujPU~bLK&?=kySdO zK-5vp_!y}l&#a7Xb7n28+&#$p&52;M@i9_w#1B3pl>|9$zi{i}JO(m#K>-S6{$Njv*$h?N7HS&i;DTz#~6Uw3Z zqjr`I3+u>k@%#95;zEU%iXrWf>voAT){xatY0i|@G-D5ed~9ElL)1~s_!z0p-rXI( zcFiDO#pUN${|1?C!r`;`$~F4PJsRA(*?4mvYnf83>zJFZnUNYGqeU^~pICr*c+c;G zp6makKxVLiA&2Ux0ow%f?V!^D_ouZSO}Vg++}*N*Bu+R|l0h}%Ot?mlofRxVG2QEm zod&@g5{{AnDH@EJQwX5Kk$7c@!^IoC!M`A zS(j_v(PAV&?a&&hWS6$ZIs zllcPW&wlSn*3v&og1_rvMg}TFQv1LCOC}d*EddxW(;aIR=344&PtD!`Rid1fS4)Ga zqnPnAl0UT_xS;+_Lsqq3$W*U7ADF*&$p`Cw$FUyc&x8)y=+DCKKP10coy?3{T@%ov zIDB)m1l8z011YN$PTyvb8!ooi^V{V$bdO$bt6@2-kU4vl4kYf_z)k7~En zpt!;0y97~3G2>$-8w|g)AkMUb`fJCQg_Hd!x={FoQ7tl~H5t$D(0;Zw~R1NZ$+T03%c98!O~bG30kX z{(v%bzTnGYVI3KMjTw?S`R#U%lm=~8$dQX&IGHSe)5z2#Qw1_UG><{lQOx)lTHd_2 zrX43nbIWsgudXfdWT7Rc*{l!jL`Gr`u${N94~q>sHX8?4w>DNdKxUOBgaqw*Xq<5e=NjRGyi2WIh;q5;;~+D z^qiSwtRbsO--jy@brdr`#_Lj#<_>$=e>p2NGwAKYX$35Jy>sRC<({nFhmzUDvfr~w zHgDbg)*8f&UQL$KqB!};egUfT!}^rf(9`7=vY>*AZq}B!8s(QG>}ANBJo$@qVI6tj zA0S{qc)Yhj&f3+Ek)zMdmZAD-P+!Vw=BZ5@`TSNd5_J?aK1M?A{U@(~ynCIM-rx7H zqZMYb;08ylU8#P6u?PLv8)UC!6W^7z`>T2+GcI*NLW|;bzj`uMXWZNcaJ^{s0Z#s? z*+n7qnzshH;$y9UV`NVDjRfSvI$13+7picV*~q!8a7+Ps~R=WEdR0rtLkRne`n=A z%=o;ij26W?ZmlJ~yOggPZc=A1V@-$8dwJ_EGq!hB(6Tsn>^ogKrMYf^yZ(+4^GBpgcwEiQ z_R@X-WDz~i%5H5#7Ojt3Ndl5O*>EmEKa34)^gKk_ON!|)^{up;l)hg=R$tH5_baye zxs#k#_kW*V_SF{$KitEuU3j@I6O7rA6VL9fY3s}CE-v-(#nlkzcYa*b;xVbr zG<^!k?^kDkYD-T;4fD5WgS}VT}gG)G`MF`Ze|$Dbo-3`^Ic4Si4j# z=%R1`M6Za?v{;OED5H>L?V|*!zL{e!Q~uHMAKkZSs%glmVIB^K1d z2+g+9Zs*q12+h}1UT0y1=A2I-f+I0PbDP!a(HL>@bUiePZ^i-3Br<4Zh)j+gXs3~R znVSLnDt-K>K#b5FsX>g;yhegWl`ykPRVb^&!u&WH?R8oqgQv8XA;xGIKiy89Q0tbW~v4F0DQ^N4{OqB_d0tU%7o{Ho{B z{T_(80{bOEhTL@9c~VWkA55Z>e>oFYRLS6!Js0MhPfcTHzVi*1G`zz?GbV0sRXc&# zJE4Y5yq?Egm$qCsXBcK`I@?b{i{i`+c@h~Y%mz4yf1&6(&rp9YxrI}0%4){exf*2E zhL)rJpGr&@&!?;g-HDOOC>sk#`eSqfD1PvrO>*O|(*&q~dNoOgsN`SHg%wpYsLnZ; zwKmQ6pwFA5s10~Yw zh5mlV*U(SXVRq|r@hOBGlmidXNXJ4B!x$OzePEtj67&w5J9n!Fc&C*2X z6u;Fw;5J4lT{=`ii|EzkuRzu7@mQI1^MW;`dg4Jiz`LutUIAt&$)(IHvz0RaP*jO& z+4~$(o%u0GhT??s19iuIOyQL8PrfKXZp_BZ5c8AVql`-a<$PFCC0SESMQ<4N>olt{ z&276j>JY1YGT+|a$dY$8yw{{_y?t!hkeNf8yt8E%`6a1EHRMo94T_QG(E$Axrg}2S z%wO3^l;2G_DwBz|N;1ew_wG#orxHVFNkbHac5@l3Q(fvwq+gjQ95REQrxYm0m#U!q zY+gJ=MkW7pMy#lky@OhXTpu&{FfS)P95c|r1l=+AxU{W5?if27R?|zv9ZD^l1;=5e z6{Fvkge}*oWi|3#zHZ3Xezw+H2(hR5Ha$K~!fi86!h6 z;?OgV^f@}0vic*xE-7YqvDC;%{X7+jO8(`XSWzYRy(>GvsJZYd`uvaTDEXbIIed5J zJFP>X|18^q;t`s|_sru863yYu&?5T$zgEadIZGljhcAQZ^WTa@pZ_lkdAs2*%9z7f z1^-it!QM3#F^lgMQk*g9ra+D{-l{?L`X9|Gf9n}5kPn_+kRU4gm$PC;mF$V!+2&ce ztA%)Nv!^Wb#Sb}TjlWco;D6SMQGRdA>Vh3*85F}_uLkH9G5s`Uwcxy$LOE-# zzIDd9W+!MlD*2c5VnvnQ$@%+W|L1Xvn3y2Ck4*{1c=?XQB1>NqpEzxu)pg5Q<`9&7 z;&y0VjPFri{|%Ad^3gD^6e{avD0;-iaPs;bRfD+CbrEHZmmf%(Sr!~t;D0JH+4WsQ zROL^l1=8=3H6uGUwAURS^;kjX#unuTsHRz&P(~&Ha%QZklHI49vWl@*ISZe ze%T`u`Q>1e44K)cB?>I6MEvFUK|&M*3#yXpl>RRjGC+C22SYUD6R zKE=yiu8`zEyb^j{IoKJDJA__W?^A7PXawRROCdR;|Gl(V-T zrOYa3y_eztRFcu3BdT)6dJ+^vtHl7ES6V$^BkxscDMR$TYEGix)h?i4|F3AU#-QY= z5|>NVvW>C%Xw$>`u(h4lM9U|ZZWbHN^UX(?T{5d-I@jdANmt_o#wCqH4Idgd!$0|N z{ufzkws+Z*#`uP0roW<9%`Nb=X?6b!5gKzSvB{-&Nw;=SoN6g}63Ck~46i!9&BLqFAm$5B=kG+ zKhG{?=lVm}=Uy#pkRczR(o0}=E)8Y;!@_6h+Pt}ILzfnAZcXv*2LGR)U6^$zKi=Z9 zZ|MI=+(_K&BPlxKN>ZU6}x<^o!{7?vm zu|}FdtRK!(sSUMX6~c%;!rD~S;}PP_x=F!u;!xMa0B{#a-?UM!Has^Rx=rCOKFI&u z^?o+;fG``n|Ig`jej%peU>rd1j_u^)=mP(70e8O2+(lSu4n}5^#t#^}u%F=S27P*~ zps&U?J7zxp$l4t^+17S$=7Au1fdH{S(oGFwee7m2tcch;9164cWaOzne?3*W3li&8 z40C$Y5fJpN?3!K3XVV4*1!V4+JPy~<#f9aunU~=SYne%tjxD zIeXEzvH@v6 zaO1-%To~(}FEm&`w?w?iIlm;!R;|Y6z`9GP5!=0ig5@3V*n}EhAARZ7NBC>_R)ss! zgsHo=l0xJG5$OH)Ju?5toeP)m;>hGU^YCdmS(&>C-M~2%cf-70zYW3#|x`OK~MnhR} zH&LbVTgmJ>ug<3(@3p(trs2t_Tiv~z5POyu1C|q8w{DdH5j!PQidp-qrf68y+P}GT zZ+!SZVTe2+ig9neSK;7;hDM#Bpg@N^7dGrWE)Tn#pv;}f@KZk3HF$wXr?Too=(PCU3|$18DXFVda6hvF`ra@C)uE_m2ppIHBDcpi!;0>bmd5eucTX4TCOC|2#w<5P?3M z^g@c*pFhF+@&oJ4VeuTDSbQg3p^R1LF2%Il6~x{%$IN9zzTk2DX@|A#Ys&k8`ihF| zCvB9)<&=lxXJ+do*6bY=1Sj_Xq$?K0T4F8)jcOMx4XPi)P`L9_*%WiWc40?A&@URV zh9FMzQyf5E7qeVAj!YLPKE7nfD03I4+mt|YUN}V`n|w{4Drje$&wOe9!)tcy)_lyi zC_pC{^*K5au|K+tAvm$)zVm`1uuiLk;{4Bxxe)y%+d3F`C{Lo$X7K3kX z7;iYQ8gUpI3~OC3u|{7Qo$3el-luRUD$LC8%GoRrh(K=@&z$%_`_6^SVmUIo zd=@_MMk#X_VYp%$Wg7SRWSAk5g6^<$>`SLCwgBrYJ|cYH3>k>$+f z@WI_krS2?x)*w?$=Jg*3T?99+LM3=>rj|~fj%RCHlZMPm1+5L{32K2HYsO8)IzTs! ziEX)TOHN=!TweJ5G}>RSf|@U%74{zszFi;`|!3mUB_pbBE+3mcu6O zO5u*b2qf&gvzdjx#5)x30^cvrW@_no1O&a~dScQ4y5r+pg|jmgcAbwBa}E)Ddxuk; zK?f%!4zU(Al>U;y#y}tNADT$&+YT4!)kVxU-IL$*b}io>aag~w8RAl0t`vwk(I;uv zS`Q0);T@$g7q@BH+@$Hb@_-2R^jW*-B8J8`F%*;^SPmQNVkb5jso@O#Suc_2<@I}Px9W93C#M;2o^m;7D3JNs( zELmFHgB_vpHoZ+0&)KWm5fF5g--IwU>V_vS(BX{>-dJ27KCWKM+=W@(5~FUSLxbLqwA5F)&_ZlVL-`;91M5I)@}##MzGn^o@-3G{9t_NZGpxZ#5A0$U;0FnSZ}@$l5wUP}tsj0)097 z;Dmc7AwMdztb2_%JXQQAx#*2HV)g1F#qQE&jdqs(PDHya=rM3*KhM9pyWxHO)wD|; z0RwM%pX{6uhvuJ+cD}O^tomesq(;v~+>+_A5zifMNr4{ba))`M>oK?69HAc37C{AZNsO8=G!7;lAMBSSbSxV7>1`YQy&x7bedqKWFP2J-Bwt za?Az9dJrpt<-~@>i$c5}E8sJQS)1yW-@uDk>@d#=lcdqhTRH-Q4*s~*A9=gIzy)Q8 zlgVN+9hpu%7QPgQD07!$aaj{}{n~TWh0Du2ZZ)9Kf4+aU%SxyBNPAsyr!Z%yWyhpe zh}B;45;=!EA?%8^_9fmJs=eS)y}sw)S|x$*EY&rm_#FX3YZ4Ai5fa72`W^20E?h@v zCes<;_mMJp5r&saD07?0p4{Tkf|^ZQ6ch#>F?Ror-4S|4Ax1&J>N}bCdmhj&=M7>H zYH!hyPMJo>|_Fdu)uN?_PYvlnE==FSqCe-&I7pgNKx`@tPHWycN zgi?2w`4z}GX`R2;Uk?S}9UR1XYsN|2?kpN#|H#biu$yhkH0@t|y!WfE8iUx?O%yCA zuKwjI0U}O63!1ezUI+QiPgJ-I`aCD_sas4(K+p~~Fn@@}C;It??8v$}u{e%=o)g`@ z7ASLCijY2EtKh9+y>Fz>g!tICbIi&}77ps}jO;Vs{NaM1k0I zMo^&9m)>D>yq3*ZxD&aUZE{gtDG!K1@9Mm$1`Qr}8W)i7J06VcFkPTF0e5an-D#+A z#y3;lJCV}`MSfeQcx#}CWu!XWq%U~Czi=NVJyZfdVRRL zzFspNI(+i+y6J6m^9GLG&=C;yiF>D$(8Paq!kjqt&6VSXuL@UX?!xpuJt@w-_vK!C zljL`5tn(Y1`CrziUYa&hNT3%VqafRvvz2>HhW>1!d>FGl!aZ6 zZITCsuhOP0+PTLMv7-XSP*CoU&vk}IlM@f$zWB=ArRYUOQXK#CySof~2p+W`S*&yW zC;t|G(xi-3$CYW^cew8|>CIjkiESSeJ1*OI9CHzZjiFi#4;$O{1qyd@>uzT2&X&jn zBGB72UJ|JP|Li*^2kJW>ecthuxfAIPJ4A6()_gEOHd}D-{W%J@270a4%pT^qbB%`8 zPfMxO*pP2AZbWP@Vt>lQZiwBo1}D4ygS>ELRYo+{{cSWcZs!!sIQ|%UKsd5f%-CKn z$n_r$PG@{%ot*g2xChHs<}Sser!(|xa@lRG?sbfFfZm-mq@a7*@i_~&g%8)zT5&z) z^3qVml7KgS-pQMs#RnOnp`X*|8uO;FV&4Tt7A#J-dDjsT^ukHMT9Dg6nw(7B&tYsOXDs~#0q9RY*fbSvLa&O+n^`OU73Kn2 zj{D**<2wR^&MTd=<$wA)Ouh@W*10S;zGPjLy3-9z#v`wup{_dxFOBw5ur;jrkDO+A zw0wM*6xWLawbkDz#n-8uAgpt@iI*0Zu+`-7J-4p*Dc|bf8 z+kIqWHR|QO20;ONh0;C}>%zv>9aEV*k?u!(it{e-VGwy$ewJW;;M}tlt=C3B5B7es zR`uX$4GZ;YV-Z_8L=3?R`GI9meT_8?>kPj<>+WxbySR_PyVxah2$%H2y>XdN1?jz^1-7vsAaK{We?YwF_~y z_DjKy;^9fRV60(SN9Sym>))z7QDNHNnb-O9fH16+c50lN{{K8Yp@(%`oM0Rq+&L+A zr{GT!*X_XW^?$sc{fOz zJCVh$lc-xlLF@B)Xt-~`4P~JZn{d8+R@#F*<9m%>MUrkRBTHwfO#A{3_xm+sSWaXw z?N6`AdL$OY|3Du$uRm&zR2wB01YH(WRDwKeFXI9-cT7I+!*ZPI^UhJJ zJFTZl$oA=++A+Nz397_v5#E~ZQ&Id90k2!Kt!eAy13AA#<}VW5zw!=o1Mf*;IdNq< zmtK!Ms)z!OK9}|!Z*@m|qO+4zk~ceH7V@4oSnRlZc}GC+S^R7HD@>7fC0_ysG!rJw z`yt>sJ3&7TKe@w|@jzwnBFuMaQ5sVwmgFC|71T5Gg-{r@&zddYGR8eRbEfdO?ZCFN z{nI5G`?U}^sZb2diEDgc1`u(%-crok!^Ll}^n0pshen5;Z*vpL10v81avNSEi}9`$ z6p*=N@|e)+VR7gt$^d2VL?&n|HTXx_#U)cd3BDE%!G%FTFrG7hSizK|ODbY79{twz zXgF$3ABni$^!Y|yHHnqLdWh?Pox-e*?VB@aTa?0`Xhp(@7YFp@0TJlB3w9%s<++<; zC?I#|!h``kCnp!Y@A@lqmtr1oL8&@Omu@`mD`+w}0%hT^0UTQUT)W19UiN-~VCbY5 z-`1J=YLyCI5a+>h1}rCTy2nX?h^u)}j9Gh7rxUcG6z+xzOZ0o>pO6PcpnK~%^+d+4 zM}$y7?vBfceFr-b*LVGtxf7YV9i|3XhAvwB?1|u4=r~*$bcoll^21hBO~zJc8&&zi*>ClPpU6@_w6_04KxtKa^n6?dJG%!#*amK{V(0rwy?XbRP3+_uY)DM zd%o@n2>Nk&-Zj)Yt-lls$aa_u7oJKu@tomtlXXTPW$r{Kz8aM4hSS|olphy-t>)pa zfj(u-GWj7`to!|rs zAQC0(QFr@y{-0`u>B7XrIWPnZ?s_Y87hxW`no=3Hr~3BkiGqfG^r!>qqbXB*uc(<4 zb*0BG_g%Hy&eav@Y&wZJm6j6huBd{3O$c%I>3xZ658U8?d^t_A?}n|vA7M(~lLthg z)u%LGL0z)wD_nlEarsWrm*CRRn|dj8C$clxhK!zHPwkzc*D-Qs^t|}=*l8TcSA{P#CmdmlM;Df1Z4?G5KB4 zu4k-4nd z)P9422Vro7fM>Ayyz8OVon2BMvKkp6UGknysJB1n1ReKoT~umSc+j}|?{}V6Oqv;1 zkbkQm;vPIN1UJNkH^U{^6;HibjQOQ&o^+JvdWE})yo=U>X^%PrM&w;6T_f3xI<@r` zLIL@{j&@}9d1zG}o%P#83mzixr{y&mqbhFRYqQEupd z+tP6BRK#hck0^2VNoNG>A?~y`24?M%JocodT!lL#MSlV1@mn4cfliuK%tjX1uHl06 z^Nz>Fw|*uU5B+yn=1ydOz@Abw+Ph@QyBUH8`@R$u27NKW|I=B|LxWAD-0ze=8d2zb z!{(Ga;yRj%VL5T-;CCqy+&!0I)*jg#pR~IB-`)Mb7WeX;JRrFHmGI)lGh}^o3NH8$ zcP>00zGm%}xeL>(%SMb7>^_eRUU#%3KyxOX$bP%vU9L3rNn2ULP0s;d?k5m$il+pE z6IZl)NP&pEC5^)Te5+E`Xg5#cZn$XC=u54Q9RWeVo$PIoY)*N?`f_*9PA<@>W46DV zz){p95Ze(wkfT2gRR4gQjQB7y?4BkP2pb28aH=)+oP*eP=-GV$p{ zxjF|&R^%QOw7hP{TLT@Y64Fz5`KYefuGa7Wal2QvPovv{^N4fzH4aX6)F@-yui4tz zXp=!lVifK~85;|?4Xcy~M4)#!>>Bca)F@1N6N}4bGNFfUtHhmsgw_fUvbjthJ>mOM z@bsf`^r^;>(Du+di^GF8D0L>v1g*z5eIx`}}Hy z7z%fzm~dRpbUFfpK7OdK9N9LM;({If4r*g(J`Z2AU6q(~h%i5|MYfe6%PqFn6x2u? z#84LWmSTga!?ihjqIog?@kIv&Q~g+$yAe0^171#Clf9b2UFwE4sy$k`dgRi}3Uh+M zp_!{j`^y6&LEwFn-039xAp`<}S>9VG*UiuISL#96LeFt?zhi zpwDN9$8Xpl`e?sR!qmfm%^$_|+O>8z;?}lTD1~dLGq5X8V<7(8Ogus;S92RP-vwI*_0d}p z6voMV;H`ts(oa58Q6~&5 zJD9y=*O6)84-M-x5>B(X?S0rH*()*J`XqY&FIJ%}s|po`yGYM8pEsp}9RWdaIT0`y zb&ZI%awzQnGTewTJ5qtT8C z3dqkpx-a3(XX5J4T$#HF^Km`MR?U9J6W@0dH2-kIg+W^#%6_!7v3s!8fB*gwxo6g^ zTAVy5U*F~AT6!URY$X4CttT~(C3mT6VL0NF8>winfxVPk({l!Pw-y3sQ z6*kSX8ovT@EjNf^IZ+zy#{hymAFRRMe&*|O0e^Gn`{uTx@h5pemhP)-t z4bPxkv6GEpIdRXKzM&!RyL(upFE$xPRw+Fd?jl}p)KPUx?+6(2a>JGgb`;X@wpRiL zWbRzJ@Y0Dhv_A3U9aClQ!VHSM$Oq|L_dG27E%=hHE{4LOT_!N1vK=iBm7bnFc-2-t z^NJbEwzneQzEB}7C$4+$rhth1lpLF( zq8P_nAEMIl{|6Zq8T8g)u2-WsTKAZ)$v@uzzoO-!xkTf=hL3u_y1rVR>H}5Q|Lynx ztul<1>@@o@ld+d5*FdeD2QC6Zom&;Shg(>n*6VRe9#M&>cEzT%duJUEFbt1C(5|Z% zQLvo2B754&_3R`D4Hey?zQyzQD|VVFRjt@>ho(Fr0-aY{5Q3~_PdjC|e-|da5bI2L zCk&OCa|knEIg`>L_pn5FCJDZ_H#$Kd)LQT?qpGPTA@1;((JSt?)!5m$E<@anDN<~X zg8#rjaQDAUv-VJ?&Gyi53Ufodr~ma^(@`D}fxc(H<1w+ijjrrIEAY=1TZx_X1RueQI+)ORtiV@mj^_kyY%xCAfx;(cmwil6*^(i*`eQ0 z&{Ja0L8P@J99qD^Yjejzi+K4sgg0h$Ipn4*%d<*wGk$u6*nClgFK_*!s(i$0p?j~y zW!dZhT#N2Py#D9w`|TgOy;7JP#q{woKkn5L5cJe(TMr;R+3Wx9yUvM^Z(vM#Fb36} zt}=HKTFHA**P&Ty2j0TyNibbUf?mSV(=xJqO+^o2y}Gw(@cyXa@}ok;Wq=nru$RBl zBTz6B)qbn^TvW2E=jVHWbN6KXx+^m@Is%S*vc1e$=QXkuEk_VQw(H=&3Ka>zy|2?* zsXKGq4a6gO7yirdL1?vMvWp57ZY^{P)IOW6zd^3#qDL&R(s`mv6VGIt`o z^fkz8^P%$2sZ2qAeX;~^&1&=hw2cOH_m|ETjMLh=Zcgce-AN-BPC+~!y5~=na7R*j zJ)T8RA?6p`VQ$MCJ1g8pEKk|7f_&Q%5cDR;@wUk16x}J8?K>CfgTZT)&h%S5os_!M zN-jfPhkRS#`^@-`VJy(jI%UFXHLH(vdKz|4a*ezm74!~nryg{Zgt*ZCo)o*{^0md7 zpNnozH14LSa5w6Sq4sm~QAfa0PYho6eYX|a3$IC_z~7Io3zJDV*tM0qv%Ou1G*Ws$ zezh}G&`{q5W#KdBDV=Ml{vdg;O4c5%mh6k`Ls)HewQ@u}*&9!{Iyp1&dc4`TG{1CN z@pRgazco4|&*hayP|quRW~ zEXVS{6#LH8X|&6lc}wI0k*CwB+*OUC$V~P~On%;R;7-AVS0}-prZRUTlgyEn`-vMp zif8TY(SfhcCOVFfFbpqy?U3~xaVt!O*d3KY6aEq+oub`|Ld$Pv z|83MxWOiTu}*q#f{k2?8d0ph&dFNNjAMcLaETs02_8U`r(?n`ETQMeoC)n9{-~8s_z_XvgPgh#`A>LVi`i+Y|6dlx+N0b}CBT*^5k^J5l@%7v7c5-$ZOD20mH@3Yyard<=UJr2(uVrA?MhDk?Zu3yM6P=74 zu)54c9uP*_BfUrdn1S>+M@XPR$C1S&?Q9rmw@P26%v_k(v1Vkm*;wqa;X>#>Ns{7? z*=%Nv%JJ058@5Pu(3nfhx>sVTA?_h}F)Sx;$#xy5aT`dejBiwVy$A?am>a>n zI{JprmyUoVc!lPJruIWd|JZd-d>Gq-zuSS=6VewdG3O9rkXS*f+#bO`%KeK_S$3a- zje(xd)88;5^x68x1I3@KJ=gBKeT47ZfVl1BoWw)fI6H6B1_m})>#}lJbvR35F50kH z-Q72n2pGe~hy`^Vm#B7d{IQJf){9b0@OP--|4BO65K81pO^DmJ?0ijHNTEAMm+vj2`ul}b-aKgtl`OB;Du=~3l;7nT$dNN z2%9Iq=i^uG)8{L7XHZr~UfBM?_vZ}@c?X5Oux+&8rp~_#BuW-a%apOgt~?+DZ8*em9A zp1Y`xy>~LSP};lL@;7(635%q(yQ7{3l?4kF?gV!-!-HeOJLcXVs!ao)YF-Nzg` z%tGAu<_}TanZ9%pH@+XN4*M>CczM?0c!fKWi*$q9zFc`gcsEwM@q!fE0CNuF0v+`U z+$oqWHXHt)D1ELHclILlb7f>(h0h6xZ_nk8G0nu%v&b^0>T~J;;XqfgN-Od^+Ef?yNp9IC7oD_ zN#7Cg>&x$XubNQO_V!sm_YHv`C?SR0!AKoCX95o@1 zJa5HOnqu3&G~?7U;7P!=>j+3RFKiHFKHETGjnt=~X8`4y5I<9U2tg9D2Aidr$vJ(+ zm^1qI#w5SZyDP-{bJxifX7Nt?Js-&AU*;~8XVPWiA=bSE1H;K@g(I$xlP1I5Oh59y z>b6t1S0%-67e=BU5x4d@J@1b6!vCVTkmh`e0b3ybS7uPa%O7_95M%v(alQmd>D*s~ zNj4YSVQVo46o=qm*h%Z87|nwV2FSAzLsSbGyhvw@^QlMUe-QXpjV zFSD1)Guc)+;=;V*ATec{ZM()EY%Em}%d;Pb$*(`h_Es z(d&1KG1G%yA4#S?irc*3u+j59GQRd@l2>LjGX3^4X5h@jNRMolLPR31 zd$L*| zjMKLfV%~adPXQ@S>#<_2>!V{aA(MZZ!%Uuul~&@-Wn&EyWjVahtgSPOP@wCh5m62- zF3fy#XDW5WjX`5`rX!8-|F1#KMzd49w7!SwbLui}yUjx~`>`=AE3^4A79tzszDL z&!k!H`Lyll8^uJeufdpDlVlXPBPa072!CW0aCzhpj|<2wb8KVY+8;>ol9?DH!aVS* z3`{iJSb{lbe19?KOAYuhJ|?3FVUkV1w~K+O*KKbpP;pFVcEXxWaKM0XY5gdr87(rV zur8@W3?QYkIXwTP1na?@>M5*u-!T>fsm{n#0wkMvoM!-G z{y--oup<+V2P#sKqx7pmAY}3{GnvUV zNoAbZ-!Z4{9%VNFqV2stgJFi1-;;a(lpyZ~R$9WtLr}tkRZRPwM5J@k2}dMm|J4-$ zPPUAs8SyL~lhVAY$^eoryMx7;H77EKK;+QlxdhmeiCV@>3M9EW=@Kz;)XvNs6HxOog|yVWspc zdk+PaXF~jbT#hBNNHjtUBwb5CQJACR>?J^=q&)~o zq?v?BFeA5*L_nn7v`dQhy8}A`@U_V&og_e{YkCN4c_v8JVU!Rg5n66a0m+mkx;7^x zY!k3XD(mM+fW#ACGSm7to7Ogf1sK?dAs%x z%62qy!EIq9>ftzS;o&rxzc;um!jaSZp8}|m(;MDkO?KuRGcaEq{*FoMT%C!qZoaZp z0z{5=ld$f{q|>oXF-XFF+cyU0a<4}SbM6_s)+Ju>q*zm0Q(gk#9$#`^2!u@jWkxf3 zCRZ&iV-r20B31%XW%0Ah93(Q1(2uCCmbkB(~ z=W7k78JT#T0*$m>w533piFl?MNLH?$&j6wU0}o1o9hqod{enP}%Z{Fw0EyVAngIAN zEsr4tB8{@M6cFYjjzpMu-4caB$mCzAe|WxPFf!tb?9P~i7EyM^#Y znc>{;##e)HpqPHoc5ABlK-%eg3>-NmPnBjmIS3@QR>#1C!I}Ex6~1uL5pIM-o0fI5a#t*o-D?kA*E|$!gtPNTGyGC3gLH3=io2|BpXl7 z#zcd9FB1YgGSM%Ll7J*u7>s2A$*k?=2s80Ps02vteQ-|w@b=2Y8qG3MvBzl1$kJ+y_FeJ;&cem^UuzGl1|-#wiq#Yz_(lz#P8= zIx+D7$iKp>v1hyxB$;OSkbxO#^BG|d`%5YW5>Xj+?M$hw$BOaqb)D=OK*+@E&zV?> z_}RWTVK-pb-P;~bQ(nU?pX+0W7CeVpK4Zgd%d}yZ&(NUVD>*RBXSoz25`H`NaAcU} z(?f_EQ|OB|%^OmQPzTup<+l^*692*7h8aNlsZXT#Ct`H~|yp z`J|VV-+$*J{QIjasTima|C*`5_y2hYary`Jq1Z966rQ=-Md921)%1d@%H z=^h99cxe#>G_j4DA_YPov475^B&BZ;X8vJwO1sy&ZVx-VMa*W8L_jjGM=;GHx=$sb329-t6!X*00|5Bu?vD>bAmkDA=RArXZ>J`w`b~s+ z+CF4@y|;vU+O9hBiZ;SLZChI|oIVTlw7tDJbon%-KA;c6kr%#pqku3^+dT

(qxJ ztRHouqd^|gf6n82*W#IWQAb*kRUfToV_DNt*rx+0v?K;-oQ^c=;n+(^RKLr7Wo0l@ zTjL~#$VA+LUH~|+3|JIqf&otqq*QhqFtGmeZ6<|z?~W7`4SXsFbmZ~eivfyEi*KWV zgn#@1OhK#OJpk0C9&0$(4cV#-)@EFOsoc|`s> zkE_>uODC_9W+K}oshy{n97Ex$0o1$~F=%X-_u<8hj8MWxb3@xtp1uO#JqZCxEP4c`g&n4%cD_xjsc2{j|&7q8~0?n1k=~N z2$O86z9PhY##l!IiLAvon2<-rpYyn4+3CX7Wv88~ZeDZY5BG1Pgo~~Y&)EYh&y}UW z%3E9!zx(vmH`Yx^?M|x@M;7@RAk6fpRRFkRgEvTkNF_*$N!AvR6Ju6Rzexd+zagc2IkwEKupLZ{Lgt%gN&;!$+%aP zCdXmTwIUA`KgntN{$exAE@&g7Zl6PT=cce+b~wR&L26QnNG#$0M402u7m9)8keL}m zAaPA~qZCM~J~zRdtlNJ@0z?MGb10xZkEvU+BIoBuBOu{7@sk8=*77zA>*9@`Qmh+q zQVh(>S*=(@9$|mZgG?Qn+vQbIXJlmic0~Sh7K%Q{C|?tG0}bU}&Uv(T1WG#Rm_P65 zOQhcCqYxs(eA@90AUMgQIh)&GiuJdcO+w7Jw}`^5sO}>HBGvVDeIU!DnRA|jmnOV9 zKmm#PPkU*O3A}|slb?3$eyM{NLHA3lx@%!C|YA&^p?J_P`y5@s)G zQr%}r;D3;=%?Jiietn!a#EJ|rT?T-88Rz$vVEQ-rKtQr#!w?4M^X4D`)Q&;Z5fJhS z`Hwsd18nca%u7VJb%VummP=9i%FlDhYgM4J@f&wP4q~E&+@q6fjMgAE|DHG^5%OgS z0P^U5j>1e>Jxl_mREIN!SXcY4mSPreU4e;={eDwGc^-Fy8KB5GzB2&srW)fVnA76v z`iX3~sf&LnpPd*c#=kE)upR*+kKjM&aru^6@?wX58py=byZUVGH54P7af?ek|^UTE3H*WIx`H51j8#;HT`S>?|$3tAn zeynAA{OqL#0%i%l`9p{~e_^x~b5>Xw14!1zw4b%j@6U>{etAk;0)#w*{+!3fQrFcP z;#M1K)`tw-538#E(~T0G2@nhWndm! z8AV~XsDy|y-*4V71wtNyf6n6qyQ0u8wyl}!_E!M=*1~fzpR&Dk=S4!w$5uFYSLY=N zxtI0wG&lnDf~|yr@Rjz2P7F-(TX$MVx?0lu=QR4xOg8==B7_BG^{WiL9Cgw0KtNd@ z&3)Ui=ec|a{*?^#V*yl1zo}jnkf?wA3zMv@v6NuGo-+v(@|gAKJkGsdbz`U4c{I}7 z{MkXP-$4{TIVbcFB&fCWNjUV><%@>IQ^=* z@UO(=r|o)@0gC~e{5bBB5dU7ORwu4L*{OjH9bR4@en^Vl2n7eZOR(Azh{(>ZbhyCB4>UJ^ir=u8({%L?jlJ-@=i} z)EV~xaD}?d?tgdF#aMq^{hY#l8(k*^BD>|EFy+_Bs)rz87A;3rDUg^ucM!t#Pw0wC z*6>tl*4xnc7xFQ)o`T;YkAOet@mE3N`yP43BcwC4wM%0}7z*DSTXkpMJ~TRY$%L7% zT~R`)_L)~p*hpPvGXo;RpNYPb0Lev0YXFePld%Z%M~~waW^30=5=`k6J^=ngVSFnC z*q#Scd7du?MJ{G%i7}TN4yN@S@(I?IO4re}u9(!D0h+urDSa`_PrI?(8qY&8Kkbl!MZ*um{ItEq;Fb_AJ}d>2 z>4_#ZMdySRXi9}^k4e^_y9j_D%GgLu#QgFUQ-0nquNQ(M7rU*cFwwnGn&X&s`+#g! z?}za3WxX{4Fh8v>2WvRq{(sKnblu?NOuHA3RJV(F)q}&f!F-Q{egrv3zaD0yxV=WTxqN1g{C(83f~9TsCo4s&B*{+y`_fEHv6U7wTHPZvr-BeP@UD4_hj z?LAZsiku_(%D{|R#g}6GWv5{h^)HtKpkA*dg_y5XLogu^zyHXivhm6B>r+MK@q(mu zi(UQEh?Pq>EdS7*A`&GZs$sT-;LB&{*{>>wcK;@XBQL8pWdMm;?~;TxkJ^hdr3GG6 z%%A4R05Ho`fBFssA3Yw#04>X-F{i!OOziFc?4U={WZJ(g{v6K6ysrga1F`X_N9{7rEn~3=PCLfUc z#UBt5T2$;dn&TS%8Spa^{rBH3^>RdErUBPBve=bt*=_H1T#C)|ln&MF_5!T@^7yy_jc-nOdkZc}0 zS&I4L$zTLT-S=u>%FnB7XDLu*s;fKxm7J=xS%^7iWFY{atW~&(v99!{?@Q#X{9_0- z-Kk9yAkoLhBD<37bDm{oq40au_7CeWpfUBVpD)Z0qlAj5wWrpuKx#e>Qiw}sZBUm6?L(ufb#5l;>wmO zYY^|<<|1TIBDG6%g%FX5TmO)Onc+w8cOq_ZEY?WnR09PhKXuwqvx29GiI~gwBcLpg z){>WEP()%ek^;$;oF@{@FuyMVsCR`-O!D)%{tV2Dq7o?(@|gJNJdRBX@1yhhW*>?X z;E?aasX&Pas~g^o8cz+sJa*gaQGX$~@<+CdTxKD)pX05{p6d7SWLyVay_GDlNCFEiavn*bcVEyjM zatc$r&lv$Bj|qRyIsITt|tqdo*vf zEk>%3Or$unPY-7bbIIDRLd=9^sQ@^y3Rj4+zB1V!0bxGpDh3d72PX@GvOK2H&J4onR;oD7|Q3fZbC%lJ2d3vAc2$G z5TxGgq!1z!^ZUF4K)b{}n8FM%_{PBeK@=e%%x|_qiuvJ-HYRdNY?T0I$D{dKS1~9u zHEk({Imu+U5OZXungmFE{T)kTU5UhE%vU=K0634a|B**v!1=)yIddqT3$FG%T8_bd z^K*;$Txz3?wqw!BIiW~N#Ivf$gQ9-XPj*o1Bp|9YXDI1Ry;tMm2Ik+ zkjI$+$fFx`ytB^9waE5xx=&4yRv16-XL4WN7{-scZt+WvMhP2UtlSr~9>$Md#SoDQ zy(+7B4NED|Fn;`w0fh16_A|S>ag7lFezizd3`Axw#uA_`k4EWtAt*9_vo-@rB z?YKT~vUFrIl{Bhx~80(fuF2a;r7zlyL z@%TQhWqGt(u4aHDCTlhUVDvL-8?A%pcNc<&=XU90%m%GjQp~68=V1+bjQVpP2N#RY zZcgxaLw27Q-T1P91B!?V$p~!9M5AMWUuX6FjS^?7y=3W}MrvK{gbO^%$@;#@s2jiFwne~F6v1TXe$RS?#u#q&?WwyC zC9>`=Gjpj#YJ+Yfh)DPucM$@~#dimZF%zf{3?QYtv;9neQx^+Tg*<5QK>?AyeT*0= zKd+iAq@W0&FJ%(Ux%v?lrtpY90+J2m9*HrfUPmRE+Xqq<5b_xDA9-kfy)g6Wy@6CW z&uc4&S+0WT5{^UN)*HcdiO5Uhj4g<7-w;u9#TlMUd=f)MLgf1#VWu}dqB%QcB1jda z5^OC7lC?!53bS&$h7^b#pI>AEWqGvZtAc>h+q4}NkW5)V9btwWO_l)Rxx^Y;*Cd{y zK*Mv1ZZo)$okd4NBItd zcKcuqV5N0BcLpg<{~Zw6qy<{ ziRP5oObY8Ut}zlIwA(*Ov99z!Cd7P|KOXS^MgQOXKlanz8IP_?v(F+%Smo8IUwe@7 zj_+aft`5l4R<*AD@>CQtN5#0jd=*kT=_kd`)7x%{F&C{^Lt%#8i)CQ83@;J`iN?qC z05EUu`nM7wGN^qg1{ z0^xWM|BpOQe9dX>I4#FvE~YnmH)_Jb~0U_dqyu zT;q65GX2Fb38tu_Rt%(6QqNOZ*R5@@AMdwp!WudGw*lmNv;?FhP-Oa)R!m|X)l-PI zi^XsNoGo|lB$#FWuZn?$#QGrxggm_doX7tAuH}`Q?o&R~^ydy(4BjdbLsk6$s#wApN$qOKV2)_s)#VUD+4 zF;I3qzHe-0fFc(17fES)#Zs6-N1en#vexg85VNVeeJd}5jPdV~hv$FfajS=l`i-1r z6hnLW^c_*&DABD8KljlJ%KO1u6W8%m5V3t~?=I(eAax@?-U2cA<0qPFCIcv}Bai6O z{CQ}G1T@*W?gj(1dO05xJ|gh z_5F32kcY>Ay)J&}6v z=@djH7K-S4j+}pku9?Z%6Da&U`7JAlX501eVyr7lf&lPVZ^RLJL)SR1d536{R96>#EcJ-V9pY|i-AyIc}uab zZ1{hOyAPnMu|E#<6nk4678H!VD>hUtIqF61*n0^oO+Y##QgUfxZ`f8>WADAMB~k3X zVnw4=MHB=)^IRnzFsgAj-0x z^A3(o4BZW|T_${%D5t$1N8`TEy= znNj2m>XDy+;=BOmF#ZF;e%1Pe3KfGZrlkDs>K&1MWMQg6)^2+S2TNWbGwTURa>%O}NgBkTaA<94^!_&Pqy{_aiWM=rwB`Q2H9toI5AN?soZ3E`CZ)fK}T zIXZBeL=LuzQpvm>2Nfu1tT+vDJ;P=o2T@1g((BlzWXbE}V~&W#`#VNZEAPF&mO@TDFo2Y^Rs^VIeBo}6JaLGZAnNE-dL7d7 zkG7|NP3)p@!u@gA;*FVi-pUI;quMAo{cE+%Y9Ge>+?)73G|a^%OVp8(ZiRYo1+<=yPxz-qWJSa-&#IuCJQFhYLg!{NT*h$RQouHn95L%Zt} z;+e6*R_d7#C>$?B+!t&I=)PbF%|9(>OO$_^)z?&C{w*TozF>h0i|W9A!CZ!f`&}uT zi|c&ST%i1t4CtuWI7)ky73{r5isxNp15RDgW?@<54l%JQKA+hgTt%BZ7P>2+)$`7okI z-Ay}Km19>@4!d_^P7B3Dhb`B$k@t^IY5B#Ud5^5H{p#d~%qVxaf*$#RsUZTCL;g-z zNzaNiIP%Y`i3*vs-A~J^?IT2N59?S4OSXr{U1OS8Cb8Aw{5^bAeYOpXGNzB9ljJKQXISWxjx1ygf3MG?l*VR9#d^ zkHtX(IlHhbz;@u@ag-B>f6X=XMLsz~VY`Ce39?@}e`|66B{`o&zwR@927w@~xG>6v|n}N56#Gl_X@;(W8ty#^i3C zwREm4pK+k|mwHkRb9vo;L5)^>n67y@@zaZ`%(vT--A(({W5(e<&>`+ebg#~!>@jsL zW!YhjrNFX zB{?X}2+;NE6&yJ{-a&!5KJ6<)T%Qi+$oK!8Rv@-_w=(LublEP##NCwDvHP+8&Fibo z@pbd4POrZ++lY;WJRfXgUU%Y;S{s-$Bjq4NkG#+C#T=CN7qug?>^1zQK)KMZi%Mp$ z-owaOwZDMFfv*YuLPg7#p$Ph#G**zvOpJ}7E4OkGzjLiFq)brTO zb4fX(!VH1*zSUcWis6|d3g!3vOeHcR^sx$An=ozr6xZ=HpaDZG@9_64`Yyj<#RG*D zmaUMWoVnt>=Et)!BIP~bj&l%obSa|_`To|lpLP6|%9k#sejQ#3_Y3S6blY}csk#37 z-PGavtZtVfM$L8N{E4-d* zwVKbIIyv6VXc}{Qu>D7mbZ@4!UDo{Qu!YRmVntlLyKR`UzJW^rQdpfR31Ud}DS>pj z@}48pXZ)o?#b9PFfa8c^BxPoKy)%Qwb>x43uOP{bR$B>VSok58^fm1&lG$~fCGtm1 zH;#4WjYff>(MKd!IIa9)maIN4-TqJ z%KmoR@lLikJ}prD58*G4+Q0XyBasQa!bONWI+b3>raiSrd`#MRPBB~6_|cO^Zp>%l z4VMqWQ`OKN`CCUk%#!!~u``ZAjqhslH zZ1nHvJ7Qf=6UK$^j1Jk$G2wFg=(KgKm@e5OvO;Vqb9M?p8+i9KGdMF|pq@i-d#hwf z{5p6*i#3xS^dzfkF;axZb!1#J5RiChpF0Za{-PTrXNR@|*qW1e zN|X~j_ScN_Qxs&>@gMWdlsrGS-hI^BO?mMSj9Zbrv(wx~IPW>zM~ynkbPvp4o7(Qe zdC#Ru0UZi(-g8Stk2vqi)f6jMVU!*0|5eGH^V22BXDx8mOdkDEMP{|i%~W7X9m;(K zY{hZG9O>avUnFP${;EJZg}F1ziNl*pWZag1fa>^R3JK^RZ_fqLYF2E7-_;FF&u(&5-TfKQ5NwyRM7 zWg5$osaDM>vrf)Q5-hIc*9=c;mGzO&1jtW0^h2}Tt~ng#tXmf(C?~k|W#s#51qwtR z|1rl*NgdIVbq!MP4co)adff`wtslvp?3ey**cRjYj_JEqHT=Q67wKona};J|AIhOe z-l?K?hA4++JQ0v(?|`!a$J_R;MarMmyBv9A%{~RPhV840u(*z&UQ8{#Qr+-xq5By8c+H*mkmwiQyqOYJK%q zUB1p!tj731IHua+Ik>jf9jd{x$;vXU4j-0WQJo?j?C*bR3HUY;fDJYngG|qVv$;t@Swc zj1fm+^e@QnU2*`fA3J_g$h_p4A{n2S&Ov5a>u*p}N4Dh-6-oA)e?@@2$E#x!IWlLb z0_Cr}4{>CAGZTS)(YY&$I{ssJng7+nc@5gqG~_HZcXjnR^8P(@I#&PpyJ{ElT#1?4 zr}8#3G#lOJ1}xabwP!4W>9B=e7ke8 z2>E_a+Pf6$_>Y-oO6pj%&~uONp^7b)ij}6n-Pv>i#`8V5kLdJHX*jV^3C-TkOs^a9 zHLnk%@q8Th9MQh$I5mFd7ztUnnmq?#JYU}llrf&~6=jU)D|)}rN;Kh+i|fegZLiXI z7tQ-5kkwoGGIHd|WfJ7CO-<5r)>%gdSxykok*MQ8=9T$h9m6NITHN2DfKPp&>@?P{ zF>}e=nrN8Qkqz4a>7Sz+yO?jkqs_9n&SJ*Pnp2N_*rGfE%JYXfYueiesg%=B&t*_C zI66hksqe!@WY+UVRRtE;QJB1&Ly~c2Iz* z<3DDVDXC+1YTMdu+R6$F*LQ2ZAp>%l_xquJrIgXiVExqg`{JfDE6-Z1-@D#%F;8MBS@qb+8muB=_bJlAd8FmzC` z(tQ8Xx8~&^u-5AE*hapK*j8^Uq@c<1x+53OTl) z3h5;`19bg(umt5;DNOU_H**dd*N^*jIP4ee;bj}CX!^5BV!6l%}gc~gYA ze!PP+t{=ChTpj=azx?_CKjy3X-|wmIQlsPbUS4DM*8NUgzWz3IeA0Zj$${?7Hhp2I zCH_gw>)z8TA5T1HMtdCTC+5ATbWori8j?&|_EMUuP%+@UDwH!e&!sFsul|&Stm~1# zDHrdzzZ^0pBsqL;3r6;y_KYJ(#?+I@oDc5+{WtT-PRm z`^|m6D`ftZTMWus1wI0q5*M#PR=2P~g(a`QwhQn);_RV<}7T59f zUmFRDcl9+?$OYxM3#8Dq8iR7?=|`01PkmRbNAulA%XF$g{Qp<-mv4IMSnXYn9BK zep=ErbY|qIuhmt^>KpDr*%~{VEg57Y_|o zh&ujbo|=+6mb&k+I&@B2C55v~a$dS%IP)>mH9J>Nr`WjX4&GSUjrHwvzm3b7>&)oa zAnKX-3;#zV^{+l@x~@zYDHr6-R-yc3-$sc{F295_tA8j;ge7%UNk+gp;Gqi{>3e6G zO3oF_E09l@0~qDxHctdnjjt?1)bSs))Rfe*WZ}1m^+KDRW{p2PH*02Jl{p`_X>EJY zijCNiwtD#4XUs=zygaB%M`rl_p+Y_X71l{J!aa;p_C7uaV86P5LZJNk(n*p(*J(At zjE$Y5Iapjr&izd)5+BkbQ6fFo+biV6%V$aX=lBC0nUK{_%cp%K88Yhlk2z{e>WFGQ zwZfk#uajA$>(-Z2Yky=;?N6U>vhM-2{n^MQk{L4Zz7f{0-A6JbziA?Rm#~oI|S(C6twcLS>^&cVnLD$<*b(; zl;zm5w;ANO)Ez>ij{lgUrlgL=w$(0ITXkk0s~uKpt;wz-%x&Ya&PIX5SpWPcr-pv* z%lxXw7%Zw7Ve~I^>Gf+X#m=CYLN_wtnu9<&M+wbJ;Rz-VCTtFVf_gU9d=jr>! zb>z-jhJd^8Pqh7p4_UdCQJ#Bhy9(vh>pF?dvZ|!{rlmiLIvSLI9V62CqxJO@4n0_- zYw0WAWIe+6Sn8FGaWb|)2k0Aaeva)iVaTyxBd|T31oX(exIUzwu|3vkdXHbpA!B>2 z6rlWd+HlJ9-TALX$eKTCA;6ODvEGnFk^}a|OOSVcFqASs!sQhyXT7Q`kg;RtNRZ#s z@CFA_NBz?4h`2Gqiv<^8HNjAa+&a;zxE|AI_JIU{T#tF%uv6DpxE^!tcGXrBaXn_Z zfF9-0JJmEpK1C^%9lCd7WX?Kc3G!JRhiWD}R#K2z%?@=qSX@WOMeS(KJO6uJKr8RQ zH;lf^&koL(pqw(fzm^j_1S!aTobLw}qK9O!_K`+tiW84hdK41 z{*DbC8|l5LdC3(aQUgg;yNe50@Uh+5?((z|kP5i-Nn4wU65v-T&PO0PpNI?6G z?YRv{Y3(0VXOXN~UNse#Y>#`VL?n4uB@b%l1DaK33WzXF*NY^defy)zW# z>iEB|fuq4+y=!d*jW_L0ckmQB)`iP`pX~B^* z&+Sqm9`E|9QqGEhERYGg`3giGRm-TO#VP;aue!VQ=~X(K_*AIHT#_bLlxGcOgMQ4k zyk5tP`SuL7Iq^Aw887aE4skhgoDGAr`;j#&X7tX5O5h_yq`oy8cY|-1^s3TWcGpgjLf|JKq5cH^`MM8s+3;GLLakb{Xbfc zkY_2lMUqJ-2n}U1l6*DWXR{Y{pqqb|1G{BJHol zGxFQELmX5LHa3+gr~b2rWIfj|CQIrVQC~!o7hI~LkdX^I7v+tcHM6hU3)G&O+KnSs z^TCu+N98i=Fnzxvrd*vU*2t{wvKdk&bKde}m+8S*Y{aIkK65JKaN^peXUy0_X81Hy zqMnBp?4+LMh_#a#=^gGbLB&w{qf%C?71bLPb(7jzPuJ0)q>fyfXGnbb^Ip`-3olbd za!|-#1YGw^#c8|07F=N-X%`~HbYY5d zfXkUvZ_2Y83p_`i>Lf7Z3)87*K0GCpk#_bs1=v5nqYp>Bzs z&e=GCHF7jbJK^@0Iq#Ukubc9SjaXAN`19L0%*V4AYfz&NGkobmJ<2w9mP=&BN>L!a zmp@UVVtBoUMEP1wB4z&C$YTt$o(9`FSX@W0TV(}_4=?{zbI!i=66GQ7FA7l3{!oRH zN#)0?Go^7utv&u|Jqa0glqafvO%no-i=fS*wOJCB5`7@(sovCNueWeL0 zhqW3gk=~{;0LQ5#`xxc#bK45!EXGMBYzKY63E;J8zl1Mm-QlK%?j_Su%wRr z<^qyD^T1dI@&UTp5;;?erHs>tX)2i!yhAf~pDBloI*iJwgWauIcT}}^tah-A{+V)> zx&7TKeaf(C)<5OXrHt=BxV)5L9&+P7Gm3LXhq&;2d=vF7i_caHr04m`B3aO2Bfybn zdtJ?UGiHhSf^~a6hl9m+=YAo1Q$dMac{jEzdp{a1UxlfO=}5s>Aq^LErOzj-j4 zL>-2u*WtVIL1*(^TrM`h&~unq%S`4RWYzCYjYDiC>z(km(Lv_pY#}AC-NKAa)-v?S zkFxkGKsmxZR3g1^O<+(lyu3-F{I2E-4)XeG+E2vVXWSQ%i|hFLxV=Q*b&DCVl70uQ z6;i1CRfKY8OnJ>u4z)B7WNcE=jyep=s3Yv*=!6W@Bxct0u=Vt}Du!iD?sIc?PYlbr z?dzBqhhZ7Rv$|h-hhZ731oX%|Rf*;xhGneL^bXiUVpxWemOp>`S4Ea@toH&KmSNeH zEUBaXZUss9X}FbIc`wgt5;@_#_BxTXxJs1eSA%Pd)V{XHS;~4{Na~^OKW<1{4t<+!R9LEYcV)Mw)d>Xs;;_vJZzOGbL~Iq#`;0y z*r2T!nnsRJW4@t_4u5plGh;rEI^@MsjQ~!&rxcz3xCLu3203SSUxnKLczKFs;>jZd zWDT!wP+)Nx*?Sry;1136CKBYQ%$}uXi#N?Vd?9~#AIHdV2}1!6zHcH#sMiISUWV7v z5RdJ3qYkkqR$;4e)WT5t_Er!1&hBhfQ|n>JtDa^)odfo+Ebq&Vx;z)qAwTtVA_wL0 z+g6liZ+S1k*4Xku%b)rcsL1l6`4R^iHy}}h#bsor9#W8a>s30D^f$ZDNMYMN%5vt& zIEhpntp^73=6;nFsMjqjql_j`5Bhrs3|F|>9nIbD-N&7xw2@Didn-0R*8e0YZDf6| zqLwwt8-+VX%OrHj`z?Dfkowq@3hBCJ8AldmUKOGI+LWc3|Saj;}-R6atj zvi%$1X>b=$vUCW0# z%P8a3gYyx4$L?g+{x*$1b8aZ+$9&W|_{u{zXzIl~4Q9Gxe$2J1=bu{`&o@n_j>Se{ z9F!M6H4{mh>=D9MWkN+RTEw;d?J{(A!j z70Tc4bYbN8S_;5xqiBB2>iD{@fzSW-^Z)J2hLx-;dRI7HKC!%cxk<*WjqVs38TP>+ z{D0s-z0SY%U!BL%*wzy>A9rK5T2yv>+^89IGi<$SU%NhRz>EO}9{CdUjcuCl`Ry7r zHq~}WKFIcgh#yQ2{B~K>!=$oAInVEuW^QsEqx|X74h6CXr~Fk|@~d;MDjv}3hWq2c#g)#e#7&mo<& z7CPQj8XF#1Yh3OyYaH9cVI!V6GQN~5p+i2RtM-=4kF~18kmdgJ7m?a;WxVB}V)*P7 zzz=dHR)owra*_&*%g7n@m6E(@Coe&Mv`Mf?dC-g%B!9!iN}-%}>6AbwxaNycuk$Ot z4A(-B_<&b|NlN9*$Itv~x*Jz!#@aky{X@J2bBX}v*6$ROf6?&}gSax&OM$pD;~_v+=iGM@mXzTnP^-Lf;}=a~=|z$9%o`C5 z@~4J(BQcP^DMu#ctrwtP=UYY@r!KYKxbb`_<`ud3;ut@S<`s!rZ1EiPinvwJ`jU-# zMY0aJ8Y9uXA__Xhydv5G0rQGHN5lsPds`7=UXg|Z#JnPlG!v#*5|J^l$Q}ijlyNYh z5^pb_1o)s6b&ZjjS44oAS7ar?&gFAVB0r_(P}b{wN-x7@+M32+@77$x8a0l2Fv-vo z*E7CsI5F%du4g=0UGdT&T+gU~=ScJxT+c{j)ba3geHAid?Q|`Bhr5!vp7D!=vQlZh z0Qs1uQ6#Qsd;uk8UG|wm*IT= z#;FahZ@g2?yEcBiFslnbc3w-Kw(lRsx>K!XwGyY{V`samHCG1UW2Zld4l(-whh~`V zQANw`XG#zsJ4cI9{yOk z!%9WsBU~yoD0{tqBa*YmTvZ{T#V;XofEXf>iP2LOsMmRxUWSwV*;Yg9oQ-4^Qo4?- zcx4?f_OyEFLUYvP#v=ZdJZV~UMrQQg%aPyiUt{EpC8tEF z*Ojr9j8jdcE@SG~{v#(JY&W^($Gyz+s5mHRSu3S=js0R$fXdo zdAIwIHD6~03CMEp6D`Zh-PcOUdR-YS$vD1V<(P6|+dKZr@8RFNt&C-^%e#9{T#1=J zk2@Cr9f6f|U9-l|I>pyO>6%hy!_Y;MlK_kzm3W(Hna)S(<@a-M^{i<^T& z+Fw4-k;!=%HUIE2MfuD`3GJ-I^-xf?N&e4Obo3jP)-6={BJ-)THrpiNykfJb;jNg)6}iABwV~~{o>_+5wM!;V z-TNE2Yg^m4Kaq{wwH3V%2R)|Swdjv-*E%yYL@-oHduftD<}~*gA#T?$Qpn_eTR4c@ zwM>M?W#D$LfW*7?udR}<&4RTo$X zTp@pdZcJIuzIBizzw`}~AZyjQp8!jaRLzE~NOH(`FUI#kb|uESncl>mB9>tdNcZN zSfTcskQYbH5h?q(UavyMz(V`ea{Q?Uh`5?yG)RJaT^Vc0%pN@d^_(e>dSMQL@zJ5@ zgJ=!_tl2gca{%1GoxA-s<^Y(cb2P>|)t5HZA?5&Z05k`{Yc0D)KhSdF?u{b8z#IT^ zDw$+{lrrW3*h&^}jcj{o2}!p0zoC-;`)_I4VOLWTndSgs5OV;GlE{xMzf#ug%2+~X zmRh*1@sm|6Fx973&U24Hmh>}7XsS;I9b&4_XClN@ zpGCWe-Mre=E+>Yq<{+lRc_UCxd)<_SmPeE!eqf1%E0 z*0)b;dW(W``24?tI^_M97XteHKTRQBqfRrjAk#vH@{gl^IWoD*8wOeZ85~)%H7bh& zl5D?sgi88)ozwDM_ZpP(`F{?Bax!;OQ~h9th%(ApKxTU2Yrfo%I~SCiBhG|JpGjt; z4_ki^9(PQ!>}g%6{cw>vut19m$u*e4?K1*8l*5{RCiy_FnWo38864%D9-B0??{N}k zb@eimbr02E1;u4#?RKJ8*?C(MBSnW(964ug0KnfGZmpT|?XKpl)4YN&^tv+EkC|3) z_`&iIjN>pKoRd3VZh-ONm8-A!e}wViuhJ^F|BmtC`Ag56;rV})P=*dM9=wJ?ik~)U zd0wAXg>wFs>jK1haDR?WiAz!-Gcjl{!ji4gHc~;tc<@+F#{-QR<*747663+=im zcT41l1sRm9+jy6%Cd`-5a%wMy;C_&cKa2A8bzq(wv zH6jupUUMP9)zK4E1#(DpMS-|Fx>_QW%13MYT#ySx*6Yexv1Y3IvkMbEjA~=#cIds& zTc5A zuPbB0nyD3*`YoLO?XqIBpyh<>`n}Ajcja19bB`%DsVB{Qt&3sHU{Ls((bvtatbq z1r~3OU;8*}l_PpbP%FP+)*FFzsIZijzvZn`DW`WEPMJ@L>cF60SH^NRlkZ%9_o%|2 z->hL-a2}EZM1EFQp{wg-;Cf(_;Fouz&6|O9^69(%k^t)aR23 zu^gbbZS}e`R;!t`&B$l9)$($zW`N1Nno$p#>*EL3H;$*V0U49M?{yl(e6OWuCSENvQQn*`;+zv`=`dt}jGoO`{l)}FWa7f0>!EAfE>WmJjDP+>_K{j^T64jaB-pQ2seJULh0rRbf!Cn_l`$e2jWsc2nlJ z=X)}!*Ojr5jIC4ifle{5dsxj_mzD>2eP^y;I(NR2RFe(Zl~-styBG6)U>-I4O*LlR z#FaYagDy=Gp&S_2LLlAW9%E!4SHNlZ;Q@Z19abrj)oEEnf+c0x_NG?ZeQB5gd0%(! zfXYv}kfS2YsgF*JWcr4V0{KERCH1;8){(KT^`4FW+ua!lTz#8?o#xU)>EffIN&x{(IF1FDbz6zxO$HCt|3Se2V7pI{5@x>=HmxdNgQxbga4KB zI+$8z{jBFAb; zI_4)tR9DIHbiG7+$@M7XsoFszlw)@cqAb7s=EESKs_iSnl9zab_T%xZtKE~(%8%At ztK{IZn>F)x1}eySs&=nPCfQg3Y>hG&kr`Gke?hrnuck7y0~S8lcJ*V9V(yS*UpCcV{>ZL zpTp(v9U?4#ZR8mI!y)ktdsWl)Hoc38ovT(q1V-PYzT)A`JhvUn!zKO zNZFxJhGyRMZYr|;YrzDGObs5uL1xkaF9{Zx@$2SD0Z9(6YM?@Xc8hk59JqNNNB-=# zAE1owFGTXaXM_axx-yoK8M5o7*Y45ltQ6B`b<0^D`OLgy`i?4aKcZOA(BW>f%&fjd z7cN#8VeD8k`ooX?yF5pR573TCyyudOl$qh|&j3HbIBh4wXxlVK?LD(@i?FziKc^QU z;CgYlKncp8a?wkEhRrd`eAfC60=1_)^9D6E(5ozEFwyKc-oMx!`(e{dntk;#X zfQiq3nLlETe08pBLwF+w`~AaYT>r z=oc@}Ur7p-+46QIe`BYo1Qo-VQIxT8*K`iD$`?gYQbx`wTZSYr5-XDY=xId*^|~@vtQmBFynHrymA_JB?&iN*MjvJ#H(TyApYJ2b{aISE?d8o_iQ{+A*%67% zcz1t=Iu1L3TOvJt=WCAIldPGV+>lfZT&s#uPIvvpAhUGcqrl=a3UuiL67Sfi0fSgi zX^~3$$_F{}5L+dXc}cOFsV6n{x-u558E`x`CHqsee2jqnJM!huYZw99@cWz3w=e>7 zbJu@=Ctw6*qs;23IE;YYuHYZwNA(@WAVxqI?Q-6?CQ8T{0jV8U%myaTz}!{UajrZlPbPmG@g?tYyI@k0HyM*AA-Wr@4PL52SiCWWBD8^=jDig`1ny z$eW;W0}rN8xz~hwe>A`O#mZM19DQt{&$=zlYR&#xu3fM`YRU`|9r8Y@oJ59>+ANT6 zXN)xqF9ecu##~V%lMK6YkTqYI1xm`u_U%kb?rV99k^br08?)>X+eSs^)0cMUs6FMP zERi4gED@nzSH^NR{pSv?H8#tOQ_PZ2C!Pzgj`O%)%xrcaoX1`CdWC5ld!=j$cgSso z^SHm!A9kkMCX5W5`$p9Akhv;k1{YH(%c+;22xNQ%ZLelspB+amDdSh_L=K7fT{D4N z<$wj}6w^fl=W4%^pk6nu z^t0+d<7f2rywg~()XnsY&FoqZ57&+y{j6aq9QzUZA zOAD3!b8-}khie@fnV5e_fqGpTYt>k7e4X9YCZ-Z*6Pz1pc)=NSOfKtfG`At%k<W z*2qseOs!a(aSy;JB~LRg`!+Dr%&vWvGeG-~YdS5*-OLq`^|~^alCfI(!q)EE-!-x1 z$Ez>x&tl4*J_Qc9#^hkhkII%acV;nuvwGpA-Bm34VL=^Y$&d09sjnTVkn@h*6UqF8 zi7LdBA2l>nqO>DFtGUOTLwj)<*;Suv-}NTBPAB{;Fxq<_~5@HPWd=e#oXO059>*>jl#LW=#?c0T%6CZ%is7;_{IG zA8Kcv+z&{wqzuz_2>9`1=WF`@7_3sBWge$MIrFE=$WO@u0(tnesS5SFGS-pldD?Ah z$NJ?h8F!+F{ZwmfCam4tv+?bDOeaihcJs}2=4{jW{yo2Zd@%kZQpce~PBK!Qk*{UP z4)+Agzee>|A)lGmm?M+2JQ-y5E|ybZaTysOf;lAKZu1j~^vXWX$hogxaZpbFR&?M^ zd|n&TfPdfjwF>pRGM16C?CP-do#D@CSmx~dzy6P&Vwp3RmJ)N4jdXW-da$r2^O^aj z@o>vDEOXYLI^@TEXfHr4b2geIy>HD|A(lCtuTXwp*^IL6?=yyEZKt#nU~w5gR}`JX zIKJD@(8@39`&lFf-vuhfGG``?{P^)Nk=&aUszSZ4j8$a1UzpR58#)>ne>NpwJ9cFQ zb9`y$tZObX+m)aGtX}_+dEE}2wCEvbYrk@zI_7;IoTQFr{oQt&UgK9H;;dRfMBzxQx8$9@HvFPVGXivU9f%B54=0P$Bcm{T9d~??V5)(2zt4b3LxT=)2YZ3;|cWPOqdzFpMIOvW<9fzjJ zh*0))3}B>V=S2#c{;YtqVlZf+mcLv$Dj_pV+vx%y3HS$a6+cpea$t=DKv(gD8JQbXsFHc>stRQM%I^SI@vr)D zu(*uueg_03T*ZGWLEh_+2W5V&vyTXI75^SbrcVwB*eZ5+RH)an(#z;FaZw}teM3E% z`Mp8QyAB6RP~ilvgIP*WagL%BwRB9rCvMhXsf!uVj_6cSgP0*>W&dYKUn;y$Sd;cBYwnfP7uOwZju|r6 zX&py>YpbAR-Y(IK(>$_SQ$06HpqyXtx=4O*KL}u_;(n{JxQwiS8YoD-PrIukxv<_u zO^*hXRmv$>trYT`u8By#Y2m^_y{?P}WGu1+t@1+7eZoYzZT(+GSz;nwou97HZcK!0 z?&kjR5F0eIQu|TU-eBN-ynqh*Al-103_d9`(&1QhmCQRcjI#Wzkp;lXh?PBM*2$$% zg2iS0n(M7VY;`F?tJ4(fINO5dq0R>jqAetbkl z#Wd~G_MU4qnfHaoXFJS(rC4{^`u^_v^UTVs=Cm%(>C9+DRqBv;u02aLyn*%-=RJq1 zBHEduu8@=y6B{s6-Rel0^|aKsUU3<}U(Z1!&dmPA2vGK!`AgHzOyH3D%xjE6Io&5k zAroh(fd8L9|MxB9SErx!@a^(BBQeP9m+!kPXEDg@)t9Qb<1ol;jP;2>voOeORR46Z z^BCl{m5RY2ucn$})e0hIC*5Yv-=AAb$Qb0c9nc^zUj|vLMn4o-{Hx1pHi%l~kZ%v>xnxJaQ$a!EeYW-$pd7GenM68H+YfL$I_spC(?-Q>#^11FXxHn?Sg@uu zm(lX}cZ=Pu@s~Euyt3b5f|P-aCflB7BQ{puH-F|=Ovqto*+04$CT@GDqC>gH)tMq0 z;qEGs-p2(6VnU9S9OcK#cLDO}THRD3>)t^EC1vE?J1-*fA;!)U>9Kr?LQafnt3vtb zxCDVr$aLk%)82zcsMnRTUQLIawKo=|jh`w%Dc^q6^;Su2SXASICSrAk-+Zim#M(c)i)Un69RiqsHvl=744L(s;3|4GX zDW?x@sQJV(fg!WbQO$L(u1ih(5xoMpC9wuimy&+P5KGc%rCEJ?} z>UCwTR?{}s_Ka=eFGr=xqDk|oPVUQmMmrqrcs))TI=GTUWX;ix1t*WZo_~)SH&~}q z#|sthw4E1k4yAv9cg{7|%x`;wLzaK`JF1YMYkUNFhDaqxtc&Z?qnwk>UCwTRnu1NI^mJO zsV{4E?Cz>8tjpxQ=j_N2KfdB>OtnYLw_D(9%&HafPY&X0%qNjL9@M?)r5>@g38U;C z9VtOvjk!h2H+&2jm>ylS7QAH6}7a zB86tMNRFMRogwhQ_eCMnNTs5~c2e6@3bI~T#!5163}5b>_xR^`)}+FTlx6Q5wH2rs9?VfGtGnI< z19{up1QD_pZ*4hPyfuDW>=2On=`rRCx!}9CL<({y2jxshOHH+~yXKCoTU4~`b!99h z)4IXSu^+|gJjPuOS@lj>#RPB1u*a$A*r2>>wJ)dRDO-zSv!W9|F@tNr6?7=;nvq|aSXfpnPPl0iB7 z+&V2MId@0I2j?UkQm-px9hp}1CsaSW{!m5sf3WwSQB8bb+;8Z;s3;bcLF@=ZAPJ~s zHVBGXv0{%1N|$QEHVF2Hz3bR}SM1K%K@mF^EF&VQND&o9Mejb5fcO7HeeR2At@|RZ zrJsD~!m|E(;)U`)*~O$7trkbQ9%L+ zL=K{mTylgDPZ+y%03ok2Hd!G0)$NW7!qb6v6cSAQasXL5H$LtG9u+roN&^Z>j|#I@c<0iYoMu)*FF~{4Vg>`!fvXz;l1m1eKp4B)o@%L7zH^ZHg_2u=9$B7y0o#!VGSF6x8@PdIwRcY&~^LpUKX3L2*%x}a_%g9K@gh_Et^ zA9c3?DlXyhRDm#Q{Z<8Ge0T&!Sk`7eOIVucrXYM)=Pd%$_{VR^*miDt_t?@;n|YI> z-tlE4PD9jLzRtVUE%bt6U#9pDoCzXk+6B3GJZK#n&SOA2IWmhO46&F_5z2esRS@PL z+0P?+joQCf{gVBQhz6(NmlU!xjj!PcQIK<*=XpYAil36v%SI`AAeqL=c#Fyue~R!m z6$8jbuYdf6Ow&cX@@5+tt6ac< zxc$Ax6Gp`cGK4+1aSFnMd_4+@C+bXy#_h4vKAxysC!ki%4X;;#iXPLvltFSc7QN&N zEx4IJYnGd$3SSuKB5SPYbh*}*YI~Enk(V@A@I7f zxB?zonZ}nVehO4v^uotVjz4`LB|ZV$PpZ_boRmjrxH4N;icHU2mm?q(z1sYJ57#E_ ziT%bgXUsr%O|0+rZh;W|`g%W&r8}VevGT!+*L5IvQ-t|KJUgsq_6pk!i zmU|&SzQgyJotgZ7vS3W zrUF@+hGytR3Ki$Jp%@4kmMv21Y0bn+ezB;poLcC>pyncGN?F3ss}hy`9gH*sQh)5$8~r zpk%S4geSTv{~-|Gdo+zgCVEZ($Ag=m{Mc-o^mVe7Z&;MRqpkgQJl-?q>wJ(2@neFX$lE4i6x7y zT>HK~>`SK7Nv|6#kX+n=RGx4oWd=z4!}%^O(Is&q0^zHGDS*U0{_)_Z?259W(}q8H za^;%A-wgcTLi9HKD@)T;_!ci5I_eZ}h8ClqOmIHe60}ZPQW%k*;jY}<`=0$~?-K_heCVv3JwRRC;nS%@B^%w-gBP-|q)V zu54Qy1m@Aokw;=4|9Egy_VvOY9==n4KrPc&hPT#xGd7kvJIBMKwZQ=(cT@YrNwaacu z_ru38xP5GyC+E!#Z5dg4uKEGOh3pH4$rh!An79>kyUGBG9 z0Gcih0VC2wvb!pfTrxj^CyYI3MKrIm%AF-TXIU6UxL>46G&EVfMnG2P@lmiNik7a) zrU)k(9OVfWbDt5-edxAAAbg&+pCMdXS66|=JpS?MrU#jO%Qpx~bD_!ZG4vDe35;~T zva}OE2BUL7orz;R>Ox5Nv;9jaE(8(Vg(s1d+iheBr#R;;^{C44&Yybr#>pbOY&fIj z>$RH&d;vD*#w@ZjkM~7wiK64X)>E?o1p|uceun;lVRe4 ziRRv3YDW>4F77Q5KKM34fy6u}|9$<)+N8O7u5Q_7UPoi7d+(8xp+~{(Zw)m%@*RA4 zENQhU4uV#ied`$B6GRK1@Feog`7c={9lNnH(RAp_OG;fbtAZl>>&A~f;j6)`Si+C) zUWBUo(ffixrp7e*M5b~fZC4@idq&|zK+-QpI`M?>>a?K<@8zT-NFM)qdQ(=1uUqT2 z9^i>rHa}z&XW;`Sn`iA?^7H{-+1zWST_WDmv906D`SV-hmCc=a61o4(&kSL5$#{V< zhSO9a@ycd*9?4~G-oX$)D6T_@S2oLev8Y(4;ld0jg2NX+9Ok8jGV@ANvTvGq`1gYNJ0%wP$0f3`bi@t|RR#~a%kUFkL* zf<{bA>){;(qU1mlk(*T?2?=>S@MIN( z#4pbE2yt)PR)K^z>@fyena7u$B%-+J&Kv^A-+CMnhQ}GOgsYnu^MqxsW-FPu#f%X1 z_{YPP?ypWtdQe`^8tRN1Wz)Lf9lQ?Ja)s}}zIYw#-R@VChF-8` z2m?Ab;0eDvJVW5Q*}Mh<5`XXKI*Y_S{_!-WyGI%{@@~HCB;IPJ*f#m30&lhQ@J-(4 z3!bS>1AI%u@K&qDy|ORg@K&qI3`WFTt==k`cq>RCI&?)t1Uh&krjYdeh_yiYs__j0 z3HD*746-th4_+N9RNSzz9)xs4&(=I)Xvlh^x%>x*l=@8*4WjW@E9F{;dHmyXN_W{| z%Y8P#58xZ7j{kh6_Ed;l`90ghEf%kt-}^n>V?K0Bx#{BBx(;45@1=}9<;W@}r*>!Y zgYbY>$Dsn@*P*j0BwjNwVh9VSYVt^MbhlC_j=4)ZzfiLmz@Z#v2D?k@ z^Z0%z6$Sn}MW9PNs8*GL(CjT>07w3ZmZ*%K{|k@s9^8-D>~fRj|h> zC*EpChRd8nJ^WqbHFqT~UgGZ(7i=5yauB$A)v}s*n+DBY-zkhpN3_{QBIEB8XCv?+ zz?B*K`tBef-67%@SBb71>*!da4r8KB5 zcRzFBeIP`H^cY-cKn{S#+4RCD-yr_`vbJ`YeL(XkLn6`%4@xK`H?7@bmM|(OOv&%< zrwK%V(e~vDg|l7?B(!C)JG0+P#}G!O`%y;rU`CiuF6Xr;(3!1R)73zON#Z;g8x9W<(K7;o^v zz)gRxUWMKW-SfJ(Iy1EMwYzHV(=yYH*0`%7!7l$d{}LNcA8%o%rKhK-tBe2XS?lTJ z|2noIrWrH2^d1kswwzEPFWH|)p$XJUiV%Zi@AQuM4ocu!I&HTMJSy(!2C}3F#3P*+ zbCA3bWcMy92#aJ}uLt{ZeVrT~<9#HtBOxr>sUf!^Ulow({%ekmggRaGNC%zRJ`=yo zrB1G{5(g(2H&$yEb>+O8L)LXZ~ zXYD4Oe?>&IBLnc9Z3hLm(^p21CIs1uy$oT|L9H1E+`4La?y2rqU+hU#1qAn0IafaE zDbzI?Mmnn8xw^PII61nxxh2-09#exm>zI18r}ElEXKo*#87UX?Px9Cq^>iI3@5#;S zpET}&{NYxoiObS)_nY?wS*>6O+v#H$7AcS*{Sc>Q(Qe<3=f?duXFdCG8g|5OswyC$ z*4n859h&bEDC}6}&Qa{(>@3F0Mc3fXHl$wda$30AT*JeMmA`&-gu$+;J3YTX{#xGj ztCP$w>!g_=AM3}Ml3@B9^$4^+GJ9CRz}?75Dr zfPgx({iTi2H0%(M9e%HzlPj(i=x%IH?qW>Sqq)rOBY%8ZtwoDE2e8g^R9R1 zbMD{nc`)22RA2UO{f^DCWjMMe2-r?%WGHJ0vMd?VINHsc@RX~c19iHH+9v@ zld?w8;Q3)=tr1__=nj2;0%YCw0o&;dkt-P_$cAiI5EezY2?)5pquO1}+!3bYwB4$J zfZAr$rZ3Q7`BK0Rs@ge=ogAFSF3LF@QIoqE1A}qATYGuvhR5CHFW!%1u`|&JzQgQ% z=T=6g2_4eExQ_eo@HD#xNMF`gU^|`ay@^MH^m;ZB7QWG~AClxp+y%P~95fQ5Cwv;4 zcw7EYWRZQWA- zaQzk1ndnkk7pHEk&YSf0X*IHJ!@@=8A(Ki#di2+v-P81zfbRq8ub=LQ!n|O9kd5|l z?)El4e0XeMRX~W|+jQUVdB?y>afy6Dl{;sNgIMOQbQe~Wy9E7Dw|UnKN0)C17@!(D zi9U6H=i-J3_U(5R_QsVaHT7J1$$KnF_Yc6>I2NDzwKhmMb>vaw#pOZ8)=%I4=8hR0 zxaMegRX~Vl21~E^918W$J_PJg@eU@Maf^7D7WpX zLwnduHhXot*>nBb13vHX3_N=D?81d0N&K~y-7ldT;!3jnt&)YWLb5iEeD#~V^pT}L zv;0*7Av%3z(ZUWrp`Ggy9y?aK6XS^DD07jZyO5gP#nkeB0WPBsO@3uxC@(bMPCBFB zKl`ZIU;uNvzg760VZ}b5W=?AFF9w-@3~8s&o$wRL`(%bGEMd{GT^bW7U#fODDrIm* z=hs$M0Y{|_X^^*XIy5;lg~1M$?!->6ZnzzDav|$%a82&U=;!$JviwKI*M6*!KlFPf zU}vJ8!cF`24Qp_D=Au;t3Y*4juT8YZb++Mfz;^l|rCfYKI#}rz3iHR*KRdc*wY%sg zIWtD<6sZCN>LuPy7J$8bdlmy!)y@UCvSKHhEV04#F*Ufei7~M2%(u7p?_AqiU$vN^ z9(bv=d&Ab#miM>aW~+EJ^{Bo0j>}Mx^`>~@?sR@Qg9KUIMJ!?A+vLCtUtE84H%>cB zx4eNWAViPTW=0ks0xGXBV23JqQryZq%Sh=)*I>>z!CK3iZ)eoAt~BkfyzHW!d>HC} zp7-==JL|I!U8y{s*PFu+>V*v*2hyh>cx)a@YxeOSzb7Oj)0v} zxBr>Ctd_~oWJ_`XO+E7uX_b2I>irlbQ=(XGr%z@#CGP{to`wuz;j7Tqj^|9P-9@K& zy%~Msv??I36kT5)UiBE7x@%C_p{jIbrEqo?EE>4b8oUfy5sGMy`lc6nn+pNtuw?K|O z_W}>t8TA&bY(xJY7kf(l%QyAPUt=_CDCLz4(#fq!aQb-WBMJ%9eLsP)xY6#^i{p%| z-HnZS<-M`XNJJ@D!VX;Hie09X#VjPLd4oX0c?({!x!!>=Nax+UgLU<`FpY z$M-p(Pfq{1C)<+A$=-kJBi@mzJT5{XKB!!AakKWEC+>=y<*Z)4dQG*v$a-&2FAvXE z1q9URmUPsC)*T)T)$W{SxYO?FirZRrr>MzYOg-)WeAA|>!X5Ws^5Q5H3Of^B`mkb3 zn|AYxjMCXr4}6*)syiA&z}095X|LY(yyc4$hvLoY=`%`xs~vSCXwKDT(^LV0x-rr} zbRMnR^{lL&SmuJe#f}nJ;%-FE?wC`&U4}=P_TWeI_tzP~Pb2zEJNI5LAxZb-h8?1B z=kMg@^2kRZ`?W(%r+>>=5Ld1Pt|(cwTiY@Jr+;3w^cSiF4-$LXnGLxEyn`}0IXl!caWtd=aP(M%KO0SRZY@Ku3AGOo^3wxt?|2edNP5-k`mhvYXH!LXJ1GeUjt{K2alla|ZZ9 z?8DMGSMELs!{~_$jEYac;|aORPfk&UK_(*^B(G(g$rFm2Ody(j-@s0R1gGmFLQ&n_ zqu2jfl>&walPHwj^y2{n;UahC=*CTSoJ=&Y7xo%~$FEc7^Mv`58AM|yf4Hwkl}T0` zhs*rV4Gic8$rZ~>#zA!H^n>?zjpo}Nx9<|podTy6_Dlx;+?Da^0Rl#(r*a)gWID_` znIcqV%ayuZw2;LYynfR@6p}05bOwQE{QNHgd7;WgXMrt`l8b0?UqLus)K4G`iTXe^ z{bsZqi{!O@4G6h+PqF}snf&3_8dWBjTVMI&dT7Q&E{`sA*KdCvqMe&G&b_~ouhUOF zKWWiGpi(tPoT+6G#vU;k5RXND45X0UsKI@J&}GvzhVa$(T%vh{9*q%rtZ*;SBEhkD zR~}iFiP4>2MAC`IqXfcvPbi8o;ofZ}@3mNm665&f5{*AsXh}3?@`t-?RGFCc(%I2( z?~q`q7t*rPv$uGQ$hebpPmR&wy+<}*JSXWM_?&(kcvVjXM$4}$Fd`lL^Z-vdHSq;Q z7|!VONYL*~1ClFmHi9C2d@)ObggQ4R0j7NE36TC zUi`aq-g9>s`jhFH$scZZQDw5fb(bjr=Zn&K%crR?f==#%fGj?9@5l`3EPp0By{ZF5 zw8?()Bf}PSSIrd2CrrG4mm*{gYXjlfl8-Es*Mj8;JfD2doFRN{ev@d>8{AVtUQ{k! zbOp1PMM+2P>7yWw-)hej4sBY z6F5)Yp6%W@7vjcixprDy%zFo{cC+ml4L;|(tZBGc3+KsKibM|0X%B=`nWF+>__WbH z67>5lXOLWZ%K!w-#6X!b&XZ6Dnf|fd zLLj_5XEck%O#X1oiz<`+lG)4JUi5RtOeC$c9q(c$n#&U{>tQA<;+WOjFcXi z3o#hRpvFv;{V)?NinQOIqg+2RlRw<`qROP#qwK;C9`i-KWySSgvx>YSAbp7}tCa#e zTXjFXrfo1pSp9t6{>N9)U1y>oktZs59dYSu9Hi7yTGJ_<&T9>>%OdHrsmjrpezf=y z(fEr1{Q#-TgszyD8sqA`;{-20-+Z%salE4$0_NYGe(S3uHl{i6xN zq2+A~smg@@IWCt)$@Mv=-15;$L*q#^9h$gLL8iZ0D^=?6A!itTLFXo2Qy?*uKivMJ z%4CnE*h|YoF=CBDqq>Yzm2){NWB5RVL|I`C9Ck=eV9+7+K2XOosqdqp+;k3E