utility_calc_1.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{'EUs': 0    0.5\n",
      "1    0.5\n",
      "dtype: float64, 'optimal': 0}\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import random\n",
    "\n",
    "def decide(probs=None, utils=None):\n",
    "    if probs is None and utils is None:\n",
    "        raise ValueError(\"Either 'probs' or 'utils' must be given\")\n",
    "\n",
    "    if utils is None:\n",
    "        if len(probs.shape) > 1:\n",
    "            raise ValueError(\"'probs' is not a vector, please specify array of utilities\")\n",
    "        utils = np.diag(probs)\n",
    "\n",
    "    elif probs is None:\n",
    "        probs = np.full(utils.shape[1], 1/utils.shape[1])\n",
    "\n",
    "    if len(probs) != utils.shape[1]:\n",
    "        raise ValueError('Mismatch between inference and utility variates')\n",
    "\n",
    "    if not isinstance(utils, pd.DataFrame):\n",
    "        utils = pd.DataFrame(utils)\n",
    "\n",
    "    # Calculate expected utilities (matrix product)\n",
    "    exputils = np.dot(utils, probs)\n",
    "    # If expected utilites are the same for every class, choose random class\n",
    "    # if(len(set(exputils)) == 1):\n",
    "    #     exputils = pd.Series(exputils, index=utils.index).sort_values(ascending=False)\n",
    "    #     return {'EUs': exputils, 'optimal': random.randint(0,len(exputils)-1)}\n",
    "\n",
    "        \n",
    "    exputils = pd.Series(exputils, index=utils.index).sort_values(ascending=False)\n",
    "    #Select one decision with max expected utility\n",
    "    optimal = exputils.idxmax()\n",
    "    return {'EUs': exputils, 'optimal': optimal}\n",
    "\n",
    "    \n",
    "\n",
    "# Example usage\n",
    "probs = np.array([0.5, 0.5])\n",
    "utils = np.array([[1, 0], [-1,2]])\n",
    "result = decide(probs=probs, utils=utils)\n",
    "print(result)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [],
   "source": [
    "def utilityMatrixCreator(minUtil, maxUtil):\n",
    "    U_AA = random.randint(minUtil if minUtil > 0 else 0, maxUtil)\n",
    "    U_BA = random.randint(minUtil-1, U_AA-1)\n",
    "\n",
    "    U_BB = random.randint(minUtil if minUtil > 0 else 0, maxUtil)\n",
    "    U_AB = random.randint(minUtil-1, U_BB-1)\n",
    "    return [[U_AA, U_AB], [U_BA, U_BB]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 7,  2],\n",
       "       [-9,  3]])"
      ]
     },
     "execution_count": 41,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.array(utilityMatrixCreator(-10, 10))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Probabilities from RNopm infer\n",
    "prob_predict_cat = [0.976039626021315, 0.0239603739786851]\n",
    "prob_predict_dog = [0.00717048738691545, 0.992829512613085]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [],
   "source": [
    "CM_df = pd.read_csv('cd_dataframe.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Unnamed: 0</th>\n",
       "      <th>True Cat</th>\n",
       "      <th>True Dog</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>Predicted Cat</td>\n",
       "      <td>1744</td>\n",
       "      <td>38</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>Predicted Dog</td>\n",
       "      <td>6</td>\n",
       "      <td>1712</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "      Unnamed: 0  True Cat  True Dog\n",
       "0  Predicted Cat      1744        38\n",
       "1  Predicted Dog         6      1712"
      ]
     },
     "execution_count": 44,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "CM_df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [],
   "source": [
    "CM = np.array([[1744, 38], [6, 1712]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [],
   "source": [
    "util = np.array(utilityMatrixCreator(-10, 10))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [],
   "source": [
    "expected_RN = np.multiply(CM, util)\n",
    "expected_RN = expected_RN.sum()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [],
   "source": [
    "test_df = pd.read_csv('test_dataframe.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [],
   "source": [
    "cm_rn = np.array([[0, 0], [0, 0]]) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [],
   "source": [
    "def func_expected_opm(df,util, cm):\n",
    "    expected_opm = 0\n",
    "    number_of_overrides = 0\n",
    "\n",
    "    for index, row in df.iterrows():\n",
    "        predicted = row['Predicted']\n",
    "        actual = row['Actual']\n",
    "\n",
    "        if predicted == 'dog':\n",
    "            optimal = decide(prob_predict_dog, util).get(\"optimal\")\n",
    "            \n",
    "            if optimal == 1:\n",
    "                if 'dog' == actual: #RPopm did not override and resnet predicted dog and it was dog\n",
    "                    expected_opm += util[1][1]\n",
    "                    cm[1][1] += 1\n",
    "                else:   #RPopm did not override and resnet predicted dog and but it was cat\n",
    "                    expected_opm += util[1][0]\n",
    "                    cm[1][0] += 1\n",
    "            else:\n",
    "                number_of_overrides += 1\n",
    "                if 'cat' == actual: #RPopm override and guessed cat and it was cat\n",
    "                    expected_opm += util[0][0] \n",
    "                    cm[0][0] += 1  \n",
    "                else: #RPopm override and guessed cat and but it was dog\n",
    "                    expected_opm += util[0][1]\n",
    "                    cm[0][1] += 1\n",
    "        else:\n",
    "            optimal = decide(prob_predict_cat, util).get(\"optimal\")\n",
    "\n",
    "            if optimal == 0:\n",
    "                if 'cat' == actual: #RPopm did not override and resnet predicted cat and it was cat\n",
    "                    expected_opm += util[0][0]\n",
    "                    cm[0][0] += 1\n",
    "                else:   #RPopm did not override and resnet predicted cat and but it was dog\n",
    "                    expected_opm += util[0][1]\n",
    "                    cm[0][1] += 1\n",
    "            else:\n",
    "                number_of_overrides += 1\n",
    "                if 'dog' == actual: #RPopm override and guessed dog and it was dog\n",
    "                    expected_opm += util[1][1]  \n",
    "                    cm[1][1] += 1\n",
    "                else: #RPopm override and guessed dog and but it was cat\n",
    "                    expected_opm += util[1][0]\n",
    "                    cm[1][0] += 1\n",
    "\n",
    "    return expected_opm, number_of_overrides"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Trivial utility-matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[5, 2],\n",
       "       [4, 9]])"
      ]
     },
     "execution_count": 51,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "util"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 71,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[1744,   38],\n",
       "       [   6, 1712]])"
      ]
     },
     "execution_count": 71,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "CM"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 72,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[1750, 1750],\n",
       "       [   0,    0]])"
      ]
     },
     "execution_count": 72,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "cm_rn"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Expected utilities with resnet alone: 24228\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Expected utilities with OPM: 24228, total number over overrides 0\n"
     ]
    }
   ],
   "source": [
    "print(f'Expected utilities with resnet alone: {np.multiply(CM, util).sum()}')\n",
    "expected_opm = func_expected_opm(test_df, util, cm_rn)\n",
    "print(f\"Expected utilities with OPM: {expected_opm[0]}, total number over overrides {expected_opm[1]}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {},
   "outputs": [],
   "source": [
    "resnetUtil = []\n",
    "opmUtil = []\n",
    "for i in range(100):\n",
    "    util = np.array(utilityMatrixCreator(-10, 10))\n",
    "    expected_opm = func_expected_opm(test_df, util, cm_rn)\n",
    "    rsnetUtil = np.multiply(CM, util).sum()\n",
    "    resnetUtil.append(rsnetUtil)\n",
    "    opmUtil.append(expected_opm[0])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\47472\\AppData\\Local\\Temp\\ipykernel_15740\\1541215544.py:12: MatplotlibDeprecationWarning: Passing the pad parameter of tight_layout() positionally is deprecated since Matplotlib 3.3; the parameter will become keyword-only two minor releases later.\n",
      "  plt.tight_layout(True)\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjoAAAI6CAYAAAAuSOLAAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAABik0lEQVR4nO3deXgUVdbH8e8hBAigBtmEAIqIKIqCouIyCiKbGxlGBcR13NdxQwVxQUF00GFm1HHcUNxARIzooHkZJTouiGhARI0gKBIQZAkIRAjJff+oG2hiEgKkU53u3+d5+kn1reqqU7eru0/uvVVlzjlERERE4lGNsAMQERERiRYlOiIiIhK3lOiIiIhI3FKiIyIiInFLiY6IiIjELSU6IiIiEreU6IjsgJldZGYfVuL63jazC0tbt5mtN7P9K2tb1Y2Z/cHMcqp4m1lmdmlVblOqlpntZ2bOzGqGHYtUPSU6UunM7Aczy/c/2sWPR0OMJ7QfMjO7x8xejCxzzvVxzo0rbXnnXH3n3EL/2ufMbERVxBkrnHP/c861CzuOiqoOSVKJz+PP/riqv5vrvMgnDreWKF9iZl0r8PqYSzz8Ps01s42+nh43s9SI+feYWYGvxzwz+9jMjo14rTOzMSXW2deXP1e1eyORlOhItJzhf7SLH9eGHZDEtlj60YtDZzjn6gMdgU7AkEpY52rgVjPboxLWFSozuxl4EBgM7AV0AfYFpplZrYhFX/H12Bj4EJhsZubnfQ+cU+I4vhD4LtrxS/mU6EiV8v8lvRbx/EEze9cCXf1/hEPNbKX/T3RQxLK1zewhM1tsZsvN7N9mlhIxv6+ZzTazdWb2vZn1NrORwB+ARyNblszsIDObZmarzSzHzM6JWE9DM5vi1zMTaFPO/nQ1syUlyn4ws1PMrDcwFOjvtz3Hzy+zFcD/93eAmV0ODCL4IVlvZm+a2eDIuvPL/9PM/lHGulqa2WQz+8XMVkXsew0zG2ZmP5rZCjN73sz28vOK/9O+2Mx+MrM1ZnalmR1lZl/6/2QfjdjGRWb2kZk9amZrzexbM+seMf9iM/vGzH41s4VmdkXJujOz28zsZ+DZkvXp5+X61+cUr9sfC383s6X+8Xczq11ivTf7/VtmZheX9R56bcxspn/P3zCzvSNi6GLBf+95ZjbHfItFaceWmQ03s0f8/GQz22Bmo/3zFDP7rXjdZa3Xz9vLzJ7xseea2QgzS4qo8w8t+CysMbNFZtZnB/sHgHPuZyCTIOEpd/8itrXQ1/8ii/g8At8AnwA3lbYtf5zdbsFncZWZTYyo1w/83zxfd8eW8vqjzewTH9cyX7+1IuY7f2zO98s8ZhYkHWaW5OtnpZktBE4rq07MbE9gOHCdc+4d51yBc+4H4BxgP+C8UuqxABgH7AM09MU/A3OBXn69ewPHAVPK2rZUEeecHnpU6gP4ATiljHl1Cf7DuYjgR2Il0MLP6wpsAf4G1AZOAjYA7fz8MQRfGnsDewBvAqP8vKOBtUAPggQ+DTjIz8sCLo2IoR7wE3AxUJPgP9yVQHs/fwIw0S93KJALfFjG/nQFlpS1/8A9wIsl5m+Nx9fDhxHzHHCAn34OGBExr5mvj1T/vCawAjiylLiSgDm+zuoBdYAT/Lw/AwuA/YH6wGTgBT9vPx/Dv/1regK/ARlAE1+vK4CTIuLfAtwIJAP9/fuwt59/GkGiaP793AgcUeL9ftC/3ymR9Qm08+9T84jY2vjpe4EZPqbGwMfAfSXWe6+P6VS/3QZlvIdZ/j0+1NfVa8Xvmd/fVX4dNQiOr1VA4zKOrZOBuX76OIL/8j+NmDengut9HXjCx9MEmAlcEVHnBcBl/n2+ClgK2I4+j0ALgh/jf+woDr/tdWz7/DUDDok8bgkSpjUR7/cSoKuf/ot/j1r49/cJYHyJ46xmOd8jRxK0rNT0y38D3FDis/IWkAq0An4Bevt5VwLfAi0Jvi+ml7U9oDfB8VLavHERMd8TcVzUBkYDi0vUx7kErT4AV/t9HgE8F/b3ciI/Qg9Aj/h7+C/W9UBexOOyiPnHEDR7/wgMjCjv6r9w6kWUTQTuJPih3ID/ofPzjgUW+ekngDFlxJPF9j9G/YH/lVjmCeBugh+OAnyS5OfdTwwkOr7s7eK6BE4Hvi4jrmP9F39pX97vAldHPG/n97n4B8UBaRHzVwH9I56/hv/B8fFv9yNL8KN8fhlxZQB/iai7zUCd0uoTOIAgqToFSC6xnu+BUyOe9wJ+iFhHfuS++/V0Kef4eCDieXsfVxJwGz4JjJifCVxYxrGVQpAYNgRuJ2jRW0KQUA4H/umXK3O9QFNgE5ASMW8gMD2izhdEzKvr37N9dvB5/NUv9y7bkuXy4qhH8Nn9U2QsJY9bgs/og346MtH5Buge8Zpm/P44KzPRKWU/bgBeL/FZOSHi+UTgdj/9HnBlxLyeZW2PoMXm5zK2+QAwLeKzvNnXyQq/jSMj68O//8sJur9mAMejRCf0h7quJFrSnXOpEY+nimc45z4FFhIkLxNLvG6Nc25DxPMfgeYE/2HWBT73zdR5wDu+HIL/3L6vYGz7AscUr8evaxBBM3Rjgi/in0rEECvGsa0p/TzghTKWawn86JzbUsq85my/Tz8S7HPTiLLlEdP5pTyPHMya6/y3fcT6mgOYWR8zm2FBF2EeQctBo4hlf3HO/VbaDjjnFhD8uN0DrDCzCWbWvJx9aB7xfFWJfd9YIuaSSr7fyT7OfYGzSxwrJxD8aJcWcz4wi6D16kTgfYLWpuN92ft+0fLWu6/f/rKIeU8QtOwU+zlimxv9ZHn7l+6c24MgCTyIbe9BmXH4z2F/gtaRZWb2HzM7qJR13wVcZWZNS5TvC7wesd5vgEK2P87KZGYHmtlbFgwMXkfwD0ejEov9HDEd+R43p+Kf4ZVAIyt9jFgzP7/YRP991sQ5d7Jz7vPIhf37/x9gGNDQOfdROduVKqJER6qcmV1D0PS7FLi1xOwGZlYv4nkrv9xKgh/YQyKSp71cMDAQgi+1ssbSuBLPfwLeL5GI1XfOXUXQCrKFIFGIjKEsGwgSsOJ9S2Jb8lXatndGaa/NAA4zs0MJWnReKuO1PwGtyvjyXkrwI1SsFcE+Ly9l2YpIKx4bEbG+pX7MzGvAQ0BT51wqMJUgwS1Wbv045152zp3g43UE3Vxl7cPSXYwffv9+FxAccz8RtHhEHiv1nHMPlBP/+wTdVJ2Az/zzXgTdq8VjU8pb708ELTqNIubt6Zw7ZDf2LwjWufcJWgofqkAcOOcynXM9CH7wvwWeKmWd3xJ0f95RYtZPQJ8S667jnMulYp+Lx/022zrn9iRoHbPyX7LVMir+Gf6EoL77RRZacGZaH4IWsJ3xPHAz8OKOFpSqoURHqpSZHUjQlHsecD7BYNuOJRYbbma1zOwPBD/mrzrnigi+ZMeYWRO/rjQz6+Vf8wxwsZl194Mg0yL++1xOMB6l2FvAgWZ2vh8wmmzBYNuDnXOFBF/a95hZXTNrT9CMX5bvgDpmdpqZJRP8J1c7Yv5yYD8z25XPWsm48a0fk4CXgZnOucVlvHYmwZf9A2ZWz8zqmNnxft544EYza+2/zO8nGFdQWutPRTQBrvf1eDZwMEFCU4ugLn4BtvgBsz0rulIza2dmJ/uE6TeCRLcoYh+GmVljM2tE0KqwOz8s55lZezOrSzC2Z5I/Fl4EzjCzXn6Aax0LBju38K/73XtEkNhcQNCtuBnfvUXQzfqLX6bM9TrnlgH/BzxsZnv647mNmZ20G/sX6e9ADzM7vLw4zKypBQP86xEkAuvZVv8lDScY85YaUfZvYKSZ7Qvg36u+ft4vfl3lXTNqD4IxQuv9Z/mqndjHiQTHZAsza0DQjVgq59xaH/8jFpzAkGxm+/l1LKHsVtOyvE8w1umRnXydRIkSHYmWN2376+i87lsXXiToz5/jnJtP8F/aC/7HDIKm6DUE/52/RNDP/q2fdxvBINoZvin7vwTjS3DOzST4oh1DMBj2fbb9x/8P4CwLzlD5p3PuV4If3AF+Oz+zbUAswLUETeA/E/z3+2xZO+m/JK8GniYY0LqB4Mux2Kv+7yoz+6JiVbfVM0B73/SfEVE+DuhAOV/A/kf6DIJxLot9TP397LH+tR8AiwiSiOt2MrZInwJtCVpARgJnOedW+Xq+nuAHYw3BQM2dOQOlNsEYiZUE70UTtp0WPYKgi+hLgsG1X/iyXfUCwXv9M8Eg7OsBnHM/AX0JjtNfCFopBrPtu3O7Y8uXfUwwVqO49eZrgjoufl6R9V5AkCh+TVB3kyiju2xn+WTreeCuHcRRg+CMqqUEY+pOooxkwzm3iKAOI1tj/0Hwfv+fmf1KMGblGL/8RoJj5SN/fHcpZbW3EBwzvxL8k/PKTuzmUwRjjeYQHBuTy1vYOfdXgjp4iCC5+pSgLro75zbtxHZxgXedc6t35nUSPbZ917pIeCw4rfVF51yLHSya0MysFUGT/j7OuXUhx3IRwWDcE8KMQ0SkLGrREalGfBfYTcCEsJMcEZHqQFciFakm/FiJ5QRnkPQOORwRkWpBXVciIiISt9R1JSIiInFLiY6IiIjErYQbo9OoUSO33377RX07GzZsoF69ejteMI6pDlQHoDoA1UGi7z+oDiC6dfD555+vdM41Lm1ewiU6++23H7NmzYr6drKysujatWvUtxPLVAeqA1AdgOog0fcfVAcQ3TowszJv86GuKxEREYlbSnREREQkbinRERERkbilREdERETilhIdERERiVtKdERERCRuKdERERGRuKVER0REROKWEh0RERGJW0p0REREJG4p0REREZG4pURHRERE4pYSHREREYlbSnREREQkbinRERERkbilREdERETilhIdERERiVtRS3TMrI6ZzTSzOWY2z8yG+/LnzGyRmc32j46+3Mzsn2a2wMy+NLMjItZ1oZnN948LI8qPNLO5/jX/NDOL1v6IiIhI9VMziuveBJzsnFtvZsnAh2b2tp832Dk3qcTyfYC2/nEM8DhwjJntDdwNdAYc8LmZTXHOrfHLXAZ8CkwFegNvIyIiIkIUW3RcYL1/muwfrpyX9AWe96+bAaSaWTOgFzDNObfaJzfTgN5+3p7OuRnOOQc8D6RHa39ERESk+onqGB0zSzKz2cAKgmTlUz9rpO+eGmNmtX1ZGvBTxMuX+LLyypeUUi4iIiKx5J13qLluXSibjmbXFc65QqCjmaUCr5vZocAQ4GegFvAkcBtwbzTjMLPLgcsBmjZtSlZWVjQ3B8D69eurZDuxTHWgOgDVAagOEn3/IbHroNH773PIvffSskcPsvbcs8q3H9VEp5hzLs/MpgO9nXMP+eJNZvYscIt/ngu0jHhZC1+WC3QtUZ7ly1uUsnxp23+SIKmic+fOrmvXrqUtVqmysrKoiu3EMtWB6gBUB6A6SPT9hwSvg48/hi5dWHz99aHUQTTPumrsW3IwsxSgB/CtH1uDP0MqHfjKv2QKcIE/+6oLsNY5twzIBHqaWQMzawD0BDL9vHVm1sWv6wLgjWjtj4iIiOyEjRuDv0OHwvTpFNatG0oY0Ryj0wyYbmZfAp8RjNF5C3jJzOYCc4FGwAi//FRgIbAAeAq4GsA5txq4z6/jM+BeX4Zf5mn/mu/RGVciIiLhe+st2H9/+PLL4HmtWqGFErWuK+fcl0CnUspPLmN5B1xTxryxwNhSymcBh+5epCIiIrK7MrJzGZ2Zw4Gff8ATr9/PhnbtadCqVdhh6crIIiIisnsysnMZMnkuB3zxIf9+/X5yGrWid587yFi0IezQlOiIiIjI7hmdmUPbxd/w5OQRLGjYkvP6j2B5zbqMzswJO7SqOetKRERE4tfSvHx+adya5484jceOPYe1KXtsLQ+bEh0RERHZdR9/zEHJm/mGWow8+dLtZjVPTQkpqG3UdSUiIiK7JisLTjmFJ2eNIyU5abtZKclJDO7VLpy4IijRERERkZ33wQdw2mnQujUtxz3JqH4dSEtNwYC01BRG9etAeqfw78ykrisRERHZOR9+CKeeCq1awXvvQZMmpDchJhKbktSiIyIiIhVXVATXXgstWgRJTtOmYUdULrXoiIiISMXVqAFvvglJSdCsWdjR7JBadERERGTHZs4MWnIKC6FlS2jePOyIKkSJjoiIiJRv1izo2RPefhtWrQo7mp2iREdERETK9sUX0KMH7L03TJ8OTZqEHdFOUaIjIiIipcvOhlNOgdTUIMmJgZt07iwlOiIiIlK6vDzYZ58gydl337Cj2SU660pERES2t3Yt7LUXdOsGX34JNatvuqAWHREREdnmq6/gwAPh+eeD59U4yQElOiIiIlLs66/h5JOD5ObYY8OOplIo0RERERH49tsgyUlKCq543LZt2BFVCiU6IiIiiW7NmmA8DgQDj9uFf9fxyqJER0REJNE1aAB33BG05Bx0UNjRVKrqPcJIREREdt2CBUFrzlFHBbd3iENKdERERBJIRnYuozNzqLloIRNfGUq9PepSf9GCan92VVnUdSUiIpIgMrJzGTJ5LvbDIl4eP5Ram37jvD63kjF3edihRY0SHRERkQQxOjOHhiuXMmH8UOoV5DNowEhm770vozNzwg4tapToiIiIJIileflcOvN19ti0gUH9R/B10/23lser+OyQExEREWDbmJylefnUMGPEyZfy/BGns7Bhi63LNE9NCTHC6FKLjoiISJwqHpOz5aefeGLyCFLXr2FLUs3tkpyU5CQG94qf6+aUpBYdERGRODU6M4c9Vy9n/PihNN6whubrfmFVvVSSzChyjuapKQzu1Y70Tmlhhxo1SnRERETiVMGSXMZPuIMmG9ZwwTn3MrdZcFuHIudY9MBpIUdXNZToiIiIxKOff+bViXfQ6NdVXHDOvXyRdvDWWfE8JqckJToiIiLxqKiIPRs14Io+1/N5s21JTryPySlJg5FFRETiyerVsGULNG9Og6+yOeuGgaSlpmBAWmoKo/p1iOsxOSWpRUdERCRe/PJLcBfyo4+GsWOhRg3SO6UlVGJTklp0RERE4sHKldC9O3z/PZx3XtjRxAwlOiIiItXdqlVwyikwfz68+SacfHLYEcUMdV2JiIhUZ87BH/8I334LU6YECY9spURHRESkOjODkSNh/Xro2TPsaGKOEh0REZHqKC8P3nkHBgyAP/wh7GhilsboiIiIVDdr10KvXnDhhfDDD2FHE9PUoiMiIlINFN+FfN3yVYx/7W7aL5tPjUmTYL/9wg4tpqlFR0REJMYV34U8b/kqnn31btrlfsdf0m8no1XnsEOLeUp0REREYtzozBzyCwrpuvBzDl/2HdedeStvtunC6MycsEOLeeq6EhERiXFL8/IB+M/Bf2DuPgewuEGz7cqlbGrRERERiWUbN/LClJF0XjIPYGuSA4l1F/JdpURHREQkVm3cCGecwfE5n7LvhtXbzUq0u5DvKiU6IiIisSg/H/r2henTsXHj+MNd1yX0Xch3lcboiIiIxJrffoP0dHj3XXj2WTjvPNJBic0uUIuOiIhIrKlZExo1gmeeCS4KKLtMLToiIiKxYtMmWLcOGjeGF18M7mMlu0UtOiIiIrFg0yY46yw46aSg60pJTqVQi46IiEjYNm+Gc86Bt96Cxx+HOnXCjihuqEVHREQkTAUFwR3Ip0yBRx+FK68MO6K4okRHREQkTHfcAa+/Dv/8J1xzTdjRxB11XYmIiITpllvg0EPhggvCjiQuqUVHRESkqm3ZErTgFBRAkyZKcqJIiY6IiEhVKiwMro3zl78Eg48lqpToiIiIVJXCQrjoInj5ZRg1Cv74x7AjintRS3TMrI6ZzTSzOWY2z8yG+/LWZvapmS0ws1fMrJYvr+2fL/Dz94tY1xBfnmNmvSLKe/uyBWZ2e7T2RUREZLcVFsIllwQXAhwxAm7Xz1ZViGaLzibgZOfc4UBHoLeZdQEeBMY45w4A1gCX+OUvAdb48jF+OcysPTAAOAToDfzLzJLMLAl4DOgDtAcG+mVFRERiz8KFwdlVw4cHZ1pJlYhaouMC6/3TZP9wwMnAJF8+Dkj30339c/z87mZmvnyCc26Tc24RsAA42j8WOOcWOuc2AxP8siIiIrHDueBv27bw9ddw113hxpNgojpGx7e8zAZWANOA74E859wWv8gSoPhWrGnATwB+/lqgYWR5ideUVS4iIhIbiorgqqto+corwfM0/UxVtaheR8c5Vwh0NLNU4HXgoGhuryxmdjlwOUDTpk3JysqK+jbXr19fJduJZaoD1QGoDkB1kLD77xxt//530qZMoejssxOzDiKEdRxUyQUDnXN5ZjYdOBZINbOavtWmBZDrF8sFWgJLzKwmsBewKqK8WORryiovuf0ngScBOnfu7Lp27VoZu1WurKwsqmI7sUx1oDoA1QGoDhJp/zOycxmdmcPSNRt56IOn6TpjCtx6K7m9eydMHZQlrOMgmmddNfYtOZhZCtAD+AaYDpzlF7sQeMNPT/HP8fPfc845Xz7An5XVGmgLzAQ+A9r6s7hqEQxYnhKt/RERESlPj79lccMrs8nNy2fYe0/zpxlvMLZLPzL6X6c7kYcomi06zYBx/uyoGsBE59xbZvY1MMHMRgDZwDN++WeAF8xsAbCaIHHBOTfPzCYCXwNbgGt8lxhmdi2QCSQBY51z86K4PyIiIqUa9NQnzF+xYevz7xu24Kmj0hl54sWk/d93jOyiy9aFJWqJjnPuS6BTKeULCc6YKln+G3B2GesaCYwspXwqMHW3gxUREdlFGdm5fPT9anCOffOW8WOD5rzcsc/W+Uvz8oF64QWY4JRiioiI7KKM7FyGTJ4LznH7+8/xztjraLPqp+2WaZ6aElJ0Akp0REREdtnozBzyN29h8AfPc+WnrzGpQ3e+37vFdssM7tUupOgElOiIiIjssqVrNnLz/17kmhmv8lLH3tzV48rtBh63bVKP9E66dk6YlOiIiIjsooG5n3PdJ68w/rCeDOt5Nc62/awe32Zvpt3UNbzgBKii6+iIiIjEo6NvuJiha9cyvn23rUlOSnISo/p1UEtOjFCLjoiIyM56/HFYsoT0zq04+t6bad6gHgakpaYoyYkxatERERHZGaNGwdCh8OOP8MADpHdKU2ITw9SiIyIiUlF//WuQ5AwaBCN/d3k3iUFKdERERCri4YfhtttgwAB47jlISgo7IqkAJToiIiI7smkTvPACnHNO8LemRn5UF3qnREREylNUBLVrw/TpUL++kpxqRi06IiIiZXnsMUhPD1p0GjSA5OSwI5KdpERHRESkNI8/DtdeCzVqbHe1Y6lelOiIiIiU9OSTcPXVcMYZMHEi1KoVdkSyi5ToiIiIRHruObjiCjjtNHj1VSU51ZwSHRERkUgdOsC558KkScEgZKnWlOiIiIgAfPVV8PfII+Gll6BOnXDjkUqhREdEROTFF+Gww4IER+KKEh0REUlsL78MF14I3brBH/8YdjRSyXTVIxERSSgZ2bmMzsxhaV4+5//4CcMnjsJOPBHefBPq1g07PKlkSnRERCRhZGTnMmTyXPILCklbu4JhEx9kVov2/DzySc5QkhOXlOiIiEjCGJ2ZQ35BIQC5ezXhqj8O4ZNWh9Hgg58447i2IUcn0aBER0RE4l5xd1VuXj69vvuY/Jq1+WD/I3n3gGMAyM/LDzlCiRYlOiIiEteGZczlpRmLcUCP+TN49I0H+azFIXzQ+oitt3ZonpoSbpASNTrrSkRE4lZGdu7WJKf7gk95LOMBvmp6AJf3G7Y1yUlJTmJwr3bhBipRoxYdERGJW6Mzc3BAt+8/4/HXR/F109Zc0P9e1tcOBh6npaYwuFc70julhRuoRI0SHRERiVtL/dibk7//jG+b7McF59zHr7XrAUGS89HtJ4cZnlQBJToiIhK3Wu5Ri8W/buauHldSd/NvbPAtOQbqrkoQGqMjIiLxado0/vP01bTZsBJnNbZLcgZ1aaXuqgShREdEROLPu+/CmWeyx171uKnfkaSlpmAE3VVj+ndkRHqHsCOUKqKuKxERiS/vvQdnnAFt28J//8tpjRpxWrewg5KwqEVHRETix4wZQZKz//5Bq06jRmFHJCFToiMiIvGjXTvo1y9Icho3DjsaiQFKdEREpPqbMwd++w0aNIAXXoCmTcOOSGKEEh0REanePvkETjgBbrgh7EgkBinRERGR6mvGDOjVC5o1g7vuCjsaiUFKdEREpHqaOTNIcpo0genToXnzsCOSGKRER0REqp8tW2DQoOCsqunTIU0X/5PS6To6IiJS/dSsCZMnw157QcuWYUcjMUwtOiIiUn3Mng2jRoFz0KEDtGoVdkQS45ToiIhI9fDll3DKKfD447BmTdjRSDWhREdERGLf3LnQvTukpARjcvbeO+yIpJpQoiMiIrFt3rwgyaldO0hy2rQJOyKpRjQYWUREYkpGdi6jM3NYmpdP89QUxth3HF2nTnBbhwMOCDs8qWbUoiMiIjEjIzuXwa/OITcvn+QtBeTm5TNobSvenPBecDdykZ2kREdERGLGPVPmUVDk2H/VEt576gq6fv8ZBUWOO6ctDDs0qabUdSUiIjEjL7+A/VbnMn7CUJKKClmyV9Ot5SK7QomOiIiEKnJMzr5rljJ+/FBqFm5h4MD7WdBI18mR3aNER0REQlM8JqegyNFowxrGjx9K7cICBg68n+8a77d1uQZ1k8MLUqo1JToiIhKKjOxcbnxlNs4/X1k3lYxDuvLmwSeSE5HkJCcZd59xSCgxSvWnREdERKrcsIy5vDRjMQ5IW7uCpKJCFjdoxl9PumjrMgY0T01hcK92pHfSTTtl1yjRERGRKpWRnbs1yWm+bgUTxg9hfa0UTr34nzjbdjLwogdOCy9IiRs6vVxERKrU6MwcHNBs3S+MHz+UvX5bz619/rJdkqMxOVJZlOiIiEiVWpqXzz7rVjJ+/FAabFzH+efcy9xm2y4GqDE5UpmU6IiISJVqnprCbe8/R8ONeVzQ/z7mNG+3dV6DusmMPutwjcmRSqMxOiIiUqUG92rHiNXX8vSKdObtE9y7yoBBXVoxIr1DuMFJ3FGLjoiIVI3ly+GKK0g/MJVhA7uQd1AHDEhLTWFM/45KciQqotaiY2YtgeeBpoADnnTO/cPM7gEuA37xiw51zk31rxkCXAIUAtc75zJ9eW/gH0AS8LRz7gFf3hqYADQEPgfOd85tjtY+iYjILlqxAk4+GX74AS6/nPQjj1T3lFSJaLbobAFuds61B7oA15hZez9vjHOuo38UJzntgQHAIUBv4F9mlmRmScBjQB+gPTAwYj0P+nUdAKwhSJJERCSGJOflBUnOokXwn//AkUeGHZIkkKglOs65Zc65L/z0r8A3QHnpe19ggnNuk3NuEbAAONo/FjjnFvrWmglAXzMz4GRgkn/9OCA9KjsjIiK7ZuVKDr/5Zli4EN56C7p2DTsiSTBVMhjZzPYDOgGfAscD15rZBcAsglafNQRJ0IyIly1hW2L0U4nyYwi6q/Kcc1tKWb7k9i8HLgdo2rQpWVlZu79TO7B+/foq2U4sUx2oDkB1AIldB3UXL+bQ9euZfd995NWoAQlaD4l8DBQLqw6inuiYWX3gNeAG59w6M3scuI9g3M59wMPAn6MZg3PuSeBJgM6dO7uuVfAfRVZWFlWxnVimOlAdgOoAErQONmyAunXBjPebNeOkHj3CjihUCXkMlBBWHUT1rCszSyZIcl5yzk0GcM4td84VOueKgKcIuqYAcoGWES9v4cvKKl8FpJpZzRLlIiISpjVr4MQTYdgwAFyyrnIs4YlaouPH0DwDfOOc+1tEebOIxf4IfOWnpwADzKy2P5uqLTAT+Axoa2atzawWwYDlKc45B0wHzvKvvxB4I1r7IyIiFZCXBz17wldfwfHHhx2NSFS7ro4HzgfmmtlsXzaU4KypjgRdVz8AVwA45+aZ2UTga4Iztq5xzhUCmNm1QCbB6eVjnXPz/PpuAyaY2QggmyCxEhGRMKxdC716wZw5MHkynHpq2BGJRC/Rcc59SHCxy5KmlvOakcDIUsqnlvY659xCtnV9iYhIFRv01Cd89P1qzBXxysu3c8Sy76g5+TU4/fSwQxMBdAsIERHZRT3+lsX8FRsAcFaD5448k6drJLFheWNeCjk2kWJKdEREZKdlZOcyf8UG6m3ayGE/L+CTfQ9j6kEnBDO/Xx1ucCIRdK8rERHZaaMzc6i7OZ+xk4bz7KR7aLx+TdghiZRKiY6IiOy0NSvWMHbScDrnfsMtp97AL/UbhB2SSKnUdSUiIjtn40ZeeGMEHZd8zY2n38xbB5+43ezj2+wdUmAiv6dER0REdigjO5fRmTkszcvn6pz/cvMPX3LbmTcz5aCTtluubZN6vHTZsSFFKfJ7SnRERKRcGdm5DJk8l/yCQgAeO7A7n17UmoPO6Ebat7+wNC+f5qkpDO7VjvRO5d27WaTqKdEREZFyjc7MoSg/n4cyH+NfXc5mYcMWzGrchmXf/sJHt58cdngi5dJgZBERKdfKlet4/PX7Oeurdzns5/lby5fm5YcYlUjFqEVHRETKtmkTz/znQU5YOIvbe11LxiHdts5qnpoSYmAiFaNER0RESrd5M5xzDid8O4O7+1zHhMN6bZ2VkpzE4F7tQgxOpGLUdSUiIqXbvDm4G/ljj9Fp5G2kpaZgQFpqCqP6ddDAY6kW1KIjIiLbKyiATZugfn147z1ISiIdlNhItaQWHRER2aagAM49F/r0gS1bICkp7IhEdosSHRERCWzZAuedB5MmQb9+UFON/lL9KdEREZEgybngApg4EUaPhhtvDDsikUqhREdERGDwYBg/Hh54AG65JexoRCqN2iVFRASuugpat4brrw87EpFKpRYdEZFEVVQUtOI4BwceqCRH4pISHRGRRFRUBJddFpxh9c47YUcjEjVKdEREEk1REVx5JYwdC3feGZxKLhKnlOiIiCQS5+Caa+Cpp2DoUBg+POyIRKJKiY6ISCKZMweeeQZuvx1GjACzsCMSiSqddSUiEseGZcxl/Kc/UegcSWYMPKYlI7KzoX17JTmSEJToiIjEqR5/y2L+ig3gHHdMf4Y5zQ7kRXci0IoRhyjJkcSgrisRkTg06KlPtiY5w957mss+y+CwZfMBGP/pTyFHJ1J1lOiIiMSZYRlz+ej71eAcQ7Ke5dJZb/DskWdwf7c/A1DoXMgRilQdJToiInFkWMZcXpyxGJzj9vef44qZkxl3xGkM73751jE5SRqbIwlEY3REROLEoKc+CVpyvCIzXuh0KnefcuV2A48HHtMyjPBEQqFER0QkDkR2VzXcuJZV9VL564kXBjMjkpy2TeoxIr1DSFGKVD11XYmIxIHiAcY3fPQymWOvofm6FUGCE5HkHN9mb6bd1DWkCEXCoURHRCQOFDrHdR+N54aPxvNem6NYtkej7eaf16UVL112bEjRiYRHiY6ISBy49pOJ3PzhS7x26Mnc3vs6nG37ej+vSyt1V0nCUqIjIlLdjR/PLR88z+RDujG4z18oqpG0ddbxbfZWkiMJTYORRUSqu7594cEHyW7TE5u1FCJv96AkRxKcEh0Rkerq5ZfhtNNgr73g1lu5D7jvTx3DjkokpqjrSkSkOvr732HQIHj44bAjEYlpSnRERKqbRx6BG2+EP/0J7rwz7GhEYpoSHRGR6uRf/4Lrr4f0dBg/HpKTw45IJKYp0RERqS42bIAHH4Qzz4RXXlGSI1IBGowsIlJd1KsHH30EjRtDrVphRyNSLahFR0Qk1o0dC1dfDUVF0KIF1K4ddkQi1YYSHRGRWPbcc3DppbBwIRQUhB2NSLWzw0THzHS1KRGRMDz/PPz5z9CjB2RkqCVHZBdUZIzOv8ysNvAc8JJzbm10QxIRSUwZ2bmMzsxhaV4+Fy36kLtefRDr3j1IcurUCTs8kWpphy06zrk/AIOAlsDnZvaymfWIemQiIgkkIzuXIZPnkpuXjwMWFdXmgzadeXP445CSEnZ4ItVWhc66cs7NN7NhwCzgn0AnMzNgqHNucjQDFBFJBKMzc8gvKKTZul9YtmdjstocRdb+nUn7YDFnHHdA2OGJVFsVGaNzmJmNAb4BTgbOcM4d7KfHRDk+EZGEsDQvn1O//ZCsJy+j6/efBYVmLM3LDzcwkWquIi06jwBPE7TebP3EOeeW+lYeERHZBZFjcvrkfMw/p/yVL9IOYmbLQ7cu0zxV3VYiu6Mip5e/7px7ITLJMbO/ADjnXohaZCIicSwjO5fBr84hNy+fnt99zD+mPMjs5u24+Kx72FgrSG5SkpMY3KtdyJGKVG8VSXQuKKXsokqOQ0QkYWRk53LjK7MpKHK0WfUTj77xIHP3OYCLzh7Oxtp1MSAtNYVR/TqQ3ikt7HBFqrUyu67MbCBwLtDazKZEzNoDWB3twERE4lHx2VXOP/++YUuG9byaqQedwPradQH44YHTwgtQJM6UN0bnY2AZ0Ah4OKL8V+DLaAYlIhKPiltyHHDSws9ZUb8B3zTZn1cO7xV2aCJxq8xExzn3I/AjcGzVhSMiEp+Kx+QUJzlPTr6PT1t24IL+9223XIO6uiO5SGUqr+vqQ+fcCWb2K2xtZQUwwDnn9ox6dCIicSAvv4AbM2fjHPxh0Rc8OXkECxq24rozb91uueQk4+4zDgkpSpH4VF6Lzgn+7x5VF46ISHzJyM5lyZp8nEvi+B9m89TkEXzfsAWDBoxgbcq2r9cGdZO5+4xDNPhYpJKVedaVme1d3mNHKzazlmY23cy+NrN5xaek+9dPM7P5/m8DX25m9k8zW2BmX5rZERHrutAvP9/MLowoP9LM5vrX/NNfrVlEJGaMzszBuaBR/LzsqSxq0JxB/UeQl7KtUbxB3WSy7+qpJEckCsobjPw5QZdVacmDA/bfwbq3ADc7574wsz0I7pM1jeDU9Hedcw+Y2e3A7cBtQB+grX8cAzwOHOOTqruBzn67n5vZFOfcGr/MZcCnwFSgN/D2DvdaRKSKLM3LhxZBonPDGbdQd3M+a+rutXW+uqtEoqvMFh3nXGvn3P7+b8nHjpIcnHPLnHNf+OlfCW4hkQb0Bcb5xcYB6X66L/C8C8wAUs2sGdALmOacW+2Tm2lAbz9vT+fcDBf8u/R8xLpERGJCnzXz6TdyGHvl/8qmmrW2S3JqGIw+63C15IhEUXmDkQ9yzn0b2YUUqTiJqQgz2w/oRNDy0tQ5t8zP+hlo6qfTgJ8iXrbEl5VXvqSUchGR2PDhh/zzhTtY12Bvkou2bDcruYYx+mwlOSLRVl7X1U3A5Wx/DZ1ijuCmnjtkZvWB14AbnHPrIofROOecmbkyX1xJzOxygn2hadOmZGVlRXuTrF+/vkq2E8tUB6oDSNw62POrrzjs1lvZ3KgRn957H5fU3pMtPtlJqmE0T00hde18srLmhxxp9CXqMRBJdRBeHZR31tXlfrKPc+63yHlmVqciKzezZIIk5yXn3GRfvNzMmjnnlvnupxW+PBdoGfHyFr4sF+haojzLl7coZfnS9uVJ4EmAzp07u65du5a2WKXKysqiKrYTy1QHqgNI0DqYOROGDoUWLaiZlUXKd99xVaLVQYSEPAZKUB2EVwcVudfVxxUs244/A+oZ4Bvn3N8iZk0Bis+cuhB4I6L8An/2VRdgre/iygR6mlkDf4ZWTyDTz1tnZl38ti6IWJeISHiaNoVjj4Xp06F587CjEUlo5Y3R2YdgzEuKmXVi29lXewJ1K7Du44HzgblmNtuXDQUeACaa2SUEV14+x8+bCpwKLAA2AhcDOOdWm9l9wGd+uXudc8X32roaeA5IITjbSmdciUh4vv8eWreGffeFzMywoxERyh+j04vgVPAWQGSLzDqChKVczrkPKf3UdIDupSzvgGvKWNdYYGwp5bOAQ3cUi4hI1H3+OZxyClx5JYwaFXY0IuKVN0ZnHDDOzP7knHutCmMSEalesrOhRw/Ya68g0RGRmFFe19VNJaYdsBL40Dm3qApiExGJfXPmBC05e+wRjMnZd9+wIxKRCOUNRt6jxGNPgqsTv21mA6ogNhGR2LZpE5x+OtSrFyQ5rVuHHZGIlFBe19Xw0sr9LRn+C0yIVlAiIrEoIzuX4W/OY83GAgBSU5J59M6HOaH7kbD/Di8YLyIhqMjp5dvxZzzp5pkiklCGZczlhldms2ZjAW1/+ZF+X71LXn4BF/9Qn4x1Fbq0mIiEYKcTHTPrBqyJQiwiIjFpWMZcXpyxGIADVi7m5Ql3cOv746i7OZ+CIsfozJyQIxSRspQ3GHkuwQDkSHsDSwkuziciEvcysnN5ySc5bVb9xPgJQ8Fg0ICRbKyVAvg7lItITCrvOjqnl3jugFXOuQ1RjEdEJKaMzszBAa1X5/LyhDvAwYCBo/i+4bY71jRPTQkvQBEpV3mDkX+sykBERGJNRnYuub615sRFX5BUVMjAAffzfaNtSU5yDWNwr3ZhhSgiO1Bei46ISMLKyM5lyOS5mCvCWQ3GHXkGUw4+kTV199q6TEpyDUb1O4z0TmkhRioi5dnpwcgiIolgdGYOjX7JZeqz13PYsu8AtiY5BpzXpRXf3NdHSY5IjFOLjohIKezHHxg/fgj1Nv/Glhrbf1WO6d9RCY5INVHeWVe/8vuzrrZyzu0ZlYhERML244+8OmEoKZvzGTRgJF833XYxwLTUFCU5ItVIeYOR9wAws/uAZcALBC22g4BmVRKdiEhVW7YMunWjYeFvDBx0P/MabUtyUpKTNPBYpJqpyBidM51z/3LO/eqcW+ecexzoG+3ARERC0bAh/OEP1Hp3Gudf04+01BSMoCVnVL8Oas0RqWYqMkZng5kNIri3lQMGArqWjojEl6VLoVYtaNQIxo0DIB2U2IhUcxVp0TkXOAdY7h9n+zIRkfiwdCl06wbp6eDKHJooItXQDlt0nHM/oK4qEYlXy5bByScHyc7YsWC6Z7FIPNlhi46ZHWhm75rZV/75YWY2LPqhiYhE2fLlQZKzZAm8/TYcf3zYEYlIJatI19VTwBCgAMA59yUwIJpBiYhUiSuugMWLYepUOOGEsKMRkSioyGDkus65mbZ9c+6WKMUjIlJ1/vUvWLhQSY5IHKtIi85KM2uDv3igmZ1FcF0dEZHqZ+VKuOsu2LIFmjdXkiMS5yrSonMN8CRwkJnlAosILhooIlK9rFoFp5wCOTnQrx907Bh2RCISZRVJdJxz7hQzqwfUcM79amatox2YiEilWr0aevSAb7+FKVOU5IgkiIp0Xb0G4Jzb4Jz71ZdNil5IIiKVbM0a6NkT5s2DjIxgWkQSQnk39TwIOATYy8z6RczaE6gT7cBERCrNd9/BokXw+uvQu3fY0YhIFSqv66odcDqQCpwRUf4rcFkUYxIR2WUZ2bmMzsxhaV4+LfeoxU2ntif9mGOCRGfPPcMOT0SqWHl3L38DeMPMjnXOfVKFMYmI7JKM7FwGT5pDQaGj/qaNjHnhFqb8rxuMukP3rBJJUBUZo3OlmaUWPzGzBmY2NnohiYjsvIzsXG58ZTYFhY56mzby3Kt3c9jP8/m5bgOGvzkv7PBEJCQVOevqMOdcXvET59waM+sUvZBERHbOoKc+4aPvVwNQd3M+z066h45Lc7i2721kHngcbCwIOUIRCUtFWnRqmFmD4idmtjcVS5BERKJuWMbcrUlOUlEhz04azhG533L9mbfyTjvdu0ok0VUkYXkY+MTMXvXPzwZGRi8kEZGKe2nG4q3ThTWS+L+2XXih06lMPWjbFY9TU5LDCE1EYsAOEx3n3PNmNgs42Rf1c859Hd2wRER2bFjGXBxQp+A3Wq9ZyjdN9ueZo9J/t9w9Zx5S5bGJSGyoSNcVwN7ABufco8AvujKyiIRt0FOf8OKMxdQu2MRTr43glZeHsFf+r79b7rwurXTGlUgC22GLjpndDXQmuK7Os0Ay8CKgzm8RCUXxuJzaWzbz1OQRHP/jHAafegNrU/bYbrm/9++oJEckwVWkReePwJnABgDn3FJgj3JfISISReM//YnaWzbzxOSRnPDDbG7rcz2vdei+3TJKckQEKpbobHbOOcAB+Jt7ioiEptA5Lvp8Cl0Xfc7tva/j1cN6bDdf3VUiUqwiZ11NNLMngFQzuwz4M/BUdMMSESlbkhnPdE7n6yb787/WR2w377wurRiR3iGkyEQk1uywRcc59xDB3cpfAw4E7nLOPRLtwEREfmfzZrjpJi5tV5ctSTV/l+Qc32ZvJTkisp2KnnU1F/gf8IGfFhGpWgUFMGAAjBnDkKSfOK9LK5LMgKCF57wurXjpsmNDDlJEYk1Fzrq6FLgLeA8w4BEzu9c5p/tdiUjVKCiAgQPh9dfhH/+A889nBKj1RkR2qCJjdAYDnZxzqwDMrCHwMaBER0Sib8sWGDQIXnsNxoyB668POyIRqUYq0nW1Coi8CtevvkxEJPry8uCrr+Chh+CGG8KORkSqmYq06CwAPjWzNwhOMe8LfGlmNwE45/4WxfhEJFEVFoJz0KgRfP45pKSEHZGIVEMVSXS+949ib/i/umigiERHYSFcdBFs2gQTJijJEZFdVpFE50Hn3G+RBWbWyDm3MkoxiUgiKyyESy6BF1+EESOgRkVPDhUR+b2KfIPMNLMuxU/M7E8Eg5FFRCpXURFcdhmMGwf33gt33BF2RCJSzVWkRWcQMNbMsoDmQEPg5GgGJSIJ6sYb4dln4e674c47w45GROLADhMd59xcMxsJvEBwxtWJzrklUY9MRBJP//7QpAkMHRp2JCISJypywcBngDbAYQS3gHjLzB5xzj0W7eBEJH5lZOcyOjOHpWs2cvrqHLpfPYD0446D444LOzQRiSMV6bqaC1zq72C+yMyOAXRKuYjssozsXIZMnkv+5i0M/++/ufCL/3DeZgc3DNRdx0WkUpU5GNnM9gRwzv3dJzn452uB4VUQm4jEoYzsXG6eOIf8zVu4+90nufCL//DkUX/kw30OYnRmTtjhiUicKe+sq6ziCTN7t8S8jGgEIyLxbVjGXG58ZTaFRUXc+d7TXPz5mzzduS/3d/szmLE0Lz/sEEUkzpSX6FjE9N7lzBMR2aGM7FxemrEYB3TO/ZpLZr3B2CPPZMTJl4K/C3nzVF0YUEQqV3ljdFwZ06U9FxEpU3F3VfEXx6wWh3DOuQ8ws8UhW5OclOQkBvdqF16QIhKXykt0mvj7WVnENP5546hHJiJxIS+/gCHvzqWwqIgbPnqZj/c9nJktD2Vmy0O3LpNkxqh+HTQQWUQqXXmJzlNsu59V5DTA01GLSETiRkZ2LktW55O/uQY3/+9FrvvkFfb8bcN2SY4BD59zuJIcEYmKMhMd59xunVllZmOB04EVzrlDfdk9wGXAL36xoc65qX7eEOASoBC43jmX6ct7A/8AkoCnnXMP+PLWwASCKzV/DpzvnNu8OzGLSOUpPoX86oMcN374Mtd98grjD+vJfd0v3bqMAYO6tFKSIyJRE8275T0H9C6lfIxzrqN/FCc57YEBwCH+Nf8ysyQzSwIeA/oA7YGBflmAB/26DgDWECRJIhIjhr85j/yCQo56/RX+8vF4JnY4haG9r8VZ8LWTZMaY/h0Zkd4h5EhFJJ5FLdFxzn0ArK7g4n2BCc65Tc65RcAC4Gj/WOCcW+hbayYAfc3MCO63Ncm/fhyQXpnxi8iuy8jOZc3GAswV0WBZLpMO7c5tfa7fmuSkJCepu0pEqkQ0W3TKcq2ZfWlmY82sgS9LA36KWGaJLyurvCGQ55zbUqJcRGLA6MwcUjb/hrMaTLvyL9wakeRo4LGIVCWLuOjx9jO2nWVVKufcDm8DYWb7AW9FjNFpCqwkOD39PqCZc+7PZvYoMMM596Jf7hngbb+a3s65S335+cAxwD1++QN8eUvg7eLtlBLH5cDlAE2bNj1ywoQJOwp9t61fv5769etHfTuxTHWQuHWQ/MSztM/6L6/deT977pPK8ojrALbcuy6pKcnhBReCRD0OiiX6/oPqAKJbB926dfvcOde5tHnlnXX1EDCbIOHYRCVcJNA5t7x42syeAt7yT3OBlhGLtvBllFG+Ckg1s5q+VSdy+dK2+yTwJEDnzp1d165dd29HKiArK4uq2E4sUx0kaB089BBMeJ4pB5/IqB8acMNejofnBl81qSnJzL67R8gBVr2EPA4iJPr+g+oAwquD8hKdTsBA4DSCs5rGA++6spqAKsDMmjnnlvmnfwS+8tNTgJfN7G9Ac6AtMJMguWrrz7DKJRiwfK5zzpnZdOAsgnE7FwJv7GpcIlJJxoyBwYNZ0uMMhh51BYWFAEEPc0pyEveceUio4YlI4ilzjI5zbo5z7nbnXEfgGYIBw1+b2ZkVWbGZjQc+AdqZ2RIzuwT4q5nNNbMvgW7AjX5b84CJwNfAO8A1zrlC31pzLZAJfANM9MsC3AbcZGYLCMbsPLOT+y4ilemFF+Cmm+Css2gxdTIjzupImr+lQ1pqisbliEgoymvRAcDMGhO07nQgGPS7oiIrds4NLKW4zGTEOTcSGFlK+VRgainlCwnOyhKRWNC7N9x6K4wYATVrkt4pjfROaWRlZXHdoK5hRyciCarMFh0z+7OZvQO8StCFdI5zrodzbkaVRScisW/qVNi8GRo3hgcfhOTEGmgsIrGtvNPLnyYYL/Mr0At42symFD+qJDoRiW1PPAGnnRaMzRERiUHldV11q7IoRKT6eeopuPLKING54YawoxERKVV597p638w6AgcA85xz31RZVCIS28aOhcsvh1NPhddeg9q1w45IRKRU5Y3RuYvgTKg/Af8xs8uqLCoRiV1r1sAttwSDj5XkiEiMK6/rqj/Q0Tm30cwaEpz2/VTVhCUiMatBA/jf/2D//aFOnbCjEREpV3mJzibn3EYA59wqMwvjvlgiErKM7FxGZ+Zw1EdTOXDLWprffw/pnXThPxGpHspLdPaPOLvKgDaRZ1s55yp04UARqZ4ysnO54/W5bNhcyJlfv8/D/xnDpy0P5fJJswF08T8RqRbKS3T6lnj+UDQDEZHYkZGdy+BJcygodJz2zf8Y89bDfNaiPZf86S7yC4O7kyvREZHqoNyzrgDMrA7BmVcAC5xzv1VFYCISnqGTv6Sg0NHn2w/5x5ujmZV2MH8+627yawVjcpbm5e9gDSIisaG8s65qmtlfCW77MA54HvjJzP5qZrr0qUgcysjO5eA732ZjQREAe27awOdpB3Px2fewsVbK1uWap6aUtQoRkZhSXtfVaGAPoLVz7lcAM9uToAvrIeAv0Q9PRKrKsIy5vDRjMQ7Y87f1rKtTn1cO78WrHU6hqEbS1uUMGNyrXWhxiojsjPLOpDoduKw4yQFwzq0DrgJOjXZgIlJ1MrJzedEnOd0XfMqHj/+ZI5d8DbBdkgMwqEsrjc8RkWqjvBYd55xzpRQWmtnvykWk+ho6+UsAun3/GY+/Poqvm7bmu8b7/m65v/fvqCRHRKqV8lp0vjazC0oWmtl5wLfRC0lEqlJGdi4bC4ro+v0s/v36SL5tsh8XnHMfv9aut91y56klR0SqofJadK4BJpvZn4HPfVlnIAX4Y7QDE5GqMTozh4NXLOSJ10fyXaN9Oa//CNbVqb91fg2Dc49pxYj0DiFGKSKya8o7vTwXOMbMTgaKL4M61Tn3bpVEJiJVYmlePssa7csjx/XnhU6nbZfkpCTX4Jv7+oQYnYjI7imvRQcA59x7wHtVEIuIVLX336cT6/iixp48etyA380e1e+wEIISEak8un+VSKKaPh369OGxT54lJXn7M6sMjckRkfigREckEb3/Ppx+Ouy/P81efZFR/TqQlpqCAWmpKYzp31FjckQkLuyw60pE4sz//genngr77gvvvgtNmpDeRDfpFJH4pBYdkUTiHAwbBq1awXvvQdOmYUckIhJVatERSSRm8PrrsGkT7LNP2NGIiESdWnREEsGMGTBwIPz2G+y9NzRrFnZEIiJVQi06IvFu5kzo1QsaN4Y1a5TkiEhCUYuOSDybNQt69oRGjYLTyZXkiEiCUaIjEq+++AJ69IAGDYIkp2XLsCMSEalySnRE4pUZtGkTJDmtWoUdjYhIKDRGRyTerFgBTZpAp07w2WdBwiMikqDUoiMST+bOhfbtYcyY4LmSHBFJcGrREanGMrJzGZ2Zw9K8fI7ftJxnxt1K7bp14Iwzwg5NRCQmKNERqaaGZczlxRmLAThg5WLGjB9KXo0afPnYeHoccEDI0YmIxAZ1XYlUQxnZuVuTnLqb83nxlWE4MwYOuJ/Bs/NDjk5EJHaoRUekGrpnyryt0xtrpTCy2yV83WR/FjZsAfkFIUYmIhJblOiIVDPDMuaSl1/AfqtzabF2BR+27sSb7U8KOywRkZikREekGikel7PvmqWMHz8UZ0a3y59kU81aW5dpUDc5xAhFRGKLEh2RaiIjO5eXZiymZd7PjB8/lNqFBQwceP92SU4Ng7vPOCTEKEVEYosSHZFqYnRmDmlrlzN+/BBSCjZx7sCR5DTeb+t8M/jbOR1J75QWXpAiIjFGZ12JVBNL8/I5Z87/UX9zPucNGME3TfbfOs+AMUpyRER+Ry06ItVE89QUxvxhEBMP78mSvZpuN29Ql1ZKckRESqEWHZFYt2QJdO/O3YfUoU6t5O2SHAPO69KKEekdwotPRCSGKdERiWW5udC1K8yaRc99khnVrwNpqSkYkJaawpj+HZXkiIiUQ11XIrFq6VLo1i24G/n//R8cdRTpoC4qEZGdoERHJBb9/HOQ5CxbBpmZ0KVL2BGJiFRL6roSiUV16kDLlvDOO3DccWFHIyJSbalFRySW/PIL1K8PqakwbVpwcRwREdllatERiRUrVgTdVQMGBM+V5IiI7Da16IjEgpUr4ZRTYOFCeOSRsKMREYkbSnREwrZqFXTvDvPnw1tvBa06IiJSKZToiITt3HMhJwemTAkSHhERqTRKdETC9vDDwTVzevYMOxIRkbijwcgiYcjLg3/9C5yDQw9VkiMiEiVKdESqWl5ekNjccAN8+23Y0YiIxDV1XYlUgYzsXEZn5rBu+SomvHY3By9bQI3Jr8HBB4cdmohIXFOiIxJlwzLm8tKMxdTdtJHnJ97FgT/P5/o/DeWUFkeQHnZwIiJxTl1XIlGUl1/ASzMW44Ajc7/h0OXfc92Zt/LW/scwOjMn7PBEROKeWnREomj52t9wzsCMD/Y/kj9c8TQr9mgIwNK8/JCjExGJf1Fr0TGzsWa2wsy+iijb28ymmdl8/7eBLzcz+6eZLTCzL83siIjXXOiXn29mF0aUH2lmc/1r/mmm6+VL7MjIzuX4B96jaEPQXdVj/gyArUkOQPPUlLDCExFJGNHsunoO6F2i7HbgXedcW+Bd/xygD9DWPy4HHocgMQLuBo4BjgbuLk6O/DKXRbyu5LZEQpGRncuQyXNZ9csaTv/bSI7/cQ4pBZu2W8aAwb3ahROgiEgCiVqi45z7AFhdorgvMM5Pj4OtYzH7As+7wAwg1cyaAb2Aac651c65NcA0oLeft6dzboZzzgHPR6xLJFSjM3Mo2riRp1+7j7Rv5nHTaTcypf1JW+cbMKhLK9I7pYUXpIhIgqjqMTpNnXPL/PTPQFM/nQb8FLHcEl9WXvmSUspLZWaXE7QU0bRpU7KysnZ9Dypo/fr1VbKdWJaodXDuPqs546URtFz8JTOvvZ4DjvkDN7MFgFpJNWi6Vx1SU1YlTN0k6nEQKdHrINH3H1QHEF4dhDYY2TnnzMxV0baeBJ4E6Ny5s+vatWvUt5mVlUVVbCeWJWod3DHqXWrX2Y9v+3Rj32O68vDc4GOWlprCR7efHHJ0VS9Rj4NIiV4Hib7/oDqA8OqgqhOd5WbWzDm3zHc/rfDluUDLiOVa+LJcoGuJ8ixf3qKU5UXCs2kT/Pwzg3sfxJCNV5JfULi1JSclOUljckREQlDV19GZAhSfOXUh8EZE+QX+7KsuwFrfxZUJ9DSzBn4Qck8g089bZ2Zd/NlWF0SsS6TqbdoEf/oTHH886W33YlS/DqT5s6rSUlMY1a+DxuSIiIQgai06ZjaeoDWmkZktITh76gFgopldAvwInOMXnwqcCiwANgIXAzjnVpvZfcBnfrl7nXPFA5yvJjizKwV42z9Eqt7mzXD22fCf/8ATT0D9+qR3qk96pzSysrK4blDXsCMUEUlYUUt0nHMDy5jVvZRlHXBNGesZC4wtpXwWcOjuxCiy2zZvhnPOgTffDO5GfvnlYUckIiIRdAsIkd0xciS88QY8+ihcdVXY0YiISAm6BYTI7rjlFjjkkKBVR0REYo5adER21pYtQUvO+vWwxx5KckREYpgSHZGdsWULnH8+DBsWjMsREZGYpkRHpKIKC+HCC2HCBHjwQRhY1nh7ERGJFUp0RCqisBAuvhhefhnuvx9uvTXsiEREpAKU6IhUxNKlMG0a3HcfDBkSdjQiIlJBOutKpDxFRWAGLVvCV19Bw4ZhRyQiIjtBLToiZSkqgiuvhJtvBueU5IiIVENKdEQiDMuYS5shU2l925u8fOTp8NRTULdu2GGJiMguUqIj4g3LmMuLMxZTWFTEvdP+zbmz3+ZfXc5i2JHnBN1XIiJS7SjREQEysnN5ccZiAIa99zTnZ0/l30f3468nXsj4mUtCjk5ERHaVBiNLwsvIzmXI5Llbn3/W8hC2JNXkgZMuAjMKnQsvOBER2S1KdCShZWTncvPEORQWFXHwL4v4psn+ZB54HJkHHrd1mSR1W4mIVFvqupKENSxjLje+MpvCoiKGZD3Lm8/dwCE/L/jdcgOPaRlCdCIiUhmU6EhCysjO5aUZi3HOcdv747hi5mRe6tSHeU3bbLfceV1aMSK9Q0hRiojI7lLXlSSk0Zk5OOcY/MHzXPXpJF7s2Ie7T7ly69lVKclJjOrXgfROaSFHKiIiu0MtOpKQlublc9KiL7hmxqu8fHhv7ux51dYkJ8lMSY6ISJxQi44kpOapKbzf+giuTB9C5oHH4izI+Q14+JzDleSIiMQJtehI4hkzhuFta5BSqybvtDt+uyRnUJdWSnJEROKIEh1JLCNHwk03ccqHbzCqXwfSUlMwIC01hTH9O2rgsYhInFHXlSSOBx6AYcPgvPPgoYdIT0pS642ISJxTi44khtGjYcgQOPdceO45SEoKOyIREakCSnQk/m3ZAm+/DQMGwLhxSnJERBKIuq4kvm3ZAjVrwltvQa1awbSIiCQMtehI/HrkETjpJPj1V6hbV0mOiEgCUqIj8emxx+D666FpU6hTJ+xoREQkJEp0JP78+99w7bXQty9MmADJyWFHJCIiIVGiI/HlhRfgqqvg9NNh4sRgXI6IiCQsJToSX449Fi65BCZNUpIjIiJKdCROzJgBzsEBB8DTT0Pt2mFHJCIiMUCJjlR/zz8Pxx0Hjz8ediQiIhJjlOhI9fbii3DRRdC9O1x8cdjRiIhIjNGFRaRaycjO5Z4p88jLL+DMr7MY89bfWN35OBq/8QakpIQdnoiIxBi16Ei1MSxjLje8Mpu8/AIabVjDg28/wsyWh9D9pJvIyFkTdngiIhKD1KIj1cKwjLm8OGPx1ucr6zXgwnOG81XTNmxMqs3ozBzdiVxERH5HiY7EvIzs3K1JTu+cj6hVWMCU9l2Z2fLQrcsszcsPKzwREYlh6rqSmDf8zXkA9PruYx6Z8lfOy56KuaLtlmmeqvE5IiLye0p0JKZlZOeyZmMBPebP4NE3HmTuPgfw57Puwdm2Qzc5yRjcq12IUYqISKxSoiMxKyM7lyGT59J9wac8lvEA85q24cJz7mV97bpblzFg9FmHa3yOiIiUSomOxKzRmTnkFxTS4ecFfNOkNReccy+/1q633TJj+ndUkiMiImXSYGSJWStXroWatfj78efy+DFnsSl5+9s6nNellZIcEREpl1p0JDZlZvLB01fQ9pcfwWy7JCfJjL/378iI9A4hBigiItWBEh2JPf/9L6SnU6dpY9anNtpuVkpyEg+fozE5IiJSMUp0JLa8+y6ccQa0bcteH73PbecdT1pqCgakpaYwql8HJTkiIlJhGqMjsWPWrCDJOeCAIOFp1Ij0RiixERGRXaYWHYkdhxwCf/5zkOQ0bhx2NCIiEgeU6Ej4PvsM8vKCu48/+ig0aRJ2RCIiEieU6Ei4PvoITj4Zrr467EhERCQOKdGR8MyYAX36QLNm8NBDYUcjIiJxSImOhOPTT6FXL2jaFKZPh+bNw45IRETikBIdqXpFRXDFFcGA4+nTIU1nVYmISHTo9HKpejVqwBtvBH9btAg7GhERiWNq0ZGqk50NN98ctOjsuy+0bBl2RCIiEueU6EjVmDMHTjkFJk2CX34JOxoREUkQSnQk+r78Erp3h3r1gjE5TZuGHZGIiCSIUBIdM/vBzOaa2Wwzm+XL9jazaWY23/9t4MvNzP5pZgvM7EszOyJiPRf65eeb2YVh7IvswFdfBUlOnTpBkrP//mFHJCIiCSTMFp1uzrmOzrnO/vntwLvOubbAu/45QB+grX9cDjwOQWIE3A0cAxwN3F2cHEkMWbYM9toLsrKgTZuwoxERkQQTS2dd9QW6+ulxQBZwmy9/3jnngBlmlmpmzfyy05xzqwHMbBrQGxhftWFLsYzsXO6ZMo+8/ALqbs7nuqNqkXfQkaR/8w0kJ4cdnoiIJKCwWnQc8H9m9rmZXe7Lmjrnlvnpn4HigRxpwE8Rr13iy8oqlxBkZOcy+NU55OUXsP+qJbz31BW0+WA6gyfNIeOrFWGHJyIiCSqsFp0TnHO5ZtYEmGZm30bOdM45M3OVtTGfTF0O0LRpU7Kysipr1WVav359lWwnViz/+VeuP6SI1GW5/PGJYdSoUUThwW24vmEBy3O+IGvt/LBDDEWiHQelUR2oDhJ9/0F1AOHVQSiJjnMu1/9dYWavE4yxWW5mzZxzy3zXVHEzQC4QecGVFr4sl21dXcXlWWVs70ngSYDOnTu7rl27lrZYpcrKyqIqthMrLr79P+y7ehkTxt/Jb0VFDBg4ir4Nm/Pw3JoYsOiBrmGHGIpEOw5KozpQHST6/oPqAMKrgyrvujKzema2R/E00BP4CpgCFJ85dSHwhp+eAlzgz77qAqz1XVyZQE8za+AHIff0ZRKCdrW3MH78UJILt3DugJEsaNRq67zmqSkhRiYiIoksjBadpsDrZla8/Zedc++Y2WfARDO7BPgROMcvPxU4FVgAbAQuBnDOrTaz+4DP/HL3Fg9Mlqp3ZXpnnvswnff37ch3jffbWp6cZAzu1S68wEREJKFVeaLjnFsIHF5K+SqgeynlDrimjHWNBcZWdoyyE374AdasIb1TJ3hoOK9MmQf5BQDUrGGMPutw0jtpjLiIiIQjlk4vl+rmxx+ha1eoWRO+/Zb0TmnbJTVZWVl0VZIjIiIhUqIju2bxYujWDdauhf/+N0h2REREYox+nWTnLVkSJDmrV8O0aXDkkWFHJCIiUiolOrLzRo2ClSuDJOeoo8KORkREpEy6e7nsvL/9DT76CI4+OuxIREREyqVERyrm55+hf39YtQpq14ZDDw07IhERkR1S15Xs2PLlwZicn36CBQugYcOwIxIREakQJTpSvhUr4OSTg7Os3n4bjjkm7IhEREQqTImOlO2XX6B7d1i0CKZOhRNPDDsiERGRnaIxOlK2334DM3jrreDCgCIiItWMWnTk99auhfr1oWVLyM6GpKSwIxIREdklatGR7a1eHQw8vuqq4LmSHBERqcaU6AgZ2bkc/8B7HH7DK3x72LEUfjUP+vULOywREZHdpkQnwWVk5zJk8lx+/fkXxk28k/2XLeSafneQ0bRD2KGJiIjsNiU6CW50Zg75m7fw79fvp/3yRVz1xyG8s9+RjM7MCTs0ERGR3abByAluaV4+mPHIcQOov3kj7x5wzLZyERGRak6JToLJyM5l+JvzWLOxgHqbNtLzx9lkHngcn+x72HbLNU9NCSlCERGRyqOuqwSSkZ3Lza/OYc3GAupuzufZSffwyBt/pcXa5dstl5KcxOBe7UKKUkREpPKoRSeBDH9zHoVFjpTNv/HspOEckfst1595K0v2akqSGUXO0Tw1hcG92pHeKS3scEVERHabEp0EsmZjAXUKfmPsa8PpvORr/nLGLUw96AQAipxj0QOnhRyhiIhI5VLXVYLpMf9Tjv5pHjeefhNvHbzt3lUakyMiIvFILToJJDUlmTfbn8TXTfbn+0Ytt5unMTkiIhKP1KKTCH77Dc45h3+03kRyDftdknNel1YakyMiInFJiU68++03+OMfYdIkTmINo88+nLTUFAxIS03h7/07MiJdV0EWEZH4pK6reLZpE/zpT/DOO/D003DBBaSDWm9ERCRhqEUnXm3eDGefDVOnwhNPwCWXhB2RiIhIlVOiE8+Sk+Hxx+Hyy8OOREREJBTquoo3BQWwbh00bAiTJoFZ2BGJiIiERi068aSgAAYMgK5dg0HISnJERCTBKdGJF1u2wKBBMHlyMB6nTp2wIxIREQmdEp14sGULnH8+vPoqPPQQ3HBD2BGJiIjEBCU68eDOO2HCBHjwQbj55rCjERERiRkajBwPrr8eWrfW2VUiIiIlqEWnuiosDK6Ps2ULNGumJEdERKQUSnSqo6IiuPRSuPJKeOONsKMRERGJWUp0qpuioqD15rnn4O67g1s8iIiISKmU6FQnRUVw1VXwzDMwbFiQ6IiIiEiZNBg5xmVk5zI6M4elefkcU7CSF158ieQhQ+Dee3VBQBERkR1Qi04My8jOZcjkueTm5eOAGcmN6HPxo2ScdbWSHBERkQpQohPDRmfmkL95C3f/9wkumjUFgAX1GzP6/74LOTIREZHqQYlODFu6ZiN3vfsUF3/+JmnrVmwrz8sPMSoREZHqQ2N0YpVz3P/xOAZ+PoWxR57JyG6XbJ3VPDUlxMBERESqDyU6ser22xn44SReOOpM7u122dYxOSnJSQzu1S7k4ERERKoHdV3FqrQ0uOYa9vj3Y6Q1qIsBaakpjOrXgfROaWFHJyIiUi2oRSeWOAeLF8O++wb3r3KOdDPSj2gRdmQiIiLVklp0Ysndd8Ohh0JOTvBcp5CLiIjsFiU6sWL4cLjvPujfH9q2DTsaERGRuKBEJxaMGAH33AMXXQRPPgk19LaIiIhUBv2ihi0jA+68E84/H55+WkmOiIhIJdKvathOPx0efxyefRaSksKORkREJK4o0QnLM8/AsmVQsyZceaWSHBERkShQohOGMWPg0kvhb38LOxIREZG4pkSnqv3jH3DTTXDWWTBqVNjRiIiIxDUlOlXp0UfhhhugXz94+eWg20pERESiRolOVdm0CZ54Avr2hfHjITk57IhERETinpoUqoJzULs2ZGXBHntArVphRyQiIpIQ1KITbU8/DWefDZs3Q8OGSnJERESqULVPdMyst5nlmNkCM7s97Hi2M3YsXHYZbNwYtOqIiIhIlarWiY6ZJQGPAX2A9sBAM2sfblSBfd55JziFvFcvmDw56LoSERGRKlWtEx3gaGCBc26hc24zMAHoG3JM8NJLtPvrX+GUU+D116FOnbAjEhERSUjmqnGXipmdBfR2zl3qn58PHOOcu7bEcpcDlwM0bdr0yAkTJkQ1rj2+/prmL73E/LvuoiiBW3LWr19P/fr1ww4jVKoD1QGoDhJ9/0F1ANGtg27dun3unOtc2ryEOOvKOfck8CRA586dXdeuXSt1/RnZuYzOzGFpXj7NU1MY3Ksvv7ZvT2Vvp7rJyspSHagOVAeoDhJ9/0F1AOHVQXXvusoFWkY8b+HLqkxGdi5DJs8lNy8fB+Tm5TNk8lzy8guqMgwREREpRXVPdD4D2ppZazOrBQwAplRlAKMzc8gvKNyuLL+gkOVrf6vKMERERKQU1brryjm3xcyuBTKBJGCsc25eVcawNC+/1PLNhUVVGYaIiIiUolonOgDOuanA1LC23zw1hdxSkp1aSdW9sUxERKT606/xbhrcqx0pyUnblaUkJ9F0L51SLiIiErZq36ITtvROaQAlzrpqR+ra+SFHJiIiIkp0KkF6p7StCU+xrCwlOiIiImFT15WIiIjELSU6IiIiEreU6IiIiEjcUqIjIiIicUuJjoiIiMQtJToiIiISt5ToiIiISNxSoiMiIiJxS4mOiIiIxC0lOiIiIhK3lOiIiIhI3FKiIyIiInFLiY6IiIjELSU6IiIiEreU6IiIiEjcUqIjIiIicUuJjoiIiMQtJToiIiISt8w5F3YMVcrMfgF+rIJNNQJWVsF2YpnqQHUAqgNQHST6/oPqAKJbB/s65xqXNiPhEp2qYmaznHOdw44jTKoD1QGoDkB1kOj7D6oDCK8O1HUlIiIicUuJjoiIiMQtJTrR82TYAcQA1YHqAFQHoDpI9P0H1QGEVAcaoyMiIiJxSy06IiIiEreU6FQyM+ttZjlmtsDMbg87nspmZj+Y2Vwzm21ms3zZ3mY2zczm+78NfLmZ2T99XXxpZkdErOdCv/x8M7swrP2pCDMba2YrzOyriLJK22czO9LX6QL/WqvaPdyxMurgHjPL9cfCbDM7NWLeEL8/OWbWK6K81M+HmbU2s099+StmVqvq9q5izKylmU03s6/NbJ6Z/cWXJ8SxUM7+J8xxYGZ1zGymmc3xdTDcl5cat5nV9s8X+Pn7Raxrp+omVpRTB8+Z2aKI46CjLw//c+Cc06OSHkAS8D2wP1ALmAO0DzuuSt7HH4BGJcr+Ctzup28HHvTTpwJvAwZ0AT715XsDC/3fBn66Qdj7Vs4+nwgcAXwVjX0GZvplzb+2T9j7XME6uAe4pZRl2/tjvzbQ2n8mksr7fAATgQF++t/AVWHvcyn71Qw4wk/vAXzn9zUhjoVy9j9hjgP/vtT308nAp/79KjVu4Grg3356APDKrtZNrDzKqYPngLNKWT70z4FadCrX0cAC59xC59xmYALQN+SYqkJfYJyfHgekR5Q/7wIzgFQzawb0AqY551Y759YA04DeVRxzhTnnPgBWlyiulH328/Z0zs1wwSf8+Yh1xYwy6qAsfYEJzrlNzrlFwAKCz0apnw//39rJwCT/+sj6jBnOuWXOuS/89K/AN0AaCXIslLP/ZYm748C/l+v902T/cJQdd+SxMQno7vdzp+omunu1c8qpg7KE/jlQolO50oCfIp4vofwvgurIAf9nZp+b2eW+rKlzbpmf/hlo6qfLqo94qKfK2uc0P12yvLq41jdHjy3usmHn66AhkOec21KiPGb5LohOBP/NJtyxUGL/IYGOAzNLMrPZwAqCH+fvKTvurfvq568l2M9q/d1Ysg6cc8XHwUh/HIwxs9q+LPTPgRId2VknOOeOAPoA15jZiZEzfQaeUKfyJeI+e48DbYCOwDLg4VCjqSJmVh94DbjBObcucl4iHAul7H9CHQfOuULnXEegBUELzEHhRlT1StaBmR0KDCGoi6MIuqNuCy/C7SnRqVy5QMuI5y18WdxwzuX6vyuA1wk+6Mt9cyP+7wq/eFn1EQ/1VFn7nOunS5bHPOfccv+FVwQ8RXAswM7XwSqC5uyaJcpjjpklE/zIv+Scm+yLE+ZYKG3/E/E4AHDO5QHTgWMpO+6t++rn70Wwn3Hx3RhRB71916Zzzm0CnmXXj4NK/xwo0alcnwFt/Qj8WgSDz6aEHFOlMbN6ZrZH8TTQE/iKYB+LR8xfCLzhp6cAF/hR912Atb6JPxPoaWYNfDN3T19WnVTKPvt568ysi++7vyBiXTGt+Mfd+yPBsQBBHQzwZ5y0BtoSDC4s9fPhW0GmA2f510fWZ8zw788zwDfOub9FzEqIY6Gs/U+k48DMGptZqp9OAXoQjFUqK+7IY+Ms4D2/nztVN1HfsZ1QRh18G5HsG8GYmsjjINzPQWkjlPXYrRHppxKcjfA9cEfY8VTyvu1PcBbAHGBe8f4R9Dm/C8wH/gvs7csNeMzXxVygc8S6/kwwAG8BcHHY+7aD/R5P0CRfQNBffEll7jPQ2X8pfA88ir+QZyw9yqiDF/w+fknwZdYsYvk7/P7kEHHGRFmfD39szfR18ypQO+x9LqUOTiDolvoSmO0fpybKsVDO/ifMcQAcBmT7ff0KuKu8uIE6/vkCP3//Xa2bWHmUUwfv+ePgK+BFtp2ZFfrnQFdGFhERkbilrisRERGJW0p0REREJG4p0REREZG4pURHRERE4pYSHREREYlbSnREEpCZFVpwh+GvzOzN4uti7OQ6upqZM7MzIsreMrOuO3jdRWbWvIx5Je+A/PHOxrUryoupnNfsZxF3c48o72pmb5Uoe87MzvLTN5hZ3Yh5UyOuS7Le/21uZpP8dEeLuCO4iOwcJToiiSnfOdfROXcowc06r9nF9SwhuB7IzrgIKC+pGOxj6+icO24X49pZO4qpMt0AbE10nHOnuuAKs0SULXXOFV+AriPBtVVEZBco0RGRT/A3zTOzNmb2jgU3bf2fmR3ky8/2rT9zzOyDiNfOAdaaWY+SKzWzI83sfb+uTDNr5ls1OgMv+RablIoEaGb/MLO7/HQvM/vAzGr4lpJ/m9ksM/vOzE73yySZ2Wgz+8yCmwxeEbGu28xsrt+XB0qLqbTYI/ZpjpnNYReSQzO7niChmm5m033ZD2bWqMRy+/n6rgXcC/T3sfU3s/lm1tgvV8PMFhQ/F5FShH2VRT300KPqH8B6/zeJ4Mqtvf3zd4G2fvoYgkvWQ3BF0zQ/ner/dgXeAk4E3vdlb/nyZOBjoLEv7w+M9dNZRFwdtURczwGL2Hbl3Zd8eV2Cq3F3I7iSbJuI5d8h+KetLUELUx3gcmCYX6Y2MAtoTXAz2o+Bun7e3iVj2kHsXwIn+unRwFel7ENX4K1S9ussP/0D0Chi3tbnEe/LfsXrJmhtejRi+bsJbqgJwWXzXwv7eNJDj1h+FN+ETEQSS4qZzSZoyfkGmGbBXamPA14NbjEDBEkCwEfAc2Y2EZgcuSLn3AdmhpmdEFHcDjjUrxeChGpZBWMb7JybVGIbG83sMuAD4Ebn3PcRsye64IaS881sIcEdlHsChxWPiyG4mWJb4BTgWefcRr/e1aVsv9TY/TiaVOdccYvWCwSJU0llXW6+si5DP5bg3j9/J7iE/rOVtF6RuKRERyQx5TvnOvpBsZkE3TDPAXnOuY4lF3bOXWlmxwCnAZ+b2ZElFhkJDAO2+OcGzHPOHVuJMXcguPNzybE0JRMI57d/nXNuu5vFmlmvCmyn1Nh3YsD2KqBBibK9gZUVfH25nHM/mdlyMzuZ4A7RgypjvSLxSmN0RBKYb9m4HrgZ2AgsMrOzIbgLsZkd7qfbOOc+dc7dBfwCtCyxnv8j+HE/zBflAI3N7Fj/+mQzO8TP+xXYY2fiNLN9fYydgD4+6Sp2th+r0obg5oo5BMnbVWaW7F9/oJnVA6YBFxef9WRme5cSU6mxu2DAcF5Ey1VZCcZ8oLmZHRwR++EEXXG7sv+lLf80wY0TX3XOFe7EukQSjhIdkQTnnCu+E/FAgh/vS/xg23lAX7/YaD+A9yuC8StzSlnVSHwC5JzbDJwFPOjXNZugWwyClqN/lzMYebRtO718tpnVBp4BbnHOLSW4c/rTZlbHL7+Y4M7QbwNXOud+I0gEvga+8DE/AdR0zr1DcIftWb7r7paSMRF0VZUV+8XAY365rf17JepzE3Ae8KxfbhJwqXNurV/kSeCd4sHIFTAdaF88GNmXTQHqo24rkR3S3ctFpNoys+cIBv5O2tGy8cTMOgNjnHN/CDsWkVinMToiItWImd0OXIXG5ohUiFp0REREJG5pjI6IiIjELSU6IiIiEreU6IiIiEjcUqIjIiIicUuJjoiIiMQtJToiIiISt/4fGEVNTiB1QaEAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 576x576 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(8,8))\n",
    "plt.scatter(resnetUtil, opmUtil)\n",
    "x = [0, 35000]  # Adjust the range based on your needs\n",
    "y = [0, 35000]  # Adjust the range based on your needs\n",
    "\n",
    "# Plot the red dotted diagonal line\n",
    "plt.plot(x, y, 'r--')\n",
    "plt.grid(True)\n",
    "plt.title(\"Expected utility comparison between ResNet and OPM\")\n",
    "plt.xlabel(\"ResNet Expected Utility\")\n",
    "plt.ylabel(\"OPM Expected Utility\")\n",
    "plt.tight_layout(True)\n",
    "# plt.legend([\"Diagonal dotted line\", \"Scatter plot\"])\n",
    "plt.savefig('demo.JPEG', transparent=False)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[1744,   38],\n",
       "       [   6, 1712]])"
      ]
     },
     "execution_count": 56,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "CM"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[176144,   3838],\n",
       "       [   606, 172912]])"
      ]
     },
     "execution_count": 57,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "cm_rn"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Slightly above threshold(p)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {},
   "outputs": [],
   "source": [
    "cm_rn = np.array([[0, 0], [0, 0]]) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {},
   "outputs": [],
   "source": [
    "def p_threshold(U_AA, U_AB, U_BA, p):\n",
    "    p = U_AB + (U_AA - U_BA)*((p)/(1-p))\n",
    "    return [[U_AA, U_AB], [U_BA, p]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {},
   "outputs": [],
   "source": [
    "u_above = np.array(p_threshold(1.,0., 0, 0.98)) #utility matrix using slightly above threshold(p)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 1.,  0.],\n",
       "       [ 0., 49.]])"
      ]
     },
     "execution_count": 61,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "u_above"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Expected utilities with resnet alone: 85631.99999999993\n",
      "Expected utilities with OPM: 85749.99999999999, total number over overrides 1782\n"
     ]
    }
   ],
   "source": [
    "print(f'Expected utilities with resnet alone: {np.multiply(CM, u_above).sum()}')\n",
    "expected_opm = func_expected_opm(test_df, u_above, cm_rn)\n",
    "print(f\"Expected utilities with OPM: {expected_opm[0]}, total number over overrides {expected_opm[1]}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[1744,   38],\n",
       "       [   6, 1712]])"
      ]
     },
     "execution_count": 63,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "CM"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[   0,    0],\n",
       "       [1750, 1750]])"
      ]
     },
     "execution_count": 64,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "cm_rn #confusion matrix "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Slightly below threshold(p)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "metadata": {},
   "outputs": [],
   "source": [
    "cm_rn = np.array([[0, 0], [0, 0]]) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "metadata": {},
   "outputs": [],
   "source": [
    "u_below = np.array(p_threshold(1.,-1 , 1, 0.95)) #utility matrix using sligthly below threshold(p)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 1., -1.],\n",
       "       [ 1., -1.]])"
      ]
     },
     "execution_count": 67,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "u_below"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 68,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Expected utilities with resnet alone: 0.0\n",
      "Expected utilities with OPM: 0.0, total number over overrides 1718\n"
     ]
    }
   ],
   "source": [
    "print(f'Expected utilities with resnet alone: {np.multiply(CM, u_below).sum()}')\n",
    "expected_opm = func_expected_opm(test_df, u_below, cm_rn)\n",
    "print(f\"Expected utilities with OPM: {expected_opm[0]}, total number over overrides {expected_opm[1]}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 69,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[1744,   38],\n",
       "       [   6, 1712]])"
      ]
     },
     "execution_count": 69,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "CM"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 70,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[1750, 1750],\n",
       "       [   0,    0]])"
      ]
     },
     "execution_count": 70,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "cm_rn"
   ]
  }
 ],
 "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.9.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}