diff --git a/src/scores/categorical/binary_impl.py b/src/scores/categorical/binary_impl.py index 3e0fc5765..7ce89f033 100644 --- a/src/scores/categorical/binary_impl.py +++ b/src/scores/categorical/binary_impl.py @@ -1,24 +1,25 @@ """ This module contains methods for binary categories """ -from typing import Optional, Sequence +from typing import Optional import numpy as np import xarray as xr from scores.functions import apply_weights from scores.processing import check_binary +from scores.typing import FlexibleDimensionTypes, XarrayLike from scores.utils import gather_dimensions def probability_of_detection( - fcst: xr.DataArray, - obs: xr.DataArray, - reduce_dims: Optional[Sequence[str]] = None, - preserve_dims: Optional[Sequence[str]] = None, + fcst: XarrayLike, + obs: XarrayLike, + reduce_dims: FlexibleDimensionTypes = None, + preserve_dims: FlexibleDimensionTypes = None, weights: Optional[xr.DataArray] = None, check_args: Optional[bool] = True, -) -> xr.DataArray: +) -> XarrayLike: """ Calculates the Probability of Detection (POD), also known as the Hit Rate. This is the proportion of observed events (obs = 1) that were correctly @@ -75,19 +76,17 @@ def probability_of_detection( hits = hits.sum(dim=dims_to_sum) pod = hits / (hits + misses) - - pod.name = "ctable_probability_of_detection" return pod def probability_of_false_detection( - fcst: xr.DataArray, - obs: xr.DataArray, - reduce_dims: Optional[Sequence[str]] = None, - preserve_dims: Optional[Sequence[str]] = None, + fcst: XarrayLike, + obs: XarrayLike, + reduce_dims: FlexibleDimensionTypes = None, + preserve_dims: FlexibleDimensionTypes = None, weights: Optional[xr.DataArray] = None, check_args: Optional[bool] = True, -) -> xr.DataArray: +) -> XarrayLike: """ Calculates the Probability of False Detection (POFD), also known as False Alarm Rate (not to be confused with the False Alarm Ratio). The POFD is @@ -144,5 +143,4 @@ def probability_of_false_detection( correct_negatives = correct_negatives.sum(dim=dims_to_sum) pofd = false_alarms / (false_alarms + correct_negatives) - pofd.name = "ctable_probability_of_false_detection" return pofd diff --git a/src/scores/probability/roc_impl.py b/src/scores/probability/roc_impl.py index e02ce6ac5..25099f627 100644 --- a/src/scores/probability/roc_impl.py +++ b/src/scores/probability/roc_impl.py @@ -1,7 +1,8 @@ """ Implementation of Reciever Operating Characteristic (ROC) calculations """ -from typing import Iterable, Optional, Sequence +from collections.abc import Iterable, Sequence +from typing import Optional import numpy as np import xarray as xr diff --git a/src/scores/processing.py b/src/scores/processing.py index f09af80ee..a91e64dd9 100644 --- a/src/scores/processing.py +++ b/src/scores/processing.py @@ -1,11 +1,13 @@ """Tools for processing data for verification""" import operator -from typing import Union +from typing import Optional, Union import numpy as np import pandas as pd import xarray as xr +from scores.typing import FlexibleDimensionTypes, XarrayLike + INEQUALITY_MODES = { ">=": (operator.ge, -1), ">": (operator.gt, 1), @@ -17,7 +19,7 @@ EQUALITY_MODES = {"==": (operator.le), "!=": (operator.gt)} -def check_binary(data: Union[xr.DataArray, xr.Dataset], name: str): +def check_binary(data: XarrayLike, name: str): """ Checks that data does not have any non-NaN values out of the set {0, 1} @@ -27,7 +29,10 @@ def check_binary(data: Union[xr.DataArray, xr.Dataset], name: str): ValueError: if there are values in `fcst` and `obs` that are not in the set {0, 1, np.nan} and `check_args` is true. """ - unique_values = pd.unique(data.values.flatten()) + if isinstance(data, xr.DataArray): + unique_values = pd.unique(data.values.flatten()) + else: + unique_values = pd.unique(data.to_array().values.flatten()) unique_values = unique_values[~np.isnan(unique_values)] binary_set = {0, 1} @@ -35,17 +40,19 @@ def check_binary(data: Union[xr.DataArray, xr.Dataset], name: str): raise ValueError(f"`{name}` contains values that are not in the set {{0, 1, np.nan}}") -def comparative_discretise(data, comparison, mode, abs_tolerance=None): +def comparative_discretise( + data: XarrayLike, comparison: Union[xr.DataArray, float, int], mode: str, abs_tolerance: Optional[float] = None +) -> XarrayLike: """ Converts the values of `data` to 0 or 1 based on how they relate to the specified values in `comparison` via the `mode` operator. Args: - data (xarray.DataArray or xarray.Dataset): The data to convert to + data: The data to convert to discrete values. - comparison (xarray.DataArray, float or int): The values to which + comparison: The values to which to compare `data`. - mode (str): Specifies the required relation of `data` to `thresholds` + mode: Specifies the required relation of `data` to `thresholds` for a value to fall in the 'event' category (i.e. assigned to 1). Allowed modes are: - '>=' values in `data` greater than or equal to the @@ -60,7 +67,7 @@ def comparative_discretise(data, comparison, mode, abs_tolerance=None): are assigned as 1 - '!=' values in `data` not equal to the corresponding threshold are assigned as 1. - abs_tolerance (Optional[float]): If supplied, values in data that are + abs_tolerance: If supplied, values in data that are within abs_tolerance of a threshold are considered to be equal to that threshold. This is generally used to correct for floating point rounding, e.g. we may want to consider 1.0000000000000002 as @@ -108,17 +115,23 @@ def comparative_discretise(data, comparison, mode, abs_tolerance=None): return discrete_data -def binary_discretise(data, thresholds, mode, abs_tolerance=None, autosqueeze=False): +def binary_discretise( + data: XarrayLike, + thresholds: FlexibleDimensionTypes, + mode: str, + abs_tolerance: Optional[float] = None, + autosqueeze: Optional[bool] = False, +): """ Converts the values of `data` to 0 or 1 for each threshold in `thresholds` according to the operation defined by `mode`. Args: - data (xarray.DataArray or xarray.Dataset): The data to convert to + data: The data to convert to discrete values. - thresholds (float or iterable): Threshold(s) at which to convert the + thresholds: Threshold(s) at which to convert the values of `data` to 0 or 1. - mode (str): Specifies the required relation of `data` to `thresholds` + mode: Specifies the required relation of `data` to `thresholds` for a value to fall in the 'event' category (i.e. assigned to 1). Allowed modes are: @@ -135,13 +148,13 @@ def binary_discretise(data, thresholds, mode, abs_tolerance=None, autosqueeze=Fa - '!=' values in `data` not equal to the corresponding threshold are assigned as 1. - abs_tolerance (Optional[float]): If supplied, values in data that are + abs_tolerance: If supplied, values in data that are within abs_tolerance of a threshold are considered to be equal to that threshold. This is generally used to correct for floating point rounding, e.g. we may want to consider 1.0000000000000002 as equal to 1 - autosqueeze (Optional[bool]): If True and only one threshold is + autosqueeze: If True and only one threshold is supplied, then the dimension 'threshold' is squeezed out of the output. If `thresholds` is float-like, then this is forced to True, otherwise defaults to False. @@ -184,7 +197,7 @@ def binary_discretise(data, thresholds, mode, abs_tolerance=None, autosqueeze=Fa return discrete_data -def broadcast_and_match_nan(*args: Union[xr.DataArray, xr.Dataset]): +def broadcast_and_match_nan(*args: XarrayLike) -> XarrayLike: """ Input xarray data objects are 'matched' - they are broadcast against each other (forced to have the same dimensions), and the position of nans are diff --git a/tests/categorical/test_binary.py b/tests/categorical/test_binary.py index f93c93e77..6f97cb1aa 100644 --- a/tests/categorical/test_binary.py +++ b/tests/categorical/test_binary.py @@ -48,6 +48,14 @@ (fcst_mix, obs1, "a", True, None, expected_poda), # Fcst mix, obs ones, only reduce one dim (fcst_bad, obs0, None, False, None, expected_pod0), # Don't check for bad data (fcst_mix, obs1, None, True, weight_array, expected_pod_weighted), # Fcst mixed, obs ones, with weights + ( + xr.Dataset({"array1": fcst0, "array2": fcst1}), + xr.Dataset({"array1": obs0, "array2": obs1}), + None, + True, + None, + xr.Dataset({"array1": expected_pod0, "array2": expected_pod1}), + ), # Test with DataSet for inputs ], ) def test_probability_of_detection(fcst, obs, reduce_dims, check_args, weights, expected): @@ -90,6 +98,14 @@ def test_probability_of_detection_raises(fcst, obs, error_msg): (fcst_mix, obs0, "a", True, None, expected_poda), # Fcst mix, obs ones, only reduce one dim (fcst_bad, obs0, None, False, None, expected_pofd0), # Don't check for bad data (fcst_mix, obs0, None, True, weight_array, expected_pofd_weighted), # Fcst mixed, obs ones, with weights + ( + xr.Dataset({"array1": fcst0, "array2": fcst1}), + xr.Dataset({"array1": obs0, "array2": obs1}), + None, + True, + None, + xr.Dataset({"array1": expected_pofd0, "array2": expected_pofd1}), + ), # Test with DataSet for inputs ], ) def test_probability_of_false_detection(fcst, obs, reduce_dims, check_args, weights, expected): diff --git a/tests/probabilty/test_roc.py b/tests/probabilty/test_roc.py index 1a626a66a..7702ac02c 100644 --- a/tests/probabilty/test_roc.py +++ b/tests/probabilty/test_roc.py @@ -4,10 +4,10 @@ import dask import numpy as np import pytest -import roc_test_data as rtd import xarray as xr from scores.probability import roc_curve_data +from tests.probabilty import roc_test_data as rtd @pytest.mark.parametrize( diff --git a/tests/test_processing.py b/tests/test_processing.py index 9c5fbbf8e..55b76bd2f 100644 --- a/tests/test_processing.py +++ b/tests/test_processing.py @@ -3,13 +3,14 @@ import pytest import xarray as xr -from tests import test_processing_data as xtd from scores.processing import ( binary_discretise, broadcast_and_match_nan, check_binary, comparative_discretise, ) +from tests import test_processing_data as xtd + @pytest.mark.parametrize( ("args", "expected"), diff --git a/tutorials/ROC.ipynb b/tutorials/ROC.ipynb new file mode 100644 index 000000000..04cde611f --- /dev/null +++ b/tutorials/ROC.ipynb @@ -0,0 +1,1397 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Reciever Operating Characteristic (ROC)\n", + "\n", + "ROC curves are used to understand the ability of a probabilitic forecast of a binary event (e.g., chance of rainfall) to discriminate between events and non-events. It is a plot that compares the hit rate or probability of detection (POD) against the false alarm rate or probability of false detection (POFD) using a set of increasing probability thresholds (e.g., 0.1, 0.2, 0.3, ...) that convert the forecast into a binary forecast.\n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [], + "source": [ + "from scores.categorical import probability_of_detection, probability_of_false_detection\n", + "from scores.probability import roc_curve_data\n", + "import matplotlib.pyplot as plt\n", + "import numpy as np\n", + "import xarray as xr\n", + "from scipy.stats import beta, binom" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Imagine we have a sequence of binary forecasts and a sequence of binary observations, where 1 indicates a forecast or observed event and 0 indicates a forecast or observed non-event" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [], + "source": [ + "fcst = xr.DataArray([1, 0, 1, 0, 1])\n", + "obs = xr.DataArray([1, 0, 0, 1, 1])" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We can then use `scores` to calculate the POD ($\\frac{H}{H + M}$) and POFD ($\\frac{F}{C + F}$), where $H, M, F, C$ are the number of hits, misses, false alarms, and correct negatives respectively." + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "
<xarray.DataArray 'ctable_probability_of_detection' ()>\n",
+       "array(0.66666667)
" + ], + "text/plain": [ + "\n", + "array(0.66666667)" + ] + }, + "execution_count": 3, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "probability_of_detection(fcst, obs)" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "
<xarray.DataArray 'ctable_probability_of_false_detection' ()>\n",
+       "array(0.5)
" + ], + "text/plain": [ + "\n", + "array(0.5)" + ] + }, + "execution_count": 4, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "probability_of_false_detection(fcst, obs)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Now suppose that we have a forecast that is the probability of rain occuring, and corresponding observations where 1 means it rained, and 2 means that it didn't rain. We generate some synthetic forecasts using a Beta distribution and generate some corresponding observations." + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [], + "source": [ + "fcst = beta.rvs(2, 1, size=1000)\n", + "obs = binom.rvs(1, fcst)\n", + "fcst = xr.DataArray(data=fcst, dims=\"time\", coords={\"time\": np.arange(0, 1000)})\n", + "obs = xr.DataArray(data=obs, dims=\"time\", coords={\"time\": np.arange(0, 1000)})\n", + "\n", + "thresholds = np.arange(0, 1.01, 0.01)" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "
<xarray.Dataset>\n",
+       "Dimensions:    (threshold: 101)\n",
+       "Coordinates:\n",
+       "  * threshold  (threshold) float64 0.0 0.01 0.02 0.03 ... 0.97 0.98 0.99 1.0\n",
+       "Data variables:\n",
+       "    POD        (threshold) float64 1.0 1.0 1.0 1.0 ... 0.1125 0.0875 0.04063 0.0\n",
+       "    POFD       (threshold) float64 1.0 1.0 1.0 1.0 1.0 ... 0.0 0.0 0.0 0.0 0.0\n",
+       "    AUC        float64 0.8074
" + ], + "text/plain": [ + "\n", + "Dimensions: (threshold: 101)\n", + "Coordinates:\n", + " * threshold (threshold) float64 0.0 0.01 0.02 0.03 ... 0.97 0.98 0.99 1.0\n", + "Data variables:\n", + " POD (threshold) float64 1.0 1.0 1.0 1.0 ... 0.1125 0.0875 0.04063 0.0\n", + " POFD (threshold) float64 1.0 1.0 1.0 1.0 1.0 ... 0.0 0.0 0.0 0.0 0.0\n", + " AUC float64 0.8074" + ] + }, + "execution_count": 6, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "roc = roc_curve_data(fcst, obs, thresholds)\n", + "roc" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We can see that the output that we have calculated POD and POFD for each threshold. We also calculated the AUC which is the area under the ROC curve which we will discuss shortly. Let's now plot the ROC curve." + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkIAAAHFCAYAAAAe+pb9AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAB0NElEQVR4nO3dd3xT1f/H8Ve6dxmFsqHsDaUIAiLiAEEQRAQEFRAVvoCIDBUHKvIVJyrTxdCfgMgUEUFUVJayCrKRWUahtEBbupvc3x/322ptQQptb9u8n49HHklubpJPkrZ595xzz7EZhmEgIiIi4oRcrC5ARERExCoKQiIiIuK0FIRERETEaSkIiYiIiNNSEBIRERGnpSAkIiIiTktBSERERJyWgpCIiIg4LQUhERERcVoKQlLozZ07F5vNlnlyc3OjfPny9OnThz///NPq8gCoVq0aAwYMsLqMbBISEnjjjTcIDQ3Fz88PX19fmjZtyuuvv05CQoLV5V2z119/neXLl2fb/vPPP2Oz2fj5558LvKYMR48eZfjw4dSuXRtvb298fHxo0KABL774IqdPn87c77bbbqNhw4aW1Xkj5s+fz/vvv59vj389vz+bNm3ilVde4dKlS9luu+2227jtttvypDYp/mxaYkMKu7lz5zJw4EDmzJlD3bp1SU5OZuPGjfz3v//F39+fAwcOULJkSUtrDA8PJyAggBo1alhax9+dO3eOO++8kyNHjjBixAjuuOMOAH766Sc++OADatSowQ8//EBwcLDFlf47Pz8/evbsydy5c7Nsj4uLY9++fdSvX5+AgIACr2vlypX06dOHoKAghg8fTmhoKDabjd27dzN79mxcXFwIDw8HzC/n6Oho9uzZU+B13qguXbqwZ88ejh8/ni+Pfz2/P++88w5jx47l2LFjVKtWLctt+/btA6B+/fp5WaYUU25WFyByrRo2bEjz5s0B80vFbrfz8ssvs3z5cgYOHGhpbaGhoQX+nHa7nfT0dDw9PXO8/ZFHHuHAgQOsW7eOW265JXP7XXfdxT333EP79u3p378/q1evLqiSgX+vOzcCAgK4+eab86Cq3Dt27Bh9+vShdu3arFu3jsDAwMzbbr/9dkaMGMGyZcsKtCbDMEhOTsbb27tAn/d6JSUl4e3tnee/PwpAkhvqGpMiKyMUnTt3Lsv2bdu2ce+991KqVCm8vLwIDQ3lq6++ynb/06dP88QTT1C5cmU8PDyoUKECPXv2zPJ4cXFxjBkzhpCQEDw8PKhYsSIjR47M1q3096b98+fP4+HhwUsvvZTtOQ8cOIDNZmPKlCmZ286ePcvgwYOpVKkSHh4ehISE8Oqrr5Kenp65z/Hjx7HZbLz11ltMnDiRkJAQPD09WbduXY7vzbZt2/j+++8ZNGhQlhCU4ZZbbuHRRx9lzZo1bN++PXO7zWZj+PDhfPTRR9SuXRtPT0/q16/Pl19+me0xbrTu5ORkRo8eTdOmTQkMDKRUqVK0atWKr7/+Osvz2Gw2EhIS+OyzzzK7RzO6PXLqGhswYAB+fn4cPnyYzp074+fnR+XKlRk9ejQpKSlZHvvUqVP07NkTf39/SpQoQb9+/di6dSs2my1b69M/TZ48mYSEBGbMmJElBP297h49emTbvnXrVtq2bYuPjw/Vq1fnjTfewOFwZN5+re9LxnMMHz6cDz/8kHr16uHp6clnn30GwKuvvkrLli0pVaoUAQEBNGvWjFmzZpFTJ8D8+fNp1aoVfn5++Pn50bRpU2bNmgWY/3R8++23nDhxIksXdYbU1FQmTpxI3bp18fT0pEyZMgwcOJDz589neY5q1arRpUsXli5dSmhoKF5eXrz66quZt/29a8zhcDBx4kTq1KmDt7c3JUqUoHHjxnzwwQcAvPLKK4wdOxaAkJCQzJoyfg5y6hpLSUlhwoQJ1KtXDy8vL0qXLk379u3ZtGlTtvdDnItahKTIOnbsGAC1a9fO3LZu3TruvvtuWrZsyYcffkhgYCBffvklvXv3JjExMfOP7enTp7nppptIS0vj+eefp3HjxsTExLBmzRouXrxIcHAwiYmJtGvXjlOnTmXus3fvXsaPH8/u3bv54YcfsnwhZChTpgxdunThs88+49VXX8XF5a//N+bMmYOHhwf9+vUDzDDRokULXFxcGD9+PDVq1GDz5s1MnDiR48ePM2fOnCyPPWXKFGrXrs0777xDQEAAtWrVyvG9Wbt2LQDdu3e/4vvXvXt3Pv74Y9auXUtYWFjm9hUrVrBu3TomTJiAr68vM2bM4MEHH8TNzY2ePXvmWd0pKSlcuHCBMWPGULFiRVJTU/nhhx/o0aMHc+bM4ZFHHgFg8+bN3H777bRv3z4zXP5bN1haWhr33nsvgwYNYvTo0fz666+89tprBAYGMn78eMAcP9W+fXsuXLjAm2++Sc2aNVm9ejW9e/e+6mNn+P777wkODs5Vi9TZs2fp168fo0eP5uWXX2bZsmWMGzeOChUqZL7ea31fMixfvpz169czfvx4ypUrR9myZQEzhA4ePJgqVaoA8Ntvv/Hkk09y+vTpzPcAYPz48bz22mv06NGD0aNHExgYyJ49ezhx4gQAM2bM4IknnuDIkSPZWrgcDgfdunVj/fr1PPPMM7Ru3ZoTJ07w8ssvc9ttt7Ft27YsrVM7duxg//79vPjii4SEhODr65vj+/TWW2/xyiuv8OKLL3LrrbeSlpbGgQMHMscDPfbYY1y4cIGpU6eydOlSypcvD1y5JSg9PZ1OnTqxfv16Ro4cye233056ejq//fYbERERtG7d+po+PymmDJFCbs6cOQZg/Pbbb0ZaWpoRHx9vrF692ihXrpxx6623GmlpaZn71q1b1wgNDc2yzTAMo0uXLkb58uUNu91uGIZhPProo4a7u7uxb9++Kz7vpEmTDBcXF2Pr1q1Zti9evNgAjFWrVmVuq1q1qtG/f//M6ytWrDAA4/vvv8/clp6eblSoUMG4//77M7cNHjzY8PPzM06cOJHlOd555x0DMPbu3WsYhmEcO3bMAIwaNWoYqamp//aWGUOGDDEA48CBA1fcZ//+/QZg/Oc//8ncBhje3t7G2bNns9Rdt25do2bNmvlad3p6upGWlmYMGjTICA0NzXKbr69vlvc3w7p16wzAWLduXea2/v37G4Dx1VdfZdm3c+fORp06dTKvT58+3QCM7777Lst+gwcPNgBjzpw5V63Xy8vLuPnmm6+6z9+1a9fOAIzff/89y/b69esbHTt2vOL9rva+AEZgYKBx4cKFqz633W430tLSjAkTJhilS5c2HA6HYRiGcfToUcPV1dXo16/fVe9/zz33GFWrVs22fcGCBQZgLFmyJMv2rVu3GoAxY8aMzG1Vq1Y1XF1djYMHD2Z7nH/+/nTp0sVo2rTpVWt6++23DcA4duxYttvatWtntGvXLvP6559/bgDGJ598ctXHFOekrjEpMm6++Wbc3d3x9/fn7rvvpmTJknz99de4uZkNm4cPH+bAgQOZrS3p6emZp86dOxMZGcnBgwcB+O6772jfvj316tW74vOtXLmShg0b0rRp0yyP1bFjx389UqlTp06UK1cuS8vImjVrOHPmDI8++miW52jfvj0VKlTI8hydOnUC4JdffsnyuPfeey/u7u65e+OuwPhfF8k/W7XuuOOOLAOoXV1d6d27N4cPH+bUqVN5WveiRYto06YNfn5+uLm54e7uzqxZs9i/f/8NvTabzUbXrl2zbGvcuHFmK0dGjRk/S3/34IMP3tBzX025cuVo0aLFVeuC3L0vt99+e44HC/z000/ceeedBAYG4urqiru7O+PHjycmJoaoqCjAbDm02+0MGzbsul7PypUrKVGiBF27ds3yc9C0aVPKlSuX7XekcePGWVpwr6RFixbs2rWLoUOHsmbNGuLi4q6rvgzfffcdXl5eWX73RDIoCEmR8fnnn7N161Z++uknBg8ezP79+7N8aWWM7RkzZgzu7u5ZTkOHDgUgOjoaMMfxVKpU6arPd+7cOf74449sj+Xv749hGJmPlRM3Nzcefvhhli1bltmcP3fuXMqXL0/Hjh2zPMc333yT7TkaNGiQpd4MGV0A/yajOySj+zAnGUcAVa5cOcv2cuXKZds3Y1tMTEye1b106VJ69epFxYoV+eKLL9i8eTNbt27l0UcfJTk5+Zpe55X4+Pjg5eWVZZunp2eWx42JicnxiLlrPYquSpUqV31/c1K6dOls2zw9PUlKSsq8ntv3Jaf3dsuWLXTo0AGATz75hI0bN7J161ZeeOEFgMznyxjH82+/C1dy7tw5Ll26hIeHR7afhbNnz173z++4ceN45513+O233+jUqROlS5fmjjvuYNu2bddV5/nz56lQoUKWbmqRDBojJEVGvXr1MgdIt2/fHrvdzqeffsrixYvp2bMnQUFBgPlHNKdBqgB16tQBzHE8Ga0bVxIUFIS3tzezZ8++4u1XM3DgQN5+++3MMUorVqxg5MiRuLq6ZnmMxo0b89///jfHx6hQoUKW6zmNScrJXXfdxfPPP8/y5cuztXhkyJiX56677sqy/ezZs9n2zdiW8UWeF3V/8cUXhISEsHDhwiy3/3NAc34pXbo0W7ZsybY9p9efk44dOzJ16lR+++23PD1yLbfvS07v7Zdffom7uzsrV67MEgj/ORdTmTJlAHPQ+D8D8bUICgqidOnSVzzy0N/f/19rzYmbmxujRo1i1KhRXLp0iR9++IHnn3+ejh07cvLkSXx8fHJVZ5kyZdiwYQMOh0NhSLJREJIi66233mLJkiWMHz+eHj16UKdOHWrVqsWuXbt4/fXXr3rfTp068X//938cPHgwMxz9U5cuXXj99dcpXbo0ISEhua6vXr16tGzZkjlz5mC320lJScl2mH+XLl1YtWoVNWrUyNO5kJo3b06HDh2YNWsWDz/8MG3atMly+4YNG5g9ezZ33313loHSAD/++CPnzp3LbBmx2+0sXLiQGjVqZLYc5EXdNpsNDw+PLF+OZ8+ezfHoqH+2muSFdu3a8dVXX/Hdd99ldukBOR4hl5Onn36a2bNnM3To0GyHz4PZ9bh8+XLuu+++XNWVm/flao/h5uaWJXQnJSXxf//3f1n269ChA66ursycOZNWrVpd8fGu9P536dKFL7/8ErvdTsuWLa+5vtwoUaIEPXv25PTp04wcOZLjx49Tv379zOkXruXnolOnTixYsIC5c+eqe0yyURCSIqtkyZKMGzeOZ555hvnz5/PQQw/x0Ucf0alTJzp27MiAAQOoWLEiFy5cYP/+/ezYsYNFixYBMGHCBL777jtuvfVWnn/+eRo1asSlS5dYvXo1o0aNom7duowcOZIlS5Zw66238vTTT9O4cWMcDgcRERF8//33jB49+l//+D/66KMMHjyYM2fO0Lp162yha8KECaxdu5bWrVszYsQI6tSpQ3JyMsePH2fVqlV8+OGH191t8fnnn3PnnXfSoUOHHCdUrFu3bo6HiAcFBXH77bfz0ksvZR41duDAgSwBIS/qzjiUeujQofTs2ZOTJ0/y2muvUb58+Wwzhjdq1Iiff/6Zb775hvLly+Pv73/FAHut+vfvz3vvvcdDDz3ExIkTqVmzJt999x1r1qwB+NeWg5CQkMzWvqZNm2ZOqAjmhH6zZ8/GMIxcB6HcvC9Xcs899zB58mT69u3LE088QUxMDO+88062uZuqVavG888/z2uvvUZSUhIPPvgggYGB7Nu3j+jo6MzD2xs1asTSpUuZOXMmYWFhuLi40Lx5c/r06cO8efPo3LkzTz31FC1atMDd3Z1Tp06xbt06unXrluvXD9C1a9fMecPKlCnDiRMneP/996latWrmkZKNGjUC4IMPPqB///64u7tTp06dbK1QYI77mjNnDkOGDOHgwYO0b98eh8PB77//Tr169ejTp0+ua5RixNqx2iL/LuOosX8evWUYhpGUlGRUqVLFqFWrlpGenm4YhmHs2rXL6NWrl1G2bFnD3d3dKFeunHH77bcbH374YZb7njx50nj00UeNcuXKGe7u7kaFChWMXr16GefOncvc5/Lly8aLL75o1KlTx/Dw8DACAwONRo0aGU8//XSWI6v+edRLhtjYWMPb2/uqR6ycP3/eGDFihBESEmK4u7sbpUqVMsLCwowXXnjBuHz5smEYfx199fbbb+fqvbt8+bLx+uuvG02bNjV8fHwMHx8fo3HjxsbEiRMzH/vvAGPYsGHGjBkzjBo1ahju7u5G3bp1jXnz5uVL3W+88YZRrVo1w9PT06hXr57xySefGC+//LLxzz9NO3fuNNq0aWP4+PgYQOYRQVc6aszX1zfbc+X0uBEREUaPHj0MPz8/w9/f37j//vuNVatWGYDx9ddfX/W9zXDkyBFj6NChRs2aNQ1PT0/D29vbqF+/vjFq1KgsRzS1a9fOaNCgQbb79+/fP9sRWdf6vmR8XjmZPXu2UadOHcPT09OoXr26MWnSJGPWrFk5Hmn1+eefGzfddJPh5eVl+Pn5GaGhoVmOmrtw4YLRs2dPo0SJEobNZstSR1pamvHOO+8YTZo0ybx/3bp1jcGDBxt//vln5n5Vq1Y17rnnnhxr/efvz7vvvmu0bt3aCAoKMjw8PIwqVaoYgwYNMo4fP57lfuPGjTMqVKhguLi4ZPk5+OdRY4Zh/q0YP368UatWLcPDw8MoXbq0cfvttxubNm3KsSZxHlpiQ0Qy2Ww2hg0bxrRp06wuxTKvv/46L774IhEREdfdGiciRYe6xkTEaWUEvrp165KWlsZPP/3ElClTeOihhxSCRJyEgpCIOC0fHx/ee+89jh8/TkpKClWqVOHZZ5/lxRdftLo0ESkg6hoTERERp2XphAq//vorXbt2pUKFCthstmxzXOTkl19+ISwsDC8vL6pXr86HH36Y/4WKiIhIsWRpEEpISKBJkybXPDDz2LFjdO7cmbZt2xIeHs7zzz/PiBEjWLJkST5XKiIiIsVRoekas9lsLFu27KqrZT/77LOsWLEiy3o7Q4YMYdeuXWzevLkAqhQREZHipEgNlt68eXPm+jkZOnbsyKxZs0hLS8txUceUlJQsU9M7HA4uXLhA6dKlr3m6dxEREbGWYRjEx8fn+bpxRSoInT17NtuCiMHBwaSnpxMdHZ3jgn6TJk3KnB1VREREiraTJ0/m6fQWRSoIQfZF+zJ69q7UujNu3DhGjRqVeT02NpYqVapw8uRJAgIC8q9QERGRPGAYkJ9rEScnw+bN8OOP8MMPcOzYtd3P3x+qV4d/fpXWrAl9+kCJEjde267IaKqV8CfQ25PLl+O4447KOS6jciOKVBAqV65ctpWho6KicHNzy1wV+588PT2zra8DEBAQoCAkIiKWSk2FyEg4ffqv05kzWa+fPg15vObwVbm5QdmyWbcFB0OLFnDTTdCggRl2SpeG/BphkmZ38O73h/hw8xHa1gris4EtuHzZvC2vh7UUqSDUqlUrvvnmmyzbvv/+e5o3b57j+CAREZGC5HCYwebwYYiLy3mflBRYsgR++gmiogq2viupXh3uvhs6doT27c3WHqucuZTEkwvC2X7iIgDVSvuS7si/47osDUKXL1/m8OHDmdePHTvGzp07KVWqFFWqVGHcuHGcPn2azz//HDCPEJs2bRqjRo3i8ccfZ/PmzcyaNYsFCxZY9RJERMTJJCebrTZHj5qBJ+N05Ih5ym3rjYcHVKgAFSv+dfrn9aAgyMPxwVnYbODnlz+PnVs/7j/H6EW7uJSYhr+nG2/c35h7Gpvjf5Pz6TktDULbtm2jffv2mdczxvL079+fuXPnEhkZSURERObtISEhrFq1iqeffprp06dToUIFpkyZwv3331/gtYuISPFz8iT8/DOEh4Pdbm5LTMzaTRUTc/XHcHWFatWgVKkrdx01aQIDBkDt2vnbxVRUpNkdvLX6AJ+sNwcoNaoYyLS+oVQt7Zvvz11o5hEqKHFxcQQGBhIbG6sxQiIiTu7UKTP4ZJyOHLm2+3l5QUiIOVamRg3zPONy1aqg0Rq5czklnS5T1nM8JpGBbarxXKe6eLq5Ztknv76/i9QYIRERkRsRHw/r1sGaNfD992aX1t+5uEDz5tC6Nfj+rzHC0zN7V1XJkmrFyUt+nm5M69uM05eS6NigXIE+t4KQiIgUOenpfw00NgxISIBLlyA21jxlXP77thMnYNMmSEv763FcXKBZM3OA8G23wS23ZD8cXPJeSrqdSasOULW0DwPbhADQsGIgDSsGFngtCkIiIlKkhIdDjx5w/Pj13T8kBDp1Mo+QatcOAgv+u9epnYhJYPj8cHafjsXDzYXOjcoTHOBlWT0KQiIiUmQsWwYPPWQOYHZx+etIKl9fM9AEBpoT+WVc/vv1oCBo29Ycy6NuLWt8+0ckzy35g/iUdEr6uPNuryaWhiBQEBIRkUIsIQF27ICtW83ZjxcvNrd36AALF+bN7MWS/5LT7Ez8dh9f/GYeCd68akmm9g2lfKC3xZUpCImISCGRlgZ79pihZ8sW83zPHnOSwr8bPhzee8+cAVkKv3S7g94fbWbXqVgAht5Wg1F31cbNNZ8mRsol/RiJiIglEhPh999h/Xr49VezxScxMft+FSqYyzu0aAG33gpt2hR8rXL93Fxd6NSoPKcuJjG5d1Pa1S5jdUlZKAiJiEiBuHgRNm40g8/69bBtW9YjuMAcy9O8+V/B56abzMPVpWhJSrUTfTmFyqV8AHiibXV6hlUiyC/72p9WUxASEZE8FxMDf/zx12nbNti92zzU/e8qVjQHMGecGjTIv6UkpGAcjopn2Lxw7IbBiuFt8PFww8XFVihDECgIiYjIDUhNhYMHs4aeP/4w1+LKSa1aZvdWRvAJCdERXMXJ4u2neGn5HpLS7AT5eXIiJpF65Qv3xEwKQiIiclXJyXDsWPYFRg8fNufyyViT65+qV4fGjc1TkybmbM3lCnbSYCkgianpvLR8L0t2nAKgTc3SvNe7KWX9rT00/looCImISBZRUebyE6tXm2N5Tp7M3qX1dwEBfwWejFPDhuDvX3A1i3UOno1n2PwdHI66jIsNRt5Zm2Hta+LqUjSa+hSEREScXFqaecTW6tXmGlw7dmTfx9//r4VF/77QaI0a5jgfdW85rze+28/hqMsEB3jyQZ9Qbq5e2uqSckVBSETECR07ZoaeNWvgxx/NxUj/LjTUXILirrvM1p0yZRR2JGdv3t+Y11ft56Uu9SldSAdEX42CkIiIE0hLM1dd//Zbs+Xn0KGstwcFmcGnY0dz1ubgYGvqlMJv75lYfj54nmHtawJQNsCL9/uEWlzV9VMQEhEpptLTzfCzaBEsXWoe0p7B1RVatYK77zZPoaE6bF2uzjAMvvg9gtdW7iM13UHNsn50bFD0R78rCImIFCPp6fDLL/DVV2b4iY7+67YyZaBbN3Pl9Tvu0Krrcu3iktMYt2Q33+6OBOCOumVpUa2UxVXlDQUhEZEizG6H8HD4+Wez9Wf9+qzjfYKCoEcP6NUL2rXT+lySe3+cusTw+eFEXEjEzcXGc53qMuiWEGzFZNCYfiVERIoIux1OnTLn8Nm50ww+v/4KcXFZ9ytd+q/wc9ttCj9y/RZsiWD813tIsxtULOHNtL6hhFYpaXVZeUq/HiIihYjdDidOmIOZ/z5x4eHDcPSoOZPzPwUEmLM1t29vBp8mTcwxQCI3qpSvB2l2gw71g3m7ZxMCfdytLinPKQiJiFggIcEMOwcOwP795vmBA+a2lJQr38/d3ZyxuU4dc4mK9u2haVMFH8k7ianp+HiY8aBjg3J8+cTNtAwpVWy6wv5JQUhEJJ+lpMCyZfDbb38Fn4iIK+/v6WmuyfX3CQwzJi+sXFmhR/KHYRh8uv4Yn6w/ytfD21A+0BugyE2QmFsKQiIi+SQiAj76CD75BM6fz357UBDUrQv16pnnGZerVFHYkYJ1MSGVMYt28eOBKAAWbTvFiDtqWVxVwVAQEhHJQ4YBP/0E06bBihXgcJjbK1aEBx6ABg3+Cj1BQdbWKgKw7fgFnlwQTmRsMh5uLrzUpT4PtaxidVkFRkFIRCQPJCfDnDkwZYrZ/ZWhfXsYNsycv0dHb0lh4nAYfPjrEd79/hB2h0FIkC/T+obSoIJzTTClX0sRkRuQmGh2fb31Fpw5Y27z84NHHoGhQ80WIJHCaO6m47y1+iAA3ZpW4L/3NcLP0/ligfO9YhGRPHD5MsycCe+8A1HmsAoqVYKxY2HAAPOQdpHC7MEWVfh652n6tqxCr+aVi+1RYf9GQUhEJBdiY83xP++999faXdWqwbhx0L+/ecSXSGFkdxis2HWabk0q4uJiw9vDlWVD2+Di4pwBKIOCkIjINYiMhA8/NMcAXbpkbqtVC55/Hvr1M+f3ESmsouKTeXrhTjYejiEyNpmht5krxzt7CAIFIRGRKzIM2LjRbAFassRc0BTMQ9xffNFcwkIDoKWw23g4mqe+3En05RS83V0pF+BldUmFin6FRUT+ISEB5s2D6dPhjz/+2t66NYwcCfffDy4ulpUnck3sDoMPfjjE1HWHMQyoE+zP9H6h1Czrb3VphYqCkIgIEB0Nv/wCP/4I8+ebY4EAvL3Nrq9hw8ylLESKgnNxyYxYEM7vxy4A0OemyrzctQHeHpqp858UhETEKcXEmCu3r1sHP/8Mu3dnvb1GDfPw94EDoWTxWmxbnMD5+BTCIy7h6+HK6z0a0a1pRatLKrQUhETEKZw5A+vX/3XavdscA/R3DRuaq7d36QJ33aXuLym6GlYMZHLvJtQvH0D1Mn5Wl1OoKQiJSLFjGHD4cNbgc+RI9v3q1zeDT/v2cOutULZsgZcqkifOXEpizKJdPNepLo0rlQCgS+MK1hZVRCgIiUiRFhNjhp7Dh82ws2ePGXzOns26n81mjvFp2/avU3CwJSWL5KmfDpxj1Fe7uJSYxrilu1n55C1OOzni9VAQEpEi5fJlc1zPmjWwenXOLT0AHh5w001mS0/btuYRX4HOtYSSFHNpdgdvrT7AJ+uPAdCoYiDT+oYqBOWSgpCIFGqGYR7CnhF8NmyAtLSs+1SsCDVrmqdataBVKzMEeXtbU7NIfjt5IZEnF4Sz8+QlAAa0rsa4znXxdNNRYbmlICQihU5MDKxdawaf7783Z3X+u2rV4O67zVP79lrXS5zL0fOX6T59I3HJ6QR4ufFWzybc3bCc1WUVWQpCIlIoJCXBBx/AsmWwdWvWI7q8vc3Ac/fd0LGj2eqj1n9xVtVK+9KsakkuJqYx7cFQKpfysbqkIk1BSEQst3YtDBkCR4/+ta1hw7+Czy23gJdWBRAnFhGTSJC/Bz4ebri42PigTyje7q54uGmOhxulICQilomKglGjzOUswBzr89JLcM89UKmStbWJFBbf/hHJc0v+4O6G5Xj7gSYABHprld+8oiAkIgXOMGDOHBg7Fi5cMLu5nnwSJk4Efy2DJAJAcpqdid/u44vfIgA4Gp1AUqpdy2TkMQUhESlQBw/C4MHmul5gzu3z8cfmUV4iYjoWncCweTvYFxkHwH9uq8Gou2rj7qqusLymICQi+cpuN5e3OH7cHAv05puQmgo+PvDqq+Zq7m76SySS6eudp3l+6W4SUu2U8vVgcq8m3FZH057nF/35EZEb4nCYh7cfPw7Hjpnnfz9FRGSf96dTJ5gxwzwMXkT+EpecxoRv9pGQaqdFSCmm9AmlXKCOFMhPCkIics2Sk81Znb/7Dvbv/yvopKZe/X5ublC1KoSEwOOPwwMP6PB3kZwEeLkzuXdTth2/wFN31MJNXWH5TkFIRK7q7Fn49lv45huzaysxMfs+rq5QpYrZwhMSYp7//VShgrmPiGS3ZPspfD3dMidFbFe7DO1ql7G4KuehICQiWRgG7NwJK1ea4Wfr1qy3V6oEXbpAy5Zm6AkJMYOOxvmI5E5iajrjv97L4u2n8Pdyo0nlQMoHal2YgqY/XSKS6eOPYcIEOH066/abboKuXc1Tkybq1hK5UQfPxjNs/g4OR13GxQaPt61OWX+NBbKCgpCIAOZipkOGmC1CPj5w111m8LnnHiinZYxE8oRhGHy17STjv95LSrqDsv6eTHkwlJurl7a6NKelICQiXL4M/fubIeihh+CTT7SkhUheszsMRn+1k+U7zwBwa+0yTO7VhCA/T4src24KQiLCM8+Y63xVqQLTpikEieQHVxcbJXw8cHWxMbpDbYbcWgMXF/UzW01BSMQJGQacPAlbtphdYjNnmttnz4bAQGtrEylODMMgMdWOr6f5dTuuc13uC61Ik8olrC1MMikIiTgJhwM++ghWrTKPBDt3Luvtw4fDHXdYU5tIcRSXnMa4pbs5H5/C/Mda4ubqgqebq0JQIaMgJOIEEhLg4Ydh2bK/trm5QePG5hFhbdtCnz7W1SdS3Ow+FcvwBTs4EZOIm4uN8JOXuKlaKavLkhwoCIkUc6dOwb33Qng4eHjA+PFw++3mYqfemrJEJE8ZhsFnm47z+qoDpNodVCzhzdS+oTSrUtLq0uQKFIREirGtW80QdPYslCkDy5dD69ZWVyVSPMUmpvHMkl2s2Wv2O3eoH8zbPZsQ6ONucWVyNQpCIsXUV1+Zh8QnJ0PDhuYs0VrkVCT/PP3VTn46EIW7q43nO9djQOtq2DT7aKGnICRSzBiGOTv0K6+Y1++5BxYsAH9/S8sSKfae61SX0xeTePuBxjSuVMLqcuQaaVlbkWJk/364776/QtCoUfD11wpBIvnhUmIqq/dEZl6vHezPd0+1VQgqYtQiJFIM7N4NEyfCokVmi5Cbmzk30GOPWV2ZSPG0/cQFnpwfzrn4FBY+4Unz/x0RpgkSix4FIZEibMcOeO01cxB0hu7d4eWXzaPCRCRvORwGH/16lHe+P4jdYRAS5Iu3h6vVZckNsLxrbMaMGYSEhODl5UVYWBjr16+/6v7z5s2jSZMm+Pj4UL58eQYOHEhMTEwBVStSOPz+O3TpAmFhZgiy2aBXL9i1y5wrSCFIJO/FXE7h0c+28ubqA9gdBvc2qcA3T95Cgwqajr0oszQILVy4kJEjR/LCCy8QHh5O27Zt6dSpExERETnuv2HDBh555BEGDRrE3r17WbRoEVu3buUxtf+LE7hwAebMgTvvhJtvhm+/BRcX6NcP9u6FhQvNCRJFJO/9fjSGzlPW8/PB83i6ufBGj0Z80Kcpfp7qWCnqLP0EJ0+ezKBBgzKDzPvvv8+aNWuYOXMmkyZNyrb/b7/9RrVq1RgxYgQAISEhDB48mLfeeqtA6xYpKFFRZovP4sWwbh2kp5vb3dzMmaLHjYNatSwtUcQp7I+M41xcCjXK+DK9XzPqlguwuiTJI5YFodTUVLZv385zzz2XZXuHDh3YtGlTjvdp3bo1L7zwAqtWraJTp05ERUWxePFi7rnnnis+T0pKCikpKZnX4+Li8uYFiOSTM2dg6VJYsgR+/dVcIyxD48bQs6cZgjQnkEj+Mgwjcx6g/v+bE6hnWKXMBVSleLDs04yOjsZutxMcHJxle3BwMGfPns3xPq1bt2bevHn07t2b5ORk0tPTuffee5k6deoVn2fSpEm8+uqreVq7SF47ccIMPkuWwD//D2je3Aw/998PNWtaU5+Is9l0OJr3fjjE7AE34e/ljs1mo3/ralaXJfnA8sHS/5x18+8J/J/27dvHiBEjGD9+PNu3b2f16tUcO3aMIUOGXPHxx40bR2xsbObp5MmTeVq/yI3YsgVatDBbd0aP/isEtW4N774Lx46Zy2Q8+6xCkEhBsDsMJq89RL9Zv7P1+EWm/XTY6pIkn1nWIhQUFISrq2u21p+oqKhsrUQZJk2aRJs2bRg7diwAjRs3xtfXl7Zt2zJx4kTKly+f7T6enp54enrm/QsQuUHp6dC3Lxw5Yg56vvVWs9XnvvugYkWrqxNxPufiknnqy3B+O3oBgN7NKzPyztoWVyX5zbIg5OHhQVhYGGvXruW+++7L3L527Vq6deuW430SExNxc8tasqurOX+DYRj5V6xIPli0yAxBpUubEyLmkONFpID8eug8Ty/cSUxCKj4errx+XyO6h+o/Emdg6YivUaNG8fDDD9O8eXNatWrFxx9/TERERGZX17hx4zh9+jSff/45AF27duXxxx9n5syZdOzYkcjISEaOHEmLFi2oUKGClS9FJFccDnj9dfPyyJEKQSJWWrrjFKO+2gVAvfIBTO8bSvUyfhZXJQXF0iDUu3dvYmJimDBhApGRkTRs2JBVq1ZRtWpVACIjI7PMKTRgwADi4+OZNm0ao0ePpkSJEtx+++28+eabVr0Ekevy7bewZ4+5BtiwYVZXI+Lc2tUuQ3CAJ3fWC+alLvXxctdM0c7EZjhZn1JcXByBgYHExsYSEKB5IKTgGQa0amXODv3ss/DGG1ZXJOJ89kfGUa/8X98BFxNSKenrYWFF8m/y6/vb8qPGRJyJYcD06WYI8vKCp5+2uiIR55Jmd/D6qv10+mA9y8JPZW5XCHJemhVKpICcOmWuBr9mjXl92DC4wgGSIpIPTl1M5MkF4YRHXALg4NnL1hYkhYKCkEg+MwyYO9ccFB0XB56e8N//mtdFpGB8v/csYxbtIi45HX8vN97u2Zi7G+ooBVEQEslXp0/DE0/AqlXm9ZYtzVBUt66lZYk4jdR0B5O+28+cjccBaFIpkGl9m1G5lI+1hUmhoTFCIvnAMODzz6FhQzMEeXrCm2/Cxo0KQSIFaUfExcwQ9NgtISwa0lohSLJQi5BIHvvtt6zLZdx0k9kKVL++pWWJOKWbq5dmbMc61An25876GpQn2alFSCSPHD8OffqYh8Zv2gQ+PuakiZs2KQSJFJTkNDv//XYfJy8kZm4b1r6mQpBckVqERG7QpUtm4PngA0hNBZsNBgyA117TmmEiBelYdALD5+9g75k4tp24yJIhrXFxyXkRb5EMCkIi1yktDT76CF55BWJizG133AHvvANNm1pZmYjz+XrnaZ5fupuEVDulfD146o5aCkFyTRSERK7Dr7+aR4MdPGher1vXDECdO5stQiJSMJLT7Lz6zV4WbDkJQIuQUkzpE0q5QC+LK5OiQkFIJBfsdrPL67XXzIVTg4Lg1Vfh8cfB3d3q6kScS2RsEgPnbOXA2XhsNhjeviZP3VELN1cNf5VrpyAkco1OnYJ+/czWIDDHAb3/PgQGWlmViPMq6eOBq4uNID8P3u8dyi21gqwuSYogBSGRa/DNN2bwuXAB/Pzgww/NUCQiBSsp1Y6HmwuuLja83F358KEwPN1cKBugrjC5Pmo/FLkKh8OcE+jee80QFBYG4eEKQSJWOHQunnunbeCDH//M3Fa5lI9CkNwQBSGRq3jxRZg82bw8apQ5J1DNmtbWJOJsDMPgq60nuXfaBv6MusxXW09yOSXd6rKkmFDXmMgVfPIJTJpkXp49GwYOtLYeEWeUkJLOC8t2s3znGQDa1grivd5N8fPU15fkDf0kieRg9Wr4z3/My6+8ohAkYoV9Z+IYPn8HR6MTcLHB6A51+E+7GpofSPKUgpDIP+zaBQ88YB4q/8gjMH681RWJOJ+ElHT6fvoblxLTKBfgxZQHQ2kRUsrqsqQYUhAS+ZtTp+Cee+DyZWjf3uwe0wSJIgXP19ON5zvV47s9kbzbqymlfD2sLkmKKQUhkf+JizND0OnT5iKpS5eCh/72ihSYPadjSXcYNK1cAoAHmlfigeaVsOm/EclHOmpMnJ5hmKGncWP44w8IDoZvv4USJayuTMQ5GIbBZ5uO02PGJoZ+sZ1LiakA2Gw2hSDJd2oREqe2dy889RT8+KN5vVIlWL4cqlWzsioR5xGblMazi/9g9d6zADSoGIgNhR8pOApC4pQuXTLXCJs61RwU7ekJY8fCc8+Br6/V1Yk4h50nLzF8/g5OXUzC3dXGuE71GNimmlqBpEApCIlTcThgzhwYNw7Onze3de8O774L1atbWpqI0zAMg1kbjvHm6gOk2Q0ql/Jm2oPNaPK/sUEiBUlBSJxGWpq5Xtj8+eb1OnVgyhTo0MHSskSc0pZjF0izG3RqWI437m9MoLe71SWJk1IQEqeQnAy9epmLp7q5mTNGjxiho8JECpJhGJkDoN/u2YQ79kbSq3lldYWJpRSEpNi7fBm6dYOffgIvL1i82DxMXkQKhsNh8PH6oxw6G8+7vZpgs9kI9HGn901VrC5NREFIireLF6FzZ/jtN/Dzg5UroV07q6sScR4xl1MYvWgXPx80B+Xd16wibWuVsbgqkb8oCEmxde6cOf7njz+gVCn47jto0cLqqkScx5ZjF3hywQ7OxaXg6ebCy10bcEvNIKvLEslCQUiKpdhYuO02OHAAypWDtWuhYUOrqxJxDg6HwYyfDzN57SEcBlQv48v0vs2oVz7A6tJEslEQkmLH4YCHHjJDUKVKsG4d1KxpdVUizmPM4l0s3XEagB6hFXmte0N8PfV1I4WTltiQYmfCBHMskKenOUu0QpBIwXogrDK+Hq683bMxk3s3VQiSQk0/nVKsLFpkzhgN8NFHEBZmbT0izsDuMDh0Lj6z66tVjdJsfO52Svhofgop/NQiJMVCYiI8+aQ5VxDA8OHQv7+1NYk4g6i4ZPp9+hs9Z27i6PnLmdsVgqSoUIuQFHnbt/81JgjMEDR5srU1iTiDXw+d5+mFO4lJSMXHw5Vj0QlUL+NndVkiuaIgJEWW3Q5vvAGvvALp6VC+vLmOWMeOVlcmUryl2x2898MhZvx8BMOAuuX8mda3GTXLKgRJ0aMgJEXS0aPw8MOwaZN5/f77zTFBpUtbW5dIcRcZm8RTC3ay5fgFAPq2rML4LvXxcne1uDKR66MgJEWKYZitPk89ZS6d4e8PU6fCI4+AlisSyX8Ltpxky/EL+Hm68XqPRtzbpILVJYncEAUhKTLOn4cnnjAPiQe45Rb4v/+DatWsrErEuTx5e03Oxycz+NYaVAvytbockRumo8akSFi1Cho1MkOQu7s5NujnnxWCRPLb6UtJvLR8D2l2BwDuri5M6tFYIUiKDbUISaGWkABjxsCHH5rX69eHL76A0FBr6xJxBmv3nWPMol3EJqUR6O3OmI51rC5JJM8pCEmhFR4OffrAoUPm9ZEj4fXXwdvb0rJEir3UdAdvfHeA2RuPAdCkUiC9b6pscVUi+UNBSAodwzAHQI8dC6mpUKECfPYZ3Hmn1ZWJFH8nLyQyfP4Odp2KBWDQLSE8e3ddPNw0kkKKJwUhKVRiYmDgQPjmG/N6t24wa5YOixcpCOv/PM/QeTuIT04n0Nuddx5owl31g60uSyRfKQhJofHrr9C3L5w+DR4e8O67MGyYDosXKSiVSvrgcBg0q1KCqX2bUbGE+qGl+FMQkkJh1izz0HiHA2rXhoULoWlTq6sSKf7iktMI8HIHICTIl4WDW1GnnD/uruoKE+egn3Sx3OHD5vpgDoc5W/T27QpBIgVhxa4z3PLGT2w6Ep25rWHFQIUgcSpqERJLORzw2GOQnAx33GEOilZXmEj+Sk6z8+o3+1iwJQKAeb9F0LpGkMVViVhDQUgs9ckn8Msv4ONjXlYIEslfR85fZti8HRw4G4/NBsNuq8nIO2tZXZaIZRSExDInT5qHyIM5P1BIiLX1iBR3y8JP8cKyPSSm2int68H7fZrStlYZq8sSsZSCkFjCMOA//4H4eGjVyhwjJCL557ejMTy9cBcAraqX5oM+TSkb4GVxVSLWUxASS8yfD99+ax4m/+mn4OpqdUUixVvLkFL0CK1I5VI+jLijFq4u6ocWAQUhscDBgzB0qHn5pZfM9cNEJG8ZhsGKXWdoV7sMJXw8sNlsvNurCTYNxBPJQsdISoGKj4f77oO4OGjbFp591uqKRIqfhJR0Rn+1i6e+3MnYxX9gGAaAQpBIDtQiJAXG4YD+/WH/fnP9sK++And3q6sSKV72R8YxbP4Ojp5PwMUGTSuXwDB0RKbIlSgISYF54w1YtswcF7RkCZQrZ3VFIsWHYRgs2HKSV7/ZS0q6g3IBXkx5MJQWIaWsLk2kUFMQkgKxejW8+KJ5edo0uPlma+sRKU7ik9N4ftkevtl1BoDb6pRhcq+mlPL1sLgykcJPQUjy3ZEj8OCD5iHzTzwBjz9udUUixYvdYbDjxEVcXWw807EOj7etjouOChO5JgpCkq+SkszB0Zcuma1AU6ZYXZFI8fD3AdAlfDyY3q8ZdodBWNWSFlcmUrToqDHJVytWwO7dULYsLF4Mnp5WVyRS9MUmpTF03g4Wbj2Zua1p5RIKQSLXQS1Ckq/27DHPu3WDihWtrUWkONh18hLDF+zg5IUkNvwZTadG5Qn01uGXItdLQUjy1f795nm9etbWIVLUGYbB7I3HeeO7/aTZDSqV9GZa32YKQSI3SEFI8lVGEKpb19o6RIqyS4mpjFn0Bz/sPwfA3Q3K8WbPxgpBInnA8jFCM2bMICQkBC8vL8LCwli/fv1V909JSeGFF16gatWqeHp6UqNGDWbPnl1A1UpupKfDn3+al9UiJHJ9klLtdJ22gR/2n8PD1YUJ3Row8yG1BInkFUtbhBYuXMjIkSOZMWMGbdq04aOPPqJTp07s27ePKlWq5HifXr16ce7cOWbNmkXNmjWJiooiPT29gCuXa3H0KKSlgY8PXOHjFJF/4e3hSo/QSizfeZrpfZvRsGKg1SWJFCs2I+MYTAu0bNmSZs2aMXPmzMxt9erVo3v37kyaNCnb/qtXr6ZPnz4cPXqUUqWub7bUuLg4AgMDiY2NJSAg4Lprl3/39dfQvTuEhsKOHVZXI1J0XEhIJSElncqlfABznqCkNDt+nhrNIM4rv76/LesaS01NZfv27XTo0CHL9g4dOrBp06Yc77NixQqaN2/OW2+9RcWKFalduzZjxowhKSnpis+TkpJCXFxclpMUjAMHzHN1i4lcuy3HLtD5g/UM/r/tJKfZAXB1sSkEieQTy36zoqOjsdvtBAcHZ9keHBzM2bNnc7zP0aNH2bBhA15eXixbtozo6GiGDh3KhQsXrjhOaNKkSbz66qt5Xr9cXWIifP+9eVkDpUX+ncNhMPOXI0xeewi7w8DHw5Xz8SmZrUIikj8sHyxt+8eSyIZhZNuWweFwYLPZmDdvHi1atKBz585MnjyZuXPnXrFVaNy4ccTGxmaeTp48meN+kne2b4ewMPjpJ/P6Pxr9ROQfoi+n0H/OFt5ecxC7w+C+0Ip88+QtCkEiBcCyFqGgoCBcXV2ztf5ERUVlayXKUL58eSpWrEhg4F+DBevVq4dhGJw6dYpatWplu4+npyeems64QKSnw5tvwiuvmJfLl4c5c6BlS6srEym8Nh2J5qkvd3I+PgUvdxcmdGvIA2GVrvgPoYjkLctahDw8PAgLC2Pt2rVZtq9du5bWrVvneJ82bdpw5swZLl++nLnt0KFDuLi4UKlSpXytV67u6FFo185cYT49He6/31xao2NHqysTKbwMw2Dy94c4H59CrbJ+rBh+C72aV1YIEilAlnaNjRo1ik8//ZTZs2ezf/9+nn76aSIiIhgyZAhgdms98sgjmfv37duX0qVLM3DgQPbt28evv/7K2LFjefTRR/H29rbqZTg1w4DZs6FJE9i0Cfz94bPPYNEiKF3a6upECjebzcb7fZrSv1VVvh7ehtrB/laXJOJ0LD0MoXfv3sTExDBhwgQiIyNp2LAhq1atomrVqgBERkYSERGRub+fnx9r167lySefpHnz5pQuXZpevXoxceJEq16CUzt/Hp54ApYvN6+3bQuffw7VqllZlUjhtuHPaHadusSw9jUBqFTSh1e7NbS4KhHnZek8QlbQPEI3Li0NPvzQHAt04QK4u8PEiTB6NLi6Wl2dSOGUbnfw/g9/Mv3nwxgGzHusJW1qBlldlkiRkV/f35qYQq6ZYcCKFfDMM3DokLmtUSOzFahpU0tLEynUzsYmM+LLcLYcuwBA35ZVCKta0uKqRAQUhOQapKWZs0S//z5s3GhuK1sWJkyAQYPATT9FIle07mAUo7/axYWEVPw83Xi9RyPubVLB6rJE5H/0FSZXdPYsfPwxfPQRnDljbvPygqefhueeA/UsilzdBz/8yXs/mM2nDSoEMK1vM0KCfC2uSkT+TkFIsjAMs9Vn+nRYssRsDQKzBejxx2HIENBMBSLXplqQOSHiI62q8nzneni5axCdSGGjICSZ7Hbo2hW+++6vba1awfDh5rxAmpdS5N/FJqUR6O0OQLemFQkJ8qVxpRLWFiUiV2T5EhtSeHz0kRmCPD3NsT87dphzA/XtqxAk8m9S0x28tnIfHd77hejLKZnbFYJECje1CAlgzgn0wgvm5XffhWHDrK1HpCg5eSGR4QvC2XXyEgA/7DtHnxZVrC1KRK6JgpAA8PzzcOmSeRj8/yb2FpFrsHpPJGMX/0F8cjqB3u6880AT7qqf83qJIlL4KAgJW7bArFnm5WnTNCmiyLVISbfz+rf7+WzzCQBCq5Rg6oOhVCqpFeNFihIFIeG118yjxR55BNq0sboakaJh2k+HM0PQ4FurM6ZjHdxdNexSpKhREBK2bTPPhw61tg6RouSJW6uz4XA0T95ek9vrqitMpKjSvy9OLibGnDgRoEEDa2sRKcyS0+zM+/0EGcsz+nu5s/Q/rRWCRIo4tQg5ub17zfNq1cDPz9JSRAqtI+cvM2zeDg6cjSc13cHANiEA2Gw2iysTkRulIOTk9uwxz9UaJJKz5eGneX7ZbhJT7QT5eVCzrP5jEClOFIScXEaLUMOG1tYhUtgkpdp5ZcVeFm47CUCr6qX5oE9TygZ4WVyZiOQlBSEnl9EipCAk8pc/z8UzbP4ODp27jM0GI26vxYg7auHqoq4wkeJGQciJGYa6xkRyEpecxpHzCZTx9+SD3k1pXTPI6pJEJJ8oCDkpwzCX0rhwAVxcoG5dqysSsZZhGJmDn8OqlmJKn1BahJSijL8W2hMpzq778Pno6Gi2bdvG9u3biYmJycuaJJ9dvAjdu8PYseb1YcPA29vSkkQsdeBsHF2nbeDg2fjMbfc0Lq8QJOIEch2E9u7dy6233kpwcDAtW7akRYsWlC1blttvv52DBw/mR42Sh7ZsgdBQWLECPDxgxgz44AOrqxKxhmEYLNgSQbdpG9lzOo7XVu6zuiQRKWC56ho7e/Ys7dq1o0yZMkyePJm6detiGAb79u3jk08+oW3btuzZs4eyZcvmV71ynQwDpk6FMWMgLQ2qV4dFi6BZM6srE7FGfHIazy/bwze7zgBwW50yTO7V1NqiRKTA2YyMaVKvwbPPPssPP/zAxo0b8fLKeghpUlISt9xyCx06dGDSpEl5XmheiYuLIzAwkNjYWAICAqwup0DExsKgQbBkiXm9Rw+YPRsCA62tS8Qqe07HMnz+Do7HJOLqYmNsxzo80bY6LjoqTKTQyq/v71x1ja1du5Znn302WwgC8Pb2ZuzYsaxZsybPipMbl54Ot95qhiB3d7MbbPFihSBxXjtPXqLHjE0cj0mkQqAXXw2+mSHtaigEiTipXHWNHT16lGZX6Utp3rw5R48eveGiJO989x388QeUKAFr1kCLFlZXJGKtRhUDCa1SAn8vN955oAklfDysLklELJSrIBQfH3/V5ih/f38uX758w0VJ3vn4Y/N80CCFIHFe+87EUb2ML17urri62Pi0f3P8PN20VpiI5H4eofj4+By7xsDsv8vFkCPJZ6dOwapV5uXHH7e2FhErGIbBnI3HmfTdfh5sUYUJ3cwp1P293C2uTEQKi1wFIcMwqF279lVv139Yhcfs2eBwQLt2UKeO1dWIFKzYxDTGLt7F9/vOAXA+PoV0uwM31+uePk1EiqFcBaF169blVx2Sx+x2+PRT87Jag8TZ7Ii4yJPzwzl9KQkPVxde7FKPh2+uqn/URCSbXAWhdu3a5VcdksfWrIGTJ6FkSbj/fqurESkYDofBpxuO8tbqg6Q7DKqW9mF632Y0rKjDJEUkZ9e11tjp06dZsmQJhw4dwmazUbt2bXr06EHFihXzuj65DoYB771nXu7fH64wpEuk2Dl/OYVpPx0m3WHQpXF5JvVopPFAInJVuZpQEWDGjBmMGjWK1NRUAgMDMQyDuLg4PDw8mDx5MkOHDs2vWvOEM0yoOGsWPPYYeHrC7t1Qq5bVFYkUnO/3nuX85RT6tqiirjCRYqRQTKj47bffMmLECIYPH87p06e5ePEily5d4vTp0wwdOpSnnnqKVRmHKYkljh+HkSPNyxMnKgRJ8eZwGExfd5h1B6Myt3VoUI5+LTUeSESuTa5ahNq1a0fbtm2ZOHFijre/+OKLrF+/nl9++SXPCsxrxblFyOGAO+6An3+GW24xz11dra5KJH9EX07h6YU7Wf9nNCV93Plp9G2U9NXkiCLFVaFoEQoPD+fhhx++4u0PP/wwO3bsuOGi5PpMnmyGH19fmDtXIUiKr81HYuj8wXrW/xmNl7sL4zrVo4SPxgKJSO7larC0w+HA3f3Kf2zc3d01oaJFfvoJnn3WvDx5MtSoYW09IvnB7jCY9tNhPvjxEA4DapX1Y3q/ZtQO9re6NBEponLVItSgQQO+/vrrK96+fPlyGjRocMNFSe6cOAG9epldYwMGaN4gKZ6S0+w8POt33vvBDEEPhFXi6+FtFIJE5IbkqkVo6NCh/Oc//8HT05MnnngCNzfz7unp6Xz00Ue8+OKLzJgxI18KlZwlJcF990FMDISFwcyZoDGiUhx5ubtSqaQ3Ph6uTOzekB7NKlldkogUA7k+fH7MmDFMnjwZf39/avyv/+XIkSNcvnyZESNG8F7GBDaFVHEaLG0YZgvQ559DUBBs3w5VqlhdlUjeSbc7SEyzE/C/uYCSUu1ExiZRvYyfxZWJSEHLr+/vXAchgN9//50FCxZw6NAhAGrXrk2fPn24+eab86yw/FKcgtDUqTBihDkoeu1aaN/e6opE8s7Z2GRGfBmOl7srcwfchIuLmjpFnFl+fX/nqmssMTGRsWPHsnz5ctLS0rjjjjuYOnUqQUFBeVaQXJtffoGnnzYvv/22QpAULz8fjGLUV7u4kJCKr4crh89f1lggEckXuRos/fLLLzN37lzuueceHnzwQX744Qf+85//5FdtcgWnTpmDo+126Nv3rwkURYq6NLuDN747wIA5W7mQkEqDCgGsHNFWIUhE8k2uWoSWLl3KrFmz6NOnDwD9+vWjTZs22O12XDVpTYG4fNlcRDUqCpo0gU8+0eBoKR7OXEriyQXhbD9xEYBHWlXl+c718HLX3xYRyT+5ahE6efIkbdu2zbzeokUL3NzcOHPmTJ4XJlnZ7TB7NtSuDVu2mKvKL10KPj5WVyZy4wzDYOi8HWw/cRF/Tzdm9GvGhG4NFYJEJN/lKgjZ7XY8PLJOYe/m5kZ6enqeFiVZ/fADNGsGgwZBZCRUrw4rVpjnIsWBzWZjYveG3FStJN+OaEvnRuWtLklEnESuusYMw2DAgAF4enpmbktOTmbIkCH4+vpmblu6dGneVejE9u2DsWMhYx3bEiXgxRdh+HBzZXmRouzkhUR2n47NDD0NKwby1eBWWixVRApUroJQ//79s2176KGH8qwYMUVFwcsvm+N/7HZwc4OhQ2H8eChd2urqRG7c6j1neWbxLpLTHFQu6UOjSoEACkEiUuByFYTmzJmTX3UI5izR778PkyZBfLy5rXt3ePNNc2yQSFGXkm5n0qoDzN10HIDQKiUo6avFUkXEOrkKQpI/HA748ksYNw4iIsxtYWHw7rvQrp21tYnklRMxCQyfH87u07EADL61OmM61sHdNVdDFUVE8pSCkMXCw2HwYNi61bxeqZLZItS3L7jo+0GKiW//iOS5JX8Qn5JOSR933u3VhNvrBltdloiIgpCVzpyBu+4yF0z18zNbhJ5+Gry9ra5MJG+duJBAfEo6N1UryZQHQykfqB9yESkcFIQs4nBA//5mCGraFFavhmD9gyzFiGEYmYOfh9xagyA/T3qEVsRNXWEiUojoL5JF3nvPnB/I2xsWLFAIkuJlefhp7puxicRUc44xFxcbvZpXVggSkUJHf5UsEB5udoOBeZRY3bqWliOSZ5JS7Ty7+A9GLtzJzpOX+L/NJ6wuSUTkqtQ1VsASE82B0GlpcN998PjjVlckkjcOR8UzbF44B8/FY7PBiNtr8VhbTX8uIoWbglABGzMGDhyAChW0YKoUH4u3n+Kl5XtISrNTxt+TD3o3pXXNIKvLEhH5VwpCBSgpCT791Lz82WeaJVqKh49/PcLrqw4AcEvNIN7r3ZQy/loDRkSKBo0RKkA7dphdYuXKwR13WF2NSN64t0lFyvh7MqZDbT57tIVCkIgUKWoRKkC//26e33yzusSk6DIMgx0RlwirWhKAcoFerBtzG36e+nMiIkWPWoQK0G+/mectW1pbh8j1upySzsiFO7l/5iZW74nM3K4QJCJFlf56FaCMIHTzzdbWIXI99p6JZfj8cI5FJ+DqYuNsbLLVJYmI3DAFoQJy5gycPGmuH9a8udXViFw7wzD44vcIXlu5j9R0BxUCvZjaN5SwqqWsLk1E5IYpCBWQDRvM8wYNzHXFRIqCuOQ0xi3Zzbe7zW6wO+uV5e2eTSjp62FxZSIieUNBqIAsWmSed+pkbR0iubHl6AW+3R2Jm4uN5zrVZdAtIZnrh4mIFAeWD5aeMWMGISEheHl5ERYWxvr166/pfhs3bsTNzY2mTZvmb4F5ID4evv3WvNy7t7W1iOTGnfWDGdOhNov/05rH2lZXCBKRYsfSILRw4UJGjhzJCy+8QHh4OG3btqVTp05ERERc9X6xsbE88sgj3FFEJuP55htzMsWaNSE01OpqRK4sNjGNZxbvyjIQevjttWhauYR1RYmI5CNLg9DkyZMZNGgQjz32GPXq1eP999+ncuXKzJw586r3Gzx4MH379qVVq1YFVOmNWbjQPO/TR/MHSeEVHnGRzlPW89W2U4xdvMvqckRECoRlQSg1NZXt27fToUOHLNs7dOjApk2brni/OXPmcOTIEV5++eVrep6UlBTi4uKynArSpUuwerV5Wd1iUhgZhsEnvx7lgQ83c/pSElVL+/BMx7pWlyUiUiAsGywdHR2N3W4nODg4y/bg4GDOnj2b433+/PNPnnvuOdavX4+b27WVPmnSJF599dUbrvd6ff01pKZC/frQsKFlZYjk6GJCKmMW7eLHA1EA3NO4PJN6NCLAy93iykRECoblg6X/OfjSMIwcB2Ta7Xb69u3Lq6++Su3ata/58ceNG0dsbGzm6eTJkzdcc25kdIupNUgKm8NR8XSesp4fD0Th4ebCf+9ryLQHQxWCRMSpWNYiFBQUhKura7bWn6ioqGytRADx8fFs27aN8PBwhg8fDoDD4cAwDNzc3Pj++++5/fbbs93P09MTT09rFoFMTIS1a83LCkJS2FQo4Y2fpxvVg3yZ1rcZ9SsEWF2SiEiBsywIeXh4EBYWxtq1a7nvvvsyt69du5Zu3bpl2z8gIIDdu3dn2TZjxgx++uknFi9eTEhISL7XnFsJCZCebl6uU8faWkQALiWmEuDljouLDR8PN2YPuImSvh5aK0xEnJalf/1GjRrFww8/TPPmzWnVqhUff/wxERERDBkyBDC7tU6fPs3nn3+Oi4sLDf8xyKZs2bJ4eXll2y4i2f12NIYRC8J59JYQhrSrAUDlUj4WVyUiYi1Lg1Dv3r2JiYlhwoQJREZG0rBhQ1atWkXVqlUBiIyM/Nc5hUTk6uwOg+nrDvP+D4dwGLA8/DSDbgnB3dXyIYIiIpazGYZhWF1EQYqLiyMwMJDY2FgCAvJ3TMT581C2rHnZud5lKSyi4pN5euFONh6OAaBnWCUmdGuAj4e6wkSkaMmv72/9NRQppjYejuapL3cSfTkFb3dXJnZvyP1hlawuS0SkUFEQEimGzsen8OjcraSkO6gT7M/0fs2oWdbP6rJERAodBaF8FBtrnmtZDSloZfw9GdepLgfPxfNy1wZ4ubtaXZKISKGkIJSPPvvMPL/1VmvrEOfwy6HzlPb1oGHFQAD6t66m1eJFRP6FDhvJJykp8PHH5uX/zf8oki/S7Q7eXH2A/rO3MGz+DuKT04Dss7aLiEh2ahHKJ0uWQFQUVKgAOcwPKZInzlxKYsSCcLaduAjArbXK6LB4EZFcUBDKJ9OmmedDhoC7lm6SfPDTgXOM+moXlxLT8Pd04437G3NP4/JWlyUiUqQoCOWD7dth82YzAD3+uNXVSHGTbnfw1pqDfPzrUQAaVwpk2oPNqFJas0SLiOSWglA+mD7dPH/gAShXztpapPhxsdk4cDYegIFtqvFcp7p4uumoMBGR66EglMdiYmDBAvOyBklLXnI4DFxcbLi42Jjcqwk7Iy5xZ/1gq8sSESnSNKoyj82eDcnJEBoKN99sdTVSHKSk23llxV7GLd2duS3Iz1MhSEQkD6hFKA/Z7TBjhnl5+HBNpCg37kRMAsPnh7P7tDk758OtqmbOEyQiIjdOQSgP/fe/cPw4lCwJDz5odTVS1H37RyTPLfmD+JR0Svi4M7lXE4UgEZE8piCUR1auhJdfNi+/9x54e1tbjxRdyWl2Jn67jy9+iwCgedWSTHkwlAol9EMlIpLXFITywKFD0K+feXnYMOjf39p6pGh7/PNtrP8zGoCht9Vg1F21cdMkiSIi+UJB6AbFx0P37hAXB7fcApMnW12RFHWP3hLC/sg43u3VlHa1y1hdjohIsaYgdAMMAwYMgP37zaU0Fi0CDw+rq5KiJinVzp9R8TSuVAKA9nXK8usz7fHx0K+niEh+U3v7DVi1CpYuNWeQXrJEkydK7h2Oiqf79I089OnvnLyQmLldIUhEpGAoCN2ADRvM8/79NWeQ5N7i7afoOnUjB8/F4+HmSlR8itUliYg4Hf3beQO2bzfPb7rJ2jqkaElMTeel5XtZsuMUAG1qlua93k0p6+9lcWUiIs5HQeg6GcZfQSgszNpapOg4eDaeYfN3cDjqMi42ePrO2gxtXxNXF82+KSJiBQWh63TiBFy4YI4PatjQ6mqkqPhyawSHoy4THODJB31Cubl6aatLEhFxagpC1ymjNahRI/D0tLYWKTqevbsuAMPb16S0n35wRESspsHS10ndYnIt9p6J5ZnFu7A7DAC83F15uWsDhSARkUJCLULXads287xZM2vrkMLJMAy++D2C11buIzXdQfUyfgxpV8PqskRE5B8UhK5Dejps3mxebt3a2lqk8IlLTmPc0t18+0ckAHfULUvv5pUtrkpERHKiIHQddu6Ey5ehRAkNlJasdp+KZdj8HURcSMTNxcZzneoy6JYQbDYdFSYiUhgpCF2H9evN8zZtwEWjrOR/vt55mrGL/iDV7qBiCW+m9Q0ltEpJq8sSEZGrUBC6DhlBqG1ba+uQwqVuuQBcXKBj3WDeur8JgT7uVpckIiL/QkEolwzjr6U1FIQk+nIKQf87AqxOOX++GX4LNcv6qStMRKSIUMdOLh08COfPg5cXNG9udTViFYfD4JNfj3LLmz+x/cTFzO21gv0VgkREihC1COVSRrdYy5bg4WFtLWKNiwmpjF60i58ORAGw8o8zhFXVWCARkaJIQSiXfv3VPFe3mHPadvwCTy4IJzI2GQ83F8Z3qU+/llWsLktERK6TglAuaaC0c3I4DD789Qjvfn8Iu8MgJMiXaX1DaVAh0OrSRETkBigI5cLJk+Ziqy4u0KqV1dVIQfp+31neWn0QgG5NK/Df+xrh56lfHxGRok5/yXMhozUoNBT8/a2tRQpWxwbl6Na0Aq2ql6b3TZU1IFpEpJhQEMqFjIVW1RpU/NkdBp9vPk7PsEr4e7ljs9n4oE+o1WWJiEgeUxDKhQMHzPMGDaytQ/JXVHwyTy/cycbDMeyIuMSUPk3VAiQiUkwpCOXC/v3meb161tYh+Wfj4Wie+nIn0ZdT8HZ3pV3tMgpBIiLFmILQNUpKguPHzct161paiuQDu8Pggx//ZOpPf2IYUCfYn2l9Q6kVrMFgIiLFmYLQNTp0yFxeo2RJKFvW6mokL0XFJ/Pk/HB+P3YBgD43Veblrg3w9nC1uDIREclvCkLX6O/dYuopKV5cbDaORSfg6+HK6z0a0a1pRatLEhGRAqIgdI00Pqh4cTgMXFzMRBvk58nMh8Io6eNO9TJ+FlcmIiIFSYuuXqOMI8Y0PqjoO3MpiV4fbebrnaczt4VVLakQJCLihNQidI3++MM8V4tQ0fbj/nOMXrSLS4lpnLqYxN0Ny+HpprFAIiLOSkHoGvzxh9ki5O4ON99sdTVyPVLTHby95gCfrD8GQKOKgUzrG6oQJCLi5BSErsHcueZ5165QurSlpch1OHkhkScXhLPz5CUABrSuxrjOdRWCREREQejfpKXBF1+YlwcOtLYWyb2LCal0nbaBS4lpBHi58VbPJtzdsJzVZYmISCGhIPQvVq2C8+chOBjuvtvqaiS3Svp60Lt5ZX47doFpD4ZSuZSP1SWJiEghoiD0LzK6xR5+GNz0bhUJETGJuLraqFjCG4AxHetgGODhpoMkRUQkK30zXEVUFKxcaV4eMMDSUuQardodyT1T1jN8/g7S7A4A3F1dFIJERCRHauO4ivnzIT0dbrpJK84XdslpdiZ+u48vfosAzNmi45PTKeXrYXFlIiJSmCkIXYFhwJw55mW1BhVux6ITGDZvB/si4wD4z201GHVXbdxd1QokIiJXpyB0BTt3mvMHeXhAnz5WVyNX8vXO0zy/dDcJqXZK+XowuVcTbqujVXFFROTaKAhdQUZrUPfuUKqUpaXIFaTbHXz0y1ESUu20CCnFlD6hlAv0srosEREpQhSEcpCUBPPmmZfVLVZ4ubm6ML1fM77eeZrh7Wvipq4wERHJJX1z5OD//g8uXICqVeGuu6yuRv5uyfZTzPz5SOb1kCBfRt5ZWyFIRESui1qE/sHhgPfeMy8/9ZTmDiosElPTGf/1XhZvP4XNBq1rlKZJ5RJWlyUiIkWcvub/4bvvzAVWAwJg0CCrqxGAg2fjGTZ/B4ejLuNig5F31qZhxUCryxIRkWJAQegf3n3XPH/8cTMMiXUMw+CrbSd5ecVektMclPX35IM+obSqoZVvRUQkbygI/U14OKxbB66uMGKE1dXI88t2s2DLSQBurV2Gyb2aEOTnaXFVIiJSnGiE6d9MnmyeP/AAVKlibS0CTSqVwNXFxjN312HugJsUgkREJM+pReh/Tp+GL780L48ebW0tzsowDKIvp1LG3ww8vW+qTPNqpahZ1s/iykREpLiyvEVoxowZhISE4OXlRVhYGOvXr7/ivkuXLuWuu+6iTJkyBAQE0KpVK9asWZMndUydaq4r1rYtNG+eJw8puRCfnMbwBeF0n76R2MQ0AGw2m0KQiIjkK0uD0MKFCxk5ciQvvPAC4eHhtG3blk6dOhEREZHj/r/++it33XUXq1atYvv27bRv356uXbsSHh5+Q3UkJMBHH5mX1RpU8HafiqXL1A18+0ck5+KS2XL8gtUliYiIk7AZhmFY9eQtW7akWbNmzJw5M3NbvXr16N69O5MmTbqmx2jQoAG9e/dm/Pjx17R/XFwcgYGBxMbGEvC/w8Lmz4d+/aBaNTh82BwsLfnPMAw+23Sc11cdINXuoGIJb6b2DaVZlZJWlyYiIoVMTt/fecGyMUKpqals376d5557Lsv2Dh06sGnTpmt6DIfDQXx8PKWushhYSkoKKSkpmdfj4uKy7ZOxnMbDDysEFZTYxDSeWbKLNXvPAdChfjBv92xCoI+7xZWJiIgzsaxrLDo6GrvdTnBwcJbtwcHBnD179poe49133yUhIYFevXpdcZ9JkyYRGBiYeapcuXKW26OiIGOYUb9+uXsNcv3eWnOANXvP4e5q4+Wu9fno4TCFIBERKXCWD5a22WxZrhuGkW1bThYsWMArr7zCwoULKVu27BX3GzduHLGxsZmnkydPZrn9q6/AbjcHSNepc32vQXJvbMc6tKpemiX/ac3ANiHX9JmLiIjkNcu6xoKCgnB1dc3W+hMVFZWtleifFi5cyKBBg1i0aBF33nnnVff19PTE0/PK88988YV5/tBD11a3XJ9Liaks2XGaR9tUw2azUcLHgwVP3Gx1WSIi4uQsaxHy8PAgLCyMtWvXZtm+du1aWrdufcX7LViwgAEDBjB//nzuueeeG6rhzz/h99/BxQV6976hh5Kr2H7iAp0/WM9rK/excOvJf7+DiIhIAbF0QsVRo0bx8MMP07x5c1q1asXHH39MREQEQ4YMAcxurdOnT/P5558DZgh65JFH+OCDD7j55pszW5O8vb0JDMz9Ipzz55vnd90F5crlzWuSvzgcBh/9epR3vj+I3WEQEuRLo0paLFVERAoPS4NQ7969iYmJYcKECURGRtKwYUNWrVpF1apVAYiMjMwyp9BHH31Eeno6w4YNY9iwYZnb+/fvz9y5c3P13IahbrH8FHM5hdGLdvHzwfMA3NukAq/3aISfpyYzFxGRwsPSeYSskDEPwY8/xnLHHQH4+MC5c+CnCYzzzNbjFxg+fwfn4lLwdHPh1Xsb0PumyhoQLSIi163YzSNkta++Ms+7d1cIymtpdgdR8SnUKOPL9H7NqFsu735gRURE8pLTBqHFi81zzR2UN+wOA1cXs8WndY0gPnwojFtqBuGrrjARESnELJ9HyCoxMRAUZA6Ulhuz8XA0d07+hWPRCZnbOjYopxAkIiKFntMGIYDGjcFdkxlfN7vDYPLaQzw063eORSfw3tpDVpckIiKSK079L3vFilZXUHSdi0vmqS/D+e2ouVJ87+aVeeXeBhZXJSIikjsKQpJrvxw6z6iFO4lJSMXHw5XX72tE91C9mSIiUvQ4dRCqUMHqCoqedQejGDhnKwD1ygcwvW8o1cvosDsRESmanDoIqUUo926pGURolRLULx/AS13q4+XuanVJIiIi182pg5BahK7Nb0djCKtaEndXF9xdXVjw+M0KQCIiUiw49VFjahG6ujS7g0mr9tPn49945/uDmdsVgkREpLhw6hYhLbR6ZacuJvLkgnDCIy4BkJruwDAMLZMhIiLFitMGoTJlNIfQlazZe5axi3YRl5yOv5cbb/dszN0Ny1tdloiISJ5z2iBUXt/r2aSmO5j03X7mbDwOQJPKJZj2YCiVS/lYW5iIiEg+URCSTJGxSSzcehKAx24J4Zm76+Lh5tTDyEREpJhz2iCkI8ayq1ral7d7NsHTzYU76wdbXY6IiEi+c9ogpIHSkJxmZ9Kq/dzdsDytapQG4J7GaioTERHn4bT9Hn5OPhnysegE7p+5ic82n2DkwnCS0+xWlyQiIlLgnLZFyJmt2HWGcUv+ICHVTilfD964v7HmBhIREaekIOREktPsvPrNPhZsiQCgRbVSTHkwlHKBXhZXJiIiYg0FIScRm5RG7482c+BsPDYbDG9fk6fuqIWbq9P2joqIiCgIOYsALzdqBfsTfTmF93o3pW2tMlaXJCIiYjkFoWIsMTWddIdBgJc7NpuN1+9rSFKqnbIB6goTEREBJz5qrLg7dC6ebtM2MuarXRiGAYC/l7tCkIiIyN+oRaiYMQyDRdtOMX7FHpLTHMQmpREZm0yFEt5WlyYiIlLoKAgVIwkp6bywbDfLd54BoG2tIN7r3ZQgP0+LKxMRESmcFISKiX1n4hg+fwdHoxNwdbEx6q7a/KddDVxcbFaXJiIiUmgpCBUDdoeRGYLKBXgxtW8oN1UrZXVZIiIihZ4GSxcDri423n6gMXfVD2bVU20VgkRERK6RWoSKqD2nYzkRk5i5SGpY1VJ88ogCkIiISG4oCBUxhmHw+eYT/Pfb/bi4QK1gP2oH+1tdloiISJGkIFSExCal8eziP1i99ywAd9YOpqy/jggTERG5XgpCRcTOk5cYPn8Hpy4m4e5qY1ynegxsUw2bTUeFiYiIXC8FoSJg9oZjTPpuP2l2g8qlvJn2YDOaVC5hdVkiIiJFnoJQEXApKY00u0GnhuV44/7GBHq7W12SiIhIsaAgVEil2x24uZqzGzx1Ry3qlvOnU8Ny6goTkWLLMAzS09Ox2+1WlyIWcXd3x9XVtUCfU0GokHE4DD5ef5TVe86ycPDNeLq54upio3Oj8laXJiKSb1JTU4mMjCQxMdHqUsRCNpuNSpUq4efnV2DPqSBUiMRcTmH0ol38fPA8ACt2nuGB5pUtrkpEJH85HA6OHTuGq6srFSpUwMPDQ63fTsgwDM6fP8+pU6eoVatWgbUMKQgVEr8fjWHEl+Gci0vB082FV+5tQM+wSlaXJSKS71JTU3E4HFSuXBkfHx+ryxELlSlThuPHj5OWlqYg5CwcDoMZPx9m8tpDOAyoXsaX6X2bUa98gNWliYgUKBcXrfrk7KxoCVQQstik7/bzyfpjAPQIrchr3Rvi66mPRUREpCDoG9dij7SqxopdZxjdoQ4PhFVSv7iIiEgBUjtkAbM7DDb8GZ15vXIpH34Z255ezSsrBImIFFGbNm3C1dWVu+++O9ttP//8MzabjUuXLmW7rWnTprzyyitZtoWHh/PAAw8QHByMl5cXtWvX5vHHH+fQoUP5VL1pxowZhISE4OXlRVhYGOvXr//X+8ybN48mTZrg4+ND+fLlGThwIDExMZm37927l/vvv59q1cyVEN5///18fAXXR0GoAEXFJfPQp7/z0KzfWXcwKnO7l3vBzpkgIiJ5a/bs2Tz55JNs2LCBiIiI636clStXcvPNN5OSksK8efPYv38///d//0dgYCAvvfRSHlac1cKFCxk5ciQvvPAC4eHhtG3blk6dOl31tWzYsIFHHnmEQYMGsXfvXhYtWsTWrVt57LHHMvdJTEykevXqvPHGG5QrVy7f6r8R6horIOv/PM/TC3cSfTkVHw9XElLSrS5JRETyQEJCAl999RVbt27l7NmzzJ07l/Hjx+f6cRITExk4cCCdO3dm2bJlmdtDQkJo2bJlji1KeWXy5MkMGjQoM8S8//77rFmzhpkzZzJp0qQc7/Pbb79RrVo1RowYkVnn4MGDeeuttzL3uemmm7jpppsAeO655/Kt/huhFqF8lm538M6agzwyewvRl1OpW86fFcNvoUvjClaXJiJSaBkGJCRYczKM3NW6cOFC6tSpQ506dXjooYeYM2cORm4fBFizZg3R0dE888wzOd5eokSJK953yJAh+Pn5XfV0pdad1NRUtm/fTocOHbJs79ChA5s2bbric7Zu3ZpTp06xatUqDMPg3LlzLF68mHvuueffX2whohahfBQZm8RTC3ay5fgFAPq2rML4LvXVFSYi8i8SE6EAJxfO4vJl8PW99v1nzZrFQw89BMDdd9/N5cuX+fHHH7nzzjtz9bx//vknAHXr1s3V/QAmTJjAmDFjrrpPhQo5/wMeHR2N3W4nODg4y/bg4GDOnj17xcdr3bo18+bNo3fv3iQnJ5Oens69997L1KlTc12/lRSE8tGWYxfYcvwCfp5uTOrRiK5N1AokIlKcHDx4kC1btrB06VIA3Nzc6N27N7Nnz851ELqeVqQMZcuWpWzZstd9f8g+h49hGFc9iGffvn2MGDGC8ePH07FjRyIjIxk7dixDhgxh1qxZN1RLQVIQykfdmlbk1MUk7mlUnmpBufj3QkTEyfn4mC0zVj33tZo1axbp6elUrFgxc5thGLi7u3Px4kVKlixJQIA5QW5sbGy27q1Lly4RGBgIQO3atQE4cOAArVq1ylXNQ4YM4YsvvrjqPvv27aNKlSrZtgcFBeHq6pqt9ScqKipbK9HfTZo0iTZt2jB27FgAGjdujK+vL23btmXixImUL1801shUEMpDpy8l8do3+/jvfQ0p7ecJwLD2NS2uSkSk6LHZctc9ZYX09HQ+//xz3n333Wzja+6//37mzZvH8OHDqVWrFi4uLmzdupWqVatm7hMZGcnp06epU6cOYI7JCQoK4q233soyWDrDpUuXrjhO6Ea6xjw8PAgLC2Pt2rXcd999mdvXrl1Lt27drvh4iYmJuLlljREZy2LcSOtWQVMQyiNr951jzKJdxCal4epqY3rfZlaXJCIi+WjlypVcvHiRQYMGZbbqZOjZsyezZs1i+PDh+Pv7M3jwYEaPHo2bmxtNmjThzJkzvPDCC9SrVy8zRPn6+vLpp5/ywAMPcO+99zJixAhq1qxJdHQ0X331FREREXz55Zc51nKjXWOjRo3i4Ycfpnnz5rRq1YqPP/6YiIgIhgwZkrnPuHHjOH36NJ9//jkAXbt25fHHH2fmzJmZXWMjR46kRYsWmaErNTWVffv2ZV4+ffo0O3fuxM/Pj5o1C0lDgeFkYmNjDcD4739j8+TxUtLsxqsr9hpVn11pVH12pXHv1PVGRExCnjy2iIgzSEpKMvbt22ckJSVZXUqudOnSxejcuXOOt23fvt0AjO3btxuGYRjJycnGhAkTjHr16hne3t5G1apVjQEDBhiRkZHZ7rt161ajR48eRpkyZQxPT0+jZs2axhNPPGH8+eef+fp6pk+fblStWtXw8PAwmjVrZvzyyy9Zbu/fv7/Rrl27LNumTJli1K9f3/D29jbKly9v9OvXzzh16lTm7ceOHTOAbKd/Pk6Gq/0sZHx/x8bmzfd3BpthFKH2qzwQFxdHYGAg//1vLM8/f2MLm568kMjw+TvYdSoWgEG3hPDs3XXxcNOsBCIi1yo5OZljx45lzmoszutqPwsZ39+xsbGZ467ygrrGrtP2ExcZMGcL8cnpBHq7884DTbir/pUHlYmIiEjhoyB0nWoH+1HK14NaZf2Y8mAolUrm4jADERERKRQUhHLhbGwywQGe2Gw2/L3cmfdYS4IDvHB3VVeYiIhIUaRv8Gv0za4z3Dn5Fz7ffCJzW6WSPgpBIiIiRZi+xf9FcpqdcUt38+SCcC6npLN237kiNT+CiIiIXJm6xq7iyPnLDJu3gwNn47HZYNhtNRl5Z62rTjkuIiLXR/9kihU/AwpCV7As/BQvLNtDYqqdID8P3uvdlLa1ylhdlohIsePu7g6YMxV7e3tbXI1YKTU1FfhrhuqCoCCUg2PRCYxZ9Ad2h0Gr6qX5oE9TygZobgsRkfzg6upKiRIliIqKAsDHx0ct707I4XBw/vx5fHx8si3dkZ8UhHIQEuTLMx3rkJRm58nba+Hqol9IEZH8VK5cOYDMMCTOycXFhSpVqhRoEFYQwuyTXLz9FE0ql6B2sD8Ag9vVsLgqERHnYbPZKF++PGXLliUtLc3qcsQiHh4euLgU7HFclgehGTNm8PbbbxMZGUmDBg14//33adu27RX3/+WXXxg1ahR79+6lQoUKPPPMM1kWhcuthJR0Xlq+h6Xhp6lV1o8Vw2/B26Pg+iZFROQvrq6uBTo+RMTSw+cXLlzIyJEjeeGFFwgPD6dt27Z06tSJiIiIHPc/duwYnTt3pm3btoSHh/P8888zYsQIlixZcl3Pvz8yjq7TNrA0/DQuNugeWhFPrRMmIiLiNCxddLVly5Y0a9aMmTNnZm6rV68e3bt3Z9KkSdn2f/bZZ1mxYgX79+/P3DZkyBB27drF5s2br+k5MxZt6/fiHrYYEaSmOygX4MWUB0NpEVLqxl+UiIiI5Ln8WnTVsuaP1NRUtm/fTocOHbJs79ChA5s2bcrxPps3b862f8eOHdm2bVuu+5R/vbyP1HQHt9Upw6qn2ioEiYiIOCHLxghFR0djt9sJDs66YntwcDBnz57N8T5nz57Ncf/09HSio6MpX758tvukpKSQkpKSeT02NhYAIyWJkXfVYkDrEFzsycTFJd/oSxIREZF8EhcXB+T9pIuWD5b+5yFyhmFc9bC5nPbPaXuGSZMm8eqrr2bbfmpmf56eCU/ntmARERGxTExMDIGBgXn2eJYFoaCgIFxdXbO1/kRFRWVr9clQrly5HPd3c3OjdOnSOd5n3LhxjBo1KvP6pUuXqFq1KhEREXn6Rsr1iYuLo3Llypw8eTJP+3wl9/RZFB76LAoPfRaFR2xsLFWqVKFUqbwdymJZEPLw8CAsLIy1a9dy3333ZW5fu3Yt3bp1y/E+rVq14ptvvsmy7fvvv6d58+aZU7T/k6enJ56entm2BwYG6oe6EAkICNDnUUjosyg89FkUHvosCo+8nmfI0mPFR40axaeffsrs2bPZv38/Tz/9NBEREZnzAo0bN45HHnkkc/8hQ4Zw4sQJRo0axf79+5k9ezazZs1izJgxVr0EERERKcIsHSPUu3dvYmJimDBhApGRkTRs2JBVq1ZRtWpVACIjI7PMKRQSEsKqVat4+umnmT59OhUqVGDKlCncf//9Vr0EERERKcIsHyw9dOhQhg4dmuNtc+fOzbatXbt27Nix47qfz9PTk5dffjnH7jIpePo8Cg99FoWHPovCQ59F4ZFfn4WlEyqKiIiIWEnrSYiIiIjTUhASERERp6UgJCIiIk5LQUhEREScVrEMQjNmzCAkJAQvLy/CwsJYv379Vff/5ZdfCAsLw8vLi+rVq/Phhx8WUKXFX24+i6VLl3LXXXdRpkwZAgICaNWqFWvWrCnAaou/3P5uZNi4cSNubm40bdo0fwt0Irn9LFJSUnjhhReoWrUqnp6e1KhRg9mzZxdQtcVbbj+LefPm0aRJE3x8fChfvjwDBw4kJiamgKotvn799Ve6du1KhQoVsNlsLF++/F/vkyff30Yx8+WXXxru7u7GJ598Yuzbt8946qmnDF9fX+PEiRM57n/06FHDx8fHeOqpp4x9+/YZn3zyieHu7m4sXry4gCsvfnL7WTz11FPGm2++aWzZssU4dOiQMW7cOMPd3d3YsWNHAVdePOX288hw6dIlo3r16kaHDh2MJk2aFEyxxdz1fBb33nuv0bJlS2Pt2rXGsWPHjN9//93YuHFjAVZdPOX2s1i/fr3h4uJifPDBB8bRo0eN9evXGw0aNDC6d+9ewJUXP6tWrTJeeOEFY8mSJQZgLFu27Kr759X3d7ELQi1atDCGDBmSZVvdunWN5557Lsf9n3nmGaNu3bpZtg0ePNi4+eab861GZ5HbzyIn9evXN1599dW8Ls0pXe/n0bt3b+PFF180Xn75ZQWhPJLbz+K7774zAgMDjZiYmIIoz6nk9rN4++23jerVq2fZNmXKFKNSpUr5VqMzupYglFff38Wqayw1NZXt27fToUOHLNs7dOjApk2bcrzP5s2bs+3fsWNHtm3bRlpaWr7VWtxdz2fxTw6Hg/j4+DxfYM8ZXe/nMWfOHI4cOcLLL7+c3yU6jev5LFasWEHz5s156623qFixIrVr12bMmDEkJSUVRMnF1vV8Fq1bt+bUqVOsWrUKwzA4d+4cixcv5p577imIkuVv8ur72/KZpfNSdHQ0drs92+r1wcHB2Vatz3D27Nkc909PTyc6Opry5cvnW73F2fV8Fv/07rvvkpCQQK9evfKjRKdyPZ/Hn3/+yXPPPcf69etxcytWfyosdT2fxdGjR9mwYQNeXl4sW7aM6Ohohg4dyoULFzRO6AZcz2fRunVr5s2bR+/evUlOTiY9PZ17772XqVOnFkTJ8jd59f1drFqEMthstizXDcPItu3f9s9pu+Rebj+LDAsWLOCVV15h4cKFlC1bNr/KczrX+nnY7Xb69u3Lq6++Su3atQuqPKeSm98Nh8OBzWZj3rx5tGjRgs6dOzN58mTmzp2rVqE8kJvPYt++fYwYMYLx48ezfft2Vq9ezbFjxzIXC5eClRff38Xq37ygoCBcXV2zJfmoqKhsqTFDuXLlctzfzc2N0qVL51utxd31fBYZFi5cyKBBg1i0aBF33nlnfpbpNHL7ecTHx7Nt2zbCw8MZPnw4YH4ZG4aBm5sb33//PbfffnuB1F7cXM/vRvny5alYsSKBgYGZ2+rVq4dhGJw6dYpatWrla83F1fV8FpMmTaJNmzaMHTsWgMaNG+Pr60vbtm2ZOHGiehEKUF59fxerFiEPDw/CwsJYu3Ztlu1r166ldevWOd6nVatW2fb//vvvad68Oe7u7vlWa3F3PZ8FmC1BAwYMYP78+epzz0O5/TwCAgLYvXs3O3fuzDwNGTKEOnXqsHPnTlq2bFlQpRc71/O70aZNG86cOcPly5cztx06dAgXFxcqVaqUr/UWZ9fzWSQmJuLikvWr09XVFfirNUIKRp59f+dqaHURkHEo5KxZs4x9+/YZI0eONHx9fY3jx48bhmEYzz33nPHwww9n7p9x+N3TTz9t7Nu3z5g1a5YOn88juf0s5s+fb7i5uRnTp083IiMjM0+XLl2y6iUUK7n9PP5JR43lndx+FvHx8UalSpWMnj17Gnv37jV++eUXo1atWsZjjz1m1UsoNnL7WcyZM8dwc3MzZsyYYRw5csTYsGGD0bx5c6NFixZWvYRiIz4+3ggPDzfCw8MNwJg8ebIRHh6eOZVBfn1/F7sgZBiGMX36dKNq1aqGh4eH0axZM+OXX37JvK1///5Gu3btsuz/888/G6GhoYaHh4dRrVo1Y+bMmQVccfGVm8+iXbt2BpDt1L9//4IvvJjK7e/G3ykI5a3cfhb79+837rzzTsPb29uoVKmSMWrUKCMxMbGAqy6ecvtZTJkyxahfv77h7e1tlC9f3ujXr59x6tSpAq66+Fm3bt1VvwPy6/vbZhhqyxMRERHnVKzGCImIiIjkhoKQiIiIOC0FIREREXFaCkIiIiLitBSERERExGkpCImIiIjTUhASERERp6UgJCIiIk5LQUhECo0BAwZgs9mw2Wy4u7tTvXp1xowZQ0JCQuY+n332GS1atMDX1xd/f39uvfVWVq5cmeVxfv7558zH+fvpxRdfzHa7i4sLgYGBhIaG8swzzxAZGVmgr1lErKUgJCKFyt13301kZCRHjx5l4sSJzJgxgzFjxgAwZswYBg8eTK9evdi1axdbtmyhbdu2dOvWjWnTpmV7rIMHDxIZGZl5eu6557LdfubMGbZu3cqzzz7LDz/8QMOGDdm9e3eBvFYRsZ6W2BCRQmPAgAFcunSJ5cuXZ257/PHHWblyJcuWLaNVq1ZMmTKFJ598Msv9Ro8ezdSpUzly5AiVK1fm559/pn379ly8eJESJUpke54r3Z6UlERoaChBQUFs2LAhn16liBQmahESkULN29ubtLQ0FixYgJ+fH4MHD862z+jRo0lLS2PJkiU3/FxDhgxh48aNREVF3dBjiUjRoCAkIoXWli1bmD9/PnfccQeHDh2iRo0aeHh4ZNuvQoUKBAYGcujQoSzbK1WqhJ+fX+YpJibmX5+zbt26ABw/fjxPXoOIFG5uVhcgIvJ3K1euxM/Pj/T0dNLS0ujWrRtTp06lf//+V72fYRjYbLYs29avX4+/v3/m9ZIlS/7r82eMFvjnY4lI8aQgJCKFSvv27Zk5cybu7u5UqFABd3d3AGrXrs2GDRtITU3N1ip05swZ4uLiqFWrVpbtISEhOY4Rupr9+/cDUK1atet+DSJSdKhrTEQKFV9fX2rWrEnVqlUzQxBAnz59uHz5Mh999FG2+7zzzju4u7tz//3339BzJyUl8fHHH3PrrbdSpkyZG3osESka1CIkIkVCq1ateOqppxg7diypqal0796dtLQ0vvjiCz744APef/99KleunKvHjIqKIjk5mfj4eLZv385bb71FdHQ0S5cuzadXISKFjYKQiBQZ77//Po0bN2bmzJm89NJL2Gw2mjVrxvLly+natWuuH69OnTrYbDb8/PyoXr06HTp0YNSoUZQrVy4fqheRwkjzCImIiIjT0hghERERcVoKQiIiIuK0FIRERETEaSkIiYiIiNNSEBIRERGnpSAkIiIiTktBSERERJyWgpCIiIg4LQUhERERcVoKQiIiIuK0FIRERETEaSkIiYiIiNP6fzK55zCHzq83AAAAAElFTkSuQmCC", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "plt.title('Receiver Operating Characteristic')\n", + "plt.plot(roc.POFD, roc.POD, 'b', label = 'AUC = %0.2f' % roc.AUC.item())\n", + "plt.legend(loc = 'lower right')\n", + "plt.plot([0, 1], [0, 1],'--')\n", + "plt.xlim([0, 1])\n", + "plt.ylim([0, 1])\n", + "plt.ylabel('POD')\n", + "plt.xlabel('POFD')\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "- A perfect ROC curve would go from the lower left corner, up to the top left corner, and across to the top right corner.\n", + "- If the curve follows the dotted line, it means the forecast has no discrimination ability. This synthetic forecast clearly does have discrimination ability.\n", + "- The AUC is the area under the curve. 1 is the maximum value that can be achieved. If we picked a random forecast value from out dataset where the event occured, and a random forecast value when the event did not occur, the AUC tells us the probability that the former forecast value will be greater than the latter." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Cautionary notes\n", + "- ROC or AUC is not a suitable overall performance measure for a probabilistic classifier as it is not a proper scoring rule and is ignorant of calibration (Hand 2009; Hand and Anagnostopoulos 2013, 2023). It is however, a useful tool for understanding potential predictive performance subject to recalibration.\n", + "- Non-concave ROC curves like the synthetic example above have several problems (Pesce et al. 2010; Geniting and Vogel 2022). To get around these problems, you can apply isotonic regression to the forecast data before you calculate the ROC curves (Fawcett and Niculescu-Mizil 2007).\n", + "- Results may be misleading if two ROC curves cross (Hand 2009)." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### References\n", + "Fawcett, T. and Niculescu-Mizil, A., 2007. PAV and the ROC convex hull. Machine Learning, 68, pp.97-106.\n", + "\n", + "Gneiting, T. and Vogel, P., 2022. Receiver operating characteristic (ROC) curves: equivalences, beta model, and minimum distance estimation. Machine learning, 111(6), pp.2147-2159.\n", + "\n", + "Hand, D.J., 2009. Measuring classifier performance: a coherent alternative to the area under the ROC curve. Machine learning, 77(1), pp.103-123.\n", + "\n", + "Hand, D.J. and Anagnostopoulos, C., 2013. When is the area under the receiver operating characteristic curve an appropriate measure of classifier performance?. Pattern Recognition Letters, 34(5), pp.492-495.\n", + "\n", + "Hand, D.J. and Anagnostopoulos, C., 2023. Notes on the H-measure of classifier performance. Advances in Data Analysis and Classification, 17(1), pp.109-124.\n", + "\n", + "Pesce, L.L., Metz, C.E. and Berbaum, K.S., 2010. On the convexity of ROC curves estimated from radiological test results. Academic radiology, 17(8), pp.960-968." + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "scoresenv", + "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.11.4" + }, + "orig_nbformat": 4 + }, + "nbformat": 4, + "nbformat_minor": 2 +}