diff --git a/.gitignore b/.gitignore index 6ec1572..2b95587 100644 --- a/.gitignore +++ b/.gitignore @@ -158,4 +158,7 @@ Thumbs.db *.swp # if examples are run -examples/*.eqx \ No newline at end of file +examples/*.eqx + +test/ +.vscode/ diff --git a/examples/model.py b/examples/model.py index f478eee..1fad03b 100644 --- a/examples/model.py +++ b/examples/model.py @@ -6,7 +6,7 @@ from dilax.model import Model, Result from dilax.optimizer import JaxOptimizer -from dilax.parameter import Parameter, lnN, modifier, unconstrained +from dilax.parameter import Parameter, compose, lnN, modifier, shape, unconstrained from dilax.util import HistDB @@ -18,29 +18,44 @@ def __call__( ) -> Result: res = Result() + mu_modifier = modifier( + name="mu", parameter=parameters["mu"], effect=unconstrained() + ) res.add( process="signal", - expectation=modifier( - name="mu", - parameter=parameters["mu"], - effect=unconstrained(), - )(processes["signal"]), + expectation=mu_modifier(processes["signal", "nominal"]), + ) + + bkg1_modifier = compose( + modifier(name="lnN1", parameter=parameters["norm1"], effect=lnN(0.1)), + modifier( + name="shape1_bkg1", + parameter=parameters["shape1"], + effect=shape( + up=processes["background1", "shape_up"], + down=processes["background1", "shape_down"], + ), + ), ) res.add( - process="background", - expectation=modifier( - name="lnN1", - parameter=parameters["norm1"], - effect=lnN(0.1), - )(processes["background1"]), + process="background1", + expectation=bkg1_modifier(processes["background1", "nominal"]), + ) + + bkg2_modifier = compose( + modifier(name="lnN2", parameter=parameters["norm2"], effect=lnN(0.05)), + modifier( + name="shape1_bkg2", + parameter=parameters["shape1"], + effect=shape( + up=processes["background2", "shape_up"], + down=processes["background2", "shape_down"], + ), + ), ) res.add( process="background2", - expectation=modifier( - name="lnN2", - parameter=parameters["norm1"], - effect=lnN(0.05), - )(processes["background2"]), + expectation=bkg2_modifier(processes["background2", "nominal"]), ) return res @@ -48,15 +63,20 @@ def __call__( def create_model(): processes = HistDB( { - "signal": jnp.array([3]), - "background1": jnp.array([10]), - "background2": jnp.array([20]), + ("signal", "nominal"): jnp.array([3]), + ("background1", "nominal"): jnp.array([10]), + ("background2", "nominal"): jnp.array([20]), + ("background1", "shape_up"): jnp.array([12]), + ("background1", "shape_down"): jnp.array([8]), + ("background2", "shape_up"): jnp.array([23]), + ("background2", "shape_down"): jnp.array([19]), } ) parameters = { "mu": Parameter(value=jnp.array([1.0]), bounds=(0.0, jnp.inf)), - "norm1": Parameter(value=jnp.array([0.0]), bounds=(-jnp.inf, jnp.inf)), - "norm2": Parameter(value=jnp.array([0.0]), bounds=(-jnp.inf, jnp.inf)), + "norm1": Parameter(value=jnp.array([0.0])), + "norm2": Parameter(value=jnp.array([0.0])), + "shape1": Parameter(value=jnp.array([0.0])), } # return model @@ -69,6 +89,7 @@ def create_model(): init_values = model.parameter_values observation = jnp.array([37]) +asimov = model.evaluate().expectation() # create optimizer (from `jaxopt`) diff --git a/examples/nll_profiling.py b/examples/nll_profiling.py index 879501c..170eaca 100644 --- a/examples/nll_profiling.py +++ b/examples/nll_profiling.py @@ -5,12 +5,12 @@ import equinox as eqx import jax import jax.numpy as jnp +from examples.model import asimov, model, optimizer from jax.config import config from dilax.likelihood import NLL from dilax.model import Model from dilax.optimizer import JaxOptimizer -from model import model, observation, optimizer config.update("jax_enable_x64", True) @@ -21,15 +21,17 @@ def nll_profiling( model: Model, observation: jax.Array, optimizer: JaxOptimizer, + fit: bool, ) -> jax.Array: # define single fit for a fixed parameter of interest (poi) - @partial(jax.jit, static_argnames=("value_name", "optimizer")) + @partial(jax.jit, static_argnames=("value_name", "optimizer", "fit")) def fixed_poi_fit( value_name: str, scan_point: jax.Array, model: Model, observation: jax.Array, optimizer: JaxOptimizer, + fit: bool, ) -> jax.Array: # fix theta into the model model = model.update(values={value_name: scan_point}) @@ -37,19 +39,29 @@ def fixed_poi_fit( init_values.pop(value_name, 1) # minimize nll = eqx.filter_jit(NLL(model=model, observation=observation)) - values, _ = optimizer.fit(fun=nll, init_values=init_values) + if fit: + values, _ = optimizer.fit(fun=nll, init_values=init_values) + else: + values = model.parameter_values return nll(values=values) # vectorise for multiple fixed values (scan points) - fixed_poi_fit_vec = jax.vmap(fixed_poi_fit, in_axes=(None, 0, None, None, None)) - return fixed_poi_fit_vec(value_name, scan_points, model, observation, optimizer) + fixed_poi_fit_vec = jax.vmap( + fixed_poi_fit, in_axes=(None, 0, None, None, None, None) + ) + return fixed_poi_fit_vec( + value_name, scan_points, model, observation, optimizer, fit + ) # profile the NLL around starting point of `0` -profile = nll_profiling( - value_name="norm2", - scan_points=jnp.array([-0.3, -0.2, -0.1, 0.0, 0.1, 0.2, 0.3]), +scan_points = jnp.r_[-1.9:2.0:0.1] + +profile_postfit = nll_profiling( + value_name="norm1", + scan_points=scan_points, model=model, - observation=observation, + observation=asimov, optimizer=optimizer, + fit=True, ) diff --git a/examples/nuisance_parameter.ipynb b/examples/nuisance_parameter.ipynb new file mode 100644 index 0000000..5dab6d5 --- /dev/null +++ b/examples/nuisance_parameter.ipynb @@ -0,0 +1,181 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [], + "source": [ + "%matplotlib widget\n", + "\n", + "import ipywidgets as widgets\n", + "import matplotlib.pyplot as plt\n", + "\n", + "import numpy as np\n", + "\n", + "from dilax.pdf import Gauss" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [ + { + "data": { + "application/vnd.jupyter.widget-view+json": { + "model_id": "a5864df1a5ce4280bddc2c8d5331b147", + "version_major": 2, + "version_minor": 0 + }, + "text/plain": [ + "interactive(children=(FloatSlider(value=0.0, description='nuisance', max=4.0, min=-4.0, step=0.01), Output()),…" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "application/vnd.jupyter.widget-view+json": { + "model_id": "95ee9277233a4fa287b148c5fd42137d", + "version_major": 2, + "version_minor": 0 + }, + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAoAAAAHgCAYAAAA10dzkAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAACmsElEQVR4nOzdd1xV9f/A8de9l703iKCguJCluM2RUZRlatowc2VmmpWRZZajMnNkpplpWZmapg1H3zTTKHPmQHFvQRRFhrLHhXvv7w/y/iJBQYED3Pfz8bgPuJ97zue8D14vbz5TZTAYDAghhBBCCJOhVjoAIYQQQghRvSQBFEIIIYQwMZIACiGEEEKYGEkAhRBCCCFMjCSAQgghhBAmRhJAIYQQQggTIwmgEEIIIYSJkQRQCCGEEMLESAIohBBCCGFiJAEUQgghhDAxkgAKIYQQQpgYSQCFEEIIIUyMJIBCCCGEECbGTOkAagu9Xs/ly5ext7dHpVIpHY4QQgghTJjBYCArKwtvb2/U6oq350kCWE6XL1/G19dX6TCEEEIIIYwuXryIj49Phc+TBLCc7O3tgeIftIODg8LRCCGEEMKUZWZm4uvra8xPKkoSwHK60e3r4OAgCaAQQgghaoQ7HZZWayeBLFiwAD8/P6ysrGjfvj179+4t13mrVq1CpVLRp0+fqg1QCCGEEKKGqpUJ4OrVq4mKimLKlCkcOHCA0NBQIiMjSU5OvuV58fHxjBs3ji5dulRTpEII8R/790OPHsVfhRBCIbUyAZwzZw4jRoxg2LBhBAYGsmjRImxsbPj666/LPEen0zFw4EDeffddGjVqVI3RCiHEvyxbBn/+CcuXKx2JEMKE1boxgFqtlpiYGCZMmGAsU6vVREREsHv37jLPe++99/Dw8GD48OFs3769OkIVQohiFy5AaiqoVLB6dXHZqlUwZAgYDODmBg0bKhtjLaDT6SgsLFQ6DCGqhbm5ORqNpsrqr3UJYGpqKjqdDk9PzxLlnp6enDx5stRzduzYwVdffUVsbGy5r1NQUEBBQYHxeWZm5h3FK4QwTdoiPSeTMjl0KYNBHf2M5QZABRiSk1GFhxvLDyVcp3k9eyzNqu4Dv7YyGAwkJSWRnp6udChCVCsnJye8vLyqZP3hWpcAVlRWVhaDBg1i8eLFuLm5lfu86dOn8+6771ZhZEKIuiYzv5Atx66y+XgS206nkleoA2D/I68xe+NczPU6bnyM3/haqNYwrudY1i/YiZW5mi5N3Hkg0JMHAr1wtDFX5D5qmhvJn4eHBzY2NrIYv6jzDAYDubm5xrkN9erVq/RrqAwGg6HSa61CWq0WGxsbfvzxxxIzeYcMGUJ6ejrr168vcXxsbCytWrUq0Yyq1+uB4q7jU6dO0bhx45uuU1oLoK+vLxkZGbIMjBCihOOXM1n+dzzrDl42Jn0ANmbmNLB1or6tAyHJZ3llWq+bzl06/0f+sPHl8KV0ruf+f/emlbma3qH1GdSxIUH1HavlPmoinU7H6dOn8fDwwNXVVelwhKhWaWlpJCcn07Rp05u6gzMzM3F0dLzjvKTWtQBaWFgQHh5OdHS0MQHU6/VER0czZsyYm45v3rw5R44cKVE2ceJEsrKymDdvXpm7e1haWmJpaVnp8Qsh6o4TVzL5aPNpfj9x1VjmaW1LqIs37bw9ae7pgLV1cWuV9clcAAwqNSqDHoNajUqvZ0gnf4a0bo3BYOD4lUy2HL/KxiNXOH01m9X7L7J6/0V6NPcg6v6mJpkI3hjzZ2Njo3AkQlS/G+/7wsLCSh8PWOsSQICoqCiGDBlCmzZtaNeuHXPnziUnJ4dhw4YBMHjwYOrXr8/06dOxsrIiKCioxPlOTk4AN5ULIUR5pGQVMH3jCdYcTASKu3NDXerR3achbRu6YGV1cxdlkbMHha5e5Lv7knD/cJrv/ArN5Yvg4VFch0pFS29HWno78sp9Tdh/4TrLdl9gw+HL/HEymT9OJvNoqDcTH26Bh4NVdd5ujSDdvsIUVeX7vlYmgE8++SQpKSlMnjyZpKQkwsLC2LRpk3FiSEJCwh1tjCyEELei1xtYsTeBWZtOkpVfBEArl3r0CWhKCx87bvUHeqGnD0f+F0+ezoKsbBV+HzyPrbkWSulpUKlUtPVzoa2fC69GNGFe9Bl+PnSZnw9d5s+Tybz2QFMGd/RDrZakSAhxZ2rdGECl3G1fuxCidkvOzOe1Hw6x/UwqAL62DjzVJJg2/k6YVeBP6fx8yMqCLl3A1rb85x1NzODttUc4dCkDgI6NXPnoiVC8nawrchu1Tn5+PnFxcfj7+2NlVTtbPpOSkpg+fTobNmzg0qVLODo6EhAQwDPPPMOQIUOke1uU6Vbvf5MbAyiEENUt+sRVXv/xMNdytFio1fT0aU6/ED9srKuvBS6oviNrRndm5Z4LfLDxJLvPp/Hg3G3M7BfCQ8GVP0NQVI7z58/TuXNnnJyc+OCDDwgODsbS0pIjR47wxRdfUL9+fR599FGlwxQmSPpJhRCiDHq9gdm/nWL40v1cy9FS39aBN1vfwzPt/Ks1+btBo1YxqKMfG1/pQqiPI5n5RYxacYAPfzuJTi+dOTXR6NGjMTMzY//+/TzxxBO0aNGCRo0a0bt3bzZs2ECvXsUzw+fMmUNwcDC2trb4+voyevRosrOzjfW88847hIWFlah77ty5+Pn5GZ9v3bqVdu3aYWtri5OTE507d+bChQsAHDp0iHvvvRd7e3scHBwIDw9nv2xHaNIkARRCiFLkFBTxwrcxfPrnWQC6ePrx3j2dCPW3R+n5CP5utvw4qhPP3eMPwII/zzFi2X4y82WXjJokLS2NzZs38+KLL2JbRn//jUH+arWaTz75hGPHjrF06VL++OMP3njjjXJfq6ioiD59+tCtWzcOHz7M7t27ef755431Dxw4EB8fH/bt20dMTAxvvvkm5uayzqQpky5gIYT4jysZeQxbso+TSVmYqdQ80SiY3iE+FRrrV9XMNWomPhJIy/oOvPnTEf44mUy/z3axbHg76jnW7XGBBoOhxHqL1cnaXFPumZlnz57FYDDQrFmzEuVubm7k5+cD8OKLLzJz5kzGjh1rfN3Pz4/333+fF154gc8++6xc18rMzCQjI4NHHnnEuLZtixYtjK8nJCTw+uuv07x5cwCaNGlSrnpF3VWDPs6EEEJ5cak5PPPlHhLT83CwsOS5ZuF0auaseKtfWfq28iHA3Z7nlu3jTHL2P0lgewI87JQOrcrkFeoInPybItc+/l4kNhZ396tz79696PV6Bg4caNxw4Pfff2f69OmcPHmSzMxMioqKyM/PJzc3t1yTRFxcXBg6dCiRkZHcf//9RERE8MQTTxh3kIiKiuK5555j+fLlRERE8Pjjj5e6CYIwHdIFLIQQ/zh2OYPHF+0iMT0PD2tb3gzvROfmNTf5uyHYx5GfRnWikbstlzPyeXzRLmIvpisdlskLCAhApVJx6tSpEuWNGjUiICAAa+viltr4+HgeeeQRQkJC+Omnn4iJiWHBggVA8e5XUNxF/N9FO24skn3DkiVL2L17N506dWL16tU0bdqUv//+GygeQ3js2DEefvhh/vjjDwIDA1m7dm2V3LeoHaQFUAghgNiL6Qz6ag9Z+UX42DrwWpt2NPCoPbsB+Tjb8OMLnRi2ZC+HLmXwzJd7WDa8Ha0bOCsdWqWzNtdw/L1Ixa5dXq6urtx///18+umnvPTSS2WOA4yJiUGv1/PRRx8Z17D9/vvvSxzj7u5OUlISBoPB2AUdGxt7U12tWrWiVatWTJgwgY4dO7Jy5Uo6dOgAQNOmTWnatCmvvvoqAwYMYMmSJfTt27fc9yPqFmkBFEKYvKOJGcbkr5G9M2+261Crkr8bXGwtWDmiAx0auZBdUMSQr/bWyZZAlUqFjYWZIo+K7szw2WefUVRURJs2bVi9ejUnTpzg1KlTfPvtt5w8eRKNRkNAQACFhYXMnz+f8+fPs3z5chYtWlSinu7du5OSksKsWbM4d+4cCxYs4NdffzW+HhcXx4QJE9i9ezcXLlxg8+bNnDlzhhYtWpCXl8eYMWPYunUrFy5cYOfOnezbt6/EGEFheiQBFEKYtBNXMnnmy/9P/t5o3456brV3dqStpRlfD21LO38XsgqKGPTVHg7VwSSwtmjcuDEHDx4kIiKCCRMmEBoaSps2bZg/fz7jxo1j6tSphIaGMmfOHGbOnElQUBArVqxg+vTpJepp0aIFn332GQsWLCA0NJS9e/cybtw44+s2NjacPHmSfv360bRpU55//nlefPFFRo4ciUajIS0tjcGDB9O0aVOeeOIJHnroId59993q/nGIGkR2Aikn2QlEiLrnXEo2TyzaTVqOloZ2Toxv3w5Pl6pN/u50J5CKyikoYuiSveyLv46TjTk/vtCpVk4MqQs7gQhxp6pyJxBpARRCmKTkrHyGfL2XtBwtvrYOvN626pO/6mRracaSYe0I9XEkPbeQIV/vJSkjX+mwhBA1hCSAQgiTk11QxLAl+7h0PQ93KxteDa/d3b5lsfunO9jfzZbE9DyGfL2XjDxZLFoIIQmgEMLEaIv0jPo2hmOXM7E3t+CVsHY09Kx9Ez7Ky9XOkmXPtsPd3pJTV7MYsXQ/+QotoiyEqDkkARRCmAyDwcBba4+w/UwqlhoNLwS2pYVvFQ7EqyF8XWxYOqwd9pZm7I2/xoQ1R25aU04IYVokARRCmIwvt8fxY8wl1CoVQwJa076Jk9IhVZtAbwcWDQpHo1ax9mAii/46r3RIQggFSQIohDAJf55KZvqvJwDo7duCiJYeNX6Hj8rWOcCNd3oFAjDrt5NsOX5V4YiEEEqRBFAIUeedTc7m5ZUH0RugvbsvT7TyQ1P+DR3qlEEd/RjUoSEGA4xddZCTSZlKhySEUIAkgEKIOi0jr5Dnl+0nq6B4oecX2gRhYWFiTX//MblXIJ0au5Kj1TFi2X4ycmVmsBCmRhJAIUSdZTAYeP2HQ5xPzcHZ0ooxrcJxsJOPPXONms8GtsbXxZqL1/J47YdY9HqZFCKEKZFPQiFEnfXl9jg2H7+KmVrNc83C6/RyLxXlZGPBwoHhWJip+f1EMou2nVM6JHELQ4cOpU+fPkqHIeoQSQCFEHXSvvhrzNh0EoA+DQJpZ0IzfssrqL4j7z3aEoDZv51i17lUhSOqW4YOHYpKpWLGjBklytetW4eqgjOQ5s2bxzfffFOJ0YnK4ufnx9y5c5UOo8IkARRC1Dmp2QWMWXkAnd5Aa1dvHgtpgFo+7Ur1ZFtf+of7oDfAy98dlO3iKpmVlRUzZ87k+vXrd1WPo6MjTk5OlRNUHWIwGCgqKlI6jEqh1Wqr9XrykSiEqFN0egNjV8VyNbMAL2s7RrYOxtLStCd93IpKpWJq7yCae9mTmq3l5e8OUqTTKx1W1dm/H3r0KP5aDSIiIvDy8mL69OllHvPOO+8QFhZWomzu3Ln4+fkZn/+3C/jHH38kODgYa2trXF1diYiIICcnB4B9+/Zx//334+bmhqOjI926dePAgQMl6lepVHz55Zf07dsXGxsbmjRpws8//1zimGPHjvHII4/g4OCAvb09Xbp04dy5/x8q8OWXX9KiRQusrKxo3rw5n3322S1/Ft27d2fMmDGMGTMGR0dH3NzcmDRpUolFyZcvX06bNm2wt7fHy8uLp59+muTkZOPrW7duRaVS8euvvxIeHo6lpSU7duzg3Llz9O7dG09PT+zs7Gjbti2///57iev7+fnx/vvvM3jwYOzs7GjYsCE///wzKSkp9O7dGzs7O0JCQtj/n/fGjh076NKlC9bW1vj6+vLyyy8bf9bdu3fnwoULvPrqq6hUqhItu7c670Y8U6dOZfDgwTg4OPD888/f8udX2SQBFELUKQu3nmXH2VQsNBpGBLbGzclM6ZBqPGsLDYueCcfun51CFvxZh8cDLlsGf/4Jy5dXy+U0Gg0ffPAB8+fP59KlS5VS55UrVxgwYADPPvssJ06cYOvWrTz22GPGRCorK4shQ4awY8cO/v77b5o0aULPnj3JysoqUc+7777LE088weHDh+nZsycDBw7k2rVrACQmJtK1a1csLS35448/iImJ4dlnnzW2tq1YsYLJkyczbdo0Tpw4wQcffMCkSZNYunTpLWNfunQpZmZm7N27l3nz5jFnzhy+/PJL4+uFhYVMnTqVQ4cOsW7dOuLj4xk6dOhN9bz55pvMmDGDEydOEBISQnZ2Nj179iQ6OpqDBw/y4IMP0qtXLxISEkqc9/HHH9O5c2cOHjzIww8/zKBBgxg8eDDPPPMMBw4coHHjxgwePNj4szx37hwPPvgg/fr14/Dhw6xevZodO3YwZswYANasWYOPjw/vvfceV65c4cqVK+U674bZs2cTGhrKwYMHmTRp0i1/dpVNZZD9gMolMzMTR0dHMjIycHBwUDocIUQpDiZcp/+i3ej0BgY0CqVfa58at9hzfj5kZUGXLmBbw3ahW3cwkbGrY1Gr4PuRHWnj56J0SOTn5xMXF4e/vz9WVlbFhQYD5OaWv5KEBEhLA5UK+vaFlBRwd4e1a4vrcnWFBg3KV5eNDeV9Uw0dOpT09HTWrVtHx44dCQwM5KuvvmLdunX07dvXmGS88847rFu3jtjYWOO5c+fOZe7cucTHx99U14EDBwgPDyc+Pp6GDRveNg69Xo+TkxMrV67kkUceAYpbACdOnMjUqVMByMnJwc7Ojl9//ZUHH3yQt956i1WrVnHq1CnMzc1vqjMgIICpU6cyYMAAY9n777/Pxo0b2bVrV6lxdO/eneTkZI4dO2ZsKXvzzTf5+eefOX78eKnn7N+/n7Zt25KVlYWdnR1bt27l3nvvZd26dfTu3fuW9x0UFMQLL7xgTLr8/Pzo0qULy/9J/pOSkqhXrx6TJk3ivffeA+Dvv/+mY8eOXLlyBS8vL5577jk0Gg2ff/65sd4dO3bQrVs3cnJysLKyws/Pj7FjxzJ27FjjMeU9r1WrVqxdu7bMeyj1/f+Pu81L5E9jIUSdkF1QxCurYtHpDbRy9aZXUP0al/zVdH1a1eev0ymsPZjIK6ti2fhKFxytb/7lr7jcXLCzu7s6UlLgnnsqfl529h1l7jNnzqRHjx6MGzeu4tf8j9DQUO677z6Cg4OJjIzkgQceoH///jg7OwNw9epVJk6cyNatW0lOTkan05Gbm3tTa1hISIjxe1tbWxwcHIzdrbGxsXTp0qXU5C8nJ4dz584xfPhwRowYYSwvKirC0dHxlrF36NChRDdpx44d+eijj9DpdGg0GmJiYnjnnXc4dOgQ169fR68vHo6QkJBAYGCg8bw2bdqUqDc7O5t33nmHDRs2cOXKFYqKisjLy7vlPXt6egIQHBx8U1lycjJeXl4cOnSIw4cPs2LFCuMxBoMBvV5PXFwcLVq0KPU+y3vef++jOkkCKISoE6asP0bCtVxcLK0ZHhIk4/7u0Hu9WxJz4ToJ13J5e+0R5g9oVeEZq+JmXbt2JTIykgkTJtzUpalWq/lvZ1xhYdmLc2s0GrZs2cKuXbvYvHkz8+fP5+2332bPnj34+/szZMgQ0tLSmDdvHg0bNsTS0pKOHTveNMngv8mdSqUyJlzW1tZlXj87OxuAxYsX0759+5tiu1M5OTlERkYSGRnJihUrcHd3JyEhgcjIyJtit/1PEj5u3Di2bNnC7NmzCQgIwNramv79+9/ynm+8r0sru/FzyM7OZuTIkbz88ss3xdvgFq3G5T3vv/dRnSQBFELUej8fusxPBy6hAgY3DaOeWw1staol7K3MmfdUGP0X7eaXw1fo1tSdx9v4Kh1WSTY2xS1xFREbW3qL344d8J8JGLe99h2aMWMGYWFhNGvWrES5u7s7SUlJGAwGYwLy7+7g0qhUKjp37kznzp2ZPHkyDRs2ZO3atURFRbFz504+++wzevbsCcDFixdJTa3YEj8hISEsXbqUwsLCmxJFT09PvL29OX/+PAMHDqxQvXv27Cnx/MYYRY1Gw8mTJ0lLS2PGjBn4+ha/5/47IaMsO3fuZOjQofTt2xcoTsBudJ/fjdatW3P8+HECAgLKPMbCwgKdTlfh85Qmk0CEELXapevFLVUAEd4BdGqi/Li12q5VA2ei7m8KwJSfjxGfmnObM6qZSlXcDVuRx40WrRvrAd34am1dsXruojU0ODiYgQMH8sknn5Qo7969OykpKcyaNYtz586xYMECfv311zLr2bNnDx988AH79+8nISGBNWvWkJKSYuxWbNKkCcuXL+fEiRPs2bOHgQMH3rJFrzRjxowhMzOTp556iv3793PmzBmWL1/OqVOngOIJJNOnT+eTTz7h9OnTHDlyhCVLljBnzpxb1puQkEBUVBSnTp3iu+++Y/78+bzyyitAccuYhYUF8+fP5/z58/z888/GMYq306RJE9asWUNsbCyHDh3i6aefNrbi3Y3x48eza9cuxowZQ2xsLGfOnGH9+vUlJnP4+fmxbds2EhMTjYl2ec5TmiSAQohaS683ELX6EFn5RfjZO/FMWBPuogdK/MsL3RrToZELuVodr/1wCF1t3yrOwwO8vCA8HBYtKv7q5VVcXo3ee++9mxKTFi1a8Nlnn7FgwQJCQ0PZu3fvLccKOjg4sG3bNnr27EnTpk2ZOHEiH330EQ899BAAX331FdevX6d169YMGjSIl19+GY8K3qerqyt//PEH2dnZdOvWjfDwcBYvXmxsDXzuuef48ssvWbJkCcHBwXTr1o1vvvkGf3//W9Y7ePBg8vLyaNeuHS+++CKvvPKKcfkTd3d3vvnmG3744QcCAwOZMWMGs2fPLle8c+bMwdnZmU6dOtGrVy8iIyNp3bp1he65NCEhIfz111+cPn2aLl260KpVKyZPnoy3t7fxmPfee4/4+HgaN26Mu7t7uc9TWq2dBbxgwQI+/PBDkpKSCA0NZf78+bRr167UY9esWcMHH3zA2bNnKSwspEmTJrz22msMGjSo3NeTWcBC1Dxfbj/P+xtOYKnRMLFNVwIb3Hn3XHWpybOA/+vS9VwenLud7IIiJjzUnJHdGld7DLeaBVlhBQVgYVHcimcwgFYLlrI9YHXp3r07YWFhtXLXDKVU5SzgWtkCuHr1aqKiopgyZQoHDhwgNDSUyMjIEotF/puLiwtvv/02u3fv5vDhwwwbNoxhw4bx22+/VXPkQojKci4lmw9/K+6O6uXbgha+NT/5q218nG2Y/EjxzMuPNp/m9NWs25xRw1la/n8XrkolyZ8wabUyAZwzZw4jRoxg2LBhBAYGsmjRImxsbPj6669LPb579+707duXFi1a0LhxY1555RVCQkLYsWNHNUcuhKgMRTo9r31/iIIiPc0d3egb3ECWfKkij7fxoUdzD7Q6PVHfx1JYl3cJEcKE1LoEUKvVEhMTQ0REhLFMrVYTERHB7t27b3u+wWAgOjqaU6dO0bVr16oMVQhRRb7Yfp7Yi+lYa8wYFhSClZVkf1VFpVIx47FgnGzMOZqYyYI/zyodkqiltm7dKt2/NUitSwBTU1PR6XTGxRpv8PT0JCkpqczzMjIysLOzw8LCgocffpj58+dz//33l3l8QUEBmZmZJR5CCOWdSspi7pYzAPRpGEiAd8VmN4qK83CwYmrvIAA+/eMsRy5lKByREOJu1boE8E7Z29sTGxvLvn37mDZtGlFRUWzdurXM46dPn46jo6PxcWNNIiGEcgr/6YbU6vS0dPKgV5CP0iGZjF6h3jwcUo8ivYGo72PJL9Td/iQhRI1V6xJANzc3NBoNV69eLVF+9epVvLy8yjxPrVYTEBBAWFgYr732Gv3792f69OllHj9hwgQyMjKMj4sXL1baPQgh7sxnf57j2OVMbMzMGRYcLLt9VLOpvYNws7PkTHI2n0SfqdZr19IFK4S4K1X5vq91CaCFhQXh4eFER0cby/R6PdHR0XTs2LHc9ej1egoKCsp83dLSEgcHhxIPIYRyjl/OZP4fxUlHP7+W+Hvd5ZIgosJcbC2Y1re4K/jzbec5mlj1XcE31p3Lzc2t8msJUdPceN+Xtifz3aqVW8FFRUUxZMgQ2rRpQ7t27Zg7dy45OTkMGzYMKF5osn79+sYWvunTp9OmTRsaN25MQUEBGzduZPny5SxcuFDJ2xBClJNOb+DNNYcp0hsIcfbkocCas5iqqYls6cXDIfXYcPgKb/x4mPVjOmOuqbq2BI1Gg5OTk3GZLxsbG9mbWNR5BoOB3NxckpOTcXJyuqs9lstSKxPAJ598kpSUFCZPnkxSUhJhYWFs2rTJODEkISEBtfr/P5BycnIYPXo0ly5dwtramubNm/Ptt9/y5JNPKnULQogKWLIzjsOXMrA2M2NQyyDp+lXYO71asvNsKsevZPLFtvO8eG/V7nd6Y3hPWWu9ClFXOTk53XJ4292otTuBVDfZCUQIZVy8lssDH28jr1DH4/7BPNWmgdIh3ZXatBPIraw9eIlXVx/CQqNm4ytdCPCwq/Jr6nQ6CgsLq/w6QtQE5ubmt2z5u9u8pFa2AAohTIPBYGDS+qPkFepoZO9C7yCZjV9T9Amrz/rYy2w9lcL4nw7z/ciOaNRV2zKr0WiqpCtMCFNU6yaBCCFMx8+HihMMM5WaIS2CsZYFn2sMlUrFB32DsbXQEHPhOst2xysdkhCiAiQBFELUSNdztLz3v+MA3F8/gJYNqr6LUVSMt5M1b/ZsAcCsTae4eE1m6gpRW0gCKISokaZtPEFajhYvazueCm0se/3WUAPbNaCdvwt5hTomrDki6/UJUUtIAiiEqHF2nU3lx5hLqIBnmgVjZyMfVTWVWq1iZr8QLM3U7Dibyg8xl5QOSQhRDvKpKoSoUfILdUxYewSATh4Nad/YReGIxO34u9ny6v1NAfhg4wnSssteZF8IUTNIAiiEqFHmRZ/hQlouThZWDA5thlo+pWqF4ff406KeA+m5hby/4YTS4QghbkM+WoUQNcbxy8ULCwM83rglbk6Vv/2RqBrmGjUzHgtGpYK1BxPZdjpF6ZCEELcgCaAQokbQ6Q1MWHMYnd5AsLMXEc2rZvV7UXVCfZ0Y0tEPgLfXHSFPq1M2ICFEmSQBFELUCEt3xXPon+3ehga3xEyWqa+VxkU2o56jFRev5TEv+ozS4QghyiAJoBBCcZeu5zJ78ykAejVojp+nlcIRiTtlZ2nGe72DAFi8/TzHL2cqHJEQojSSAAohFGUwGJi07ii5Wh2N7J3pE1y79/oVcH+gJw8FeRV36689gk4vawMKUdNIAiiEUNT/Dl/hzxvbvQUGY2khKz7XBe882hJ7SzMOXUxnuWwTJ0SNIwmgEEIx6bla3v3fMQDu825MS197hSMSlcXTwYo3HmoOwIe/neJyep7CEQkh/k0SQCGEYooXDdbiaW3HANnurc4Z2K4B4Q2dydHqmLz+mGwTJ0QNIgmgEEIRu86l8v3+4m3DBjYJxt5Wo3BEorKp1So+6BuMmVrF7yeu8tuxJKVDEkL8QxJAIUS1yy/U8daaG9u9NaBjE9nura5q5mXPC90aAzDl52Nk5hcqHJEQAiQBFEIoYP4fZ4hPy8XRwpJBIc1lu7c6bkyPAPxcbbiaWcCHm04pHY4QAkkAhRDV7MSVTD7/q3i7t35+QXg4y3ZvdZ2VuYYP+gYD8O2eC8RcuK5wREKIakkACwsLuXjxIqdOneLatWvVcUkhRA2k0xt4c80RivQGgp09iWwp272Zik4BbvQP98FggLfWHEFbpFc6JCFMWpUlgFlZWSxcuJBu3brh4OCAn58fLVq0wN3dnYYNGzJixAj27dtXVZcXQtRAy3bHc+hiOtYaMwa3DJLt3kzM2z1b4GJrwamrWXyx7ZzS4Qhh0qokAZwzZw5+fn4sWbKEiIgI1q1bR2xsLKdPn2b37t1MmTKFoqIiHnjgAR588EHOnJH9IoWo6xLT8/jwt+LxXw/7NsffS7Z7MzXOthZMeqQFAJ/8cZbzKdkKRySE6aqSv7/37dvHtm3baNmyZamvt2vXjmeffZZFixaxZMkStm/fTpMmTaoiFCFEDfDf7d56BzWQNf9MVJ+w+qw5kMj2M6m8tfYI343ogEreDEJUO5VBVuYsl8zMTBwdHcnIyMDBwUHpcISoVX45fJkxKw+iUal4O7wLof6mu+NHfj5kZUGXLmBrq3Q0ykhIy+WBuX+RX6hnZr9gnmwr+z8LUVF3m5dUSRdwVlYWr732mnHMX0BAAD179mTatGmcPHmyKi4phKihMnILeefn4wDcVy+A4Iamm/yJYg1cbYi6vykA0zacIDkrX+GIhDA9VZIADh48mB9++IGnn36aadOm8dJLL/HHH3+wfPlyWrZsSe/evbl8+XJVXFoIUcN8sPEEqdkFeFrb8kRQY1nzTwDwbGd/guo7kJlfxLv/O650OEKYnCr5KN68eTPr169n0qRJPP/887zyyiuYm5uzceNGzp8/j6enJ23btiUuLq4qLi+EqCF2nU1l9f6LADzVOARnR9nuTRQz06iZ8VgIGrWKDYevEH3iqtIhCWFSqiQB9PT0JDc3t9TXGjZsyBdffMGoUaN45ZVXquLyQogaIL9Qx4S1xdu9dfZoQCfZ7k38R1B9R4bf4w/AxHVHyS4oUjgiIUxHlSSAY8aM4dlnn+XQoUNlHvPMM8/wxx9/VMXlhRA1wNzfz3AhLRcnSyueCmwua/6JUr0a0ZQGLjZcychn9m+yTZwQ1aVKEsCoqCh69epF69atefDBB1m0aBF6vb7EVP9Vq1bh5uZWFZcXQijsaGIGi7cXb/fWv2EQ3u6y3ZsonbWFhml9gwBYujueAwmyTZwQ1aHKhmPPnj2bXbt2YW9vz2uvvUZeXh6hoaE0atQIV1dXpk6dyocffnjH9S9YsAA/Pz+srKxo3749e/fuLfPYxYsX06VLF5ydnXF2diYiIuKWxwsh7lyRTs+baw6j0xsIc63HfS08lQ5J1HBdmrjzWOv6GAww4SfZJk6I6lClnTLt27fnhx9+QKvVcuDAAU6fPk1mZiZubm706NEDDw+PO6p39erVREVFsWjRItq3b8/cuXOJjIzk1KlTpda5detWBgwYQKdOnbCysmLmzJk88MADHDt2jPr169/tbQoh/mXJzniOJmZiY2bG080CsbBQOiJRG0x8OJCtp1KM28SN6SGbAwhRlWrlQtDt27enbdu2fPrppwDo9Xp8fX156aWXePPNN297vk6nw9nZmU8//ZTBgweX65qyELQQt/fvBX6f8g+hf7iv7PjxH7IQdNnWHUxk7OpYLDRqfh3bhcbudkqHJESNVSMXgk5ISKjQ8YmJieU+VqvVEhMTQ0REhLFMrVYTERHB7t27y1VHbm4uhYWFuLjIrEQhKovBYOCttUfIL9TT1NGVR1r6SPInKqR3mDddm7qj1emZsOYIen2ta58QotaokgSwbdu2jBw5kn379pV5TEZGBosXLyYoKIiffvqp3HWnpqai0+nw9Cw5rsjT05OkpKRy1TF+/Hi8vb1LJJH/VVBQQGZmZomHEKJsPx1IZMfZVMzVap5uGoy1tWR/omJUKhXT+gRhba5hb9w14xqSQojKVyVjAI8fP860adO4//77sbKyIjw8HG9vb6ysrLh+/TrHjx/n2LFjtG7dmlmzZtGzZ8+qCKNUM2bMYNWqVWzduhUrK6syj5s+fTrvvvtutcUlRG2WklXA1F+Kd3OIrN+Ulg2kb1PcGV8XG157oCnvbzjBBxtPcF9zDzwcyv6sFkLcmSppAXR1dWXOnDlcuXKFTz/9lCZNmpCamsqZM2cAGDhwIDExMezevbvCyZ+bmxsajYarV0uuGn/16lW8vLxuee7s2bOZMWMGmzdvJiQk5JbHTpgwgYyMDOPj4kX5S1SIsrz7v2Nk5BXiY+tAv5b+st2buCtDO/kR4uNIVn4Rb687Si0cqi5EjVels4Ctra3p378//fv3N47zu9tZtxYWFoSHhxMdHU2fPn2A4kkg0dHRjBkzpszzZs2axbRp0/jtt99o06bNba9jaWmJpaXlXcUqhCnYdPQKvxy+glql4ukmITjYS/Yn7o6ZRs3MfiE8+ukOthy/ys+HLtM7TFZsEKIyVfkn9c6dO/H396dBgwY0aNAAT09Pxo8ff1dj6qKioli8eDFLly7lxIkTjBo1ipycHIYNGwbA4MGDmTBhgvH4mTNnMmnSJL7++mv8/PxISkoiKSmJ7Ozsu74/IUzZ9RwtE9cdBeBer0a0aeSocESirmhRz4GX/lkKZsrPx0jJKlA4IiHqlipPAEeOHEmLFi3Yt28fp06d4sMPP+T333+ndevWFZr9+29PPvkks2fPZvLkyYSFhREbG8umTZuME0MSEhK4cuWK8fiFCxei1Wrp378/9erVMz5mz55dKfcohKl653/HSM3WUs/GjgHBTdBolI5I1CWjujcmsJ4D6bmFTFx3RLqChahEVb4OoLW1NYcOHaJp06bGMoPBwBNPPAHADz/8UJWXrzSyDqAQJf12LImRy2NQAVHBnenU3EnpkGoFWQewYo5fzuTRT3dQpDcwf0AreoV6Kx2SEDVCjVwH8N9atGhBcnJyiTKVSsV7773Hpk2bqvryQogqcD1Hy9tr/+n6rdeYdgFOygYk6qxAbwfG9AgAYPL6o6RmS1ewEJWhyhPAoUOH8tJLL900i1Za0oSovd793zFSswvwsrHjqaAmmFXpdDJh6kZ3D6BFPQeu5xYyef1RpcMRok6o8o/tsWPHAtCkSRMee+wxwsLC0Ol0fPvtt8yaNauqLy+EqGSbjyWxLvYyKuDpxiG4OsnAP1G1LMzUzH48hN6f7mTjkSQ2HL7CwyH1lA5LiFqtyhPAK1euEBsby6FDh4iNjeWbb77hzJkzqFQqZs2axa+//kpISAghISE8+OCDVR2OEOIupOdqefufWb/d6zWifYCzwhEJU9HS25HR9wbwSfQZJq0/SodGLrjayVJdQtypKp8EUpr8/HyOHDlSIjE8evQo6enp1R1KuckkECFg7KqDrIu9jJeNLe917iKtf3dAJoHcOW2Rnkc/3cHJpCx6Bnux4OnWqGTDaWGi7jYvUWTkjpWVFW3btqVt27ZKXF4IcQd+OXz5/7t+A0Il+RPVrrgrOJQ+C4q7gtfHXqZPK1kgWog7IUv2CyFuKykj3zjr9756AXSQrl+hkKD6jrx8X/EC0ZPWHyUxPU/hiISonSQBFELcksFg4PUfD5GRV0gDO0cGhMiCz0JZo7s3plUDJ7Lyixj3/SH0elkgWoiKkgRQCHFLy3ZfYPuZVMzVaoY0C8XJQT42hLLMNGo+fiIMa3MNu8+n8fXOOKVDEqLWkU9yIUSZziZn88HGEwA84tuCED97hSMSopifmy2THgkEYNamU5xKylI4IiFqF0kAhRClKtTpeXV1LAVFepo7udE3qCFq+cQQNciAdr7c19wDrU7P2NWxFBTplA5JiFpDPs6FEKX6JPoMRxIzsDEzZ2hgKLY2styGqFlUKhUz+oXgYmvBiSuZzNlyWumQhKg1JAEUQtxkf/w1Fvx5FoAn/IMJ8LZSOCIhSudub8n0x4IB+GLbeXafS1M4IiFqB0kAhRAlpOdqeWVVLHoDtHGrT2RgPWStXVGTRbb04sk2vhgMMHb1Qa7laJUOSYgaTxJAIYSRwWDgjR8Pk5ieh7u1DcNCgrCwUDoqIW5vyqOBNHa35WpmAeN+OIQCm1wJUatIAiiEMFr+9wU2H7+KRqViWNPWeLkqslmQEBVmY2HGp0+3xsJMzR8nk/lqhywNI8StSAIohADg2OUM3v+leMmXXr4taNPYUeGIhKiYFvUcjEvDzNx0ksOX0pUNSIgaTBJAIQQ5BUW8tPIgWp2eIGcP+of4yW4folZ6pn0DHgryolBnYMzKg2TlFyodkhA1kiSAQggmrT/K+dQcnCytGB4cirW1zPoQtZNKpWLGYyHUd7Im4Voub609KuMBhSiFJIBCmLjv919kzYFEVMDggFY08JRZH6J2c7Qx55MBrdCoVfzv0GVW7k1QOiQhahxJAIUwYUcTM5i47igAkfWbck8zF4UjEqJyhDd05vXIZgC8+/NxYi+mKxuQEDWMJIBCmKj0XC0vfBuDtkhPoJMHA1sFyLg/UaeM7NqIyJaeaHV6Rn8bQ1p2gdIhCVFjSAIohAnS6w28siqWS9fzcLOyYWRoGDYy7k/UMSqVig8fD6WRmy2XM/J5edVBdHoZDygESAIohEmaF32Gv06nYK5WM7x5OD4e5kqHJESVcLAyZ9GgcKzNNew8m8ZHm08pHZIQNYIkgEKYmD9PJjMv+gwAj/sF06axg8IRCVG1mnraM7N/CACfbT3Hb8eSFI5ICOVJAiiECTmbnM3Lqw4C0NmjIb2CfVDLp4AwAY+GevNsZ38AXvv+EKevZikckRDKko9+IUxEeq6W55buIyu/iEb2zjzbKlD2+RUmZULP5nRo5EJ2QRHDl+7jWo5W6ZCEUIwkgEKYgEKdntErDhCflourlTWjQ8NxcpD//sK0mGvULBwYTkNXGy5eyzPOghfCFMlvACFMwLv/O8auc2lYaTQ837wN/vUslQ5JCEU421rw1ZA22FuasTfuGpPWyU4hwjRJAihEHbdsdzzf/p2ACngmoBXhATLpQ5i2AA975j/dCrUKVu+/yFc74pQOSYhqJwmgEHXY1lPJvPu/4wD09GnO/YGeqGS5PyHo3syDtx8OBOCDjSf4/fhVhSMSonrV2gRwwYIF+Pn5YWVlRfv27dm7d2+Zxx47dox+/frh5+eHSqVi7ty51ReoEAo5cimD0SsOoNMbaOvmw4CwRpiZKR2VEDXHs539GNDOF70Bxnx3gAMJ15UOSYhqUysTwNWrVxMVFcWUKVM4cOAAoaGhREZGkpycXOrxubm5NGrUiBkzZuDl5VXN0QpR/RLSchn2zV5ytTqaObrxQngw1rLThxAlqFQq3usdRPdm7uQX6hn+zT7OpWQrHZYQ1aJWJoBz5sxhxIgRDBs2jMDAQBYtWoSNjQ1ff/11qce3bduWDz/8kKeeegpLSxn8Luq2tOwChizZS2q2Fh9bB15q3Vpm/ApRBnONms8GtibUx5HruYUM+XovyZn5SoclRJWrdb8VtFotMTExREREGMvUajURERHs3r270q5TUFBAZmZmiYcQNV2utojhS/cTl5qDq5U1L4e2pZ6bbPMmxK3YWJjx1dC2+LnacOl6HkOX7CMrv1DpsISoUrUuAUxNTUWn0+Hp6Vmi3NPTk6SkytveZ/r06Tg6Ohofvr6+lVa3EFWhoEjHyOUxxF5Mx9bcnNEt29G4vpXSYQlRK7jZWbL02Xa42Vlw/EomI5fHkF+oUzosIapMrUsAq8uECRPIyMgwPi5evKh0SEKUqVCnZ8zKg2w/k4qlRsPzzdsS6m+ndFhC1CoNXW35emhbbCw07DqXxugVB2ShaFFn1boE0M3NDY1Gw9WrJafsX716tVIneFhaWuLg4FDiIURNpNMbGPfDIbYcv4q5Ws1zzdrQuZmzLPcixB0I8XHiqyFtsTRT88fJZF5ZdZAinSSBou6pdQmghYUF4eHhREdHG8v0ej3R0dF07NhRwciEqH4Gg4G31x5hfexlNCoVQwJa072FmyR/QtyFjo1d+WJwGyw0an49msS4Hw6h08tuIaJuqXUJIEBUVBSLFy9m6dKlnDhxglGjRpGTk8OwYcMAGDx4MBMmTDAer9VqiY2NJTY2Fq1WS2JiIrGxsZw9e1apWxDirun1BiauO8qqfRdRAQMbt+KBIE/UtfJ/tRA1S7em7nz6dCs0ahXrYi/z1poj6CUJFHVIrVwW9sknnyQlJYXJkyeTlJREWFgYmzZtMk4MSUhIQP2v34KXL1+mVatWxuezZ89m9uzZdOvWja1bt1Z3+ELcNZ3ewIQ1h/l+/yVUwFONQnkkpB4ajdKRCVF3PNDSi7lPhvHKqoOs3n+RQr2eWf1CMNPIX1mi9lMZZBfscsnMzMTR0ZGMjAwZDygUVaTT8/qPh1l7MPGflr8weoXUl10+aon8fMjKgi5dwNZW6WhEefx86DKvro5FpzfwcEg95j4ZhrkkgUJhd5uXyK8MIWqRQp2eV1fH8svhK6hVKp5p3IqHg+tJ8idEFXo01BsLjZqXvjvAhsNX0Bbp+fTpVliaSZO7qL3kTxghaolcbREjlu3nl8NX0KhUDA1ozSMhkvwJUR0eDPLii0FtsDBTs+X4VZ5bup+cgiKlwxLijkkCKEQtkJZdwIAv/mbrqRQs1GqGN23DQyFeMuZPiGp0b3MPlgxti7W5hu1nUnnqi79JySpQOiwh7ogkgELUcBev5dJ/0W4OXcrA1tycFwM7cH+Qh8z2FUIBnQPc+O75DrjYWnAkMYPHFu4kLjVH6bCEqDD5FSJEDXb4UjqPLdxFXGoOLpbWvBrcic7NnSX5E0JBYb5OrBnViQYuNly8lke/hbs4mHBd6bCEqBD5NSJEDfW/Q5d5fNFuUrIK8LaxZ3x4J1o1tpNFnoWoAfzcbPlpVCdCfBy5lqPlqS/+Zn1sotJhCVFukgAKUcPo9QbmbD7FS98dpKBIT6CTO2+170hAfSulQxNC/Iu7vSXfjejAfc09KCjS88qqWGb8elJ2DRG1giSAQtQgOQVFvLjyAJ/8UbxLTXevRozv1JZ6buYKRyaEKI2tpRlfDG7DqO6NAVj01zlGLNtPVn6hwpEJcWuSAApRQ5xMyqTXpzv49WgSZmo1A/xDGdm+BXa20ucrRE2mUasY/2Bz5j0VhqWZmj9OJtP7052cuJKpdGhClEkSQCFqgO/3X6TPgp2cT8nBydKKMYHteay1DxYWSkcmhCiv3mH1+eGFjtRztOJ8ag69F+xk5Z4EZMMtURNJAiiEgnIKinjt+0O88eNh8gv1NHdyZ1Lbe+jSwkVm+gpRC4X4OLHh5S7c28wdbZGet9Ye4eVVsdIlLGoc+RUjhEL2x1/joXnb+enAJVRAT59mvNW5LX71LJUOTQhxF1xsLfhqSFsmPNQcjVrF/w5d5uFPdrA37prSoQlhJAmgENWsoEjHzE0neeLz3SRcy8XZ0ooxgR0Y0jYAWxsZ7ydEXaBWqxjZrTHfj+xIfSdrEq7l8uQXu5m24Tj5hTqlwxNCEkAhqtPhS+n0/nQnC7eeQ2+Atm4+vNexK91busqevkLUQeENnfl1bBeeaOODwQCLt8fx8Cfbib2YrnRowsSpDDI6tVwyMzNxdHQkIyMDBwcHpcMRtUxWfiEfbT7Nst3x6A1gb27B437B3B/oJRM9TEx+PmRlQZcuYGurdDSiOkWfuMqba46QklWASgXPtG/IuMhmOFrLMk+i4u42L5EEsJwkARR3wmAw8OvRJN793zGuZhZvGh/u6s3TgYH4eclYP1MkCaBpu56j5b1fjrP2YPGuIW52lkx8uAW9w7xRyTY/ogIkAawmkgCKijqamMH7G47z9/nigd/u1jY87hdM12ZumMsf/CZLEkABsOtsKhPXH+V8Sg4A7f1dmPhwIME+jgpHJmoLSQCriSSAorySMvL58LdTrDl4CYMBzNVquns14vGWAbg6aZQOTyhMEkBxQ0GRjsXbzjP/j7MUFOkB6BPmzbjIZvg42ygcnajpJAGsJpIAittJzS5g8bbzLN0dT35h8Yd5uKs3/Zo2o4m3jazrJwBJAMXNLl3P5aPNp43dwhZmagZ3aMjIbo1xt5ehIqJ0kgBWE0kARVluJH7Ldl8g75/lHRrZO/NYo0DaNnKS2b2iBEkARVmOXMrgg40n2H0+DQArczVPt2vIC90a4eFgpXB0oqaRBLCaSAIo/isuNYclO+P4Yf8lY+LX0M6Rh3yb0iXAHSsrGdAtbiYJoLgVg8HAX6dTmPv7GeNSMRZmap5s48uwzn40crdTNkBRY0gCWE0kARRQ/OH89/lrfLUjjuiTV7nxv6eBnSMP+jSla4A71taS+ImySQIoysNgMLD9TCrzos8Qc+G6sbxHcw+e7exP5wBXmTVs4u42L5HOKSHK4XqOlrUHE/l+/0VOJmUZy1s6e3Cvtz8dG7lKi58QotKoVCq6NnWnSxM3dp9L4+udcUSfTOaPfx5NPe14qm0D+rSqj4utLCYqKk5aAMtJWgBNT5FOz65zaazef5Etx66i1RVP7LBQq2nj5sP9Df0J9LGTMX6iQqQFUNypuNQclu6K54f9F8nRFg87MdeoiGjhyeNtfOjaxB0zjcw2MxXSBVxNJAE0DTq9gT1xaWw4fIVNR5NIy9EaX/O1daCduy/d/OpTz9VcZvWKOyIJoLhbmfmFrD+YyA8xlzh8KcNY7mprQWSQFz2D6tGhkYskg3WcJIDVRBLAuiunoIidZ1P581QKW45fJTW7wPiarZk5YS7edPPxJcjHEUtZkUHcJUkARWU6cSWTH/ZfYl1sItf+9Qers405kS29uLe5B50D3LCzlK6KukYSwGoiCWDdodcbOJ2cxfbTqWw9nczeuGsU6v7/v4GNmTnBzl608ahHuK8rDnZqZKy1qCySAIqqUKjTs/tcGr8eLe69uJ5baHzNTK2ijZ8z3Zp60K2pO8297FGr5UOttpMEsJpIAlh7Fer0HLucyd64NPbGXWNf/HUy8gpLHONmZUNzR3dC3Dxo4+smSZ+oMpIAiqpWpNPz9/lrbD6exLbTKcSn5ZZ43dHanDYNnWnr70JbPxeC6ztiYSbdxbWNzAIW4l8KdXrOXM3mSGI6RxIzOJKYyYkrmWj/2WbpBkuNBj87ZwKdPGhdzx1/N1tZvkUIUSeYadTc08SNe5q4ARCfmsO2MylsPZXC3+fTyMgrJPpkMtEnkwGwNFPTop4DQfUdCPJ2JKi+I0097SUprOMkARS1UkGRjvjUXM4kZ3E2Odv4OJ+ac1OyB8Xduo3sXWhk70wLN1eaeThgb6uWiRxCiDrPz80WPzdbBnf0o1Cn5/jlTPbFX2Nv3DX2X7jOtRwtsRfTjQtPQ/Hs4sbudgR42JX42sjdFitz2dO8Lqi1XcALFizgww8/JCkpidDQUObPn0+7du3KPP6HH35g0qRJxMfH06RJE2bOnEnPnj3LfT3pAq5eudoiUrIKSLyex6XreVy6nvvP1+LvkzLz0ZfxzrXWmOFj64iPrSMN7R1p4uZIAxcbrK1U0q0rFCddwKImMRgMxKflcjQxg6OXM4q/JmbeNEzmBpUKvB2tqe9kjY9z8aO+szU+zjbUd7LGw8ESGwtpW6oOJtkFvHr1aqKioli0aBHt27dn7ty5REZGcurUKTw8PG46fteuXQwYMIDp06fzyCOPsHLlSvr06cOBAwcICgpS4A5MS6FOT2ZeIRl5hWTmF5Fx4/u8Qq7naEnJLiAlq/iR+s/3N9a4uhUrjRle1nZ4WNvhZW1HfXs7GjrZU9/ZGitLlbTuCSHEbahUKvzdbPF3s6VXqDdQnBReup7H6atZnEsp7l05l5LD2eRsMvIKSUzPIzE9j73xpddpa6HBzd4SdztL3O0tcfvnq7OtBY7W5iUeDlZmOFibYy5L1lS7WtkC2L59e9q2bcunn34KgF6vx9fXl5deeok333zzpuOffPJJcnJy+OWXX4xlHTp0ICwsjEWLFpXrmrWlBdBgMKA3gN5gQG8wYDB+X/zVoIcivZ5CnYFCnf6fR8nvi3R6tDo9Rf+U//v7giI9eYU6crU68gt15GqL/vW9jjytjrzC4q/ZBcXJXm45krnSmKvVuFha42xhg4ulNa5WNrjbWONha42XgzVutpZYSqInahlpARS1lcFgIDVby0Vjj0xuiV6ay+n5xn3RK8rWQoODtTk2FhpsLMywttD8870Ga3Oz///+n68WGjXmZmrMNWos//lqrlFjYabGXKMqft34XI2FpnjIj0atQqNSob7xVaUylqtVqhKv13Qm1wKo1WqJiYlhwoQJxjK1Wk1ERAS7d+8u9Zzdu3cTFRVVoiwyMpJ169aVeZ2CggIKCv5/PbiMjOLFNjMzM+8i+rIV6vT0nLfdmMCV/PpPAkfJhM5gfPz/MTWZpUaDtZk51mpzrDRmWJuZYWtmgZ25JQ4WFjhaWeJoaYGztSXONpbYWWowN1ehKXW4iRadTktubmmvCVFzFRRATg5kZoLuzn5XCqEYSyDASUOAkx3425V4zWAwkKPVkZZd3JuTlq01fk3JKiA9T0tWfhGZ+UVk5hWSmV9ITkHxf4KsguI/jGqSG0mhSgUaNcaEUaWi+Cv8M6zoRtn/f1/8mor1YzpXWZf4jXzkTtvxal0CmJqaik6nw9PTs0S5p6cnJ0+eLPWcpKSkUo9PSkoq8zrTp0/n3Xffvanc19f3DqIWQgghhKmpd3MaUemysrJwdHSs8Hm1LgGsLhMmTCjRaqjX67l27Rqurq6oqmgmQWZmJr6+vly8eLFGdzNXBVO9d1O9bzDdezfV+wbTvXdTvW8w3Xuvjvs2GAxkZWXh7e19R+fXugTQzc0NjUbD1atXS5RfvXoVLy+vUs/x8vKq0PEAlpaWWP5n3y8nJ6c7C7qCHBwcTOo/yr+Z6r2b6n2D6d67qd43mO69m+p9g+nee1Xf9520/N1Q64bPW1hYEB4eTnR0tLFMr9cTHR1Nx44dSz2nY8eOJY4H2LJlS5nHCyGEEELUZbWuBRAgKiqKIUOG0KZNG9q1a8fcuXPJyclh2LBhAAwePJj69eszffp0AF555RW6devGRx99xMMPP8yqVavYv38/X3zxhZK3IYQQQgihiFqZAD755JOkpKQwefJkkpKSCAsLY9OmTcaJHgkJCaj/tTZIp06dWLlyJRMnTuStt96iSZMmrFu3rsatAWhpacmUKVNu6no2BaZ676Z632C6926q9w2me++met9guvdeG+67Vq4DKIQQQggh7lytGwMohBBCCCHujiSAQgghhBAmRhJAIYQQQggTIwmgEEIIIYSJkQSwhisoKCAsLAyVSkVsbKzS4VSLRx99lAYNGmBlZUW9evUYNGgQly9fVjqsKhUfH8/w4cPx9/fH2tqaxo0bM2XKFLRardKhVYtp06bRqVMnbGxsqm3BdaUsWLAAPz8/rKysaN++PXv37lU6pCq3bds2evXqhbe3NyqV6pb7sNcl06dPp23bttjb2+Ph4UGfPn04deqU0mFVuYULFxISEmJcBLljx478+uuvSodV7WbMmIFKpWLs2LFKh1IqSQBruDfeeOOOt3mpre69916+//57Tp06xU8//cS5c+fo37+/0mFVqZMnT6LX6/n88885duwYH3/8MYsWLeKtt95SOrRqodVqefzxxxk1apTSoVSp1atXExUVxZQpUzhw4AChoaFERkaSnJysdGhVKicnh9DQUBYsWKB0KNXqr7/+4sUXX+Tvv/9my5YtFBYW8sADD5CTk6N0aFXKx8eHGTNmEBMTw/79++nRowe9e/fm2LFjSodWbfbt28fnn39OSEiI0qGUzSBqrI0bNxqaN29uOHbsmAEwHDx4UOmQFLF+/XqDSqUyaLVapUOpVrNmzTL4+/srHUa1WrJkicHR0VHpMKpMu3btDC+++KLxuU6nM3h7exumT5+uYFTVCzCsXbtW6TAUkZycbAAMf/31l9KhVDtnZ2fDl19+qXQY1SIrK8vQpEkTw5YtWwzdunUzvPLKK0qHVCppAayhrl69yogRI1i+fDk2NjZKh6OYa9eusWLFCjp16oS5ubnS4VSrjIwMXFxclA5DVBKtVktMTAwRERHGMrVaTUREBLt371YwMlFdMjIyAEzq/7VOp2PVqlXk5OSYzParL774Ig8//HCJ/+s1kSSANZDBYGDo0KG88MILtGnTRulwFDF+/HhsbW1xdXUlISGB9evXKx1StTp79izz589n5MiRSociKklqaio6nc64Y9ENnp6eJCUlKRSVqC56vZ6xY8fSuXPnGrcLVVU4cuQIdnZ2WFpa8sILL7B27VoCAwOVDqvKrVq1igMHDhi3oq3JJAGsRm+++SYqleqWj5MnTzJ//nyysrKYMGGC0iFXmvLe+w2vv/46Bw8eZPPmzWg0GgYPHoyhFm5aU9H7BkhMTOTBBx/k8ccfZ8SIEQpFfvfu5N6FqKtefPFFjh49yqpVq5QOpVo0a9aM2NhY9uzZw6hRoxgyZAjHjx9XOqwqdfHiRV555RVWrFiBlZWV0uHclmwFV41SUlJIS0u75TGNGjXiiSee4H//+x8qlcpYrtPp0Gg0DBw4kKVLl1Z1qJWuvPduYWFxU/mlS5fw9fVl165dta4LoaL3ffnyZbp3706HDh345ptvSuxpXdvcyb/5N998w9ixY0lPT6/i6KqfVqvFxsaGH3/8kT59+hjLhwwZQnp6usm0cqtUKtauXVviZ1DXjRkzhvXr17Nt2zb8/f2VDkcRERERNG7cmM8//1zpUKrMunXr6Nu3LxqNxlim0+lQqVSo1WoKCgpKvKY0M6UDMCXu7u64u7vf9rhPPvmE999/3/j88uXLREZGsnr1atq3b1+VIVaZ8t57afR6PVC8JE5tU5H7TkxM5N577yU8PJwlS5bU6uQP7u7fvC6ysLAgPDyc6OhoY/Kj1+uJjo5mzJgxygYnqoTBYOCll15i7dq1bN261WSTPyh+r9fGz/CKuO+++zhy5EiJsmHDhtG8eXPGjx9fo5I/kASwRmrQoEGJ53Z2dgA0btwYHx8fJUKqNnv27GHfvn3cc889ODs7c+7cOSZNmkTjxo1rXetfRSQmJtK9e3caNmzI7NmzSUlJMb7m5eWlYGTVIyEhgWvXrpGQkIBOpzOueRkQEGB8/9cFUVFRDBkyhDZt2tCuXTvmzp1LTk4Ow4YNUzq0KpWdnc3Zs2eNz+Pi4oiNjcXFxeWmz7u65MUXX2TlypWsX78ee3t741hPR0dHrK2tFY6u6kyYMIGHHnqIBg0akJWVxcqVK9m6dSu//fab0qFVKXt7+5vGd94Yy14jx30qOgdZlEtcXJzJLANz+PBhw7333mtwcXExWFpaGvz8/AwvvPCC4dKlS0qHVqWWLFliAEp9mIIhQ4aUeu9//vmn0qFVuvnz5xsaNGhgsLCwMLRr187w999/Kx1Slfvzzz9L/fcdMmSI0qFVqbL+Ty9ZskTp0KrUs88+a2jYsKHBwsLC4O7ubrjvvvsMmzdvVjosRdTkZWBkDKAQQgghhImp3YOMhBBCCCFEhUkCKIQQQghhYiQBFEIIIYQwMZIACiGEEEKYGEkAhRBCCCFMjCSAQgghhBAmRhJAIYQQQggTIwmgEEIIIYSJkQRQCCGEEMLESAIohFBE9+7dGTt2rNJhCIWkpaXh4eFBfHy8scxgMDBnzhz8/f2xsbGhT58+ZGRkGF9/6qmn+OijjxSIVoi6R7aCE0Io4tq1a5ibm2Nvb690KOI/unfvTlhYGHPnzq2ya0RFRZGVlcXixYuNZePGjWP9+vV89dVX2Nra0qdPH/r378/HH38MwNGjR+natStxcXE4OjpWWWxCmAJpARRCKMLFxUWSv1JotVqlQ6g0Zd1Lbm4uX331FcOHDzeW7dmzhzlz5rB69Wq6du1KeHg4I0aMYOPGjcZjgoKCaNy4Md9++22Vxy5EXScJoBCiwrp3787LL7/MG2+8gYuLC15eXrzzzjvG1/38/G5qPQoLCytxzH+7gH/88UeCg4OxtrbG1dWViIgIcnJyANi0aRP33HMPTk5OuLq68sgjj3Du3LkKxaTX65k1axYBAQFYWlrSoEEDpk2bVuL16dOn4+/vj7W1NaGhofz444+3/TmMGTOGMWPG4OjoiJubG5MmTeLfHSu3i/1GHWPHjsXNzY3IyMgK3fNLL73E2LFjcXZ2xtPTk8WLF5OTk8OwYcOwt7cnICCAX3/9tdz3OXToUP766y/mzZuHSqVCpVIRHx9frp9PWffyXxs3bsTS0pIOHToYy2bPns19991H69atjWWenp6kpqaWOLdXr16sWrXqlv8uQojbkwRQCHFHli5diq2tLXv27GHWrFm89957bNmy5Y7qunLlCgMGDODZZ5/lxIkTbN26lccee8yYSOXk5BAVFcX+/fuJjo5GrVbTt29f9Hp9uWOaMGECM2bMYNKkSRw/fpyVK1fi6elpPHf69OksW7aMRYsWcezYMV599VWeeeYZ/vrrr9v+HMzMzNi7dy/z5s1jzpw5fPnll8bXyxP70qVLsbCwYOfOnSxatKjC9+zm5sbevXt56aWXGDVqFI8//jidOnXiwIEDPPDAAwwaNIjc3Nxy3ee8efPo2LEjI0aM4MqVK1y5cgVfX99y/3xKu5f/2r59O+Hh4cbnBQUFbNiwgb59+5Y4Lj8//6au3nbt2rF3714KCgpu+e8ihLgNgxBCVFC3bt0M99xzT4mytm3bGsaPH28wGAyGhg0bGj7++OMSr4eGhhqmTJlSoo5XXnnFYDAYDDExMQbAEB8fX67rp6SkGADDkSNHyhVTZmamwdLS0rB48eJS68vPzzfY2NgYdu3aVaJ8+PDhhgEDBpQZR7du3QwtWrQw6PV6Y9n48eMNLVq0KHfs3bp1M7Rq1arsmy3jvBvn/vuei4qKDLa2toZBgwYZy65cuWIADLt37y73ff7738ZgKP/Pp7z30rt3b8Ozzz5rfL5r1y4DYLCysjLY2toaHxYWFobIyMgS5x46dKhC7xUhROnMFMw9hRC1WEhISInn9erVIzk5+Y7qCg0N5b777iM4OJjIyEgeeOAB+vfvj7OzMwBnzpxh8uTJ7Nmzh9TUVGMrWEJCAkFBQbeN6cSJExQUFHDfffeVev2zZ8+Sm5vL/fffX6Jcq9XSqlWrW8beoUMHVCqV8XnHjh356KOP0Ol0aDSacsX+79awG+7knjUaDa6urgQHBxvLbrRyJicn3/F9VuS80u7lv/Ly8rCysjI+P336NLa2tsTGxpY47uGHH6Zz584lyqytrQGMLZpCiDsjCaAQ4o6Ym5uXeK5SqYxJilqtLjEODqCwsLDMujQaDVu2bGHXrl1s3ryZ+fPn8/bbb7Nnzx78/f3p1asXDRs2ZPHixXh7e6PX6wkKCrppkkFZMd1IGsqSnZ0NwIYNG6hfv36J1ywtLW957u2UJ3ZbW9s7Og9Kv+d/l91ITvV6/R3fZ0XOK+1e/svNzY3r168bn2dmZuLm5kZAQICx7MKFC5w5c4Z+/fqVOPfatWsAuLu73/Y6QoiySQIohKh07u7uXLlyxfg8MzOTuLi4W56jUqno3LkznTt3ZvLkyTRs2JC1a9cyZMgQTp06xeLFi+nSpQsAO3bsqFA8TZo0wdramujoaJ577rmbXg8MDMTS0pKEhAS6detWobr37NlT4vnff/9NkyZN0Gg0pKWl3VHsd3re7ZT3Pi0sLNDpdBU+r7xatWpVYiavm5sbGRkZGAwGY8I6bdo0evbsSWBgYIlzjx49io+PD25ubncdhxCmTBJAIUSl69GjB9988w29evXCycmJyZMno9Foyjx+z549REdH88ADD+Dh4cGePXtISUmhRYsWODs74+rqyhdffEG9evVISEjgzTffrFA8VlZWjB8/njfeeAMLCws6d+5MSkoKx44dY/jw4djb2zNu3DheffVV9Ho999xzDxkZGezcuRMHBweGDBlSZt0JCQlERUUxcuRIDhw4wPz5842LFd9p7JVxz6Up7336+fmxZ88e4uPjsbOzw8XF5Y5/PqWJjIxkwoQJXL9+HWdnZ3r06EF+fj4zZszgqaeeYsWKFfzvf/9j7969N527fft2Hnjggbv+WQhh6iQBFEJUugkTJhAXF8cjjzyCo6MjU6dOvWULoIODA9u2bWPu3LlkZmbSsGFDPvroIx566CEAVq1axcsvv0xQUBDNmjXjk08+oXv37hWKadKkSZiZmTF58mQuX75MvXr1eOGFF4yvT506FXd3d6ZPn8758+dxcnKidevWvPXWW7esd/DgweTl5dGuXTs0Gg2vvPIKzz//PFDcFX4nsd/peeVRnvscN24cQ4YMITAwkLy8POLi4u7451Oa4OBgWrduzffff8/IkSPx9PTkm2++4fXXX2fq1Kn06NGDHTt24OvrW+K8/Px81q1bx6ZNm+765yCEqZOdQIQQ4g5Vx44ZddWGDRt4/fXXOXr0KGp1+VYkW7hwIWvXrmXz5s1VHJ0QdZ+0AAohhKh2Dz/8MGfOnCExMfGmlr6ymJubM3/+/CqOTAjTIAmgEEIIRfx7J5jyKG0CjxDizkgXsBBCCCGEiZGt4IQQQgghTIwkgEIIIYQQJkYSQCGEEEIIEyMJoBBCCCGEiZEEUAghhBDCxEgCKIQQQghhYiQBFEIIIYQwMZIACiGEEEKYGEkAhRBCCCFMjCSAQgghhBAmRhJAIYQQQggTIwmgEEIIIYSJkQRQCCGEEMLESAIohBBCCGFiJAEUQgghhDAxkgAKIYQQQpgYM6UDqC30ej2XL1/G3t4elUqldDhCCCGEMGEGg4GsrCy8vb1RqyvenicJYDldvnwZX19fpcMQQgghhDC6ePEiPj4+FT5PEsBysre3B4p/0A4ODgpHI4QQQghTlpmZia+vrzE/qShJAMvpRrevg4ODJIBCCCGEqBHudFiaTAIRQgghhDAxkgAKIYSosP37oUeP4q9CiNpHuoCFEEJU2LJl8OefsHw5tGmjdDR1j06no7CwUOkwhILMzc3RaDRVVr8kgEIIIW4rT6tj75E8Ei7ryNXqWLbCEdCw7Fs9re/LxNnWkmYNLGjWpOp+YZkCg8FAUlIS6enpSociagAnJye8vLyqZPk5SQCFEEKUkJ6rZX/8dQ5evE7sxXTOJeeQlJnPhZkP/+soQ/Gx11QM7e1kLA2f+juB3g60qGdP6wbOdGjkiqO1efXeQC12I/nz8PDAxsZG1p01UQaDgdzcXJKTkwGoV69epV9DEkAhhBAkpOWy8egVok9cJebCdfSGm4+p9+hhrvwSBHo1cCMx+eerWo9rz0OkZhew7XQK206nFBerINTXifsDPekV4o2vi0213E9tpNPpjMmfq6ur0uEIhVlbWwOQnJyMh4dHpXcHSwIohBAmKr9Qxy+Hr/D9/ovsjbtW4jUPa1v87Jxp5OCMv7M9Po62OD9qzrn+KgYPvrmu/XtVNGkZxLlkP45fyeTY5Uz+PpfG+dQcDiakczAhnVmbTtG6gRNPt2/IIyH1sDKX7uJ/uzHmz8ZGkmRR7MZ7obCwUBJAIYQQd+d6jpblf1/gm13xXMvRAsXteE0cXQl29qKdjwd+7jaYlfIb4saOUyoVGAzFz/X64rXIHKzMadXAmVYNnI3HJ6bnsfVUMr8cusLfcWkcSEjnQEI6H2w8wYB2vjzb2R9XO8tquOvaQ7p9xQ1V+V6QBFAIIUxERm4hn209y7LdF8gr1AHgYmlNe3dfevj70NDdmts1Mjg7g6sruLvD/ffDzp1w+TJ4eJR+fH0nawa2b8jA9g1JzsznxwOX+Hb3BS5n5LPgz3Ms2RnPkE5+PN+lEc62FpV8x0KIsqgMBkMpIz3Ef2VmZuLo6EhGRobsBCKEqFXyC3Us332BT/88S0ZecTejj60D99VvTI8AL+xsK7YkrFYLOh1kZ8M994C5OVhWoBGvSKfn9xPJLPjzLEcSMwCwtdAw+t4Anuvij6WZaXYN5+fnExcXh7+/P1ZWVkqHU6mSkpIYNGgQu3btwtzcnPT09FLLREm3ek/cbV5SYxeCXrBgAX5+flhZWdG+fXv27t1b5rGLFy+mS5cuODs74+zsTERExE3HDx06FJVKVeLx4IMPVvVtCCGEov46ncIDH29j2sYTZOQV4m1jz8gWbZnZ4x4eDfWucPIHYGFR3AUMxV8rkvwBmGnUPBjkxc9jOvPl4Da09HYgR6vjw99OEfnxNv44ebXCMQlllfY79t+/Zz/++GOuXLlCbGwsp0+fLrPsbvn5+TF37txKqauuq5FdwKtXryYqKopFixbRvn175s6dS2RkJKdOncKjlH6GrVu3MmDAADp16oSVlRUzZ87kgQce4NixY9SvX9943IMPPsiSJUuMzy0r+qklhBC1REpWAVN/Oc7Phy4D4GRhyUM+zegZ6IONdc0YY6ZSqYgI9OS+Fh6si01k+saTxKfl8uw3+4ls6cn7fYJxt5fP6driv79j4f9/z547d47w8HCaNGlifK20MlF9amQL4Jw5cxgxYgTDhg0jMDCQRYsWYWNjw9dff13q8StWrGD06NGEhYXRvHlzvvzyS/R6PdHR0SWOs7S0xMvLy/hwdnYutT4hhKjNNh65QsScv/j50GVUQFdPf6Z36U7/cN8ak/z9m0qlom8rH/4Y152RXRthplbx27GrPPDxX/xy+LLS4Yly+u/v2Bu/Z/38/Pjpp59YtmwZKpWKoUOHlloGkJ6eznPPPYe7uzsODg706NGDQ4cOlbjO//73P9q2bYuVlRVubm707dsXgO7du3PhwgVeffVVYwukKFuNawHUarXExMQwYcIEY5larSYiIoLdu3eXq47c3FwKCwtxcXEpUb5161Y8PDxwdnamR48evP/++7LWkhCizsgpKOLd/x3j+/2XgOJxfk83CaZNI6fbTu6oCewszZjQswWPhnkz7ofDnLiSyZiVB9l0NIkPHgvGwcr0FpQ2GAzGCTvVzdpcUylJ1L59+xg8eDAODg7MmzcPa2trtFrtTWUAjz/+ONbW1vz66684Ojry+eefc99993H69GlcXFzYsGEDffv25e2332bZsmVotVo2btwIwJo1awgNDeX5559nxIgRdx13XVfjEsDU1FR0Oh2enp4lyj09PTl58mS56hg/fjze3t5EREQYyx588EEee+wx/P39OXfuHG+99RYPPfQQu3fvLnVtnYKCAgoKCozPMzMz7/COhBCi6h25lMHLqw4Sl5qDCuhRrzHPhDXFwa5GdvTcUktvR9a/2JkFf55lwZ9n+eXwFY4kZrDg6dYE1XdUOrxqlVeoI3Dyb4pc+/h7kdhYlD9N+OWXX7CzsytR9tZbb/HWW29haWmJtbU1Xl5extf+W7Zjxw727t1LcnKyset49uzZrFu3jh9//JHnn3+eadOm8dRTT/Huu+8a6wkNDQXAxcUFjUaDvb19ieuI0tW4BPBuzZgxg1WrVrF169YSM2aeeuop4/fBwcGEhITQuHFjtm7dyn333XdTPdOnTy/xBhNCiJrq+30XmbjuKFqdHmdLKwY1CeOepq61otWvLBZmal69vyndm7kzZuVBLqTl8thnu5j4SAsGdWgo3Xs10L333svChQtLlP23J+5WDh06RHZ29k09c3l5eZw7dw6A2NhYad2rJDUuAXRzc0Oj0XD1aslZYFevXr1tRj979mxmzJjB77//TkhIyC2PbdSoEW5ubpw9e7bUBHDChAlERUUZn2dmZuLr61uBOxFCiKqlLdIz9ZfjLP/7AgBBzp48HxpKffe601XaqoEzG1/uwrgfD7Hl+FUmrz/GscRMpvYJwsKs9rVuVpS1uYbj70Uqdu2KsLW1JSAg4I6vl52dTb169di6detNrzk5ORXH9E9Xsbh7NS4BtLCwIDw8nOjoaPr06QNgnNAxZsyYMs+bNWsW06ZN47fffqNNmza3vc6lS5dIS0src4NlS0tLmSUshKixUrMLGPVtDPvir6MCIus3ZWCrgBo5yeNuOdqY88WgcL7cHsf0X0+wev9F4lJzWPhM6zq/i4hKpapQN2xt1rp1a5KSkjAzM8PPz6/UY0JCQoiOjmbYsGGlvm5hYYFOp8yYydqmRv75FBUVxeLFi1m6dCknTpxg1KhR5OTkGP/BBw8eXGKSyMyZM5k0aRJff/01fn5+JCUlkZSURHZ2NlD8V8Xrr7/O33//TXx8PNHR0fTu3ZuAgAAiI5X5y0oIIe7U+ZRs+n62k33x17E2M+O55m14tn2TOpn83aBSqRjRtRFfDW2LvaUZe+Ov8einOzmZJOOza4qCggLj798bj9TU1HKfHxERQceOHenTpw+bN28mPj6eXbt28fbbb7N//34ApkyZwnfffceUKVM4ceIER44cYebMmcY6/Pz82LZtG4mJiRW6timqkQngk08+yezZs5k8eTJhYWHExsayadMm48SQhIQErly5Yjx+4cKFaLVa+vfvT7169YyP2bNnA6DRaDh8+DCPPvooTZs2Zfjw4YSHh7N9+3Zp5RNC1Cr746/x2Ge7uHgtDzcrG94I60xkkGetHu9XEfc282Dti53wc7UhMT2Pxxft5u/zaUqHJYBNmzaV+B1cr1497rnnnnKfr1Kp2LhxI127dmXYsGE0bdqUp556igsXLhh//3fv3p0ffviBn3/+mbCwMHr06FFi44f33nuP+Ph4GjdujLu7e6XfY10iW8GVk2wFJ4RQ2q9HrvDK6li0RXoa2jkxNrwNDTyU+SM2Px+ysqBLF7C1rf7rp+dqGbFsP/vir2NhpuaTp8J4MKj0IT21RV3eCk7cGZPcCk4IIcT/+37fRUavPIC2SE+Qsydvd2yvWPJXEzjZWLB8eHseCPREW6Rn1IoDfPvPZBghxO1JAiiEEDXcNzvjeOOnwxgM0NGjAW90CsfVyTQmBtyKlbmGzwa2ZkC7BhgMMHHdURZuPad0WELUCpIACiFEDfbZ1rO887/jQPGWbmPaB2FrU3cne1SUmUbNB32DeLlH8fIjMzedZH70GYWjEqLmkz8hhRCiBjIYDMzZcpr5f5wF4H7vAIa1aYqlpSR//6VSqYh6oBkWZmpmbz7NR1tOU6g38GpEE1kwWogySAugEELUQPOizxiTv4d9mjO8XTNJ/m5jTI8mvPlQcwA+iT7Dh7+dQuY5ClE6SQCFEKKG+WzrWeb+XtyN2btBCwa3bYx53dnco0q90K0xEx9uAcBnW8/x8ZbTCkckRM0kCaAQQtQgX+2IY9amUwA87NOMp8MbYSaDdSrkuS6NmNIrEIBP/jjL53/JxBAh/ksSQCGEqCG+/fsCU38pnvDxgHcTBrcNkOTvDg3r7M/rkc0AmP7rSVkiRoj/kARQCCFqgLUHLzFx3VEA7q3XiGFtm0jyd5devDeA0d0bAzBp/VHWHrykcERC1BySAAohhMK2nkrm9R8OA3CPpx/PtWmOhYVM+KgMr0c2Y2gnPwwGGPfDYX47lqR0SOIOde/enbFjx95VHUlJSdx///3Y2tri5ORUKXHVVpIACiGEgg5dTGf0igMU6Q20dvVmZNtArKwk+assKpWKyY8E0j/cB53ewMvfHWR//DWlw6pzUlJSGDVqFA0aNMDS0hIvLy8iIyPZuXOn0qGV8PHHH3PlyhViY2M5fbpyJghVRmKqBEkAhRBCIedTshn2zT5ytTqaObrxUrtQbKwl+atsarWKmf1CiGjhSUGRnuFL93M2OVvpsKrc/v3Qo0fx16rWr18/Dh48yNKlSzl9+jQ///wz3bt3Jy0treovXgHnzp0jPDycJk2a4OHhoXQ4JWi12mq9niSAQgihgOTMfAZ/vZdrOVp8bR0Z2zYcBzv5SK4qGrWK+QNa0aqBExl5hQz5ei/JmflKh1Wlli2DP/+E5cur9jrp6els376dmTNncu+999KwYUPatWvHhAkTePTRR0scN3LkSDw9PbGysiIoKIhffvkFgLS0NAYMGED9+vWxsbEhODiY77777pbXLSgoYNy4cdSvXx9bW1vat2/P1q1byzzez8+Pn376iWXLlqFSqRg6dCgAc+bMITg4GFtbW3x9fRk9ejTZ2SX/QNi5cyfdu3fHxsYGZ2dnIiMjuX79OkOHDuWvv/5i3rx5qFQqVCoV8fHxAPz111+0a9cOS0tL6tWrx5tvvklRUZGxzu7duzNmzBjGjh2Lm5sbkZGRFfip3z0ZYiyEENUsK7+QIUv2cel6Hu5WNkS1aYuHs3wcVzVrCw1fDWlLv4W7iEvNYeiSfawe2QF7q5q7yKLBALm55T8+IQHS0kClglWrisu++w6eeKK4LldXaNCgfHXZ2BTXczt2dnbY2dmxbt06OnTogKWl5U3H6PV6HnroIbKysvj2229p3Lgxx48fR6PRAJCfn094eDjjx4/HwcGBDRs2MGjQIBo3bky7du1Kve6YMWM4fvw4q1atwtvbm7Vr1/Lggw9y5MgRmjRpctPx+/btY/DgwTg4ODBv3jysra0BUKvVfPLJJ/j7+3P+/HlGjx7NG2+8wWeffQZAbGws9913H88++yzz5s3DzMyMP//8E51Ox7x58zh9+jRBQUG89957ALi7u5OYmEjPnj0ZOnQoy5Yt4+TJk4wYMQIrKyveeecdY0xLly5l1KhRinSVqwyyTHq5ZGZm4ujoSEZGBg4ODkqHI4SopYp0xV2Qf51Owd7cgjdadyKwga3SYVVYfj5kZUGXLmBby8JPSMvlsYU7Sc3W0qWJG18NaYuFmfKtr/n5+cTFxeHv74+VlRUAOTlgZ6dMPNnZ5f+3/emnnxgxYgR5eXm0bt2abt268dRTTxESEgLA5s2beeihhzhx4gRNmzYtV52PPPIIzZs3Z/bs2UBxi1lYWBhz584lISGBRo0akZCQgLe3t/GciIgI2rVrxwcffFBqnX369MHJyYlvvvmmzOv++OOPvPDCC6SmpgLw9NNPk5CQwI4dO0o9/t9x3fD222/z008/ceLECeN2hJ999hnjx48nIyMDtVpN9+7dyczM5MCBA2XGUtp74oa7zUuUf8cLIYQJmfrLcf46nYKFWs3IFm1rZfJX2zVwteHroW2xsdCw/UwqE9cdkS3j7lK/fv24fPkyP//8Mw8++CBbt26ldevWxkQrNjYWHx+fMpM/nU7H1KlTCQ4OxsXFBTs7O3777TcSEhJKPf7IkSPodDqaNm1qbIG0s7Pjr7/+4ty5ii38/fvvv3PfffdRv3597O3tGTRoEGlpaeT+0/R6owWwIk6cOEHHjh1L7EXduXNnsrOzuXTp/5cjCg8Pr1C9lUn6HIQQopp8szOOpbuLFyQeGBBGh6ZOygZkwkJ8nFjwdGuGL93H9/svEeBhx/NdGysd1k1sbIpb4ioiNhbuuefm8h07ICysYteuCCsrK+6//37uv/9+Jk2axHPPPceUKVMYOnSosbu1LB9++CHz5s1j7ty5xvF4Y8eOLXNiRHZ2NhqNhpiYGGM38g12FWgyjY+P55FHHmHUqFFMmzYNFxcXduzYwfDhw9FqtdjY2Nw29rthq2DzubQACiFENfjzZDLv/bPLx8M+zXkoqF65xleJqnNvcw8mPly8Zdz0X0/y+/GrCkd0M5WquBu2Io8b+YpaXfKrtXXF6rnb92dgYCA5OTkAhISEcOnSpTKXXtm5cye9e/fmmWeeITQ0lEaNGt1ymZZWrVqh0+lITk4mICCgxMPLy6vcMcbExKDX6/noo4/o0KEDTZs25fLlyyWOCQkJITo6usw6LCws0Ol0JcpatGjB7t27S7Qs79y5E3t7e3x8fModX1WSBFAIIarYyaRMXvruIHoDtHPzYWDrRvyn0UIoZFhnP55u3wCDAV5ZdZATVzKVDumueXiAlxeEh8OiRcVfvbyKy6tCWloaPXr04Ntvv+Xw4cPExcXxww8/MGvWLHr37g1At27d6Nq1K/369WPLli3ExcXx66+/smnTJgCaNGnCli1b2LVrFydOnGDkyJFcvVp2Qt60aVMGDhzI4MGDWbNmDXFxcezdu5fp06ezYcOGcsceEBBAYWEh8+fP5/z58yxfvpxFixaVOGbChAns27eP0aNHc/jwYU6ePMnChQuNYwT9/PzYs2cP8fHxpKamotfrGT16NBcvXuSll17i5MmTrF+/nilTphAVFYVaXTNSr5oRhRBC1FHJWfkM/2Y/2QVFBDi4MKptMJaW0vRXU6hUKt59tCWdGruSo9Xx3NL9pGQVKB3WXfHxgfh42LMHRo4s/hofX1xeFezs7Gjfvj0ff/wxXbt2JSgoiEmTJjFixAg+/fRT43E//fQTbdu2ZcCAAQQGBvLGG28YW84mTpxI69atiYyMpHv37nh5edGnT59bXnfJkiUMHjyY1157jWbNmtGnTx/27dtHg/JOcwZCQ0OZM2cOM2fOJCgoiBUrVjB9+vQSxzRt2pTNmzdz6NAh2rVrR8eOHVm/fj1m/+zVOG7cODQaDYGBgbi7u5OQkED9+vXZuHEje/fuJTQ0lBdeeIHhw4czceLEcsdW1WrsLOAFCxbw4YcfkpSURGhoKPPnzy9zKvjixYtZtmwZR48W76MZHh7OBx98UOJ4g8HAlClTWLx4Menp6XTu3JmFCxeWOlW8NDILWAhRUdoiPQMW/03Mhet4WNvydvtO+LhbKB1WpajNs4BLk56rpe9nxcvDtGrgxHcjOmBlXr3NtLea8SlMk8nNAl69ejVRUVFMmTKFAwcOEBoaSmRkJMnJyaUev3XrVgYMGMCff/7J7t278fX15YEHHiAxMdF4zKxZs/jkk09YtGgRe/bswdbWlsjISPLz6/ZCoEII5bzzv2PEXLiOtZkZo1q2qTPJX13kZGPBV0Pa4GhtzsGEdCaskZnBom6rkQngnDlzGDFiBMOGDSMwMJBFixZhY2PD119/XerxK1asYPTo0YSFhdG8eXO+/PJL9Hq9cdCmwWBg7ty5TJw4kd69exMSEsKyZcu4fPky69atq8Y7E0KYihV7LrByTwIqYFBAK4L9FFrMTZRbI3c7Fg5sjUatYu3BRL7eGa90SEJUmRqXAGq1WmJiYoiIiDCWqdVqIiIi2L17d7nqyM3NpbCwEBcXFwDi4uJISkoqUaejoyPt27cvs86CggIyMzNLPIQQojz2x1/jnZ+PAdDTpxkRgR4y47eW6BTgxls9WwDwwcYT7DqXqnBEQlSNGpcApqamotPp8PT0LFHu6elJUlJSueoYP3483t7exoTvxnkVqXP69Ok4OjoaH76+vhW9FSGECbqSkccL3x6gUGcg1KUeT7dqLDN+a5lnO/vRt1V9dHoDY1Ye5NL1CuzFJkQtUeMSwLs1Y8YMVq1axdq1a+9qEO2ECRPIyMgwPi5evFiJUQoh6qL8Qh0vLI8hNbuA+jb2vBgegpWVNP3VNiqViumPBRNU34FrOVpe+DaG/ELd7U8UohapcQmgm5sbGo3mpvV/rl69etvFHWfPns2MGTPYvHmzcf9BwHheReq0tLTEwcGhxEMIIcpiMBh4e+1RDl3KwNbcnNEhbXB1ks2Waisrcw2LngnHxdaCo4mZ1TopRCafiBuq8r1Q4xJACwsLwsPDS6y6fWNCR8eOHcs8b9asWUydOpVNmzbRpk2bEq/5+/vj5eVVos7MzEz27NlzyzqFEKK8lu6K56cDl1ABQ5q0prlvBffREjWOj7MNnz7dyjgpZEkVTwoxNzcHMO5BK8SN98KN90ZlqpF/nkZFRTFkyBDatGlDu3btmDt3Ljk5OQwbNgyAwYMHU79+feNijTNnzmTy5MmsXLkSPz8/47i+G5tDq1Qqxo4dy/vvv0+TJk3w9/dn0qRJeHt733ahSSGEuJ2YC9d4f8MJAB5t0ILuzd0UjkhUlk6NiyeFTP3lONM2nqB5PXs6Na6af1+NRoOTk5NxyTMbGxtUMnvIJBkMBnJzc0lOTsbJyemm/Y4rQ41MAJ988klSUlKYPHkySUlJhIWFsWnTJuMkjoSEhBJbqSxcuBCtVkv//v1L1DNlyhTeeecdAN544w1ycnJ4/vnnSU9P55577mHTpk2y2KYQ4q5cy9EyZuVBivQGWrl680Sov0z6qGOe7ezH0cQM1h5M5KWVB9nwche8HKvmd8eNYUllrXsrTIuTk1OF9jauiBq7E0hNIzuBCCH+S6c3MHTJXrafScXT2pZ3Ot+Dh3ON/Lu60tW1nUBuJ79QR9/PdnHiSiZtGjrz3fMdMNdU3SgqnU5HYWFhldUvaj5zc/NbtvzdbV5iGp9UQghRBT794yzbz6RioVYzomW4ySR/psjKXMPCga3pNX8H+y9cZ8avJ5n0SGCVXU+j0VRJt58QN1Tap1V6ejpr165l+/btXLhwgdzcXNzd3WnVqhWRkZF06tSpsi4lhBCK23EmlbnRpwHo7x9MmJ+9whGJqubnZsvsJ0IZuTyGr3bEEd7QmZ7B9ZQOS4g7ctft15cvX+a5556jXr16vP/+++Tl5REWFsZ9992Hj48Pf/75J/fffz+BgYGsXr26MmIWQghFJWXk88qqgxgM0N7Nl94hPrLTh4mIbOnFyK6NAHjjx8OcT8lWOCIh7sxdtwC2atWKIUOGEBMTQ2Bg6c3heXl5rFu3jrlz53Lx4kXGjRt3t5cVQghFFOr0vPTdAdJytNS3ceD5Ni0xk55fk/J6ZDMOXkxnb9w1Rn17gLUvdsLGQt4Eona560kgaWlpuLq6VtnxNYVMAhFCAEzfeILPt53HSmPGpHb30NzHBGZAlMLUJoH8V3JmPj0/2UFqdgGPtarPR0+EypItolrdbV5y113AFU3mamPyJ4QQAFuOX+XzbecBGBgQYrLJnwAPByvjItFrDiaycm+C0iEJUSGVMoe9a9eupKenG5///PPP5OXlVUbVQghRIySk5fLa97EAdPX0J7KlDP43dR0aufJ6ZDMA3v35OIcvpSsbkBAVUCkJ4I4dO9BqtcbnzzzzDFeuXKmMqoUQQnH5hTpGr4whM78IPzsnhrdpLos9CwBGdm3EA4GeaHV6Rn17gOs52tufJEQNUCWrWMra0kKIumTqL8c5mpiJrZk5L7ZqjZ1NjdtGXShEpVLx4eOhNHS1ITE9j6jvY9Hr5XegqPnkU0wIIW5hfWwiK/YkoAKGNGtFIy9rpUMSNYyjtTkLB4Zjaabmz1MpLPzrnNIhCXFblTZv/bfffsPR0REAvV5PdHQ0R48eLXHMo48+WlmXE0KIKnc2OYsJa44AEOEdwL3N3RWOSNRUgd4OvNe7JeN/OsJHm0/RqoETnRq7KR2WEGWqlL2A1erbNySqVCp0Ot3dXkoxsgyMEKYlV1tE7093ciY5myYOrkzp1h5rK1nm4wZTXwamNAaDgdd/PMyPMZdws7Nk48v34OFgpXRYoo5SfBkYKG7xu92jNid/QgjTYjAYeHvtUc4kZ+NoYcnoVq0k+RO3pVKpmNo7iOZe9qRmFzDmu4MU6fRKhyVEqWQMoBBC/Md3ey+y9mAiapWKZ5u3poGHpdIhiVrC2kLDgoGtsbXQsDfuGh9tOa10SEKUqlISwNOnT7N3794SZdHR0dx77720a9eODz74oDIuI4QQVe5oYgbv/HwMgJ4+zejc1EXhiERt09jdjpn9QwBYuPUc0SeuKhyREDerlARw/Pjx/PLLL8bncXFx9OrVCwsLCzp27Mj06dOZO3duZVxKCCGqTEZeIaNWxKDV6Wnp5MGAVo2Q3b3EnXgkxJuhnfwAiPr+EBev5SobkBD/USkJ4P79+3nooYeMz1esWEHTpk357bffmDdvHnPnzuWbb76pjEsJIUSVMBgMvP7DIS5ey8PV0ppRrcOwspTsT9y5t3q2INTXiYy8Ql5ceYCCIhkLL2qOSkkAU1NT8fHxMT7/888/6dWrl/F59+7diY+Pr4xLCSFElfhyexybj1/FTKVmeIvW1HM1VzokUctZmKlZ8HQrnGzMOXwpg2kbTigdkhBGlZIAuri4GLd+0+v17N+/nw4dOhhf12q1sjuIEKLG2h9/jRmbTgLQp2Eg7QKclA1I1Bk+zjZ8/GQYAMt2X+DnQ5eVDUiIf1RKAti9e3emTp3KxYsXmTt3Lnq9nu7duxtfP378OH5+fpVxKSGEqFRp2QWMWXkQnd5AKxdv+oU2kHF/olLd28yDF+9tDMCbPx3mbHK2whEJUUkJ4LRp0zh58iQNGzZk/PjxzJo1C9t/rQy6fPlyevToUe76FixYgJ+fH1ZWVrRv3/6mGcb/duzYMfr164efnx8qlarUySbvvPMOKpWqxKN58+YVukchRN2j0xsYuzqWpMx8PK1tGdk6GAsLyf5E5Xs1oikdG7mSq9UxekUMudoipUMSJq5SEkA/Pz9OnDjBwYMHuXDhAqNGjSrx+rvvvsvEiRPLVdfq1auJiopiypQpHDhwgNDQUCIjI0lOTi71+NzcXBo1asSMGTPw8vIqs96WLVty5coV42PHjh3lv0EhRJ30SfQZtp9JxUKjYURgOO7OlbY7phAlmGnUzBsQhru9JaevZjNx7VEZGiUUVWkLQZuZmREaGoq3t/dNr4WGhuLq6lqueubMmcOIESMYNmwYgYGBLFq0CBsbG77++utSj2/bti0ffvghTz31FJaWZS/WamZmhpeXl/Hh5iZ7NAphyradTuGTP84A0L9hEGH+9gpHJOo6D3srPh3QCo1axZqDiazad1HpkIQJu+sEcMaMGeTl5ZXr2D179rBhw4YyX9dqtcTExBAREfH/AarVREREsHv37ruK88yZM3h7e9OoUSMGDhxIQkLCLY8vKCggMzOzxEMIUTdcychj7OpYDAbo4O7LoyE+Mu5PVIv2jVwZ90AzAKb8fIyjiRkKRyRM1V0ngMePH6dBgwaMHj2aX3/9lZSUFONrRUVFHD58mM8++4xOnTrx5JNPYm9f9l/Zqamp6HQ6PD09S5R7enqSlJR0xzG2b9+eb775hk2bNrFw4ULi4uLo0qULWVlZZZ4zffp0HB0djQ9fX987vr4QouYo1OkZs/Ig13K0+Ng6MLxVS8xlxRdRjUZ2bcR9zT3QFukZveIAGXmFSockTNBdJ4DLli3j999/p7CwkKeffhovLy8sLCywt7fH0tKSVq1a8fXXXzN48GBOnjxJ165dKyPuCnnooYd4/PHHCQkJITIyko0bN5Kens73339f5jkTJkwgIyPD+Lh4UZrqhagLZvx6kpgL17E2M+P5lq1xcdQoHZIwMWq1io+eCMXH2ZqEa7m8/sMhGQ8oql2ljHgODQ1l8eLFfP755xw+fJgLFy6Ql5eHm5sbYWFh5R5v5+bmhkaj4erVkvsmXr169ZYTPCrKycmJpk2bcvbs2TKPsbS0vOWYQiFE7fPrkSt8tSMOgAH+oQQ2sL3NGUJUDScbCz4b2Jr+C3ez+fhVvtwex4iujZQOS5iQSpsEAsXj9cLCwujduzdPPfUUERERFZpsYWFhQXh4ONHR0cYyvV5PdHQ0HTt2rLQ4s7OzOXfuHPXq1au0OoUQNVtcag6v/3gYgHu9GhEZ5CXj/oSiQnycmPRICwBmbDrJ/vhrCkckTEmlJoAajabU5VrS0tLQaMrXzRIVFcXixYtZunQpJ06cYNSoUeTk5DBs2DAABg8ezIQJE4zHa7VaYmNjiY2NRavVkpiYSGxsbInWvXHjxvHXX38RHx/Prl276Nu3LxqNhgEDBtzlHQshaoM8rY5R38aQXVBEYwcXBoU1w0xWfBE1wDMdGvJoqDc6vYExKw+Sml2gdEjCRFTqR2BZYxgKCgqwsLAoVx1PPvkkKSkpTJ48maSkJMLCwti0aZNxYkhCQgJq9f/nrZcvX6ZVq1bG57Nnz2b27Nl069aNrVu3AnDp0iUGDBhAWloa7u7u3HPPPfz999+4u7vf4Z0KIWoLg8HAxHVHOZmUhYOFJc8HtcLRvlL/9hXijqlUKqY/FsyxyxmcS8lh7KpYlj7bDo1amqdF1VIZKmHk6SeffALAq6++ytSpU7GzszO+ptPp2LZtG/Hx8Rw8ePBuL6WYzMxMHB0dycjIwMHBQelwhBDltHpfAuN/OoIKGN2iPT2CZA3QypCfD1lZ0KUL2MpQyrt2+moWvT/dSV6hjlfua8Kr9zdVOiRRw91tXlIpLYAff/wxUPyX9qJFi0p091pYWODn58eiRYsq41JCCFFuxy5nMGn9MQAe8mlG1+aS/ImaqamnPdP6BhH1/SE++eMM4Q2d6dpUeqlE1amUBDAurnhW3b333suaNWtwdnaujGqFEOKOZeQVMurbA2iL9AQ6efBUaGMZ9ydqtMda+7Av/jrf7U1g7OpYNrx8D/UcrZUOS9RRlToQ5s8//5TkTwihOIPBwLgfDpFwLRdXK2ueDw3F1kbGVImab0qvQFp6O3AtR8uYlQcp1OmVDknUUZWaAPbr14+ZM2feVD5r1iwef/zxyryUEEKUafH282w5fhUztZphTVvj61G+SWhCKM3KXMNnA1tjb2VGzIXrzPz1pNIhiTqqUhPAbdu20bNnz5vKH3roIbZt21aZlxJCiFLtPpfGzE2nAOjTIJD2TZyUDUiICmroasvsx0MB+HJHHJuOXlE4IlEXVWoCmJ2dXepyL+bm5mRmZlbmpYQQ4iaX0/MYs/IAOr2BcNf69A1ugFpWfBG1UGRLL0Z08Qfg9R8OE5+ao3BEoq6p1I/G4OBgVq9efVP5qlWrCAwMrMxLCSFECfmFOl74Noa0HC0+tg483zoYKysZ9ydqrzcebE6bhs5kFRQxesUB8gt1Sock6pBKnRM3adIkHnvsMc6dO0ePHj0AiI6O5rvvvuOHH36ozEsJIYSRwWBg8vqjHL6Uga25OS8EhePmVL7dh4Soqcw1aj59ujUPf7Kd41cymbL+GDP7hygdlqgjKrUFsFevXqxbt46zZ88yevRoXnvtNS5dusTvv/9Onz59KvNSQghhtGJPAt/vv4QKGBzQiua+NkqHJESl8HK0Yt5TrVCpYPX+i3z79wWlQxJ1RKXsBGIKZCcQIWqmmAvXeeqL3RTqDDzs05zBbWW9v+ogO4FUr8+2nmXWplOYqVWsHNGBdv4uSockFHa3eUmlD49OT0/nyy+/5K233uLatWsAHDhwgMTExMq+lBDCxCVn5jPq2xgKdQbCXLwYENZIkj9RJ43q1phHQupRpDcwekUMl9PzlA5J1HKVmgAePnyYpk2bMnPmTD788EPS09MBWLNmDRMmTKjMSwkhTJy2SM/oFQdIziqgno0dz4eFYm0tkz5E3aRSqZjVP4QW9RxIzdYycnmMTAoRd6VSE8CoqCiGDh3KmTNnsLKyMpb37NlT1gEUQlSqqb8cZ/+F61ibmfF8YBs8XaXpT9RtNhZmfDEoHGcbc44kZjBhzRFkFJe4U5WaAO7bt4+RI0feVF6/fn2SkpIq81JCCBP27d8XWP7PYPhnGocR7CeD0IRp8HWxYcHA1mjUKtYeTOSrHXFKhyRqqUpNAC0tLUtd8Pn06dO4u7tX5qWEECZq17lU3vn5GAAP1W9GRKAnKun5FSakU2M3Jj7cAoAPNp5gx5lUhSMStVGlJoCPPvoo7733HoWFhUDxmIWEhATGjx9Pv379KvNSQggTdCEth9ErDlCkN9Da1ZuBrWTGrzBNQzv50a+1D3oDjPnuAAlpuUqHJGqZSk0AP/roI7Kzs/Hw8CAvL49u3boREBCAvb0906ZNq8xLCSFMTFZ+Ic8t3U96biEN7Rx5oXWITPoQJkulUjGtbxChvk6k5xby3LJ9ZOUXKh2WqEUq9W9nR0dHtmzZws6dOzl06BDZ2dm0bt2aiIiIyryMEMLE6PQGXlkVy5nkbJwsLRkd0gZX2elDmDgrcw2fPxPOo5/u4PTVbMasPMhXQ9pgppENsMXt3fW7xMXFhdTU4vEHzz77LFlZWXTu3JnRo0fzxhtvSPInhLhrszad5I+TyZir1TzXrA0B9a1uf5IQJsDL0YqvhrTFylzNX6dTeH/DCaVDErXEXSeAWq3WOPFj6dKl5Ofn33VQQghxw48xl/h823kAnvIPpX0TJ2UDEqKGCfZxZO6TYQB8syueZbvjFY1H1A53nQB27NiRPn36MGzYMAwGAy+//DLPPvtsqY/yWrBgAX5+flhZWdG+fXv27t1b5rHHjh2jX79++Pn5oVKpmDt37l3XKYSoGXafS2PCmsMARNQL4JEQb9TSuyXETR4MqscbDzYD4J2fj7H1VLLCEYma7q4/Sr/99lt69uxJdnY2ABkZGVy/fr3UR3msXr2aqKgopkyZwoEDBwgNDSUyMpLk5NLfzLm5uTRq1IgZM2bg5eVVKXUKIZR3NjmLkcv3U6gzEOrixZDwpjLjV4hbGNWtMf3Di2cGv7TyIKevZikdkqjBVIZKXEbc39+f/fv34+rqesd1tG/fnrZt2/Lpp58CoNfr8fX15aWXXuLNN9+85bl+fn6MHTuWsWPHVlqdN9ztpstCiPJLySqg72c7uXQ9D397JyZ06CCTPmqY/HzIyoIuXcBW1uGuMbRFep75ag97467h42zNuhc742ZnqXRYogrcbV5SqZNA7r33XiwsLO64Lq1WS0xMTImJI2q1moiICHbv3l2tdRYUFJCZmVniIYSoernaIp5buo9L1/PwsLZhTGhbSf6EKCcLMzWfPxOOn6sNl67n8dzS/eRpZc9gcbMaNQkkNTUVnU6Hp6dniXJPT8873kruTuucPn06jo6Oxoevr+8dXV8IUX43lns5dCkDO3NzRrVsh1+9O/+jUghT5GxrwVdD2+JkY07sxXRe+u4ARTq90mGJGuauR9TcmAQSHh5unARibW1d6rFff/313V6u2kyYMIGoqCjj88zMTEkChahi7284zpbjVzFXqxnetI3s8SvEHWrsbseXg9sw8Ms9/H4imUnrj/JB32BUsm+i+EelTgJRqVR3NQnEzc0NjUbD1atXS5RfvXq1zAkeVVWnpaUlDg4OJR5CiKrz5fbzLNkZD8AA/1Duae4ie/wKcRfa+LnwyYBWqFXw3d6LzP/jrNIhiRrkrlsAPT09mTFjBlA8CWT58uV3PAnEwsKC8PBwoqOj6dOnD1A8YSM6OpoxY8bUmDqFEJVrzYFLxgVsH/FtzsOy3IsQlSKypRfv9g5i0rqjzNlyGk8HS55s20DpsEQNUCkfsT179iQjI4O4uDhcXV2ZMWMG6enpxtfT0tIIDAwsV11RUVEsXryYpUuXcuLECUaNGkVOTg7Dhg0DYPDgwUyYMMF4vFarJTY2ltjYWLRaLYmJicTGxnL27Nly1ymEUM4fJ6/y+o/Fa/118/Tn6VaNZLkXISrRoA4NefHexgC8tfYof5y8epszhCmolGVg1Go1SUlJeHh4AODg4EBsbCyNGjUCirtbvb290enKNxPp008/5cMPPyQpKYmwsDA++eQT2rdvD0D37t3x8/Pjm2++ASA+Ph5/f/+b6ujWrdv/tXfncVGW+//HXzNsAzKAgqyyKCAugKTIZmUdSdTqaIuZdXKpo0ePmkVWYpna6RyPrZqZ6ze1PKZHKzV/Rcc4xzRB3EllcUVcANlX2Yb79wc5RaK5ADMwn+fjcT9grvu67/mME9N7rvu+r5udO3fe1D5vhkwDI0TzO3iukKdXJVNVW0+okwfPh/ehg40c920LZBqYtkVRFF7alMKXhy5ibWHGvyaE09ero6HLEnfgTnNJiwRArVZLSkrKbQdAYyQBUIjmlZFTxshliZRW1dHToTMzwkNxsJPjvm2FBMC2p1ZXz3NrD7DrRB52GnM2TIykl7v8/6ytMvg8gEIIcasuFFUy5pNkSqvq6Kp14Pm+fSX8CdHCLMzULPtTX/p5d6S0qo4xnyRzJq/c0GUJA2mWT1yVSnXNpeVyqbkQoim5pVU8vSqZ3NJq3Gxsmd63P86OctKfEK3BxtKcT8b1p5ebHfnlNfxpVTIXiioNXZYwgGb51FUUhXHjxmFl1XC7maqqKiZNmkSHn48LVFdXN8fTCCHauLyyap5auZdzBZU4aayZHhKOp7NM9CxEa7K3tuCz58J4YnkSp/Mq+NOqZP49KRJnrcbQpYlW1CznAN7s1bSrV6++06cyGDkHUIg7U1hRw+gVe8nILaOTlYbYkEh6etkYuixxm+QcwLYvu+QKjy9N4mLxFXq4atkwMQIHG/lC1lYYxUUgpkACoBC3r6SylqdW7eX4pVLsLa14ITiS4K6SGtoyCYDtw7mCCkYuS+JyWTWBHnasey5cQmAbIReBCCGMWllVLWNW7+P4pVLsLC2ZFhgh4U8II+Ht2IF1fw7HsYMlxy6W8vSqZIorawxdlmgFEgCFEC2mrKqWcav3k3K+GFsLC6b0CucuX1tDlyWE+JXuLlrWT4jAsYMlxy9JCDQVEgCFEC2ipLKWP/3fPg6eK8LG3JzJPcPp5yenTwhhjAJctXw+MQIn24YQ+NTKZIoqJAS2ZxIAhRDNrrCihqdW7dWP/E3tHUF4d3tkdighjFd3Fy2fT2gIganZDSOBEgLbLwmAQohmlVdWzegVe3855693BGH+Ev6EaAv89SHQitTsUkav3Mvl0ipDlyVagARAIUSzySmpYtSKJDJyy7C3tOL5wEhC/e0k/AnRhvi7aNkwMRxnrRXpOWWMXJ7E+UKZLLq9kQAohGgWZ/LKeXxZImfyKuhkZc0LwZFywYcQbZSfs5bNk6Lw7GTNuYJKHl+WyMncMkOXJZqRBEAhxB1LOV/M48uSuFB0hc7WNsSGyFQvQrR1Xo42bJ4URXcXW3JLqxm5PImU88WGLks0EwmAQog7svtkHqNX7qWwogYvW3te7Rcld/gQop1wsdOwcWIkfTwdKK6s5amVe/nxZL6hyxLNQAKgEOK2bT1ykWfX7KeyRkeAvROv9I+gq5uVocsSQjSjjh0s+defw4nydaSiRse41fvYfPCCocsSd0gCoBDilimKwspdZ5i+4Qi1OoW7HN15JaI/bk7mhi5NCNECbK3M+WRcfx7u405dvcKMTSl8sOMEcjfZtksCoBDiltTU1RP35VH+/k0aAHe7+PBiRAgOdvJxIkR7prEwY9GoECbf5wvAooSTzNj0EzV19QauTNwO+bouhLhpxZU1TF53iKQzBaiA4V69eCLEBysrmedFCFOgVqt4dUgPPDvaMHvrMb44dIHskissfbof9jYWhi5P3AL5yi6EuCln8sp55ONEks4UoDEz47nuoYzu11XCnxAm6KlwL1aNDcXG0ozE0wUMX/KjTBPTxhhtAFyyZAk+Pj5oNBrCw8PZt2/fDftv2rSJHj16oNFoCAoK4ptvvmm0fty4cahUqkbLkCFDWvIlCNFu7DqRxyMfJ3I2v2GOvxeDohgS7IK5HEMQwmTdH+DMpkmReDhYk1lQySMfJ7IjNdfQZYmbZJQBcOPGjcTGxjJnzhwOHTpEnz59iImJ4fLly032T0xMZPTo0Tz33HMcPnyYESNGMGLECI4dO9ao35AhQ8jOztYvn3/+eWu8HCHarPp6hY/+e5Kxq/dRcqUWH1sH4kIHyN09hBAA9Ha3Z9vUAYR37UR5dR0TPj3Ahwknqa+Xi0OMnUoxwkt4wsPD6d+/Px999BEA9fX1eHp6Mm3aNGbOnHlN/1GjRlFRUcH27dv1bREREYSEhLBs2TKgYQSwuLiYLVu23FZNpaWl2NvbU1JSgp2d3W3tQ4i2pORKLS/9O4Xv0xq+0Ud09uS5u3rTyd7MwJUJY1BVBWVlcM890EHm/DZ5tbp63tqeytqkcwAM6e3K2yODsdPIeYEt5U5zidGNANbU1HDw4EGio6P1bWq1mujoaJKSkprcJikpqVF/gJiYmGv679y5E2dnZwICApg8eTIFBQXN/wKEaAfSc0oZ/tGPfJ+Wi4VazRM+QTwfGSzhTwjRJAszNfOGB7LgsSAszdTEH8/h4cU/cvRCiaFLE9dhdAEwPz8fnU6Hi4tLo3YXFxdycnKa3CYnJ+d3+w8ZMoRPP/2UhIQEFixYwA8//MDQoUPR6XRN7rO6uprS0tJGixDtnaIofL4vixFL9pBZUEknK2umB0Yysp8XVjK/sxDid4zq78XGv0Tg4dBwD+HHliayZs9ZmS/QCJnMKdxPPvmk/vegoCCCg4Px9fVl586dDBo06Jr+8+fPZ968ea1ZohAGVVxZQ9yXR/n2WMMXpx72TkwMvgtvV0sDVyaEaEvu8urIN8/fw8ubU/hPai5zv05l75lCFjwejL21HBI2FkY3Aujk5ISZmRm5uY2vJMrNzcXV1bXJbVxdXW+pP0C3bt1wcnLi1KlTTa6Pi4ujpKREv5w/f/4WX4kQbcfeMwUMXbSbb4/lYKZS8ZBnD167J0zCnxDittjbWLD8mX7MebgXFmYq4o/nMGzRbhJPy32EjYXRBUBLS0v69etHQkKCvq2+vp6EhAQiIyOb3CYyMrJRf4AdO3Zctz/AhQsXKCgowM3Nrcn1VlZW2NnZNVqEaG+q63S88106o1fuJbukCmfrDsQGRTGmvy821nKZrxDi9qlUKsYP6MoXk6Pw6mTDxeIrPLUymTe/TqWqtunTr0TrMboACBAbG8vKlStZu3YtaWlpTJ48mYqKCsaPHw/AmDFjiIuL0/efPn068fHxvPfee6SnpzN37lwOHDjA1KlTASgvL+fll19m7969ZGZmkpCQwPDhw/Hz8yMmJsYgr1EIQ0s5X8zDi39kyf9OoygQ1rkLcyLuJiLAATO51kMI0UyCuzjwzfR7GB3mBcAne87y4Ie7STlfbNjCTJxRngM4atQo8vLyeOONN8jJySEkJIT4+Hj9hR5ZWVmo1b9k16ioKNavX8/rr7/OrFmz8Pf3Z8uWLQQGBgJgZmbGTz/9xNq1aykuLsbd3Z3Bgwfzt7/9DSs5s12YmKpaHQu/P8mKXaepV0BrYcmj3oEMDXTDQk7PEUK0AFsrc+Y/GsTgXi68+sVPnM6r4NGliUwe6MvUP/ihsZBvna3NKOcBNEYyD6BoDw5kFuo/fAH6Orrzp1695Vw/cctkHkBxu4oqapi99Rjbf8oGwMfRhr8/EsQAPycDV9a23GkukQB4kyQAirYsv7yaBd+ms+ngBQDsLK143DuQB3q5YinZT9wGCYDiTsUfy2bOtuPkllYD8OhdHrz2YE8cbeXI3M2401xilIeAhRDNQ1evsD75HO98l0FpVR3QcK7f07170qWzJD8hhOEMCXRjgJ8T736Xwad7z/Hl4Yv8N+MyL8cEMCrUE3Mzo7xMod2QEcCbJCOAoq3Zn1nIvK+Pc+xiwyTmXTrYMbJbIBF+HTGXr37iDskIoGhOR84XE/flUdKyGz6vAly0vP5QT+7x72zgyoyXHAJuJRIARVtx6nI5C+LT2ZHaMDemtbk5Qz0C+GMvb7S2MrWLaB4SAEVzq9PVs27vORYmnKS4shaAQT2cmfVgT3w72xq4OuMjAbCVSAAUxu5yWRULvz/Jxv3n0dUrqFUqwpw8GdmjO94uVqgk+4lmJAFQtJTiyhoWJZzks6Rz1NUrmKtVjAz1ZNof/HB3sDZ0eUZDAmArkQAojFV+eTUrd5/hs6RzVNY0TK4a1NGFEb4BBHpq5XCvaBESAEVLO51Xzj/+XxoJ6ZcBsDRT81S4F3+93xdnrcbA1RmeBMBWIgFQGJvLpVWs2HWGdcnnqKqtB8DH1oHhPj2J8O0kV/eKFiUBULSW/ZmFvPtdBslnCwHQWKgZE+nDn+/uirOd6QZBCYCtRAKgMBZZBZV8sucsn+/LorquIfh52zoQ08WPe/2csZZbuIlWIAFQtCZFUUg8XcC7/8ngcFYx0DAi+GhfDybc280kzxGUANhKJAAKQ1IUhf2ZRfzfj2fYkZpL/c9/tV21Dgzp0p0Bvk4S/ESrkgAoDEFRFP6XcZmP/3eaA+eKAFCpYHAvFybe60tfLwdUJnLCs8wDKEQ7VlWr49tj2azek8lPF0r07T3sOxPdpRuR3RzRaEzjw04IIVQqFX/o4cIferhwILOQZT+c5vu0y3x3PJfvjufS292OZyK8+WOIOzaWEnFuREYAb5KMAIrWlJZdysb95/ny0AX9BM4WajX9HD0Y7N2V3nJxhzAwGQEUxuJEbhkrdp1hW8olan4+LUarMefxfl14KswLfxetgStsGXIIuJVIABQtrbiyhm+P5bBh/3lSzhfr2ztZWRPu5MlgXy+6dLZCLZPjCyMgAVAYm6KKGjYdPM+6vVlkFVbq24M87HmsrwcP93FvV7eZkwDYSiQAipZQUV3HjtRcvk65xK6TedTqGv4czVQqgjq6EOXqRbiPE7Yd5DCvMC4SAIWxqq9X2HUyj3V7s9iZcZm6n0+aNleruC/AmRF3uXN/gDMdrNr2YRQ5B1CINqakspadJy7zn+O5JKTn6qdwAfCw0dLX0YP7u3ahi5MVZmYGLFQIIdog9c9B774AZwrKq/k65RJfHr7ITxdK+D4tl+/TcrEyV3Nv984M6e1KdE8X7G0sDF12q5MRwJskI4DiTpzNr+D71IYPngPnitDV//Jn11ljw12O7kR6uNPDTSvz94k2QUYARVtzMreMLw9f5Juj2Zwr+OUQsblaRaSvIwO7d2Zg9874Odu2iSuJ5RBwK5EAKG5FQXk1SWcKSDxdQOKpfDJ/9WED4GZjS097F8Ld3ejtZidTuIg2RwKgaKsURSE9p4z4YznEH8shI7es0Xo3ew0Du3fm3u6difJ1xMHGOL+VSwBsJRIAxY3klFRxKKuIA5lFJJ7OJz2n8QeKmUqFr50jvR2c6e/uQldnGxnpE22aBEDRXpzJK+e/6Zf54UQeyWcL9VcSX9XdxZb+Pp0I69qJ/j6djOZ+xBIAW4kEQHFVZU0dadllHM4q4nBWMYeyisguqbqmn4eNFl87J3p1cqSPeycc7SzknD7RbkgAFO3RlRodyWcL2HUin10n8zh1ufyaPl06WnOXV0eCPewJ9LAn0MMOrab1zyGUi0CEaCGKopBdUkVadunPSxlp2aWcLajgt1+bVIBHBzu8OjjQo6Mjwa6OuHa0wsL0zisWQog2y9rSTH8BCUB+eTUHMgvZd7aI/ZmFHL9UwoWiK1wousLXKZf023Vz6kCghz293O0IcNXSw1WLq53GqM8llBHAmyQjgKYlt7SKwR/souRKbZPrtRZWeNs64GPrgH+njgR0tqej1lwmZxYmQ0YAhSkqr67jcFYRP10o4eiFEo5eLOFi8ZUm+9ppzIl/4d4WO2TcbkcAlyxZwjvvvENOTg59+vRh8eLFhIWFXbf/pk2bmD17NpmZmfj7+7NgwQKGDRumX68oCnPmzGHlypUUFxczYMAAli5dir+/f2u8HGGkFAVqa6G6umGpqmpYSkqtqK6pR61S4WJti7uNFg8bO7zs7PB1tMPF3gqr9jOfqBBCiJtga2XOPf6duce/s76tsKKGoxdLOHaxhLTsUjJyyjiTX0F1XT0udhoDVntjRhkAN27cSGxsLMuWLSM8PJyFCxcSExNDRkYGzs7O1/RPTExk9OjRzJ8/n4ceeoj169czYsQIDh06RGBgIABvv/02H374IWvXrqVr167Mnj2bmJgYUlNT0WiM9w0Sd6auriHg1dZCTc0vP6uqoLwcKioaHl9dV1/fcGNxCwsVr4fdjZudNbY2ZjKyJ4QQokmdOljqp5C5qrpOx8WiK5ip5RDwLQkPD6d///589NFHANTX1+Pp6cm0adOYOXPmNf1HjRpFRUUF27dv17dFREQQEhLCsmXLUBQFd3d3XnrpJWbMmAFASUkJLi4urFmzhieffPJ3a5JDwIanKA2BTqdr/PPXS01Nw0jelSsNS21t4/VX/2tvCHkNi6XlLz+N+HQNIYyKHAIWwrDa3SHgmpoaDh48SFxcnL5NrVYTHR1NUlJSk9skJSURGxvbqC0mJoYtW7YAcPbsWXJycoiOjtavt7e3Jzw8nKSkpJsKgOL3KUrDCJqi/PL71cdN/f7rRadr/PPqqN2vR+/q6hrWXf2p0/3SH34Jb+bmYGbW8NPcHGxsGgKeubkEPCGEEAKMMADm5+ej0+lwcXFp1O7i4kJ6enqT2+Tk5DTZPycnR7/+atv1+vxWdXU11dXV+sclJSVAQ+JuKYoCFy5Afn6LPcV1n/fqyNhvx4N/u+63v99o+XUYhMYB8EahsKkxaTMzUKsbfv76d7X6l1B3dbuamub99xFCXKu6umHUvLS04YuYEKJ1Xc0jt3sg1+gCoLGYP38+8+bNu6bd09PTANUIIYQQQlyrrKwMe3v7W97O6AKgk5MTZmZm5ObmNmrPzc3F1dW1yW1cXV1v2P/qz9zcXNzc3Br1CQkJaXKfcXFxjQ4r19fXU1hYiKOjo1HP62MMSktL8fT05Pz583K+pBGQ98O4yPthXOT9MC7yftw8RVEoKyvD3d39trY3ugBoaWlJv379SEhIYMSIEUBD+EpISGDq1KlNbhMZGUlCQgIvvPCCvm3Hjh1ERkYC0LVrV1xdXUlISNAHvtLSUpKTk5k8eXKT+7SyssLqN/N8ODg43NFrMzV2dnbyB2xE5P0wLvJ+GBd5P4yLvB8353ZG/q4yugAIEBsby9ixYwkNDSUsLIyFCxdSUVHB+PHjARgzZgweHh7Mnz8fgOnTpzNw4EDee+89HnzwQTZs2MCBAwdYsWIFACqVihdeeIG33noLf39//TQw7u7u+pAphBBCCGEqjDIAjho1iry8PN544w1ycnIICQkhPj5efxFHVlYWarVa3z8qKor169fz+uuvM2vWLPz9/dmyZYt+DkCAV155hYqKCiZOnEhxcTF333038fHxMgegEEIIIUyOUc4DKNq26upq5s+fT1xc3DWH0UXrk/fDuMj7YVzk/TAu8n60HgmAQgghhBAmRv37XYQQQgghRHsiAVAIIYQQwsRIABRCCCGEMDESAEWrqK6uJiQkBJVKxZEjRwxdjknKzMzkueeeo2vXrlhbW+Pr68ucOXOokXvntaolS5bg4+ODRqMhPDycffv2GbokkzR//nz69++PVqvF2dmZESNGkJGRYeiyBPDPf/5TP32baDkSAEWreOWVV257tnLRPNLT06mvr2f58uUcP36cDz74gGXLljFr1ixDl2YyNm7cSGxsLHPmzOHQoUP06dOHmJgYLl++bOjSTM4PP/zAlClT2Lt3Lzt27KC2tpbBgwdTUVFh6NJM2v79+1m+fDnBwcGGLqXdk6uARYv79ttviY2N5YsvvqB3794cPnz4urfgE63rnXfeYenSpZw5c8bQpZiE8PBw+vfvz0cffQQ03OXI09OTadOmMXPmTANXZ9ry8vJwdnbmhx9+4N577zV0OSapvLycvn378vHHH/PWW28REhLCwoULDV1WuyUjgKJF5ebmMmHCBD777DNsbGwMXY74jZKSEjp16mToMkxCTU0NBw8eJDo6Wt+mVquJjo4mKSnJgJUJaPhbAOTvwYCmTJnCgw8+2OhvRLQco7wTiGgfFEVh3LhxTJo0idDQUDIzMw1dkviVU6dOsXjxYt59911Dl2IS8vPz0el0+jsaXeXi4kJ6erqBqhLQMBL7wgsvMGDAgEZ3kBKtZ8OGDRw6dIj9+/cbuhSTISOA4pbNnDkTlUp1wyU9PZ3FixdTVlZGXFycoUtu1272/fi1ixcvMmTIEEaOHMmECRMMVLkQxmHKlCkcO3aMDRs2GLoUk3T+/HmmT5/Ov/71L7k9ayuScwDFLcvLy6OgoOCGfbp168YTTzzB119/jUql0rfrdDrMzMx4+umnWbt2bUuXahJu9v2wtLQE4NKlS9x3331ERESwZs2aRvfVFi2npqYGGxsbNm/ezIgRI/TtY8eOpbi4mK1btxquOBM2depUtm7dyq5du+jatauhyzFJW7Zs4ZFHHsHMzEzfptPpUKlUqNVqqqurG60TzUMCoGgxWVlZlJaW6h9funSJmJgYNm/eTHh4OF26dDFgdabp4sWL3H///fTr149169bJh2orCw8PJywsjMWLFwMNhx69vLyYOnWqXATSyhRFYdq0aXz11Vfs3LkTf39/Q5dkssrKyjh37lyjtvHjx9OjRw9effVVOSzfQuQcQNFivLy8Gj22tbUFwNfXV8KfAVy8eJH77rsPb29v3n33XfLy8vTrXF1dDViZ6YiNjWXs2LGEhoYSFhbGwoULqaioYPz48YYuzeRMmTKF9evXs3XrVrRaLTk5OQDY29tjbW1t4OpMi1arvSbkdejQAUdHRwl/LUgCoBAmYseOHZw6dYpTp05dE8DlQEDrGDVqFHl5ebzxxhvk5OQQEhJCfHz8NReGiJa3dOlSAO67775G7atXr2bcuHGtX5AQrUwOAQshhBBCmBg5+1sIIYQQwsRIABRCCCGEMDESAIUQQgghTIwEQCGEEEIIEyMBUAghhBDCxEgAFEIIIYQwMRIAhRBCCCFMjARAIYQQQggTIwFQCCGEEMLESAAUQpiUNWvW4ODgcMf7mTt3Li4uLqhUKrZs2XLH+2sts2fPZuLEic22v5qaGnx8fDhw4ECz7VMI0fIkAAohxC1KS0tj3rx5LF++nOzsbIYOHXrH+2yuYHojOTk5LFq0iNdee03flpeXx+TJk/Hy8sLKygpXV1diYmLYs2ePvo+Pjw8qlarRcvV+0paWlsyYMYNXX321RWsXQjQvc0MXIIQQbc3p06cBGD58OCqVysDVNKbT6VCpVKjV136/X7VqFVFRUXh7e+vbHnvsMWpqali7di3dunUjNzeXhIQECgoKGm375ptvMmHCBP1jMzMz/e9PP/00L730EsePH6d3794t8KqEEM1NRgCFEEZn8+bNBAUFYW1tjaOjI9HR0VRUVOjXf/LJJ/Tu3RsrKyvc3NyYOnWqft37779PUFAQHTp0wNPTk7/+9a+Ul5ff8Pm2bt1K37590Wg0dOvWjXnz5lFXV9dk37lz5/Lwww8DoFar9QFw//79PPDAAzg5OWFvb8/AgQM5dOhQo22Li4v5y1/+gouLCxqNhsDAQLZv387OnTsZP348JSUl+hG2uXPnAlBUVMSYMWPo2LEjNjY2DB06lJMnT+r3eXXkcNu2bfTq1QsrKyuysrKarH3Dhg362q/Ws3v3bhYsWMD999+Pt7c3YWFhxMXF8cc//rHRtlqtFldXV/3SuXNn/bqOHTsyYMAANmzYcMN/ZyGE8ZAAKIQwKtnZ2YwePZpnn32WtLQ0du7cyaOPPoqiKAAsXbqUKVOmMHHiRI4ePcq2bdvw8/PTb69Wq/nwww85fvw4a9eu5b///S+vvPLKdZ9v9+7djBkzhunTp5Oamsry5ctZs2YNf//735vsP2PGDFavXq2vNTs7G4CysjLGjh3Ljz/+yN69e/H392fYsGGUlZUBUF9fz9ChQ9mzZw/r1q0jNTWVf/7zn5iZmREVFcXChQuxs7PT73PGjBkAjBs3jgMHDrBt2zaSkpJQFIVhw4ZRW1urr6myspIFCxawatUqjh8/jrOz8zV1FxYWkpqaSmhoqL7N1tYWW1tbtmzZQnV19U29P9cTFhbG7t2772gfQohWpAghhBE5ePCgAiiZmZlNrnd3d1dee+21m97fpk2bFEdHR/3j1atXK/b29vrHgwYNUv7xj3802uazzz5T3NzcrrvPr776Svm9j0+dTqdotVrl66+/VhRFUb777jtFrVYrGRkZTfb/bV2KoignTpxQAGXPnj36tvz8fMXa2lr597//rd8OUI4cOXLDeg4fPqwASlZWVqP2zZs3Kx07dlQ0Go0SFRWlxMXFKSkpKY36eHt7K5aWlkqHDh30y6JFixr1WbRokeLj43PDGoQQxkNGAIUQRqVPnz4MGjSIoKAgRo4cycqVKykqKgLg8uXLXLp0iUGDBl13+++//55Bgwbh4eGBVqvlmWeeoaCggMrKyib7p6Sk8Oabb+pHw2xtbZkwYQLZ2dnX3aYpubm5TJgwAX9/f+zt7bGzs6O8vFx/OPbIkSN06dKF7t273/Q+09LSMDc3Jzw8XN/m6OhIQEAAaWlp+jZLS0uCg4NvuK8rV64AoNFoGrU/9thjXLp0iW3btjFkyBB27txJ3759WbNmTaN+L7/8MkeOHNEvY8aMabTe2tr6lv69hBCGJQFQCGFUzMzM2LFjB99++y29evVi8eLFBAQEcPbsWaytrW+4bWZmJg899BDBwcF88cUXHDx4kCVLlgAN05U0pby8nHnz5jUKN0ePHuXkyZPXhKUbGTt2LEeOHGHRokUkJiZy5MgRHB0d9c/7e7XfCWtr69+9GMXJyQlAH6Z/TaPR8MADDzB79mwSExMZN24cc+bMuWZ7Pz8//fLbK5YLCwsbnRcohDBuEgCFEEZHpVIxYMAA5s2bx+HDh7G0tOSrr75Cq9Xi4+NDQkJCk9sdPHiQ+vp63nvvPSIiIujevTuXLl264XP17duXjIyMRuHm6tLUlbTXs2fPHp5//nmGDRumv0AlPz9fvz44OJgLFy5w4sSJJre3tLREp9M1auvZsyd1dXUkJyfr2woKCsjIyKBXr143XRuAr68vdnZ2pKam/m7fXr16Nbro5mYcO3aMu+6665a2EUIYjkwDI4QwKsnJySQkJDB48GCcnZ1JTk4mLy+Pnj17Ag1X4U6aNAlnZ2eGDh1KWVkZe/bsYdq0afj5+VFbW8vixYt5+OGH2bNnD8uWLbvh873xxhs89NBDeHl58fjjj6NWq0lJSeHYsWO89dZbN123v78/n332GaGhoZSWlvLyyy83GvUbOHAg9957L4899hjvv/8+fn5+pKeno1KpGDJkCD4+PpSXl5OQkECfPn2wsbHB39+f4cOHM2HCBJYvX45Wq2XmzJl4eHgwfPjwW/p3VavVREdH8+OPPzJixAigIUyOHDmSZ599luDgYLRaLQcOHODtt9++5f3v3r2bv/3tb7e0jRDCgAx9EqIQQvxaamqqEhMTo3Tu3FmxsrJSunfvrixevLhRn2XLlikBAQGKhYWF4ubmpkybNk2/7v3331fc3NwUa2trJSYmRvn0008VQCkqKlIUpemLLeLj45WoqCjF2tpasbOzU8LCwpQVK1Zct8amLgI5dOiQEhoaqmg0GsXf31/ZtGmT4u3trXzwwQf6PgUFBcr48eMVR0dHRaPRKIGBgcr27dv16ydNmqQ4OjoqgDJnzhxFURSlsLBQeeaZZxR7e3v9azpx4oR+m6Zez/V88803ioeHh6LT6RRFUZSqqipl5syZSt++fRV7e3vFxsZGCQgIUF5//XWlsrJSv91vX8dvJSYmKg4ODo22EUIYN5Wi/Dy3ghBCiHZNURTCw8N58cUXGT16dLPtd9SoUfTp04dZs2Y12z6FEC1LzgEUQggToVKpWLFixXUnub4dNTU1BAUF8eKLLzbbPoUQLU9GAIUQQgghTIyMAAohhBBCmBgJgEIIIYQQJkYCoBBCCCGEiZEAKIQQQghhYiQACiGEEEKYGAmAQgghhBAmRgKgEEIIIYSJkQAohBBCCGFiJAAKIYQQQpgYCYBCCCGEECZGAqAQQgghhImRACiEEEIIYWIkAAohhBBCmBgJgEIIIYQQJub/A0hUyd5epyvhAAAAAElFTkSuQmCC", + "text/html": [ + "\n", + "
\n", + "
\n", + " Figure\n", + "
\n", + " \n", + "
\n", + " " + ], + "text/plain": [ + "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "gauss = Gauss(mean=0, width=1)\n", + "\n", + "fig, axs = plt.subplots(2)\n", + "\n", + "unc = 0.5 # 50% uncertainty\n", + "effect = Gauss(mean=1, width=1 + unc)\n", + "\n", + "linsp = lambda max_x: np.linspace(-4, max_x, 1000)\n", + "\n", + "\n", + "def sf(x):\n", + " gx = Gauss(mean=1.0, width=1.0 + unc)\n", + " g1 = Gauss(mean=1.0, width=1.0)\n", + " return gx.inv_cdf(g1.cdf(x + 1))\n", + "\n", + "\n", + "x = linsp(4)\n", + "axs[0].plot(x, gauss.pdf(x), label=\"Gauss\")\n", + "axs[1].plot(x + 1, effect.pdf(x + 1), label=\"Effect\")\n", + "\n", + "param_art = axs[0].plot(\n", + " [0.0], gauss.pdf(0.0), marker=\"*\", color=\"red\", label=\"Nuisance parameter\"\n", + ")\n", + "param_cdf_art = axs[0].fill_between(\n", + " linsp(0),\n", + " gauss.pdf(linsp(0)),\n", + " color=\"b\",\n", + " alpha=0.2,\n", + " label=f\"CDF: {gauss.cdf(0):.4f}\",\n", + ")\n", + "\n", + "sf_art = axs[1].plot(\n", + " [sf(0.0)], effect.pdf(sf(0.0)), marker=\"*\", color=\"green\", label=\"Scale factor\"\n", + ")\n", + "sf_cdf_art = axs[1].fill_between(\n", + " linsp(0) + 1,\n", + " effect.pdf(linsp(0) + 1),\n", + " color=\"b\",\n", + " alpha=0.2,\n", + " label=f\"CDF: {effect.cdf(1):.4f}\",\n", + ")\n", + "\n", + "\n", + "@widgets.interact(nuisance=widgets.FloatSlider(min=-4, max=4, step=0.01, value=0.0))\n", + "def update(nuisance):\n", + " # Plot the nuisance parameter on the gauss\n", + "\n", + " print(f\"Nuisance parameter: {nuisance:.2f}\")\n", + " print(f\"Scale factor: {sf(nuisance):.4f}\")\n", + " print(f\"Constraint (logpdf): {gauss.logpdf(nuisance):.4f}\")\n", + " print(f\"Constraint CDF: {gauss.cdf(nuisance):.4f}\")\n", + " print(f\"Effect CDF: {effect.cdf(sf(nuisance)):.4f}\")\n", + "\n", + " global param_art, param_cdf_art, sf_art, sf_cdf_art\n", + " param_art[0].remove()\n", + " param_cdf_art.remove()\n", + " sf_art[0].remove()\n", + " sf_cdf_art.remove()\n", + " param_art = axs[0].plot(\n", + " [nuisance],\n", + " gauss.pdf(nuisance),\n", + " marker=\"*\",\n", + " color=\"red\",\n", + " label=\"Nuisance parameter\",\n", + " )\n", + " param_cdf_art = axs[0].fill_between(\n", + " linsp(nuisance), gauss.pdf(linsp(nuisance)), color=\"b\", alpha=0.2\n", + " )\n", + " sf_art = axs[1].plot(\n", + " [sf(nuisance)],\n", + " effect.pdf(sf(nuisance)),\n", + " marker=\"*\",\n", + " color=\"blue\",\n", + " label=\"Scale factor\",\n", + " )\n", + " sf_cdf_art = axs[1].fill_between(\n", + " sf(linsp(nuisance)), effect.pdf(sf(linsp(nuisance))), color=\"b\", alpha=0.2\n", + " )\n", + " plt.draw()\n", + "\n", + "\n", + "axs[0].legend()\n", + "axs[1].legend()\n", + "axs[0].set_xlabel(r\"nuisance parameter ($\\theta$)\")\n", + "axs[0].set_ylim(0)\n", + "axs[0].set_ylabel(r\"$p(\\theta)$\")\n", + "axs[1].set_xlabel(r\"scale factor (SF)\")\n", + "axs[1].set_ylim(0)\n", + "axs[1].set_ylabel(r\"Effect(SF)\")\n", + "plt.tight_layout()\n", + "plt.show()" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "JAX", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.10.12" + }, + "orig_nbformat": 4 + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/pyproject.toml b/pyproject.toml index 5614043..676fbc4 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -86,6 +86,7 @@ disallow_untyped_defs = false disallow_incomplete_defs = false check_untyped_defs = true strict = false +ignore_missing_imports = true [tool.ruff] @@ -93,7 +94,6 @@ select = [ "E", "F", "W", # flake8 "B", # flake8-bugbear "I", # isort - "ARG", # flake8-unused-arguments "C4", # flake8-comprehensions "EM", # flake8-errmsg "ICN", # flake8-import-conventions diff --git a/src/dilax/optimizer.py b/src/dilax/optimizer.py index bfeb6f0..3f3836d 100644 --- a/src/dilax/optimizer.py +++ b/src/dilax/optimizer.py @@ -1,7 +1,7 @@ from __future__ import annotations from collections.abc import Hashable -from typing import Callable +from typing import Any, Callable import equinox as eqx import jax @@ -15,23 +15,10 @@ class JaxOptimizer(eqx.Module): Example: ``` - optimizer = JaxOptimizer.make( - name="ScipyMinimize", - settings={"method": "trust-constr"}, - ) - - # or - - optimizer = JaxOptimizer.make( - name="LBFGS", - settings={ - "maxiter": 30, - "tol": 1e-6, - "jit": True, - "unroll": True, - }, - ) + optimizer = JaxOptimizer.make(name="GradientDescent", settings={"maxiter": 5}) + # or, e.g.: optimizer = JaxOptimizer.make(name="LBFGS", settings={"maxiter": 10}) + optimizer.fit(fun=nll, init_values=init_values) ``` """ @@ -55,6 +42,39 @@ def settings(self) -> dict[str, Hashable]: def solver_instance(self, fun: Callable) -> jaxopt._src.base.Solver: return getattr(jaxopt, self.name)(fun=fun, **self.settings) - def fit(self, fun: Callable, init_values: dict[str, float]) -> jax.Array: + def fit( + self, fun: Callable, init_values: dict[str, jax.Array] + ) -> tuple[dict[str, jax.Array], Any]: values, state = self.solver_instance(fun=fun).run(init_values) return values, state + + +class Chain(eqx.Module): + """ + Chain multiple optimizers together. + They probably should have the `maxiter` setting set to a value, + in order to have a deterministic runtime behaviour. + + Example: + ``` + opt1 = JaxOptimizer.make(name="GradientDescent", settings={"maxiter": 5}) + opt2 = JaxOptimizer.make(name="LBFGS", settings={"maxiter": 10}) + + chain = Chain(opt1, opt2) + # first 5 steps are minimized with GradientDescent, then 10 steps with LBFGS + chain.fit(fun=nll, init_values=init_values) + ``` + """ + + optimizers: tuple[JaxOptimizer, ...] + + def __init__(self, *optimizers: JaxOptimizer) -> None: + self.optimizers = optimizers + + def fit( + self, fun: Callable, init_values: dict[str, jax.Array] + ) -> tuple[dict[str, jax.Array], Any]: + values = init_values + for optimizer in self.optimizers: + values, state = optimizer.fit(fun=fun, init_values=values) + return values, state diff --git a/src/dilax/parameter.py b/src/dilax/parameter.py index a78d6f6..9792948 100644 --- a/src/dilax/parameter.py +++ b/src/dilax/parameter.py @@ -12,7 +12,7 @@ class Parameter(eqx.Module): value: jax.Array = eqx.field(converter=as1darray) - bounds: tuple[jnp.array, jnp.array] = eqx.field( + bounds: tuple[jax.Array, jax.Array] = eqx.field( static=True, converter=lambda x: tuple(map(as1darray, x)) ) constraints: set[HashablePDF] = eqx.field(static=True) @@ -20,7 +20,7 @@ class Parameter(eqx.Module): def __init__( self, value: jax.Array, - bounds: tuple[jnp.array, jnp.array], + bounds: tuple[jax.Array, jax.Array] = (-jnp.inf, jnp.inf), ) -> None: self.value = value self.bounds = bounds @@ -48,9 +48,6 @@ def constraint(self) -> HashablePDF: def scale_factor(self, parameter: Parameter, sumw: jax.Array) -> jax.Array: ... - def __call__(self, parameter: Parameter, sumw: jax.Array) -> jax.Array: - return jnp.atleast_1d(self.scale_factor(parameter=parameter, sumw=sumw)) * sumw - class unconstrained(Effect): @property @@ -61,6 +58,9 @@ def scale_factor(self, parameter: Parameter, sumw: jax.Array) -> jax.Array: return parameter.value +DEFAULT_EFFECT = unconstrained() + + class gauss(Effect): width: jax.Array = eqx.field(static=True, converter=as1darray) @@ -72,7 +72,9 @@ def constraint(self) -> HashablePDF: return Gauss(mean=0.0, width=1.0) def scale_factor(self, parameter: Parameter, sumw: jax.Array) -> jax.Array: - return parameter.value * self.width + 1 + gx = Gauss(mean=1.0, width=self.width) + g1 = Gauss(mean=1.0, width=1.0) + return gx.inv_cdf(g1.cdf(parameter.value + 1)) class shape(Effect): @@ -88,8 +90,8 @@ def __init__( self.down = down # -1 sigma @eqx.filter_jit - def vshift(self, parameter: Parameter, sumw: jax.Array) -> jax.Array: - factor = parameter.value + def vshift(self, sf: jax.Array, sumw: jax.Array) -> jax.Array: + factor = sf dx_sum = self.up + self.down - 2 * sumw dx_diff = self.up - self.down @@ -112,10 +114,9 @@ def constraint(self) -> HashablePDF: return Gauss(mean=0.0, width=1.0) def scale_factor(self, parameter: Parameter, sumw: jax.Array) -> jax.Array: - return jax.numpy.clip( - (sumw + self.vshift(parameter=parameter, sumw=sumw)) / sumw, - a_min=1e-5, - ) + sf = parameter.value + 1 + # clip, no negative values are allowed + return jnp.maximum((sumw + self.vshift(sf=sf, sumw=sumw)) / sumw, 0.0) class lnN(Effect): @@ -141,7 +142,9 @@ def constraint(self) -> HashablePDF: def scale_factor(self, parameter: Parameter, sumw: jax.Array) -> jax.Array: width = self.scale(parameter=parameter) - return jnp.exp(parameter.value * width) + g1 = Gauss(mean=1.0, width=1.0) + gx = Gauss(mean=1.0, width=width) + return g1.inv_cdf(gx.cdf(jnp.exp(parameter.value))) class poisson(Effect): @@ -209,7 +212,7 @@ class modifier(ModifierBase): effect: Effect def __init__( - self, name: str, parameter: Parameter, effect: Effect = unconstrained() + self, name: str, parameter: Parameter, effect: Effect = DEFAULT_EFFECT ) -> None: self.name = name self.parameter = parameter @@ -255,10 +258,10 @@ class compose(ModifierBase): ``` """ - modifiers: tuple[modifier] + modifiers: tuple[modifier, ...] names: list[str] = eqx.field(static=True) - def __init__(self, *modifiers: tuple[modifier]) -> None: + def __init__(self, *modifiers: modifier) -> None: self.modifiers = modifiers # check for duplicate names @@ -267,22 +270,20 @@ def __init__(self, *modifiers: tuple[modifier]) -> None: msg = f"Modifier need to have unique names, got: {duplicates}" raise ValueError(msg) - @property - def names(self) -> list[str]: - names = [] + # set names + self.names = [] for m in range(self.n_modifiers): modifier = self.modifiers[m] if isinstance(modifier, compose): - names.extend(modifier.names) + self.names.extend(modifier.names) else: - names.append(modifier.name) - return list(names) + self.names.append(modifier.name) @property def n_modifiers(self) -> int: return len(self.modifiers) - def scale_factors(self, sumw: jax.Array) -> jax.Array: + def scale_factors(self, sumw: jax.Array) -> dict[str, jax.Array]: sfs = {} for m in range(self.n_modifiers): modifier = self.modifiers[m] @@ -294,9 +295,9 @@ def scale_factors(self, sumw: jax.Array) -> jax.Array: return sfs def scale_factor(self, sumw: jax.Array) -> jax.Array: - return jnp.atleast_1d( - jnp.prod(jnp.stack(list(self.scale_factors(sumw=sumw).values())), axis=0) - ) + sfs = jnp.stack(list(self.scale_factors(sumw=sumw).values())) + # calculate the product in log-space for numerical precision + return jnp.exp(jnp.sum(jnp.log(sfs), axis=0)) - def __call__(self, sumw: jax.Array) -> tuple[jax.Array, jax.Array]: + def __call__(self, sumw: jax.Array) -> jax.Array: return jnp.atleast_1d(self.scale_factor(sumw=sumw)) * sumw diff --git a/src/dilax/pdf.py b/src/dilax/pdf.py index c91f2e1..9771c69 100644 --- a/src/dilax/pdf.py +++ b/src/dilax/pdf.py @@ -48,10 +48,10 @@ def inv_cdf(self, x: jax.Array) -> jax.Array: class Gauss(HashablePDF): - mean: float = eqx.field(static=True) - width: float = eqx.field(static=True) + mean: float | jax.Array = eqx.field(static=True) + width: float | jax.Array = eqx.field(static=True) - def __init__(self, mean: float, width: float) -> None: + def __init__(self, mean: float | jax.Array, width: float | jax.Array) -> None: self.mean = mean self.width = width @@ -76,9 +76,9 @@ def inv_cdf(self, x: jax.Array) -> jax.Array: class Poisson(HashablePDF): - lamb: float = eqx.field(static=True) + lamb: float | jax.Array = eqx.field(static=True) - def __init__(self, lamb: float) -> None: + def __init__(self, lamb: float | jax.Array) -> None: self.lamb = lamb def __hash__(self): diff --git a/src/dilax/util.py b/src/dilax/util.py index f7819de..100382c 100644 --- a/src/dilax/util.py +++ b/src/dilax/util.py @@ -2,7 +2,7 @@ import collections import pprint -from collections.abc import Hashable, Mapping +from collections.abc import Hashable, Iterable, Mapping from typing import Any, Callable, TypeVar import jax @@ -27,8 +27,7 @@ def _pretty_key(key): key = FrozenDB.keyify(key) if len(key) == 1: return next(iter(key)) - else: - return tuple([_pretty_key(k) for k in key]) + return tuple([_pretty_key(k) for k in key]) def _indent(amount: int, s: str) -> str: @@ -45,8 +44,7 @@ def _pretty_dict(x): rep += f"{_pretty_key(key)!r}: {_pretty_dict(val)},\n" if rep: return "{\n" + _indent(2, rep) + "\n}" - else: - return "{}" + return "{}" K = TypeVar("K") @@ -66,7 +64,7 @@ def _prepare_freeze(xs: Any) -> Any: return {FrozenDB.keyify(key): _prepare_freeze(val) for key, val in xs.items()} -def _check_no_duplicate_keys(keys: tuple[Hashable, ...]) -> None: +def _check_no_duplicate_keys(keys: Iterable[Hashable]) -> None: keys = list(keys) if any(keys.count(x) > 1 for x in keys): msg = f"Duplicate keys: {tuple(keys)}, this is not allowed!" @@ -134,8 +132,7 @@ def only(self, *keys) -> FrozenDB: def subset(self, *keys) -> FrozenDB: new = {} for key in keys: - key = self.keyify(key) - new.update({k: v for k, v in self.items() if key <= k}) + new.update({k: v for k, v in self.items() if self.keyify(key) <= k}) return self.__class__(new) def copy(self): @@ -180,9 +177,16 @@ def as1darray(x: jax.Array) -> jax.Array: return jnp.atleast_1d(jnp.asarray(x)) -if __name__ == "__main__": - import jax.numpy as jnp +def dump_jaxpr(fun: Callable, *args, **kwargs) -> str: + jaxpr = jax.make_jaxpr(fun)(*args, **kwargs) + return jaxpr.pretty_print(name_stack=True) + + +def dump_hlo_graph(fun: Callable, *args, **kwargs) -> str: + return jax.xla_computation(fun)(*args, **kwargs).as_hlo_dot_graph() + +if __name__ == "__main__": hists = HistDB( { # QCD