From 121524ea57fe67f8d66939b721a9ec32702fc45a Mon Sep 17 00:00:00 2001 From: frostedoyster Date: Fri, 19 Jul 2024 20:05:06 +0200 Subject: [PATCH] Revert some changes --- llpr.ipynb | 379 ------------------ src/metatrain/experimental/soap_bpnn/model.py | 1 - tests/resources/llpr.py | 221 ---------- tests/resources/options.yaml | 28 +- tests/resources/split.py | 13 - 5 files changed, 6 insertions(+), 636 deletions(-) delete mode 100644 llpr.ipynb delete mode 100644 tests/resources/llpr.py delete mode 100644 tests/resources/split.py diff --git a/llpr.ipynb b/llpr.ipynb deleted file mode 100644 index 964c8182f..000000000 --- a/llpr.ipynb +++ /dev/null @@ -1,379 +0,0 @@ -{ - "cells": [ - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "\n", - "# Computing LLPR uncertainties\n", - "\n", - "This tutorial demonstrates how to use an already trained and exported model\n", - "from Python. It involves the computation of the local prediction rigidity\n", - "([LPR](LPR_)) for every atom of a single ethanol molecule, using the\n", - "last-layer prediction rigidity ([LLPR](LLPR_)) approximation.\n", - "\n", - "\n", - "The model was trained using the following training options.\n", - "\n", - ".. literalinclude:: options.yaml\n", - " :language: yaml\n", - "\n", - "You can train the same model yourself with\n", - "\n", - ".. literalinclude:: train.sh\n", - " :language: bash\n", - "\n", - "A detailed step-by-step introduction on how to train a model is provided in\n", - "the `label_basic_usage` tutorial.\n" - ] - }, - { - "cell_type": "code", - "execution_count": 1, - "metadata": { - "collapsed": false - }, - "outputs": [], - "source": [ - "import torch\n", - "from metatensor.torch.atomistic import load_atomistic_model" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Exported models can be loaded using the `load_atomistic_model` function from the\n", - "metatensor.torch.atomistic` module. The function requires the path to the exported\n", - "model and, for many models, also the path to the respective extensions directory.\n", - "Both are produced during the training process.\n", - "\n" - ] - }, - { - "cell_type": "code", - "execution_count": 2, - "metadata": { - "collapsed": false - }, - "outputs": [], - "source": [ - "model = load_atomistic_model(\"model.pt\", extensions_directory=\"extensions/\")\n", - "model = model.to(\"cpu\")" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "In metatrain, a Dataset is composed of a list of systems and a dictionary of targets.\n", - "The following lines illustrate how to read systems and targets from xyz files, and\n", - "how to create a Dataset object from them.\n", - "\n" - ] - }, - { - "cell_type": "code", - "execution_count": 3, - "metadata": { - "collapsed": false - }, - "outputs": [], - "source": [ - "from metatrain.utils.data import Dataset, read_systems, read_targets # noqa: E402\n", - "from metatrain.utils.neighbor_lists import get_system_with_neighbor_lists # noqa: E402\n", - "import omegaconf\n", - "\n", - "\n", - "qm9_systems = read_systems(\"ethanol_reduced_100.xyz\", dtype=torch.float64)\n", - "\n", - "target_config = {\n", - " \"energy\": {\n", - " \"quantity\": \"energy\",\n", - " \"read_from\": \"ethanol_reduced_100.xyz\",\n", - " \"file_format\": \".xyz\",\n", - " \"key\": \"energy\",\n", - " \"unit\": \"kcal/mol\",\n", - " \"forces\": {\n", - " \"read_from\": \"ethanol_reduced_100.xyz\",\n", - " \"key\": \"forces\",\n", - " \"file_format\": \".xyz\",\n", - " },\n", - " \"stress\": False,\n", - " \"virial\": False,\n", - " },\n", - "}\n", - "targets, target_info = read_targets(target_config, dtype=torch.float64)\n", - "\n", - "requested_neighbor_lists = model.requested_neighbor_lists()\n", - "qm9_systems = [\n", - " get_system_with_neighbor_lists(system, requested_neighbor_lists)\n", - " for system in qm9_systems\n", - "]\n", - "dataset = Dataset({\"system\": qm9_systems, **targets})\n", - "\n", - "# We also load a single ethanol molecule on which we will compute properties.\n", - "# This system is loaded without targets, as we are only interested in the LPR\n", - "# values.\n", - "ethanol_system = read_systems(\"ethanol_reduced_100.xyz\", dtype=torch.float64)[0]\n", - "ethanol_system = get_system_with_neighbor_lists(\n", - " ethanol_system, requested_neighbor_lists\n", - ")" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "The dataset is fully compatible with torch. For example, be used to create\n", - "a DataLoader object.\n", - "\n" - ] - }, - { - "cell_type": "code", - "execution_count": 4, - "metadata": { - "collapsed": false - }, - "outputs": [], - "source": [ - "from metatrain.utils.data import collate_fn # noqa: E402\n", - "\n", - "\n", - "dataloader = torch.utils.data.DataLoader(\n", - " dataset,\n", - " batch_size=10,\n", - " shuffle=False,\n", - " collate_fn=collate_fn,\n", - ")" - ] - }, - { - "cell_type": "code", - "execution_count": 5, - "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "4\n" - ] - } - ], - "source": [ - "ll_params = []\n", - "for name, param in model.named_parameters():\n", - " if \"last_layers\" in name and \"weight\" in name:\n", - " ll_params.append(param)\n", - "print(len(ll_params))" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "We now wrap the model in a LLPRUncertaintyModel object, which will allows us\n", - "to compute prediction rigidity metrics, which are useful for uncertainty\n", - "quantification and model introspection.\n", - "\n" - ] - }, - { - "cell_type": "code", - "execution_count": 9, - "metadata": { - "collapsed": false - }, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "tensor([[ 6.1554e+10, -1.9105e+11, 1.6567e+11, ..., -3.8915e+10,\n", - " -4.0176e+10, 7.7346e+10],\n", - " [-1.9105e+11, 5.9509e+11, -5.1658e+11, ..., 1.2133e+11,\n", - " 1.2518e+11, -2.4106e+11],\n", - " [ 1.6567e+11, -5.1658e+11, 4.4873e+11, ..., -1.0540e+11,\n", - " -1.0869e+11, 2.0935e+11],\n", - " ...,\n", - " [-3.8915e+10, 1.2133e+11, -1.0540e+11, ..., 2.4764e+10,\n", - " 2.5530e+10, -4.9173e+10],\n", - " [-4.0176e+10, 1.2518e+11, -1.0869e+11, ..., 2.5530e+10,\n", - " 2.6336e+10, -5.0716e+10],\n", - " [ 7.7346e+10, -2.4106e+11, 2.0935e+11, ..., -4.9173e+10,\n", - " -5.0716e+10, 9.7679e+10]], dtype=torch.float64,\n", - " grad_fn=)\n" - ] - } - ], - "source": [ - "from metatensor.torch.atomistic import ( # noqa: E402\n", - " MetatensorAtomisticModel,\n", - " ModelMetadata,\n", - ")\n", - "\n", - "from metatrain.utils.llpr import LLPRUncertaintyModel # noqa: E402\n", - "\n", - "\n", - "llpr_model = LLPRUncertaintyModel(model)\n", - "llpr_model.eval()\n", - "# llpr_model.compute_covariance(dataloader)\n", - "\n", - "from metatrain.utils.loss import TensorMapDictLoss\n", - "loss_fn = TensorMapDictLoss(\n", - " weights={\n", - " \"energy\": 1.0,\n", - " \"energy_positions_gradients\": 0.0,\n", - " },\n", - " reduction=\"sum\"\n", - ")\n", - "\n", - "llpr_model.compute_covariance_as_pseudo_hessian(dataloader, target_info, loss_fn, ll_params)\n", - "print(llpr_model.covariance)\n", - "\n", - "llpr_model.compute_inverse_covariance()\n", - "\n", - "# calibrate on the same dataset for simplicity. In reality, a separate\n", - "# calibration/validation dataset should be used.\n", - "llpr_model.calibrate(dataloader)\n", - "\n", - "exported_model = MetatensorAtomisticModel(\n", - " llpr_model.eval(),\n", - " ModelMetadata(),\n", - " llpr_model.capabilities,\n", - ")" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "We can now use the model to compute the LPR for every atom in the ethanol molecule.\n", - "To do so, we create a ModelEvaluationOptions object, which is used to request\n", - "specific outputs from the model. In this case, we request the uncertainty in the\n", - "atomic energy predictions.\n", - "\n" - ] - }, - { - "cell_type": "code", - "execution_count": 26, - "metadata": { - "collapsed": false - }, - "outputs": [], - "source": [ - "from metatensor.torch.atomistic import ModelEvaluationOptions, ModelOutput # noqa: E402\n", - "\n", - "\n", - "evaluation_options = ModelEvaluationOptions(\n", - " length_unit=\"angstrom\",\n", - " outputs={\n", - " # request the uncertainty in the atomic energy predictions\n", - " \"mtt::aux::energy_uncertainty\": ModelOutput(per_atom=True),\n", - " # `per_atom=False` would return the total uncertainty for the system,\n", - " # or (the inverse of) the TPR (total prediction rigidity)\n", - " # you also can request other outputs from the model here, for example:\n", - " # \"energy\": ModelOutput(per_atom=True),\n", - " # \"mtt::aux::last_layer_features\": ModelOutput(per_atom=True),\n", - " },\n", - " selected_atoms=None,\n", - ")\n", - "\n", - "outputs = exported_model([ethanol_system], evaluation_options, check_consistency=True)\n", - "lpr = outputs[\"mtt::aux::energy_uncertainty\"].block().values.detach().cpu().numpy()" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "We can now visualize the LPR values using the `plot_atoms` function from\n", - "``ase.visualize.plot``.\n", - "\n" - ] - }, - { - "cell_type": "code", - "execution_count": 27, - "metadata": { - "collapsed": false - }, - "outputs": [ - { - "data": { - "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZQAAAGdCAYAAAA1yoVoAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy80BEi2AAAACXBIWXMAAA9hAAAPYQGoP6dpAAB1CklEQVR4nO3dd3QUVf/H8fdsyaZX0hu9l0BCCRY6iEixII+NIqJUQazYUFDRRxELXaogD9goooBIValJCBB6TyCdhPS2O/P7A8nPKAGS3WR2k/s6Zw4nm92dz5Jkv3vL3CspiqIgCIIgCGbSqB1AEARBqBlEQREEQRAsQhQUQRAEwSJEQREEQRAsQhQUQRAEwSJEQREEQRAsQhQUQRAEwSJEQREEQRAsQqd2gDshyzKJiYm4uLggSZLacQRBqIUURSEnJ4eAgAA0mqr7LF5YWEhxcbHZz2NnZ4e9vb0FEt05mygoiYmJBAcHqx1DEASBhIQEgoKCquS5CwsLqRfqTHKqyezn8vPz48KFC9VaVGyioLi4uADXf5Curq4qpxEEoTbKzs4mODi49P2oKhQXF5OcauJCdCiuLpVvBWXnyNQLv0RxcbEoKP90o5vL1dVVFBRBEFRVHd3uri4aswqKWmyioAiCINQmJkXGZMayvSZFtlyYChAFRRAEwcrIKMhUvqKY81hziIIiCIJgZWRkzGljmPfoyrO9TjpBEATBKokWiiAIFZaVlUVycjIFBQUUFxeXXvPg6+uLh4eH2vFsnklRMJmx96E5jzWHKCiCINxSXl4eBw4cIDo6mujoaPYfOMiF8+fKvX9wSCgd2renffsIwsPD6dixY5VOta2JxBiKIAg1yokTJ5g3bx5Lly4jNzcHnd6As0cg9q4BNGzfHoOTBxqtHkmjQ5FNyLKRovxM8jIvs+2PQ2z4aSMlxYXYOzgw9KmnGDt2LG3atFH7ZQlVSBQUQRBKmUwm1q5dy5ez57B7107sHV3xCulAg5AwHFx9kaTbDbvWwzukHQCKIlOYk0765VhWfPMtCxcupGOnTkwYP55HH30UvV5f9S/IRskomGywhSIpikqdbRWQnZ2Nm5sbWVlZ4sJGQagip06dYvjw4ezbtw93n/p414vEM7AVGq35nztl2URm4nHSLuwlM/k0rdu0YcXXX9O6dWsLJK8e1fE+dOMc50764WLGhY05OTINmiZX+3ummOUlCLWcyWRi5syZtG7ThqPHz9Oi6xiadRlLnZC2FikmABqNFq+gVjS951la9Xie8/GphIdHMH36dEpKSixyDkF9oqAIQi127tw57r77bl5++WW8QjrQosckXL0bVOk5nT1DaNFtIr4N72XqO+8Q0b49x48fr9Jz2pobs7zMOdQgxlAE4Q7l5+eTmJhIUlJS6b9JSUnk5eVhNBoxmUxotVp0Oh0ODg74+/vj7+9PQEBA6b/Ozs5qv4xShw8fpkfPXhQUKzTvMgZX7/rVdm6NVkdIq754BrbkfPS3REZ2ZvPmTURGRlZbBmsm/3WY83g1iIIiCDeRmZlJTExM6VTZg1H7uXD+Upn7GOx1ePjaY3DUoNVKSFpQZDAZFYoLZDJTCynMN5Z5THBIIO0jOhAefn1KbXh4OHXq1KnOlwbAgQMH6NWrN+hdad51JHqDOoXO2TOYZl3HcmbPMrr36MEvP/9Mt27dVMkimE8UFEHg+iZuMTExbNiwgXXr13L0SBwADk56Qpo70eheB3o82wBPfwPuPna4++hxcNbeduXZglwTWWnFZKYWk5lcTPzJPE7E7WDz1p/Jz7m+iVLT5k0YNOBBBgwYQIcOHdBqtVX6WuPi4ujduw+SwZPGd41Ep6/eTZj+Sad3oPFdIzmzbzn9+j3Ajh3b6dixo6qZ1GYyc5aXOY81h5jlJdRaiqKwY8cOvv32Wzb8tI6kxBScXO1o3dWVlne7U6+VM3517dFoLb9cuSwrpMYXcjEuj2N/XuPIzhyyrhZSx9uT/g8MZPDgwfTu3dvixSUpKYk2YWEUGu1oes9z6OwcLPr85jAZizj1x2IoziA6OooGDap2LKeiqnOW15HjPmbP8mrdPLXa3zNFQRFqnWvXrrF8+XLmzJvNmVNn8Q1xIqy7K2E9PGjUzgWdvvrnqsgmhXOHczm0LYPD27O5ci6XuvVCGDN6HE8//bRFusUURaF//wFs27GbFj1ewM7e+q5eNxYXcGz7Z7Rr05xdu3ZW6Va7FVWdBSXWAgUlTBSUmxMFRbCEY8eOMWvWLFb97xuKi4uJ6O1J9yd8aRzhUi2bJt0pRVE4fySX7d+kcGBTBhq0DB48mEmTXiA8PLzSz7ty5UqeeuopmnQejmdgSwsmtqys1LMc3zWfL774ggkTJqgdp5QoKLcnCopQ4126dIm3336LFStW4ulnT5chdbh3sA/u3nZqR7utnIwSfv8xlV2rr5ISn8cjgx/h/ffep3HjxhV6nqSkJJo2a47BrR4NOz5RRWkt50LMj1xLPETc0aNW0/VVnQUl5rgvzmYUlNwcmXbNU8SFjYJgKWlpabzwwgs0btyI9b98xxNv1eWjra0ZMDbIJooJgIunnvufCWTGllaMnNGAnX/+TPPmzXj22WdJTEy84+d57rnRFJfIhIYNqrqwFhTSuh8avRMjRjyNDXzmtThZMf9QgygoQo0jyzJffPEF9RvUY8GiOTww1o+Ptram55N+6Oxs81deo5W452EfPtjSmsEvB7P6u69p0LA+H3zwAUaj8ZaPPXz4MD/9tIGglg+gNzhVU2LzaHUGQloP5Pffd/P777+rHUe4Q7b51yUI5Th79ixdut7LxIkT6dDfmf/+1oYBY4Owd6raqbjVxc6g4b6nA/hoW2u6Pu7JW2+9ScdO7YmLiyv3MfPmzcPByQ2vYNta6dfdrynO7r7MmTNH7SjVzoRk9qEGUVCEGuFGq6RV65acvniIV1c0Z+g79XHxrJkr2jq66PjPq3V5c00LUrLOEh7e7qatlezsbL7+egVeoR3QaGyrqEqSRJ26nfjhhx9JTk5WO061EgVFEFSSlpZGj57dmThxInc97MG0DS1p1tFN7VjVon4bF6aubUGv4d689dabdO7cicuXL5d+f8WKFRQWFeJTv5OKKSvPOzQCSaNh8eLFakcR7oAoKIJNO3LkCBHt23Ho6H5eWd6cp96uV2O6t+6UnUHD4JdCeWNNCy5cOU54RFv27t0LwPz5C/AMaInBwTYLrM7OAc+gtsybN79WDc7LimT2oQZRUASb9eOPPxLZuROSYxZvf9+c5pG2+aZpKQ3auPDW981xDyqha9cuzJs3j7i4o7j7N1M7mlk8/Jtx5cplrly5onaUaiO6vAShmiiKwvTp03n44Ydpca8Tr/+vGV4BBrVjWQW3Ona8vLwpnQZ4MHbsWACc3ANUTmUeJ48gAKKjo1VOItyOKCiCTVEUheeff563336bBycGM/azhhgca1cX1+3o7TSMeL8+T7xZF0mCxFM7UWST2rEqzc7BDXtHV6KiotSOUm1MaMw+1CBWGxZshizLjBkzhoULFzJsWn26/cdX7UhWS5Ikeg31x9VLz/wXY1EUE406PIFkYzO94PprcXQPrFUFRTFzHEQRYyiCUD5FURg/fjxfffUVIz9sIIrJHerYrw7jv2hEZuJRzh78H4qi1tZL5nFwCyQ6OkbtGNWmVo6hfPjhh0iSxKRJk255v++++46mTZtib29Pq1at+OWXX8w5rVALvfLKK8ybN4/h79Xjnod81I5jU8J7ezFmViOuJhzmXPR3NjlbSm9wJjs7S+0Ywm1UuqAcPHiQBQsW0Lp161veb8+ePTz22GOMHDmSQ4cOMWjQIAYNGnTLK3sF4e8WLlzIJ598whNv1qXLYNEyqYz293kx6r8NSLtwkMRTO9SOU2EarZ6ioiKbLIaVYVI0Zh9qqNRZc3NzeeKJJ/jqq6/w8PC45X0///xz7rvvPl5++WWaNWvG9OnTadeuHbNnz65UYKF22b17N+PGjaXHE370Guqvdhyb1nmgNwPGBRIf9wuZicfVjlMh0l/7otxu3bKaQkZCRmPGYUNdXuPGjaNfv3707Nnztvfdu3fvv+7Xp0+f0guvbqaoqIjs7Owyh1D7XLp0iYceHkTjCBceez1U7Tg1wqAJwbTt7snZgyvJz7ad5UxkUwlarRa9vmYupVNTVLigrF69mpiYGGbMmHFH909OTsbXt2w3ha+v7y3X5pkxYwZubm6lR3BwcEVjCjYuLy+P/v37oXUoYuznDVXZRbEm0mgknvukIT4hek7vXUJJcb7ake6IbCrBYLBXO0a1qRWD8gkJCUycOJFvvvkGe/uq++FOmTKFrKys0iMhIaHKziVYp5deeonTZ08xYW5DnD3Ep1JLsnfS8sKCxkjkcCHmB7Xj3JGCnDRC69ZVO0a1qRVjKNHR0aSmptKuXTt0Oh06nY5du3bxxRdfoNPpMJn+ffGUn58fKSkpZW5LSUnBz8+v3PMYDAZcXV3LHELtsX37dubPn8/gl4MJbmIb+3fYGu9ge4a+U5erCYe5evmI2nFuqzA7kY4d2qsdQ7iNChWUHj16cPToUWJjY0uPiIgInnjiCWJjY9Fq/33RVGRkJNu2bStz29atW4mMjDQvuVAj5eTkMOLpYTTr6E73x8WMrqrUsZ8X7Xp6cjH2e0qKctWOUy7ZZCQ3M5Hw8HC1o1Sb64Py5h1qqNCV8i4uLrRs2bLMbU5OTnh5eZXePnToUAIDA0vHWCZOnEiXLl2YOXMm/fr1Y/Xq1URFRbFw4UILvQShJnnllVdISU1m+uJWaDTq/FHUFpIkMfTderze9wgXDq2lcaen1I50U/nZycgmIxEREWpHqTaymcunyKgzvdriHW3x8fEkJSWVft25c2dWrVrFwoULadOmDd9//z3r1q37V2EShN9///16V9dLQfiE1J4BWDW5e9uVdn1lJB5TO85N5V69hFarpU0b29pxsjYyey2vnTt33vJrgMGDBzN48GBzTyXUYIqi8PLLL1KvpSvdnyh/fE2wvI79vNj9XRoXj23Ew6+p1a33dTUhmp49e+Hg4KB2lGpj7sC6SaULQMVcTMEqrF+/nv37D/LIi4Giq6uaSZLEoy8Hk38tjbRL1rVEfG5GPNnp8YwfP07tKNXKvIsarx9qEAVFUJ3RaGTK66/SorMHLe5yVztOrVS3pTPt7/PiyonNmEwlascplXxuD0FBwfTt21ftKNXKpEhmH2oQBUVQ3YoVKzh54jSPvBikdpRa7eHJwRQV5JBy9k+1owBQUpRHxuUjjBs39qYzSAXrIwqKoCqTycS706bS/j4v6rVyVjtOreZX14Eug31IPL3NKlopV05uR6uBp59+Wu0o1c5WN9gSBUVQ1aZNm7h0MYG+I8XCj9bgvqf9KSksICPhsKo5ctIvknxmN++99x4+PrVvuwJZ0Zh9qEEUFEFVc+fOoV5LV+q1Fq0Ta+BXz4Hmke6kXFCv28tkKuFCzHeER7Rn8uTJquUQKk4UFEE158+fZ/PmLXR73BtJEjO7rEWPJ33JSU8gN/OyKudPiNtESUEmK75eXmvHTkSXlyBU0IIFC3B00dOxn5faUYS/CevmgbuPgZRze6r93BlX4kg+8zvvvfceTZs2rfbzWwsZ82Z6qbXRsygogipMJhNLly2m84OeGBxq56dQa6XVSXR7zIerCTGYjMXVdt6s1DOc3b+Shx56SHR12ShRUARV7N+/n7TUq3ToW0ftKMJNtL/PC5PRSFbqmWo537WU05zes4wePXqwatWqWtvVdYO4sFEQKmDDhg24ednToI0YjLdG/vXt8Q52qJatgtMTYjn15xJ69erBunVrsbOzq/JzWrtasR+KIFjKuvVrad3VBY1WDMZbI0mSCO/lTlbKMRSlanrkS4rzOXvgf5zZt5JHHx3M+nXratV6XTWRKChCtTtz5gynTp4mrLun2lGEW2jbw5Oi/Nwqme2VkXiMuK0zKcw4w7Jly1j1zTdiv/i/qRX7oQiCJWzatAm9QUvLu9zUjiLcQsO2Lji66LmWdAIXzxCLPGfetSskntpBenws9/Xty6KvviIwMNAiz12TmL/asDptBVFQhGp34MAB6jZzxuBYuwderZ1WJ9GgrRMJpy+Z9TyyqYSrl4+QdmEfWWkX8PPzZ+nSpQwbNkxcf1QOc68lUes6FFFQhGp3MGo/IRFiAy1bUK+VE8f3nOHMgf/h4hmCk0cQju4BaLXld0/JJiP5WYnkZl4mL/MyWcknKCrIoXv3HowfP5P+/fuj04m3nppI/FSFapWTk8OZ0+e4Z3h9taMId6BuC2dMRgVvxwLOHN2I0ViCpNHg4u6P1t4djUaHpNWhyCZkUwmmomxyryUhm4xotVqaNmvOE6Of4bnnnqNJkyZqvxybISsSshlL0JvzWHOIgiJUq9jYWBRFoW4LJ7WjCHegbsvrP6cPPnifvn37EhcXR1RUFNHR0SQlJZGXl09hYSH29vY4ODjg5+dLeHg44eHhtG7dWszaqiTz95QXXV5CLRATE4OdQUtAQ0e1owh3wMPXDndve6Kjoxk0aFBpsRCEmxEFRahW58+fxyfEEa1ODMbaAkmS8K1rz4ULF9SOUquYuwS9WsvXi4IiVKukpCTcvMXsLlviVkdLYtIVtWPUKiYkTGZcS2LOY80hLmwUqtWVKwm4eYvPMbbE3deOxERRUITbEwVFqFaJSYm4+4i1mmyJu7eepKQktWPUKmLHRkG4DUVRSE5KEQXFxrh525GdlUtBQYHaUWoNE//f7VW5Qx2ioAjVpqCggMLCIpw9RJeXLXHxvP7zSk9PVzmJYO3EX7ZQbUpKSgDQ6cXnGFui013/eRmNRpWT1B5ilpcg3MaNN6RavneSzdH89S4hCkr1EYtDCsJtKIqidgTBDOLnV30UM5egV8S0YaGmu7HfhUmtEUOhUuS/GiZivxLhdkQLRag2N96QjCVVswOgUDWMxus/L7FCcPURXV6CcBsODg7Y2xvIuyb64m1Jbub1n5eXl5fKSWoPW11tuEJlbN68ebRu3RpXV1dcXV2JjIxk06ZN5d5/2bJlSJJU5rC3F/tg1FaSJOHr50NmSrHaUYQKuJZajKubM46OYkFP4dYq1EIJCgriww8/pFGjRiiKwvLlyxk4cCCHDh2iRYsWN32Mq6srp06dKv1a7NBWuwUEBJKVdk7tGEIFXEsrwc/PT+0YtUqt2LGxf//+Zb5+//33mTdvHvv27Su3oEiSJH4ZhVKBAUGcSDypdgyhAq6lFhMQ0FjtGLVKrejy+juTycTq1avJy8sjMjKy3Pvl5uYSGhpKcHAwAwcO5NixY7d97qKiIrKzs8scQs0QEBBAdroYlLclWWkmAgOC1I4h2IAKF5SjR4/i7OyMwWBg9OjRrF27lubNm9/0vk2aNGHJkiWsX7+elStXIssynTt35vLly7c8x4wZM3Bzcys9goODKxpTsFL16tUj5VIeJqO4psEWKIpC6sUi6tatq3aUWkVGY/ahhgqftUmTJsTGxrJ//37GjBnDsGHDOH78+E3vGxkZydChQwkLC6NLly78+OOPeHt7s2DBglueY8qUKWRlZZUeCQkJFY0pWKl27dpRXGQi8Vy+2lGEO3AtpZjMtAKxS2M1MymS2YcaKjxt2M7OjoYNGwIQHh7OwYMH+fzzz29bJOD6dQht27bl7Nmzt7yfwWDAYDBUNJpgA9q2bYskSVw8lkdwE7GvvLW7eCwPgIiICJWTCLbA7HaRLMsUFRXd0X1NJhNHjx7F39/f3NMKNsrFxYWGjepzMS5P7SjCHbgYl0cdb0+CgsQYSnW6MShvzqGGCrVQpkyZQt++fQkJCSEnJ4dVq1axc+dOtmzZAsDQoUMJDAxkxowZAEybNo1OnTrRsGFDrl27xscff8ylS5d45plnLP9KBJvRPqIjB4/9onYM4Q5cPJZPeHh7Md2/milmrjas2MKV8qmpqQwdOvT6vuBubrRu3ZotW7bQq1cvAOLj49Fo/v+FZGZmMmrUKJKTk/Hw8CA8PJw9e/aUO4gv1A4dOnTgux/WUFRgwuAglh62VrJJ4cLhPAaM76B2lFrHVveUr1BBWbx48S2/v3PnzjJfz5o1i1mzZlU4lFCz3X///UyaNInje7Jo28NT7ThCOc4eyiE7s4j7779f7SiCjRCrDQvVrlGjRjRu0pBD2zLVjiLcwqHtmXj7eNGhg2ihVDdZMXccRZ3coqAIqhg08CEO78xGVus3X7itw9uzGdB/UJlubKF63Nix0ZxDDeI3RVDFgAEDyEov5PyRXLWjCDeRdL6AxPO5DBgwQO0ogg0RBUVQRadOnajj7cnBX66qHUW4iYNbrmLvYKBnz55qR6mV5L92bDTnUIMoKIIqtFotw4c9zZ9rr1JcKLZwtCaySWH3mqv8Z8hjYsl6ldjqlfKioAiqGT16NLlZxRwQrRSrcnhXJumJ+YwbN07tKIKNEQVFUE2DBg3o06c3O/6XrnYU4W92rEojIqKdWG5FRWJQXhAqYezYcZw7nMWFo2Jw3hqkxhdyZHcG48ZNUDtKrSZj5tIrYgxFqI369etHSGgQmxcnqR1FADYvScTD040hQ4aoHUWwQaKgCKrSarW89eZU9v+SzsVjopWiptT4QnZ9m8Zrr76Og4OD2nFqNcXMGV6KaKEItdXw4cNp3KQhP8y89cZrQtX68bPL+Ph4M2GC6O5Sm62uNiwKiqA6nU7HB+9/yNE/Mjm+N0vtOLXSpeN57NuYxrvvTBetEysgBuUFwQwPPfQQERHt+H7mZbEciwq+n5lAo8YNGDFihNpRBBsmCopgFSRJ4pNPPuX8kWx2rklRO06tcuCXdI7+nslHH36MTlfhTVyFKiC6vATBTF26dGHUqFF899/LpF0uVDtOrZB9tYSV0+J56OGHGDRokNpxhL+IpVcEwQI++eQT6nj5sPSNCyiK6PqqaivevYhe48i8ufPEroyC2URBEayKq6srixct5fjea+xYLbq+qtKBX9I5uDmduXPm4+Pjo3Yc4W9El5cgWEivXr0YNWoU336UwJWz+WrHqZHSrxSx4t3rXV2DBw9WO47wD6KgCIIFzZw5kwb1GvHFmLPkXitRO06NUpRv4osxZ/B082HB/AWiq0uwGFFQBKvk4uLChg0bKc7RMX/SOUxGMZ5iCbKssOjVc6QnGPlpw8/UqVNH7UjCTYgWiiBYWP369fnh+7Wc2J/F6o8uqh2nRvhp7mUObrnKNytX0apVK7XjCOUQBUUQqkC3bt344osv2bo8me2rktWOY9P2bUxn7ReXmT59upgiLFQJcRWTYPXGjBnDqVOn+OKdL7Cz13D3Q2JGUkVFb81g4ctnefKpJ3njjTfUjiPchgJmXUuiVgexKCiC1ZMkiVmzZpGfn8+S1xej0Up0HuitdiybcWh7BvMmneGhhx5m6ZKlYhDeBpjbbaVWl5coKIJN0Gg0LFiwAKPRyFevLKekWKbLYF+1Y1m9g5uvMn/yWfr3H8Cqb1aJpVVshCgoglDFNBoNixcvxt7envlvzCcnw0i/ZwPEJ+5ybF+VzDfTL/LI4EdZ8fUK9Hq92pGEGk4UFMGmaDQa5s6di7e3N9OnT+fK6QJGvF8PO3ut2tGshrFY5pv3L7LjfylMmDCBWbNmodWK/x9bIlooglBNJEli2rRptGzZkmHDh5J84SQT5jbE08+gdjTVZWeUMPf5s5w7lMtXX33FM888o3YkoRJstaCIacOCzXr00UfZ8+deSrKcmf7wcU4dzFY7kqouxuXy3iPHSb+gYfv2HaKYCNVOFBTBprVt25aogzG0aNqOD588zqoPLlJUYFI7VrUyFsv88Fk80wbHEejTiOioQ9x9991qxxLMoCiS2YcaKlRQ5s2bR+vWrXF1dcXV1ZXIyEg2bdp0y8d89913NG3aFHt7e1q1asUvv/xiVmBB+CdfX1927tjNJ598wq7VV3ln4DFOR9WO1srFY7lMe/g4mxYm8/ZbU9m/7yAhISFqxxLMVCv2QwkKCuLDDz8kOjqaqKgounfvzsCBAzl27NhN779nzx4ee+wxRo4cyaFDhxg0aBCDBg0iLi7OIuEF4QatVsvkyZM5HHuEuv6tmPHEcb557wJ5WUa1o1WJglwT338az7RH4vBwqEtUVDRTp04VM7kEVUmKmbsYeXp68vHHHzNy5Mh/fW/IkCHk5eWxcePG0ts6depEWFgY8+fPv+NzZGdn4+bmRlZWFq6urubEFWoBk8nE559/zltvv4mkNdF3lB+9hvphcLD9mU4lxTI7V6ewcV4yhbkyU6a8zuuvvy4KSTWojvehG+fouO55dE6Vn2RizCti/6Avqv09s9JjKCaTidWrV5OXl0dkZORN77N371569uxZ5rY+ffqwd+/eyp5WEG7rRmvl/LkLjBj6LOu/SGRK76PsXJ2CsURWO16lyCaFPevTeOO+o/zvg3geGvAYZ86cFa2SGspWx1AqPG346NGjREZGUlhYiLOzM2vXrqV58+Y3vW9ycjK+vmWvZvb19SU5+daL/BUVFVFUVFT6dXZ27egPFyzL19eXL7/8khdeeIE333qTZW//j58XJNPlP3W49xEfXL2s/404L8vIHz+msvN/6SRdzKNX7548/OYj6HQ6li9fTlJSEomJiSQlJZCWlkpxcQlGoxGTyYRWq0Wn06HX66jj5Y1/QDABAQH4+/uX/tukSRMaNGggLg4VLKLCBaVJkybExsaSlZXF999/z7Bhw9i1a1e5RaUyZsyYwbvvvmux5xNqt/r167Pqm1W8+sqrfPbZZ/xvzirWf3mFiPs86f64Lw3bOlvdG+qFozn8siiJQ79lIpsUfH19cHOT2Prrb2z99TcAvOsYCPDV4ecDzepp6dJei8FOQqcFrVaDySRjNBVTXFxEWkY2SSmn2fcHJKUaSU0r4kZnt5ubM+3atSMioiPh4eGEh4eLIqMyW70OxewxlJ49e9KgQQMWLFjwr++FhIQwefJkJk2aVHrb1KlTWbduHYcPHy73OW/WQgkODhZjKIJFXL16lWXLljF33mzOn7uIfz1nwnq40ra7Bw3buqDRVv8foywrnI3JZtuqFI79nk1hXgnGv+YTBPobaNdaT3hrO9q1tqdVUzv8fHTY2VU+p9GokJxq5NipYmKOFBF9pIjoI0biLxcCEBDgQ//+D9K/f3+6d++Og4ODJV6mTavOMZTwH14wewwl+uFZ1f6eafaV8rIsl3nz/7vIyEi2bdtWpqBs3bq13DGXGwwGAwaDuOpZqBpeXl68+OKLvPDCC2zdupVvv/2WnzasZ9OiY7h6GmjVxZVW97hTr5UTPiH2VfJJXVEU0q8UcfpgNr//mMb52DxMRhMmEwQF6Bn0iBu9ujgS0cYePx/LL2ih00kEBegJCtDTp5tT6e1p6UaijxTx2+58Nmz5mgULFuDoaE/v3r0ZMGAQAwcOxNPT0+J5hLIUM1soao2hVKiFMmXKFPr27UtISAg5OTmsWrWKjz76iC1bttCrVy+GDh1KYGAgM2bMAK5PG+7SpQsffvgh/fr1Y/Xq1XzwwQfExMTQsmXLOw4pZnkJVU2WZQ4cOMCGDRtYt/5HThw/BYCTqx2hLZwIbeFASDNHPP0MuPvocfe2w+B4+1ljxYUmrqWWkJVWTEZyMfEn8rh0rIDzR3LJzylBowFZhvA29gzs40T/Pk60amZnFd1NiqJw8kwJP/2ay0+/FrI3Kg+DwY7//Odxxo0bR0REhNoRq1V1tlDafT8ZrRktFFNeETGPfGrdLZTU1FSGDh1KUlISbm5utG7durSYAMTHx6PR/P/Esc6dO7Nq1SrefPNNXn/9dRo1asS6desqVEwEoTpoNBo6depEp06d+OCDD0hLSyMmJoaoqCiioqOI2nKATYvOlnmMo4sd7t4G7B01aHSg0UgoCphKFIoLZa6lFpGbXVzmMW7uruh1WvJzSvDzsePZp5wZ+bgrQQHWN0FAkiSaNbajWWNPXhkPyalGlq/JZsGK1Sxbtoz27dsxduwEhgwZIrrELEwBzBmMUGuDLbPHUKqDaKEI1iArK+uvGVVJf5tdlUR+fj5GoxGj0YhWq0Wv1+Pg4ICfnx8BAQF4eXmxe/duliz5irS0DLrd5czo4c4M7OOMXq9+S6SiTCaFX7blMW9ZDlt25OLp6cZLL73KxIkTcXR0VDtelanOFkqb719E62hGCyW/iMOPzLTuFoog1GZubm64ubnRrFmzO7q/yWRi5cqVjBnzLJcvJzL8Py5Mfi6UZo3tqjhp1dJqJfr3dqZ/b2fOXSzm86+uMXXqm3z55We8/fa7jBw5UlwbU0uJxSEFwcIURWHDhg20adOS4cOHE9Eqm6M7Q/hqpq/NF5N/alDXji/e9+H47yF071zI2LFjaNGiKd9++y2ybJsXkVoDW72wURQUQbCgixcv0qtXDwYOHIiP+2X2/hLMd4v8aNqoZhWSf6ofqufr2b5Ebw2hYUgqQ4YM4e67Izl58qTa0WzSjetQzDnUIAqKIFiALMvMmzePli2bc+bUXjauDGDrd350aGuvdrRq1aaFgY0r/dn2QyDpKUcJC2vDJ598gslUu7YUqK1EQREEM91olYwdO5YnHrLjyI5A+vZwsoqpv2rp2tmRmK0BjB3uyCuvvCxaKxWkKOYfahAFRRDMsGjRIlq2bM7Z0/vYsiaQef/1wcVZ/FkBODpq+OQdb3atC+JqahxhYW2YNWsWNjCxVHViDEUQapHi4mJGjx7NqFGj+M9APYe3B9Dz3po7ZdYcd3VwIGZrAKOHOjJ58mSefPIJCgoK1I4lVAExbVgQKigtLY3Bgx9iz549LJzpw8jH3dSOZPUcHTV8Os2bjuH2jHzhW06fPsm6dT8RGBiodjSrZG4rQ7RQBMEGHDlyhPbt23Hi+AG2fR8gikkFDRnowu51AaQkHScioi379u1TO5JVErO8BKGG27BhA507d8LTNYP9vwRwVwex3EhltGttz/5NATQIyadLl3v55ptv1I5kdcSgvCDUYKtXr+ahhx6kT1cdu9f5ExIkrgQ3h6+3jq3f+vHYgw48+eSTLFy4UO1IggWIMRRBuI3ly5fz9NNP8+Qjziz61AetCvul1EQGg4bFs3xwcZZ47rnnKCws5Pnnn1c7llW43sowZwzFgmEqQBQUQbiFlStXMmLECEY+7sq8/3qj0YhiYkmSJPHZdG/sDRomTpyIVqtl3LhxasdSna0OyouCIgjl+O677xg2bBjD/yOKSVWSJIkP3/TCaFQYP3489vb2jBw5Uu1YQiWIgiIIN7Fr1y4ef/wxhgx0ZsHHophUNUmS+OSdOhQWKYwaNQpfX18eeOABtWOpRsG8PU3UunRUDMoLwj9cuHCBhx8exD0d7Vn6uRgzqS6SJPHlB94M6OPM44//h+PHj6sdSTXiSnlBqAFyc3MZOPABXJ0KWbPQ1yY3wLJlGo3E8i99CA2UGTCgHxkZGWpHEipAFBRB+Issywwd+iQXzp9m3XIfvDxvv2e8YHkuzhrWLfflWmYiQ4YMxmg0qh2p+ikWOFQgxlAE4S/Tpk1j7dr1/LjUn5ZNK7/9am1VUqJw/HQx0UcKiTlSxMGYIi4nmigskikuUdDrJQx2Ej51tLRvZyC8tYF2re1p09wOe/uyn23rhehZs9CbPkN28OKLL/L555+r9KpUYm63lZjlJQjq2bVrF++++y7TXvFi4H3OasexGSaTwpad+cxdksX23/MpKrn+0dhV54yj0QcHnHBGi4QGpVBGRibjagFrz2WwfHU6sqKg1UCHtg6MGeHKIw84YzBcLy7d7nLks+l1mPD6F/To0YMBAwao+VKrlblXu6t1HYqk2MBa0tnZ2bi5uZGVlYWrq6vacYQaJi8vj9atWxDgnc6OH/3FjK47kH7VxNLVWcxdkk18YgnuWje8TaG44oELHuik239WNSkmcskih0zSNZdJl9PwdNfy7FOuPPuUG6HBehRFYcDQZGLiHDh27CSenp7V8Opurjreh26co97SN9A4Vn5zNjm/kAsj3q/290wxhiLUeq+99hpJSZdZPKuOKCa3UVKi8N6nVwlue4E3PsikJMmf9nQj3NSTUKkxHpL3HRUTAK2kxU3yJEhqQJjShUh643StHp/OzaFBx4uMfTWV3DyFBR97U1hwjYkTa89V9GKWlyDYoF27djF79mw+mOJBw3o1e993cx09UUTH+xJ495NM/Esac5fcjxZ0wE3yssjulE6SK02kMO4yPUAjpQ1LVubS8p54Tp4tZtY0T1au/IYNGzZY4JXYAEUy/1CBKChCrZWXl8fTTw/j7o5OjB/prnYcq1VSovD+rAwieiVw4ZQdEUo3GkmtsJOqZuKCVtIRIjWig9ybglQ3eg2+wt6oAnp3deK5554RU4mtmCgoQq313nvvia6u28jLl3ngiUTe+TiDIFMTIkw9cZWqZxzDQXIiTO5CE9qy9Js8LlwqIT8vkylTplTL+dUklq8XBBty5coVPvvsU1541lV0dZXjWpaJXoOvsOvPIsKUe2gotUQjVe+1OZIkESw1IELuwZV4DTqNwqJFX3H69OlqzVHtbPQ6FFFQhFrp3XffxckRXhrroXYUq5SXL3P/44nEHjYRJnfBU/JRNY+T5EpbUzeKcw1oJIkXX5ysah7h5kRBEWqdU6dOsWTJYqZMcMXNVVwN/09Go8JDw5OIiS2htelu3Kqpi+t2HCQn2pi6oFfs+HnjL+zYsUPtSFVGzPISBBvx1ltvEOCnZ8xwsR/8zXwyN5Ntf+TTUu6Mm+Sldpwyro+r3IsGLU888bjacaqWjXV3gSgoQi0TExPDd9/9wDsvu/9ruQ8B4k4W8c7HGYQojfGSfNWOc1NOkitNCCMpKZn3339f7TjC34i/KKFW+eyzWdQPteepR1zUjmJ1jEaF4RNScVCcqE8LtePckj+heOHL9GnTyczMVDuOxdWKLq8ZM2bQvn17XFxc8PHxYdCgQZw6deqWj1m2bBmSJJU57O0rv6SAIFRWWloaa9asYfQwZ7HHyU18MjeT2LhCmpg6oK3m2VwVJUkSzYjAWGzimWeeUTuO5dWGWV67du1i3Lhx7Nu3j61bt1JSUkLv3r3Jy8u75eNcXV1JSkoqPS5dumRWaEGojCVLliBJMsOHiPXg/inzmonpMzMIobHVDMLfjr3kQGPa8OOPP3LkyBG141iYZIGj+lVoteHNmzeX+XrZsmX4+PgQHR3NvffeW+7jJEnCz8+vcgkFwQJMJhPz589myEAnsc/JTXz9bTbFJRBCY7WjVIg/oZwjjtmzZ7Nw4UK149R6Zo2hZGVlAdx2BdDc3FxCQ0MJDg5m4MCBHDt27Jb3LyoqIjs7u8whCObYvHkzFy9eFjO7bkKWFWYvzsaHQAySbXVHayQNgdRn+bLlNet9ojZ0ef2dLMtMmjSJu+66i5YtW5Z7vyZNmrBkyRLWr1/PypUrkWWZzp07c/ny5XIfM2PGDNzc3EqP4ODgysYUBAAWLpxPu9aOtA8TG2f90/Y/CjgfX0yg0kDtKJUSSH1KSkpYsWKF2lEsp7YVlHHjxhEXF8fq1atveb/IyEiGDh1KWFgYXbp04ccff8Tb25sFCxaU+5gpU6aQlZVVeiQkJFQ2piCQl5fHli2/8tiDjhZZFbemmbv0Gq5aF9ypo3aUSrGXHKiDP5/N+gwb2N6pRqvUjo3jx49n48aN7N69m6CgoAo9Vq/X07ZtW86ePVvufQwGAwaD+CQpWMbWrVspKiqmf28ntaNYHVlW+G1XAX6mZjZdbP0J5ci5vSQmJhIYGKh2HPOZuwS9LUwbVhSF8ePHs3btWrZv3069evUqfEKTycTRo0fx9/ev8GMFoTJ++uknmjZyoFF9sQjkP505X0JegYwrtjGzqzyuXF+TLTo6WuUkllErVhseN24cK1euZNWqVbi4uJCcnExycjIFBQWl9xk6dGiZ5aWnTZvGr7/+yvnz54mJieHJJ5/k0qVLNXPuuGB1TCYTP/20jv69bWuwubpEHykEwAV3dYOYyYADOvTs379f7Si1WoUKyrx588jKyqJr1674+/uXHmvWrCm9T3x8PElJSaVfZ2ZmMmrUKJo1a8b9999PdnY2e/bsoXnz5pZ7FYJQjgMHDpCWliG6u8oRc6QIZ51DlW2WVV0kScIVD379davaUSzDRgflKzSGcicDXjt37izz9axZs5g1a1aFQgmCpfzyyy94edrRKVy0UG7m4KEiHI2eal0HZ1GueHIs7taXJNiM2jCGIgi25uDB/XQKtxNLrZTjSqIJexzVjmERDjhRUJiPLMtqR6m1REERaixFUYiKOkh4azEYX57CIgUtNWPlAM1fb2e5ubkqJzGfpJh/qEEUFKHGio+P5+rVa7RrbdvjA1WpuERBqjEF5frrOHz4sMpJLMBGx1BEQRFqrBtTSMNbi/GT8uj1Ego1o4tI/ut1xMXFqZzEAm6MoZhzqEAUFKHGio6Oxs/HQIBfpa7frRUcDBIyJrVjWMSN11HzVh62HaKgCDVWbGwM7Vrr1Y5h1fx8tRSSr3YMiygkH4NeQ2xsDbi4UXR5CYJ1ib90gXohNWN8oKpEhBnI12eqHcMicqVMAv21JCTEqx3FfKKgCIJ1SUxKJsBXdHfdSrvWBnJK8jAqJWpHMVueNpNGDfQkJ6dhMtWMbjw1XblypcKPEQVFqJGKiorIyMjCTxSUWwpvc30GXDa23UopUgrINxbRprkdJpNMenq62pHMo2ILJTk5mQkTJtCoUaMKP1YUFKFGSk5OBiDAV3R53UrThnbYGzQ2X1Bu5I9s7wBQZvknm1TFs7wyMzN57LHHqFOnDgEBAXzxxRfIsszbb79N/fr1OXjwIEuXLq1wbFFQhBopMTERAH/RQrklrVaia2cHrmpte8+hVOkyoYF62ra8fhHrjZ+/cHOvvfYae/bsYfjw4Xh5efHCCy/wwAMPEBMTw/bt29m3bx9Dhgyp8POKgiLUSDc+ofr7iIJyO6OHu5JpukaWkqF2lEopVopIJYFxI93w89EjSbbfQqnqK+U3bdrE0qVL+eSTT/jpp59QFIWwsDA2btxIp06dKp1bFBShRrqx/Iari/gVv537ezgR5KfnCufUjlIpiVxAp4PhQ1zR6yUcHHS2v/xKFY+hJCYm0qxZMwDq1q2Lvb09Tz75pNmxxV+bUCMZjUYAdKKBcltarcTYp11J1SRQrBSpHadCFEUhSXuO/zzojJfn9fEynVYq/fkLN6coCrq//XFotVocHBzMfl7x5ybUSEajEUkCjUasMnwnRjzmytT/ZpAgn6UBLdSOc8dSSCDPVMCYEcGlt+l0oqDcjqIo9OjRo7SoFBQU0L9/f+zsyi6kGhMTU6HnFQXlHxRFITExkaioKGJiYkhNTaWgoIDi4mIMBgMODg4EBgYSHh5OeHg43t7eakcWbuJO9u4R/p9PHR0vjnHnv7NP4qsE4Sy5qR3ptoqVIs5qYxnU24n2YWXXa7P1n7+EeSsG3+5j1NSpU8t8PXDgwMqf7G9EQQHS09NZtmwZ27dt5+DBg6RfvT6H3UHniL3GCY2iQVIkZElBkUzkmXIoNl3vGgjwD6RDx/bcd999PP7447i4uKj5UoS/6PV6FAVkWRGtlDv09ouerPsln5MXD9DO1AONZN094qelQ9g7mZj7UWCZ241GBb3expfcqeINtv5ZUCyl1hYURVHYv38/c+bM4ds132IyyXhKPrjiS5CuGa4aTww4IknSv8q9olUo0OaSLWeQnZrB7xv3sX79el6c/CLDRwxnzJgxtGhhO90GNdGNprzRCHZiO5Q7YjBoWD7bh8h+CVziFPVopnakcqUol0nmMt985Ievd9m3MaOp7PiAcHP79u3jp59+ori4mB49enDfffeZ/Zy18n998+bNvPrKaxw5ehhnvSuhcnMCdfWxk+5smXNJknDEBUetC36EAlCgyeNK0VmWLFjGnDlzuOfue5j56Uzat29flS/ljiiKwuXLl4mNjSUhIYHExESSkpK4kniFxMTLJCcnkZdXgMlooqTEiCRJ6HRatDotzs6O+Pn5ExgYTIB/AP7+/gQEBBAcHEzbtm0JCAi4XnStjJPT9T3kc3Ll0sFa4fYiwux5ZZwH/519HE/FFzfJU+1I/1Ko5HNGG8Og3k4MGehc5ntGo0JBgbH052+zzF2P6zaP/f777xkyZAgODg7o9Xo+/fRTPvroI1566SUzTlrLCsq1a9eYPHkyS5cuxUvnRzt9V7wkfySd+W+IDpITDXVtqK+0JFV3maP7jtOpUySvvvoKU6dOxWCovk2erly5wsGDB4mOjuZg1EGiow+Snnb9GgOdXoOHtz3u3jrcfTT4tdTTtJs9Dk5OaHQS2r/ee00mkI0KBbkmMtNSSEu7zJmDMtfSjGSmFWIyXt97wtevDuHhHYgIjyA8PJz27dvj7+9fba+1PDcyJKUaRUGpoLdf9GTHnwXEHv6dtqauVjWeUqQUcli7Cx9fhXn/9fnXh5mUNBOKglX8DpqligvKjBkzGDVqFHPmzEGr1TJjxgw++OADswuKpNjA6FV2djZubm5kZWXh6upaqefYtGkTT48YSUb6VRpKYQRqGlTpJ2tZkbloOs55+RiNGjVixcqvq6y1IssyUVFRbNiwgfXr1xIXdxwAD2976ja3p34rB+q1dKJec0c8fPVmjynIskJGcjHn4/K5EJfPhWMFXDhWQNbV6+NKYWGtGTjwQQYMGEDbtm1VacFcuHCB+vXrs3l1AL262PinVRVkXjPRZeAVzp+DNqYuOEuV+7uzpCKlkCPa3Rjc8/nz5yDqh/57nCQqtpCOfROIjo6mXbt2Fj2/Jd6H7vQcoR+8j8a+8hvDyYWFXHr9jXKzOjs7ExsbS8OGDQEoLi7GycmJK1eu4OPjU+nz1vgWiizLvPDCC3zxxRd46wLoqO2Lg1T1bzAaSUN9XUu85UBOnD9Ip06RfPzxf5k8ebJFnt9oNPLrr7+ybt06Nvy0jpTkNFzc7Qjr6sLEUfVp0s4ZD199lbyZazQSdQIM1Akw0KG3B3C9W+1qUjGnonOJ2XGFmbNm8O677xIQ6MuA/g/y4IMP0qNHD7Ta6mkt3PiEmpgiVp2tDA93LVu/D6Dnw4nEnt9Ba9M9uKrY/VWg5HFEuxsH9yK2/RBw02ICkJhyfbpwQEBAdcazOHP3hb/dY/Pz88sUGjs7O+zt7cnNzRUFpTwlJSWMGDGCVd+soqkunGBN42r/tOyi8SBC6clZ5TAvvvgiGRkZTJ8+vdI5UlJSWLRoEfMXzOVyQiIBdZ1o38+Z8O5NaNzOGa0Fuu8qQ5L+v8jc1d8LY4nMqehcorZdY90vXzN//nzq1gthzOhxPP3009SpU6dK89jb2+Ph4UpyirgeobJ8vXXsWh/I/Y8nEnNkJ/XkloTQqFr/hhRFIYXLnNXG4OML2368ecvkhqQUExqNxvan81dxlxfAokWLcHb+/zEoo9HIsmXLyvxtPv/88xU6bY3t8jKZTDz22GP8+MOPtNBE4qcNqeKUt3fReILTpkO8+uqrfPjhh3f8OEVR+OOPP5gzdw4//vADGi107u9Or8e9qd/S+rtzFEXh7OE8tq5KY98v15DQ8uijjzJ27Dg6duxYZW9QLVs2oXunVD57z8bfXFSWny/z5odX+WLRNdw1njQ1tcdJqvrp8UVKIaelGFKURB5+wJk5M7zxrnPrz8DvfnKVRf/TceVKisXzVGeXV933zO/yuvhm+V1edevWve3fnSRJnD9/vkLnrZEtFEVReO655/jh+x9orb0LH23w7R9UDerqmiEh8dFHH+Hu7s5rr712y/srisLWrVuZ8vqrxETHElDXkcde9ufeh7xwdrOdH50kSTQKc6ZRmDNPvlbCrh+u8uv/fmTFipV07NieDz/8L127drX4eYOD63IhQaw6ay5HRw2fTvPmwfudGfF8KgevbKWe3IIgGqCVLP97KCsyKSRwThuLg7PM6o/9GNz/zgrYhfgSgoLqWjxTtaviFsrFixdv+f3Lly8zbdq0Cp/Wuq9cqqS5c+eyePFimms7WE0xuSFU15T62pZMmTKFzZs3l3u/AwcO0L1HN/r06UOecp4pSxrxyZZm3D/C16aKyT+5eurpP8qPWb8145WFDblacIpu3bpx3319OHTokEXPFRbWjkNHRZeXpdzTyYEjO4MZ+7QrZ6Wj7NFu5LRymDwlxyLPX6QUcF45zj7dLxzjIP36GjjxR8gdFxOAQ0dNtG2r/lR9c1X1asO3c/XqVRYvXlzhx9W4gnL+/HleeullgrWNCNDWVzvOTTXQtsJbF8DTI0aSlZVV5nsnT57k4UcepmPHjpy/EsVL8xvyzupGtLnHrUZd8a3RSLTr5s707xsx6Yv6HD39J+3ateOxxx/j3DnLrHobHh7OlaRCklNFUbEUR0cNs6Z7c3pvXZ4f7UK263n2soVYaReJykVylaw7XvZEURTylVxSlASOspc/pV9INJzkif/YEb01hG+/8r9tF9ff5efLHD9dQHh4eGVfnmAm2/2oexOyLDNi+Ai0Rh2NtGFqxymXJEk01bRnf9omXnjhBZYsWUJRURHTpk3jo48+xNPXwJiP6nLPQC802ppTRG5GkiQ69fWkfS8Pdv6Qzo9fruOH77/nzTffYsqUKWYtoXHjjSX6SBH9etaoX3XV1Q/V8+GbdXjnJU+++ymXOYuzOHg4CgC9RouL5IaTyRMHnNCgRYMGBRkTMkXkk6vJJFdzjSLT9b3sm9Sz4/Wn6/DUYBfcXCs3EzD2WBGyrNSMglLFS69UlRr1VzZv3jx2/76bcH13dJJ1r+XjIDnRUApj6dKlhIWFsWDhPE6fPs1D4/0YMMoPvaHGNR5vSauT6DHEm3sGerF2bhLTpr3D2nU/sHzZClq3bl2p56xbty4eHq7EHCmkX0/rn7xgi+ztNTw12JWnBrtyLcvEobgiYo4UEX24kAMx8VxKKaGo+P9bLAa9hKeHjrva2hHRxpV2rQ2EtzHgU4GWSHlijhRhZ6enZcuWZj+X6qphlldVqDEFJT09nZdffoVgbSO8NH5qx7kjgZoGpGgTmDRxEnWbO/H+j00JbeqodixV2dlrGDI5kIhe7ix47TwREeG89dbbvPbaaxVurUiSRHh4BNFHDlRRWuHv3N20dLvLkW53lf0dVhSFkhLQ66nSKcfRRwpp3brFv5Zgt0VVfR3KQw89dMvvX7t2rVLnrTEfg5csWUJxURENtK3UjnLHJEmisSYMBYUBz3nX+mLydw1aOfH+j4154Blv3n13Kh06RnD69OkKP0/79h3ZF12MLFv97PgaS5Ik7OykKr9+ZW+UkfbtI6v0HDWFm5vbLY/Q0FCGDh1a4eetUAtlxowZ/Pjjj5w8eRIHBwc6d+7MRx99RJMmTW75uO+++4633nqLixcv0qhRIz766CPuv//+Coctj8lkYs7sOfhKIXe8wKO1cNF44KHz4rdVV4m830vtOFZFb7jeWmnf2505L56jfYcI1qz+tkKrovbt25cZM2Zw4FAhncLN35FOsE6nzhZz5nwBM/v2VTuKZVRxl9fSpUvNePLyVaiFsmvXLsaNG8e+ffvYunUrJSUl9O7dm7y8vHIfs2fPHh577DFGjhzJoUOHGDRoEIMGDSIuLs7s8Dds2bKF+IR4gjSNLPac1SmIxhzbn82VcwVqR7FK9Vs6Mf37RjRsp6Nfv37MnDnzjmcSRUZG4uXlzk9byv8dFWzfT7/m4eBgoEePHmpHsQxzpwyr1CCvUEHZvHkzw4cPp0WLFrRp04Zly5YRHx9PdHR0uY/5/PPPue+++3j55Zdp1qwZ06dPp127dsyePdvs8DfMnj0Hd30d3CTb/ITvqwnBoLXjt/+lqR3Fajm66HhpXn0eeMaHl156iWHDh1FYWHjbx+l0Ovr1G8DGrba1V7pQMes25eLh7sG+fftsfrdGW2bWGMqNayg8PctfNG7v3r307NmzzG19+vRh79695T6mqKiI7OzsMkd5cnNz2bJlM37K7ZcSsFYaSYsf9diz8ZraUayaRivx+MtBjP+kHmu+XUWXLveQnp5+28cNGDCAuJP5nL9UUg0pheqWftXE/uhCMlKz6dGjB40bNeaLL76o9MCyVVAscKig0gVFlmUmTZrEXXfddctpesnJyfj6+pa5zdfXl+Tk5HIfM2PGjDIDRMHB5V/tHhsbiyzLeEi2vV6Tu+RN1tViMlOL1Y5i9e4e6MXUbxpz+vxRunS9h5SUW6/b1Lt3b+zs9Gz8NbeaEgrV6ZdtecgKdND3JcLQm9xLRl6Y9AIBAYF88cUXyLKsdsSKq20FZdy4ccTFxbF69WpL5gFgypQpZGVllR4JCQnl3jcqKgqtRoeTFW0CVBmumuutvPNx+SonsQ0NWjvx1oqGJKdd5J577+LKlSvl3tfFxYUePbqzep0Yo6qJVv2Yg6e+DvYaRzy1vrS2u4d77B+kTkkwEydO5N577uXs2bNqx6wVKlVQxo8fz8aNG9mxYwdBQUG3vK+fn9+/PkGmpKTg51f+tSIGgwFXV9cyR3mio6Nx03qikWx7BrQ9Thh0dlw4JgrKnQps6MDbqxqSmZNI9x5db9lSefbZ0eyPySPmyO3HXQTbce5iMb/tzsdfalzmdoPkSDO7DkQYenH44FFatWxlU60VtdfyqqwKvQsrisL48eNZu3Yt27dvp169erd9TGRkJNu2bStz29atW4mMtMx88f37DuBksu3WCVyfq++ieHL+qJiNVBF+ofa88XUD0jMv06NHt3LHVB544AGCgvyZvzzrpt8XbNOCr7Ow0+jx04be9PueWj866O7HxxTKxIkTefzxxykpEWNpVaVCBWXcuHGsXLmSVatW4eLiQnJyMsnJyRQU/H9XwtChQ5kyZUrp1xMnTmTz5s3MnDmTkydP8s477xAVFcX48eMt8gLOXziPk8b2CwqAE25cOStmI1WUf1173ljegMvJ5xkw8AGKiv79f6jT6XjuubGsWpvHtSyxi2NNUFAgs2hlDn5So1suo6+T9DS160Abuy58/933DBgwsMx7lmA5FSoo8+bNIysri65du+Lv7196rFmzpvQ+8fHxJCUllX7duXNnVq1axcKFC2nTpg3ff/8969ats8h6O0ajEZPJiLaGrCCjRUdxoW00ya1NYEMHXppfj+joKMaMGXPTqaPPPPMMRqPE19+WP2tQsB3f/pRLVo6JYN2dXX/mqwshTN+V3379jcGDB1t3S8VGB+Ur9E58J/O7d+7c+a/bBg8ezODBgytyqjty4zoETQ1ZQUaDlpIiMYe+shqFOfPM9GDmvnJ9wc1/bl/q5+fHww8/zLzlGxg/0r1GbQdQ2yiKwuxFWXjr/HHU3PnuiV7aAFpzL5s3Xb+mbuXKlVZ5uUFVr+VVVWz6nVijsen4/6KgUMNeUrW798E69BvpywsvvMDWrVv/9f2JEydx+lwBq9dZZlMoQR2btucTc7SQYG2zCj+2jjaA5rpIVq1axfLly6sgnYXYWOsEbLyg2P+157JMzegTlzGht7fpH4lVeOLlIFrf7crgRx/hzJkzZb7XqVMnBgx4gKn/zaK4WLQGbZEsK7w27Speeh+8NAGVeg5/XT0CdPV5fsLzt5xyLlSMTb97aTQaHBwcKFZqxkB2iVKMo3PlNhcS/p9GKzFhVl2c3I088cRjGI1ld2z84IMPuRBfxKJvxIwvW/S/tTkcO11EA007s7qrmugjMBXKjBw50vqWa7HRMRSbLigArVq2JkfJVDuGReRqM6jbUqyIawlOrjrG/DeE6OgYZs6cWeZ7LVq0YOjQp5g+K4vcPDEJwpYUFyu8+UEGvrog3LXmrY6hlww00XZgy5YtVtf1VSuuQ7FGHTt1IF9n+7N2ZEUmx3SN+i3FniiW0qitM/2e9uHtt9/i+PHjZb737rvTuJal8NmCmvFhpLb4amUWCUklNNC1tcjzeWuD8NfV5+WXXqG4WCx7ZC6bLyjh4eFkF2diVKx4CuAdyFWyMMmyKCgWNnhiIN5BdgwfPrRM11doaChjx47j47nZXEqw7d+d2iI13cjU/2bgr62Ps8bdYs9bT9eC9KtprFu3zmLPaTbR5aWO8PBwALKVDJWTmCdbyUCSILSZKCiWZGevYfSHwTft+nrnnXdw9/Bm1Evp1teHLvzLuFfTKMjT0kjfzqLP66xxx0vvx+wvLbelhrlEl5dKmjZtiouzC1flpNvf2YpdVRIJauiIvaMYlLe0Rm2duX+ED1Onvl1moVE3NzcWLVrKtt25LFxh+92mNdm3G3L48ZdcGms6YpAsP84YKDXi9z9+59ixYxZ/7trE5guKTqdj+IjhJEsXkRXbnD5cpBSQKl+m+xDb3CDMFjw8PgB7Z4mpU6eWub1Pnz6MHDmSV6ZliK4vK5WabmTsK2n46ULw09WtknP4aINx0Dkxb968Knn+ChNdXuoZM2YMBcZ8UuTyl7m3ZpdNZ9HrJe59UBSUquLgrGXQGB+WL1/+rwH6Tz/9FA9PH9H1ZaVudHU10XeosnNoJC11CGLzpi1Vdo4KEQVFPc2aNaPLvV1I5JzaUSpMVmQSpbPc/aAnTq41Y00ya9XzP97UCbDnjTdfL3O7q6srX321hG27c/lswTV1wgk3tXR1VpV2df2dq8aLc+fP3nKHWOHWakRBARg/YTxXjSlky7Y1OJ8qX6bAWEDvx33UjlLj6Q0aHpnoy7q169m3b1+Z7/Xp04eXX36ZV6Zf5dedYgsBa7A3qoAxL6cRqG1YZV1df3djk7tDhw5V+bluRwzKq2zgwIE0bNCIU0QjK7ZxsVqJUswZomnV2ZW6zcXsrupwd39PQps4M+X1V//1vRkzZtCnT28eG53GmfPimgQ1XU4sYdCwZFyoQzO7quvq+jsnyQ2dRk90dHS1nO+WRJeXuvR6PStWfk2W8SqXTCfUjnNHTpsOIeuKefaDULWj1BoarcSD43zYuWM3R44cKfM9rVbL//63Bl+/UAYNTyUr2zYnedi6ggKZgUOTKci2o5W+CxqpemY+aiQNrjpPUVDMUGMKClxf+O/Fl17kvBJHrmzd6zSly4lcMZ3jqdeD8A40qB2nVono6Y6nj/1NZ/S4ubmxYcPPJKVqeWJsKiaTGKSvToqiMPKFVI6dKKGVrluVj5v8k53JgSuXxWKRlVWjCgrAu+++S7169TmhHLDaacQlSjEnlP207ORKj//UUTtOraPTa+g2xIOvVyy/6QBs48aN+fbbH9iyI58xr6Qiy6KoVAdFUXjtvausWZ9DM91dpWMa1UmLlvz8/Go/7z+JMRQr4eDgwMpvVpBDJsdM+1GsbDzFpBiJlXeCoZjnPgy1ys19aoMej3pTVFjIihUrbvr93r17s2TJEhavymbSW2I6cXV45+MMPpmbSVN9e/x06nQDS5KWwkIrWL1cdHlZj44dO7J69WpS5HhOmKKs5s3ApJg4bPqdPG0Gry9tJLq6VOTpZ0dETw/mzP2y3N+PYcOGsWDBAuYsucbkt0VRqSqKojBt5lXem5VBI307QvRNVcxiwt4g/i4rq0YWFICHH36YJUuXcEU+xzHTPtVnfhmVEmLlHWRpU3h5YUMatXVWNY8AvZ6ow4njp9i7d2+593n22WeZM2cOXyy6xvgpaaL7y8IUReGND67y7icZNNSHUU/fQtU8Jkw4OKq/hYStdnnV6Cvphg0bhoODA0888QTFciHNNR2xl6p/em6ufI04ZQ/F+hzeWNKIphEu1Z5B+LfmHV1wr2PPunXr6Ny5c7n3Gzt2LAaDgVGjRpGbp7DwE28Mhhr7WazaGI0KL7ydxtylWTTRhxOqb652JEo0hfgH+Ksdw/xuK9HlVTUeffRRNm3ahNZTYZ9pE1dM56qt60JWZM4bj7HPuBmnwELeXdNEFBMrotFItO3mzIYNa29735EjR7Jq1Sq++6mAHo8kk5xqvO1jhPJlZJq47z+JzF+WTXO7jlZRTBRFJku+WrqCuVBxNb6gAPTs2ZMTJ4/znyeGcMy4n8PybgqVqp3JkStf46D8K+fkw/R7xoePfm4ulqa3QuE93Dl16iynT5++7X3/85//sHv3H1y84kTHvolEHy6shoQ1z/FTRbTvfZk9+4y0s+tJkK6x2pEAyFOyMZpKiIiIUDuKGJS3dh4eHixfvpyNGzcieZrYY/yZY8b9Fl2qRVEUMuQUDht/Z59xE46B+Uz7tilPvBKEnegisUqtOrtgZ9Dy008/3dH9O3TowMGDMfgHtuDeQYmsXpdTxQlrlo1bc+nU9zIZKQ6019+Pp9ZP7UilbrwXtGtn2f1WKkOywKGGWvcu169fP06cPM4706Yie+ezr2Qz+4o3k2i6UOldH0uUIuJNp9gr/0xUyTY0QWkMfTOIj35uTqMwMfhuzQwOWlp2dmX9HXR73RAYGMiuXX/wyCP/4YkxybzwVhr5+dY1Pd3aFBcrvP3RVQYNS8K5JIAIXV8cNdbV/ZstX6VuaD3c3NzUjmKzavSgfHk8PDx44403ePXVV/n888956aWXiDPuBSScJBdcJc/rh8YLe8kRDRo0aJExYcJEgZJLtpxBtpJBrvYqeSV5aDQQ0cuDPk82pnlHF3F9iQ1p192VJVP3kpOTg4vLnb3JOTg48PXXKwgPj+C1115h0/YiFs/y4q4O6s8QsjYxRwoZNj6Vk2eLqa9rQ31dK6v7+5AVmXTpMk/d94TaUa6z0UH5WllQbtDpdLi4uCBJGjpKvcjmepHIljNJJh7FVP6nTjuDhrrNHWnXxon6LX1o0ckFT1+7akwvWErjtk7IssyhQ4e499577/hxkiQxadIk7rvvPkaMGEqXQQeZ9Kw701/1wsGh1jX+/6W4WOH9zzL44PMMXLTudDD0VOXq9zuRZrpMfkkuo0ePVjsKYP7UXzFtWCXR0dG46txxUdxxwZ1A6gNw1hRHkuEkL8yrj7FIwViioLOTsLPX4Omjx7+ePRqtdX3KEionsIEDBnsdUVFRFSooNzRt2pQ//tjLp59+yltvvcHPvxWx4GMv7o2sva2Vg7GFPP389VZJXW1r6utaVtsij5VxRTlDp46dCAsLUzvKdaKFYpuiDkbhaHSFf/yu55JJwzAn2twt+lNrOq1OIrSZo1mrzGq1Wl5++WUeeOABRowYSreHonigtwvvT/GgZdPac+X1mfPFvDnjKt9vzMVNb92tkhvy5CzSSxL5bMJ/1Y5i82p9uzw9PR0D//4kmavNoH5rMc23tqjbwsDBqP1mP0+zZs3Ys2c/q1at4vhZT8K6xzP8+ZQav199UoqRMa+k0vyeS/y8uYQWdpG01/Wz+mICcNF4HA8PTx555BG1o5RlY1OGQRQUCgoK0Uhl/xuKlAIKjIXUbykKSm1Rv4UTZ8+cJyfH/GnAGo2Gxx57jBMnTjN79hx+3W1H07vjeeGttBpXWBKTjbz+fjqNIuP5/mcIDg5FgwFfbei//q6s0VVTEleMZ3n//fcwWNEaXra69EqFf+K7d++mf//+BAQEIEkS69atu+X9d+7ciSRJ/zqSk5Mrm9miTCYj0j/+G3KUawBiF8VapF4LRxRF+demW+aws7Nj7NixnDt3kbffnsby74w07HSJgcOS2bw9z2bXBVMUhR1/5vPoqGTqRlxk9tJCJk9+jfPnL7F161aM2mLOlsSqHfO2jEoxJ+X9dO3Sleeee07tODVChQtKXl4ebdq0Yc6cORV63KlTp0hKSio9fHysYw91OzvDv/ZNKSQfSYI6AWLWVm3hHXT9Z3358mWLP7eTkxNvvPEGV64ks2DBQi6n1qXfE4k06XyZmfMySb9qnfv2/NO1LBOzF1+jZZcr9HzkCsfP+fH551+SmJjMe++9h5ubG40aNeLDD2cQbzxJhilF7ci3dLokBllnYsnSJWg0VtaastEr5Ss8KN+3b1/69u1b4RP5+Pjg7u5e4cdVNQd7e2TK/kEXUYizqx6d3sp+yYQq4+Csxd5BR2JiYpWdw8nJiWeeeYaRI0eyb98+5s6dw5sffsuU969ybycnHuhlT/8+TjSoaz0fZC4llLBxax4//VrAzj35KIrEgw8+yLyF4+jSpctNryeZMGEC3675lqPR+4nQ9MZOslch+a2lGhO4bDzD3C/mUq9ePbXj/IutThuutnfMsLAw/P396dWrF3/++ect71tUVER2dnaZo6rUa1CPAimvzG3FSgHuPvoqO6dgfSRJwsPHQFJSUrWcKzIykhUrVnL58hVmz56LvetdTPkgi8aRl2jV9Qqvv5/Ojj/zyc6p3tZLbp7M7/sKmPrfq7TreYX6HS4yeWoGir4DM2d+xqVL8Xz77Xd07dq13IsTtVotK1auwOCiI9a4E6NSXK2v4XYyTCnEGf9g4MBBoqvLwqp82rC/vz/z588nIiKCoqIiFi1aRNeuXdm/f3+5a+bMmDGDd999t6qjAdfXZtr/+4EyTcQiCgjwq/Uzqmsdd5+qbaHcjLe3N6NHj2b06NHk5uby66+/smHDBhav3sBHs6/vbd64gSPtWmsIb21PeGsDrZoZ8HDXmHW1uaIoZGXLxJ0sJuZIIdFHiog+YuLkmXwUBdzdXejX70GmvNWf++67r8LLkTRo0ICtv22ly71dOFy0kza6rugk9VtemaZUjhh3cW+Xe1mzZrX1dXXdIK5DubkmTZrQpEmT0q87d+7MuXPnmDVrVrnbr06ZMoXJkyeXfp2dnU1wcHCV5AsPDyevJJcibSGGv5rmJdoCPP3U/+UXqpebt4bEpCuqnd/Z2ZmHHnqIhx56CJPJxIkTJ4iKiiI6Opro6AOs33yYgoJ0AAwGDf6+BgJ8Nfj5SPj7avHx0mIwSOi0ElotmExgNCkUF0PaVRNJqUaSUhSSUmSSUospKDD99Vx2tGnTkq49OvLiK+GEh4fTsmVLdDrz3h7atm3Lr1t/pU/vPsQUbiNM103V7q900xWOlPzOXXd3Zv2G9VY1q+ufbLXLS5WP4R06dOCPP/4o9/sGg6Haftg39j7IVjLwlgIAKJYKcfd2r5bzC9bDw1vPuX2WH5SvDK1WS8uWLWnZsiXDhw8HwGg0cvLkSU6cOFFmgkti4mVO779MWno6xcUlGI1GjEYTOp0WnU6HXq+jTh1v/P2DqN8kiLu6+hMQEIC/vz9NmjShefPm6PVV08XbqVMndv++mx49enIgexNNNB3w1gZWybnKY1JMnC85zEXTce6//36+//577O2tb1ynJlCloMTGxuLvbwW7ogF169bF1dWN7JwMvLleUIyyEXtHK20KC1XG3lFLfn7V7pNjDp1OV1pkbEmbNm04dCiGkU+PZOtvWwnUNaCxPgJ9NXSBZZnSOSHvJ59spk+fzquvvmp2y6ta2GiXV4XfNXNzc4mNjSU2NhaACxcuEBsbS3x8PHC9u2ro0KGl9//ss89Yv349Z8+eJS4ujkmTJrF9+3bGjRtnmVdgJkmSGDCgP6nahNKdHBUUtDqxTldto9VLGI1iJ8aqEBwczJZft7Bo0SKyDCnsN/5MijEeRamaZf9LlCJOF8dwoHgz9ZuHcujQId544w3bKCZgs9OGK1xQoqKiaNu2LW3btgVg8uTJtG3blrfffhuApKSk0uICUFxczIsvvkirVq3o0qULhw8f5rfffqNHjx4WegnmGzduHLnGbK4q1y+2VBRFLPxYC2m1iIJShSRJYuTIkRw/cZy7ukRyuHgXe40buVByjGLFMrtfZstXOVa8l9+L1nJFOs306dM5cHC/zbXqbPVK+QqX665du95yT/Zly5aV+fqVV17hlVdeqXCw6tSxY0dat2rNlePnqIN1dMUJ6tBoxAeJqnajtXLgwAHmzJnDmtVrOF90BF9NCHU0QbhpvLCXnO5oFpusmMhVrnHNlE4KF8ksSSXAP4B3xk1l5MiR+PlZz46QtYGNtP+qliRJTHh+As+OepYCbR6SJCGbbHNZDKHyjEYFrVb8SVQHSZLo2LEjHTt25NNPP2Xp0qXMmzuPIxd3A2Cvc8AZD1zwwE5yQIMWjaTBpJiQMZGvZJOnuUZWSQayYkKr0dK1W1cmTJhAv379bKdrqzw2OoZi4//rlvPYY4/x4uQXuZB3DI1Gg8koCkptIxux/TciG1SnTh1efvllXn75ZZKTk/+aJh1NVFQUBw9EkZR1jaKiQhTlesG3NxgICQmlZ6d+hIdfn+bcpk0bHB1rztp7kqIg3aIn6E4erwbx1/MXJycnZn46k1GjRmGvMZCfaxvrKwmWU5BrwtnZVe0YtZqfnx/9+vWjX79+ZW5XFAVZltFqrXeTLkEsX1/GyJEj6dmzJ0aMXEutWcuMC7eXkVpCQED1XiMh3BlJkmpXMakts7xqMkmSWLx4MYoCGSmioNQ2WWkygQFBascQBJud5SUKyj+EhIRw7733cjXJuha0E6peVprRai64FQRbJArKTXTv3p2sdHE9Qm2iKApXUwsJCAhQO4og2GyXlxiUv4mAgADyso2UFMnoDaLm1gZ52SZKikyihSJYBVtdHFK8W95EaGgoAKmXi1ROIlSX1ITrP+uQkBCVkwiC7RIF5SbCwsIAuHDMehcKFCzrfFw+Go2GVq1aqR1FEGy2y0sUlJvw8vIiJDRIFJRa5EJcHs2aN6lRF8cJtkvM8qph2kd04EKcZRasE6zfxWPFtI/oqHYMQbhOtFBqlvDwCC4ez0eWxRIsNZ2xWCb+dG7pZmuCIFSOKCjliIiIID+3hORLYmC+pks4U0BJsSwKimBVbK27C0RBKVd4eDiSJHEqKlftKEIVOxmVi16vo02bNmpHEYTrFMX8QwWioJTD09OTTp06ELMjS+0oQhWL2ZFD165dxYC8IJhJFJRbGDjwQY7+kUNxUdVsUyqoLz/HyIkD2QwcOEjtKIJQSszyqoEGDBhAYYGRY3uz1Y4iVJHDv2djLJHp37+/2lEE4f+JWV41T9OmTWnQsC7R20S3V00VvS2LNm1aiSvkBcECREG5BUmSGND/QQ7tyBXTh2sgY4lM7M4cBgwYpHYUQShDks0/1CAKym08+OCDXE0p4Pj+HLWjCBZ2eHc2udnFPPjgg2pHEYSyRJdXzXT33XfTtFljtq5KVzuKYGFbV6UTEdGOtm3bqh1FEGoEUVBuQ5Ikxo2dQNTWTDJSxKZbNUXypUJid19j3LgJakcRhH8Rs7xqsKeeegqDvT3bvxWtlJrit/+l4e7hypAhQ9SOIgj/Ji5srLnc3Nx46smh7FiTgbFEXJNi64oLZXb9kMnIp0fh4OCgdhxB+BfRQqnhxowZw9WUQqJ+u6Z2FMFMe37OIOdaMaNHj1Y7iiDUKKKg3KE2bdrQvUc3fpydgmwSU4htlbFYZt3cVPoPeICGDRuqHUcQbk7M8qr5ZnzwIfGn8/hjw1W1owiVtO3bdFISCvjg/RlqRxGEcokur1qgQ4cOPPjQg3z/eQolYn0vm1OYZ2LtnBSeeupJWrZsqXYcQahxREGpoA/e/4D0pCK2rkpTO4pQQT8vTSE/W2batOlqRxGEWxOzvGqHpk2bMmLECNbPTyU/x6R2HOEOZWeU8PPiNMaOHUdoaKjacQThlmpNl9fu3bvp378/AQEBSJLEunXrbvuYnTt30q5dOwwGAw0bNmTZsmWViGo93n33XYoL4LvPr6gdRbhD//skEY1kx+uvv652FEGosSpcUPLy8mjTpg1z5sy5o/tfuHCBfv360a1bN2JjY5k0aRLPPPMMW7ZsqXBYaxEYGMj7789g89dpnDgo1viydrG7s9jxXRozP/kUb29vteMIwu3Z6CwvXUUf0LdvX/r27XvH958/fz716tVj5syZADRr1ow//viDWbNm0adPn4qe3mo8//zzfP/Dtyx8/QgfbmiCwUGrdiThJvJzjCx68zI9e/Zg1KhRascRhDtibreVzXR5VdTevXvp2bNnmdv69OnD3r17y31MUVER2dnZZQ5ro9VqWbpkOZnJRlbPFF1f1mrFjCsU5mpYvHgJkiSpHUcQarQqLyjJycn4+vqWuc3X15fs7GwKCgpu+pgZM2bg5uZWegQHB1d1zEpp3LgxH3zwoej6slJ/7+oSG2gJNkVWzD9UYJWzvKZMmUJWVlbpkZCQoHakcj3//PNEdu7IvFcSyM4oUTuO8JfM1GIWTkkQXV2CbbLRMZQqLyh+fn6kpKSUuS0lJQVXV9dyF+YzGAy4urqWOayVVqtl1TerkYvs+fz5i2LxSCtQXCTz6biLGHRufP31CtHVJdgcCTOnDauUu8oLSmRkJNu2bStz29atW4mMjKzqU1eb0NBQfvxhHaeic/n6/ctqx6nVFEXhqzfjSThZyPp1P+Hv7692JEGoNSpcUHJzc4mNjSU2Nha4Pi04NjaW+Ph44Hp31dChQ0vvP3r0aM6fP88rr7zCyZMnmTt3Lt9++y0vvPCCZV6BlbjnnnuYO3cev36TytZVqWrHqbV+XpLC7+vSWbx4Ke3bt1c7jiBUjo1eKV/hacNRUVF069at9OvJkycDMGzYMJYtW0ZSUlJpcQGoV68eP//8My+88AKff/45QUFBLFq0yKanDJdn1KhRHD58mPnT5xHYwIHmHV3UjlSrxO7OYtV/r/Daa6/x+OOPqx1HECrNVqcNS4qiUimrgOzsbNzc3MjKyrLq8RSAkpIS+tzXm4NRe3jj64bUa+GodqRa4fShXGaMOEf3br1Yv24DWq24LkiwrOp4H7pxjru7v4NOZ1/p5zEaC/lj+zvV/p5plbO8bJler2ftj+to2rglM0acI+H0zadGC5Zz7mgeH448R3i7DqxZ/a0oJoLtE7O8hBvc3Nz49dffqBfamPeHnRVFpQqdj8tjxohztGoexi8/b8LJyUntSIJgNklRzD7UIApKFfHw8OC3rdsJDWrM9CfPcuFYntqRapzTh3J5b+hZWjRtw+bNv+LiIsasBEFNoqBUoTp16rBj+y4aN2jB+0PPEbfX+paQsVWHdmUxY8Q52oV1YOvWbbi7u6sdSRAsR7bAoQJRUKqYh4cH27btoFOHe5gx4gxbVqZiA/MgrJaiKPz0VTL/HXWGHt17s3nTFtEyEWoc0eUllMvV1ZVNmzYzYcJElr4bz6K34jEWiyvqK6q4SGbuyxf55r+XmTLlddav2yDGTATBilT4OhShcnQ6HbNmzaJVq1aMHv0cieeLmfRlXdy89GpHswkZKcXMGneBhFNFrFq1iscee0ztSIJQdcydqSVmedUOTz/9NDt37uJqvJ43HzrN0T1iXOV2Du28xpsPnSYv3Yk//tgjiolQ89nolfKioKigc+fORB2MoXnj9rw/7DSL3r5EQa7Yn/6f8rKNzHv1Ih+NOktE2D1ERx0iPDxc7ViCUOVqzZ7ygmUEBwezfdsOZs+ezZ4Nubz6wEkxC+xvDu28xiv9ThLzWyGLFi1i8+Yt+Pn5qR1LEIRbEAVFRRqNhnHjxnH0SBzNGrXnvaHXWyu514xqR1NN1tWS0lZJeOt7OBZ3nJEjR4ol6IXaRXR5CZVVv359tm/bwZw5c9j7Ux4Texxn3fwkCvNrTzdYQa6J7764wqQexzm0rai0VWKtu3UKQlWSZPMPNYiCYiU0Gg1jx47l/LkLjBw+mh++SOGFnifYuiq1Rm/aVVIk88uyFCb1OM7GhVcZN2bi9f8D0SoRBJsjCoqV8fHx4fPPP+fUqdP0u+8RlryTwMv3n2TnD+kUF9acwlKYb2LbmjRevO8kK2dc4eFBT3DmzFk+/vhjvLy81I4nCOoSXV6CJdWrV48VX68gNjaW8FbdmP/aRcbdE8fKDxNIvlSodrxKSzxfyPL34hl3zzEWvRXPXR16ExcXx+LFi0X3liDcYKOrDYsLG61c69at+WnDRs6ePcv8+fNZvOQrNi6OI+xed3o9Xoc297qi06v3uSAn08j5uDwuHMvnfFw+WWlGSgplSkoU9HYSBgcNHn467Awa4k8Vcj4uD686HkwYO5nnnnuOunXrqpZdEATLEgXFRjRs2JBPPvmE6dOns2bNGubM+ZKPR8fg5GpH2L3OhPdwp829rji5Vu2PVFEUzh7O47f/pRH3Zy5XU4oA0Gt1uOCBneyKFi0SEoXI5GPispJDLpkoyEhIeLjVobCwkOLi4irNKgi2ytz1uNRay0vs2GjDYmNjWb9+PevXr+XQocNodRqatXehXTdXGrdzIqSpI3YGy7ReigpM/Lkxgy1fp3HpZD6OOke8TcG4Sp64SB444nzLQXRZkckjm2wlg2wlg3RtIoXGArp3686E5yfwwAMPoNOJzzeC9arOHRu7hU8xe8fGHdEzqv09UxSUGiIhIYGNGzeyfv06duzYQXFxCVqdhuBGTtRraaBeC0fqtXDEJ9iAi4cOjebOZlDJJoWfl6bw45dJFOSb8Nb4E0QjvCQ/s2ZhmRQTqUoCiZoLZBhT8ffz57PPP+PRRx+t9HMKQlUSBeX2REGpgQoLCzl69CjR0dFERUVxMGo/x4+dwGi8fl2LVqfBo44Bdx89bt4aPLz12Dtp0OokNNrrRUI2KWSllxCzI5vsjBKCpYaEaprgIDlbPG+2kslF5QQpcgIPP/Qwc+fNxcfHx+LnEQRzVGtBaTcFndaMgmIqZEdM9RcU0cdQA9nb29O+fXvat29felthYSHHjh3j8uXLJCYmkpSURGJiIolJiSSevExeXi5Go5GSkutX6efnFXAtMxsHyYn22ntxl+pUWV5XyYPWUmeSiefnDb/QdEcz5s+fJ1orQq1lq2MooqDUEvb29oSHh992ccXc3FwGDhjIjh07CJYa0VDTCq1UPb8mfpoQPBUfTmbFMGTIEHbsuL7WmVarrZbzC4LVUDDvWhIxbVhQW0ZGBvf1uY/Dhw7TVtMFL41vtWewk+xpLXXmMudYsGABmZmZrFixAr1e7BsjCNZOFBQBgJycHHr27MWJoycIowtuGk9V8wRpGqCX7fju2+9RFIVVq1aJlopQe5h7tbvo8hLUUlhYSP8H+nPsyDHa0gVXyUPtSAD4aoJBlvjuu+9wc3NjwYIFYn0voXaQAXN+1cXikIJa3nrrLf74409ac5fVFJMbfDVBNJMi+Oqrr1ixYoXacQRBuAVRUGq5vXv38unMT6lPczwkb7Xj3FSgpj7+mlAmjJ9AYmKi2nEEocrdmOVlzqEGUVBqsYKCAoY+NRQ3rRehUhO149xSE6ktJQUmRo0ahQ1cOiUI5hGrDQu25u233+bChYs0VSKQJOv+VdBLBhrLbfnll19E15cgWCnrfhcRqsy5c+f4dOan1KM5zpJtrD7gownEXxPKpImTKCgoUDuOIFSd2tRCmTNnDnXr1sXe3p6OHTty4MCBcu+7bNkyJEkqc9jbV35JAcEy5s+fj15jR4jUSO0oFVJPakHmtUy+++47taMIQtWpLQVlzZo1TJ48malTpxITE0ObNm3o06cPqamp5T7G1dWVpKSk0uPSpUtmhRbMU1BQwFcLF+Enh1bbVfCW4iS54K3158svZ6sdRRCEf6hwQfn0008ZNWoUI0aMoHnz5syfPx9HR0eWLFlS7mMkScLPz6/08PWt/iuwhf+3Zs0asrKvEaRpoHaUSglQGhAVdZDo6Gi1owhC1ZAtcKigQgWluLiY6Ohoevbs+f9PoNHQs2dP9u7dW+7jcnNzCQ0NJTg4mIEDB3Ls2LHKJxbM9uUXs/HWBuAouagdpVLqSP446VyYM2eO2lEEoUrUimnD6enpmEymf7UwfH19SU5OvuljmjRpwpIlS1i/fj0rV65ElmU6d+7M5cuXyz1PUVER2dnZZQ7BMlJSUog5FI2vEqJ2lErTSBq8TUH8tOEnMYVYqJlqyxhKRUVGRjJ06FDCwsLo0qULP/74I97e3ixYsKDcx8yYMQM3N7fSIzg4uKpj1ho3uoncJC+Vk5jHTfIk/Wo6V65cUTuKIAh/qVBBqVOnDlqtlpSUlDK3p6Sk4Ofnd0fPodfradu2LWfPni33PlOmTCErK6v0SEhIqEhM4Raio6Mx6OxxwEntKGZxla4vXinGUYQaSVbMP1RQoYJiZ2dHeHg427ZtK71NlmW2bdtGZGTkHT2HyWTi6NGj+Pv7l3sfg8GAq6trmUOwjKioKFwUd5tfZNGAAw46R1FQhJrJRru8KjxndPLkyQwbNoyIiAg6dOjAZ599Rl5eHiNGjABg6NChBAYGMmPGDACmTZtGp06daNiwIdeuXePjjz/m0qVLPPPMM5Z9JcIdObj/IM6yO9j4SvCSJOEsuxN1MErtKIIg/KXCBWXIkCGkpaXx9ttvk5ycTFhYGJs3by4dqI+Pj0ej+f+GT2ZmJqNGjSI5ORkPDw/Cw8PZs2cPzZs3t9yrEO6IyWQiKSWJ5pogtaNYhIPixLlz59SOIQhVwNxWhjotFEmxgWky2dnZuLm5kZWVJbq/zJCXl4ezszMtNZ3w14SqHcdsZ0xHMAXmEp8Qr3YUoRaojvehG+foWW8COo2h0s9jlIv47cKX1f6eKdbyqkWKiooA0Nh6f9dfNJKGwsIitWMIgvAX21p3QzCLTnf9x62odRmthSmKIvaaF2omWcGsbiuVZnmJglKLODg4ACBjUjmJZciYcBALjQo1kSJfP8x5vApEl1ctotPpcHRwpJhCtaNYRBGFeNWx7Qs0BaEmEQWlFpEkibC2bclWMtWOYhH5umzad2ivdgxBsDwbvQ5FFJRapkOH9uTrbH9tNKNSQk7JNcLDw9WOIgiWVxuulBdsX3h4ODklWZQoxWpHMUsO11BQREERaibRQhFswY03YFvv9spWMjDYGcQFsoJgRURBqWUaN26Ml4cXaYptr9J7VZNMp06dSqdCC0KNomBmC0Wd2KKg1DJarZZnRz9LiiYek2JUO06l5CpZXDUmM3rMaLWjCELVEF1egq147rnnMMolJCmX1I5SKZfls3h51uGhhx5SO4ogCH8jCkotFBoayv39+pGkuWBzOx4alRJSNPGMHvMcdnZ2ascRhKohy+YfKhAFpZYaP34c14xXyVRS1Y5SIYnKBYyykWeffVbtKIJQdUSXl2BLevXqRccOHTmtOWQzYykFSh7npWOMeHoEISEhascRBOEfREGppTQaDcu/Xk6xppCz8lG149yWoiicJJo63nWYOXOm2nEEoWqJFopga5o0acL7H7xPgnKGTCVN7Ti3dEU5T7opiaXLluDm5qZ2HEGoWuJKecEWTZo0iQ4dOnBKirbaq+fzlRzOSkd4+umn6dOnj9pxBEEohygotZxWq+XrFV8jOSgc4U+rG08pVPI5LP1BcEgwn376qdpxBKFaKIps9qEGUVAEGjduzOYtm8nX53CUPVZTVIqUAg5Lf+Dq7cyOndtFV5dQeyhmdneJMRRBTZGRkWzc+BM5+msc5g+MSomqeQqUPGKknTh6Gdi+Y7uY1SXULmJQXrB1PXr0YNu23yh2yCda2k6WclWVHKnyFaKl7dQJ9GTvvr00adJElRyCIFSMKChCGZ07d2bf/n3Ub1GPg/I2zpiOICvVs2VwsVJEnLyPw/IfdO/TnX3791G3bt1qObcgWBVxpbxQUzRv3pwDB/czffp0LmvPcFDaVqWtFUVRSJUvc0D6lXzna6xYsYKfftqAr69vlZ1TEKya6PISahK9Xs8bb7zBoUOHqN+iLgdMvxHDTpLleIu1WIxKCQnyGQ5IWzks/0n3Pt04cfIETz75JJIkWeQcgiBUH7GZhHBLLVu25GDUAdatW8fsL2eza/cuHHSO+Jnq4qsJwgk3NNKdfy4xKUZyyCRJvnR9CX1kBvUfyLjx4+jWrZsoJIIAKLKMIlW+20qtacOioAi3pdPpeOSRR3jkkUc4fvw48+bNY+mSpVzIP45Wo8NV44GT0RUXyQN7HNGgRSNJyIqMCRMF5JKtZJKvyybbmImCgo+3D6+PfZ1Ro0YRGBio9ksUBOuiKJi1S5ZKXV6ioAgV0rx5c7788ks++ugjDh06RHR0NFFRUezfd4BT52KQbzIYqNfpadGiBR069iYiIoLw8HBatWqFXq9X4RUIglBVREERKsXR0ZG77rqLu+66q/S2/Px8MjIyKCwspLi4GIPBgIODA3Xq1BF7lwhCRcgKSKKFItRijo6OODo6qh1DEGyfogBmjIOIWV6CIAiCLRMtFEEQBCujyAqKGV1eam3tXakWypw5c6hbty729vZ07NiRAwcO3PL+3333HU2bNsXe3p5WrVrxyy+/VCqsIAhCraDI5h8qqHBBWbNmDZMnT2bq1KnExMTQpk0b+vTpQ2rqzfcm37NnD4899hgjR47k0KFDDBo0iEGDBhEXF2d2eEEQhJpIkRWzDzVISgXbRh07dqR9+/bMnj0bAFmWCQ4OZsKECbz22mv/uv+QIUPIy8tj48aNpbd16tSJsLAw5s+ff0fnzM7Oxs3NjaysLFxdXSsSVxAEwSKq433oxjm6Sg+ikyo/rd6olLBTWVvt75kVGkMpLi4mOjqaKVOmlN6m0Wjo2bMne/fuvelj9u7dy+TJk8vc1qdPH9atW1fueYqKiigqKir9OisrC7j+ny0IgqCGG+8/1TE+YVSKzOq2MqLO9hMVKijp6emYTKZ/Ldrn6+vLyZMnb/qY5OTkm94/OTm53PPMmDGDd99991+3BwcHVySuIAiCxeXk5FTZZm92dnb4+fnxR7L548x+fn7Vfv2XVc7ymjJlSplWzbVr1wgNDSU+Pl7s2mem7OxsgoODSUhIEN2HZhL/l5ZjC/+XiqKQk5NDQEBAlZ3D3t6eCxcuUFxcbPZz2dnZYW9vb4FUd65CBaVOnTpotVpSUlLK3J6SkoKfn99NH+Pn51eh+wMYDAYMBsO/bndzc7PaXzZb4+rqKv4vLUT8X1qOtf9fVscHWnt7+2ovBJZSoVlednZ2hIeHs23bttLbZFlm27ZtREZG3vQxkZGRZe4PsHXr1nLvLwiCINimCnd5TZ48mWHDhhEREUGHDh347LPPyMvLY8SIEQAMHTqUwMBAZsyYAcDEiRPp0qULM2fOpF+/fqxevZqoqCgWLlxo2VciCIIgqKrCBWXIkCGkpaXx9ttvk5ycTFhYGJs3by4deI+Pj0ej+f+GT+fOnVm1ahVvvvkmr7/+Oo0aNWLdunW0bNnyjs9pMBiYOnXqTbvBhIoR/5eWI/4vLUf8X9YMFb4ORRAEQRBuRiwOKQiCIFiEKCiCIAiCRYiCIgiCIFiEKCiCIAiCRVh9QanoUvnCze3evZv+/fsTEBCAJEm3XEtNuLUZM2bQvn17XFxc8PHxYdCgQZw6dUrtWDZp3rx5tG7duvSCxsjISDZt2qR2LKGSrLqgVHSpfKF8eXl5tGnThjlz5qgdxebt2rWLcePGsW/fPrZu3UpJSQm9e/cmLy9P7Wg2JygoiA8//JDo6GiioqLo3r07AwcO5NixY2pHEyrBqqcNV3SpfOHOSJLE2rVrGTRokNpRaoS0tDR8fHzYtWsX9957r9pxbJ6npycff/wxI0eOVDuKUEFW20K5sVR+z549S2+73VL5gqCGG9sreHp6qpzEtplMJlavXk1eXp5YmslGWeVqw1C5pfIFobrJssykSZO46667KrT6g/D/jh49SmRkJIWFhTg7O7N27VqaN2+udiyhEqy2oAiCLRg3bhxxcXH88ccfakexWU2aNCE2NpasrCy+//57hg0bxq5du0RRsUFWW1Aqs1S+IFSn8ePHs3HjRnbv3k1QUJDacWyWnZ0dDRs2BCA8PJyDBw/y+eefs2DBApWTCRVltWMolVkqXxCqg6IojB8/nrVr17J9+3bq1aundqQaRZblMluAC7bDalsocPul8oU7l5uby9mzZ0u/vnDhArGxsXh6ehISEqJiMtszbtw4Vq1axfr163FxcSndztrNzQ0HBweV09mWKVOm0LdvX0JCQsjJyWHVqlXs3LmTLVu2qB1NqAzFyn355ZdKSEiIYmdnp3To0EHZt2+f2pFs0o4dOxTgX8ewYcPUjmZzbvb/CChLly5VO5rNefrpp5XQ0FDFzs5O8fb2Vnr06KH8+uuvascSKsmqr0MRBEEQbIfVjqEIgiAItkUUFEEQBMEiREERBEEQLEIUFEEQBMEiREERBEEQLEIUFEEQBMEiREERBEEQLEIUFEEQBMEiREERBEEQLEIUFEEQBMEiREERBEEQLEIUFEEQBMEi/g/7Am8AfrdMdgAAAABJRU5ErkJggg==", - "text/plain": [ - "
" - ] - }, - "metadata": {}, - "output_type": "display_data" - } - ], - "source": [ - "import ase.io # noqa: E402\n", - "import matplotlib.pyplot as plt # noqa: E402\n", - "from ase.visualize.plot import plot_atoms # noqa: E402\n", - "from matplotlib.colors import LogNorm # noqa: E402\n", - "\n", - "\n", - "structure = ase.io.read(\"ethanol_reduced_100.xyz\")\n", - "norm = LogNorm(vmin=min(lpr), vmax=max(lpr))\n", - "colormap = plt.get_cmap(\"viridis\")\n", - "colors = colormap(norm(lpr))\n", - "ax = plot_atoms(structure, colors=colors, rotation=\"180x,0y,0z\")\n", - "custom_ticks = [1e10, 2e10, 5e10, 1e11, 2e11]\n", - "cbar = plt.colorbar(\n", - " plt.cm.ScalarMappable(norm=norm, cmap=colormap),\n", - " ax=ax,\n", - " label=\"LPR\",\n", - " ticks=custom_ticks,\n", - ")\n", - "cbar.ax.set_yticklabels([f\"{tick:.0e}\" for tick in custom_ticks])\n", - "cbar.minorticks_off()\n", - "plt.show()" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# missing: test that the two covariances are the same!\n", - "# code to do force uncertainties (probably external)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [] - } - ], - "metadata": { - "kernelspec": { - "display_name": "Python 3", - "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.12.3" - } - }, - "nbformat": 4, - "nbformat_minor": 0 -} diff --git a/src/metatrain/experimental/soap_bpnn/model.py b/src/metatrain/experimental/soap_bpnn/model.py index e63df3932..54980c0b6 100644 --- a/src/metatrain/experimental/soap_bpnn/model.py +++ b/src/metatrain/experimental/soap_bpnn/model.py @@ -300,7 +300,6 @@ def load_checkpoint(cls, path: Union[str, Path]) -> "SoapBpnn": def export(self) -> MetatensorAtomisticModel: dtype = next(self.parameters()).dtype - print(dtype) if dtype not in self.__supported_dtypes__: raise ValueError(f"unsupported dtype {self.dtype} for SoapBpnn") diff --git a/tests/resources/llpr.py b/tests/resources/llpr.py deleted file mode 100644 index dbbd12a10..000000000 --- a/tests/resources/llpr.py +++ /dev/null @@ -1,221 +0,0 @@ -import matplotlib.pyplot as plt -import numpy as np -import torch -from metatensor.torch.atomistic import ( - MetatensorAtomisticModel, - ModelEvaluationOptions, - ModelMetadata, - ModelOutput, - load_atomistic_model, -) - -from metatrain.utils.data import Dataset, collate_fn, read_systems, read_targets -from metatrain.utils.llpr import LLPRUncertaintyModel -from metatrain.utils.loss import TensorMapDictLoss -from metatrain.utils.neighbor_lists import get_system_with_neighbor_lists - - -model = load_atomistic_model("model.pt", extensions_directory="extensions/") -model = model.to("cuda") - -train_systems = read_systems("train.xyz", dtype=torch.float64) -train_target_config = { - "energy": { - "quantity": "energy", - "read_from": "train.xyz", - "file_format": ".xyz", - "key": "energy", - "unit": "kcal/mol", - "forces": { - "read_from": "train.xyz", - "file_format": ".xyz", - "key": "forces", - }, - "stress": { - "read_from": "train.xyz", - "file_format": ".xyz", - "key": "stress", - }, - "virial": False, - }, -} -train_targets, _ = read_targets(train_target_config, dtype=torch.float64) - -valid_systems = read_systems("valid.xyz", dtype=torch.float64) -valid_target_config = { - "energy": { - "quantity": "energy", - "read_from": "valid.xyz", - "file_format": ".xyz", - "key": "energy", - "unit": "kcal/mol", - "forces": { - "read_from": "valid.xyz", - "file_format": ".xyz", - "key": "forces", - }, - "stress": { - "read_from": "valid.xyz", - "file_format": ".xyz", - "key": "stress", - }, - "virial": False, - }, -} -valid_targets, _ = read_targets(valid_target_config, dtype=torch.float64) - -test_systems = read_systems("test.xyz", dtype=torch.float64) -test_target_config = { - "energy": { - "quantity": "energy", - "read_from": "test.xyz", - "file_format": ".xyz", - "key": "energy", - "unit": "kcal/mol", - "forces": { - "read_from": "test.xyz", - "file_format": ".xyz", - "key": "forces", - }, - "stress": { - "read_from": "test.xyz", - "file_format": ".xyz", - "key": "stress", - }, - "virial": False, - }, -} -test_targets, target_info = read_targets(test_target_config, dtype=torch.float64) - -requested_neighbor_lists = model.requested_neighbor_lists() -train_systems = [ - get_system_with_neighbor_lists(system, requested_neighbor_lists) - for system in train_systems -] -train_dataset = Dataset({"system": train_systems, **train_targets}) -valid_systems = [ - get_system_with_neighbor_lists(system, requested_neighbor_lists) - for system in valid_systems -] -valid_dataset = Dataset({"system": valid_systems, **valid_targets}) -test_systems = [ - get_system_with_neighbor_lists(system, requested_neighbor_lists) - for system in test_systems -] -test_dataset = Dataset({"system": test_systems, **test_targets}) - -train_dataloader = torch.utils.data.DataLoader( - train_dataset, - batch_size=4, - shuffle=False, - collate_fn=collate_fn, -) -valid_dataloader = torch.utils.data.DataLoader( - valid_dataset, - batch_size=4, - shuffle=False, - collate_fn=collate_fn, -) -test_dataloader = torch.utils.data.DataLoader( - test_dataset, - batch_size=4, - shuffle=False, - collate_fn=collate_fn, -) - -loss_weight_dict = { - "energy": 1.0, - "energy_positions_grad": 1.0, - "energy_grain_grad": 1.0, -} -loss_fn = TensorMapDictLoss(loss_weight_dict) - -llpr_model = LLPRUncertaintyModel(model) - -parameters = [] -for name, param in llpr_model.named_parameters(): - if "last_layers" in name: - parameters.append(param) - print(name) - -# llpr_model.compute_covariance(train_dataloader) -llpr_model.compute_covariance_as_pseudo_hessian( - train_dataloader, target_info, loss_fn, parameters -) -llpr_model.compute_inverse_covariance() -llpr_model.calibrate(valid_dataloader) - -exported_model = MetatensorAtomisticModel( - llpr_model.eval(), - ModelMetadata(), - llpr_model.capabilities, -) - -evaluation_options = ModelEvaluationOptions( - length_unit="angstrom", - outputs={ - "mtt::aux::last_layer_features": ModelOutput(per_atom=False), - "mtt::aux::energy_uncertainty": ModelOutput(per_atom=False), - "energy": ModelOutput(per_atom=False), - }, - selected_atoms=None, -) - -force_errors = [] -force_uncertainties = [] - -for batch in test_dataloader: - print("new_batch") - systems, targets = batch - systems = [system.to("cuda") for system in systems] - for system in systems: - system.positions.requires_grad = True - targets = {name: tmap.to("cuda") for name, tmap in targets.items()} - - outputs = exported_model(systems, evaluation_options, check_consistency=True) - energy = outputs["energy"].block().values - energy_sum = torch.sum(energy) - energy_sum.backward(retain_graph=True) - - predicted_forces = -torch.concatenate( - [system.positions.grad.flatten() for system in systems] - ) - true_forces = targets["energy"].block().gradient("positions").values.flatten() - - force_error = (predicted_forces - true_forces) ** 2 - force_errors.append(force_error.detach().clone().cpu().numpy()) - - last_layer_features = outputs["mtt::aux::last_layer_features"].block().values - last_layer_features = torch.sum(last_layer_features, dim=0) - ll_feature_grads = [] - for ll_feature in last_layer_features.reshape((-1,)): - ll_feature_grad = torch.autograd.grad( - ll_feature.reshape(()), - [system.positions for system in systems], - retain_graph=True, - ) - ll_feature_grad = torch.concatenate( - [ll_feature_g.flatten() for ll_feature_g in ll_feature_grad] - ) - ll_feature_grads.append(ll_feature_grad) - ll_feature_grads = torch.stack(ll_feature_grads, dim=1) - - force_uncertainty = torch.einsum( - "if, fg, ig -> i", - ll_feature_grads, - exported_model._module.inv_covariance, - ll_feature_grads, - ) - force_uncertainties.append(force_uncertainty.detach().clone().cpu().numpy()) - -force_errors = np.concatenate(force_errors) -force_uncertainties = np.concatenate(force_uncertainties) - - -plt.scatter(force_uncertainties, force_errors, s=1) -plt.xscale("log") -plt.yscale("log") -plt.xlabel("Predicted variance") -plt.ylabel("Squared error") - -plt.savefig("figure.pdf") diff --git a/tests/resources/options.yaml b/tests/resources/options.yaml index 492cfaa0b..977c68e9f 100644 --- a/tests/resources/options.yaml +++ b/tests/resources/options.yaml @@ -3,33 +3,17 @@ seed: 42 architecture: name: experimental.soap_bpnn training: - batch_size: 8 - num_epochs: 100 - log_interval: 1 + batch_size: 2 + num_epochs: 1 training_set: systems: - read_from: train.xyz + read_from: qm9_reduced_100.xyz length_unit: angstrom targets: energy: - key: energy + key: U0 unit: eV -validation_set: - systems: - read_from: valid.xyz - length_unit: angstrom - targets: - energy: - key: energy - unit: eV - -test_set: - systems: - read_from: test.xyz - length_unit: angstrom - targets: - energy: - key: energy - unit: eV +test_set: 0.5 +validation_set: 0.1 diff --git a/tests/resources/split.py b/tests/resources/split.py deleted file mode 100644 index 4c5902b62..000000000 --- a/tests/resources/split.py +++ /dev/null @@ -1,13 +0,0 @@ -import ase.io -import numpy as np - - -structures = ase.io.read("ethanol_reduced_100.xyz", ":") -np.random.shuffle(structures) -train = structures[:50] -valid = structures[50:60] -test = structures[60:] - -ase.io.write("train.xyz", train) -ase.io.write("valid.xyz", valid) -ase.io.write("test.xyz", test)