df3_NAO_reliability_diagrams_1980-2010.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "e309f9af",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'\\nb=109\\ntime=109\\ntotal=109*51\\n'"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from __future__ import (absolute_import, division, print_function)\n",
    "from six.moves import (filter, input, map, range, zip)  # noqa\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from sklearn.linear_model import LinearRegression\n",
    "import warnings\n",
    "import copy\n",
    "\n",
    "a=0\n",
    "\n",
    "total=31*40 # 31 years, 40 ensembles\n",
    "time=31\n",
    "b=a+31\n",
    "s='1980-2010'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "09023260",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\hss84546\\AppData\\Local\\Temp\\ipykernel_14196\\469063573.py:76: RuntimeWarning: invalid value encountered in double_scalars\n",
      "  y_co[n-1]= E[1]/E[2]\n",
      "C:\\Users\\hss84546\\AppData\\Local\\Temp\\ipykernel_14196\\469063573.py:77: RuntimeWarning: invalid value encountered in double_scalars\n",
      "  x_co[n-1]=F/D\n"
     ]
    }
   ],
   "source": [
    "def lower_tercile_dp3(a,b,time,total,data,data_obs,label):\n",
    "\n",
    "    PCs_data_full = data\n",
    "    PCs_full=PCs_data_full.reshape(40, 40)\n",
    "    PCs=PCs_full[:,a:b]\n",
    "    PCs_data = PCs.reshape(40*time)\n",
    "    PCs_data_copy = copy.copy(PCs_data)\n",
    "    val= np.percentile(PCs_data, 33, axis=None, out=None, overwrite_input=False, interpolation='linear', keepdims=False)    \n",
    "\n",
    "    for m in range(total):\n",
    "        if PCs_data[m] < val :\n",
    "            PCs_data_copy[m]=1\n",
    "        else :\n",
    "            PCs_data_copy[m]=0\n",
    "\n",
    "        #print(ens_1o)\n",
    "\n",
    "    PCs_binary=PCs_data_copy.reshape(40, time)\n",
    "\n",
    "    #probabilities of warm winter for each year\n",
    "    binary_matrix=PCs_binary\n",
    "    #sum columns to give number of ensembles that say warm event for that year\n",
    "    column_sums = binary_matrix.sum(axis=0)\n",
    "\n",
    "    probabilities= column_sums/40\n",
    "\n",
    "    prob= np.zeros(time*3).reshape((time,3))\n",
    "    prob[:,1]= probabilities\n",
    "    prob[:,0]=np.arange(time)\n",
    "    prob_copy = copy.copy(prob[:])\n",
    "\n",
    "    bins=np.arange(0,1.1,0.1)\n",
    "    #to catagorize probabilities into correct bins\n",
    "    for n in range(1,11):\n",
    "        for p in range(time):\n",
    "            if bins[n-1]<=prob[p,1]< bins[n] :\n",
    "                prob_copy[p,2]=n\n",
    "    \n",
    "    #observed data:\n",
    "    #sort observed data and set lower tercile(warm events) to 1\n",
    "    obsv= np.zeros(time*3).reshape((time,3))\n",
    "    obsv[:,1]= data_obs\n",
    "    obsv[:,0]=np.arange(time)\n",
    "\n",
    "    val_o= np.percentile(data_obs, 33 , axis=None, out=None, overwrite_input=False, interpolation='linear', keepdims=False)    \n",
    "\n",
    "    for n in range(1,time+1):\n",
    "            if obsv[n-1,1] < val_o :\n",
    "                obsv[n-1,1]=1\n",
    "            else :      obsv[n-1,1]=0\n",
    "\n",
    "    #obsv_re_sort is an array of years and corresponding 1's and 0's for cold winters over UK\n",
    "    obsv_bins= np.zeros(time*4).reshape((time,4))\n",
    "\n",
    "    obsv_bins[:,0]=obsv[:,0]\n",
    "    obsv_bins[:,1]=obsv[:,1]\n",
    "    obsv_bins[:,2]=prob_copy[:,2]\n",
    "    obsv_bins[:,3]=prob_copy[:,1]\n",
    "    obsv_bins_sort= obsv_bins[obsv_bins[:, 2].argsort()]\n",
    "    #obsv_bins_sort is array with corresponding bin of year from forecast data\n",
    "\n",
    "    y_co=np.zeros(10)\n",
    "    x_co=np.zeros(10)\n",
    "    C=np.zeros(40).reshape((10,4))\n",
    "\n",
    "    for n in range(1,11):\n",
    "        A=np.zeros(time*4).reshape((time,4))\n",
    "        for p in range(1,time+1):\n",
    "            if n == obsv_bins_sort[p-1,2]:\n",
    "                A[p-1,:]=obsv_bins_sort[p-1,:]\n",
    "        B=A.sum(axis=0)#so B[1] is number of observed occurances (1's) in bin\n",
    "        D=B[2]/n #no. of years in bin\n",
    "        F=B[3] # to avg x coord\n",
    "        E=[n,B[1],D,0]\n",
    "        C[n-1,:]=E\n",
    "        y_co[n-1]= E[1]/E[2]\n",
    "        x_co[n-1]=F/D\n",
    "\n",
    "    coords=np.zeros(40).reshape((10,4))\n",
    "    coords[:,0]=x_co\n",
    "    coords[:,1]=y_co\n",
    "    coords[:,2]=C[:,2]\n",
    "    #print(coords)\n",
    "    coords_new=coords[~np.isnan(coords).any(axis=1)]\n",
    "\n",
    "    size=coords_new[:,2]\n",
    "    sh=coords_new[:,0].size\n",
    "    coords_new[:,3]=sh\n",
    "\n",
    "    a_file = open(\"NAO_\"+s+label+\"_coords_lower_tercile.txt\", \"w\")\n",
    "    for row in coords_new:\n",
    "        np.savetxt(a_file, row)\n",
    "    a_file.close()\n",
    "    \n",
    "    rego= np.zeros(1000)\n",
    "\n",
    "    for i in range(1000):\n",
    "        \n",
    "        random_indices = np.random.choice(time, size=time, replace=True)\n",
    "        rand_sample = obsv_bins_sort[random_indices,:]\n",
    "        #sum columns to give number of ensembles that say warm event for that year\n",
    "        y_co_unc=np.zeros(10)\n",
    "        x_co_unc=np.zeros(10)\n",
    "        C_unc=np.zeros(40).reshape((10,4))\n",
    "\n",
    "        for n in range(1,11):\n",
    "            A_unc=np.zeros(4*time).reshape((time,4))\n",
    "            for p in range(1,time+1):\n",
    "                if n == rand_sample[p-1,2]:\n",
    "                    A_unc[p-1,:]=rand_sample[p-1,:]\n",
    "\n",
    "            B_unc=A_unc.sum(axis=0)#so B[1] is number of observed occurances (1's) in bin\n",
    "            D_unc=B_unc[2]/n #no. of years in bin\n",
    "            F_unc=B_unc[3] #avg x coord\n",
    "            E_unc=[n,B_unc[1],D_unc,0]\n",
    "            C_unc[n-1,:]=E_unc\n",
    "            with warnings.catch_warnings():\n",
    "                warnings.simplefilter('ignore', RuntimeWarning)\n",
    "                y_co_unc[n-1]= E_unc[1]/E_unc[2]\n",
    "                x_co_unc[n-1]=F_unc/D_unc\n",
    "\n",
    "        coords_unc=np.zeros(30).reshape((10,3))\n",
    "        coords_unc[:,0]=x_co_unc\n",
    "        coords_unc[:,1]=y_co_unc\n",
    "        coords_unc[:,2]=C_unc[:,2]\n",
    "        coords_new_unc=coords_unc[~np.isnan(coords_unc).any(axis=1)]\n",
    "\n",
    "        size_unc=coords_new_unc[:,2]\n",
    "        sh_unc=coords_new_unc[:,0].size\n",
    "\n",
    "        p_unc=coords_new_unc[:,0].reshape(-1,1)\n",
    "        q_unc=coords_new_unc[:,1]\n",
    "        model_unc= LinearRegression()\n",
    "        model_unc.fit(p_unc,q_unc,size_unc[0:sh_unc])\n",
    "        p_new_unc = coords_new_unc[:,0]\n",
    "        q_new_unc = model_unc.predict(p_new_unc[:, np.newaxis])\n",
    "        rego[i]=model_unc.coef_\n",
    "\n",
    "    unc_reg_l= np.zeros(2)\n",
    "    unc_reg_l[0]= np.percentile(rego, 12, axis=None, out=None, overwrite_input=False, interpolation='linear', keepdims=False)    \n",
    "    unc_reg_l[1]= np.percentile(rego, 87, axis=None, out=None, overwrite_input=False, interpolation='linear', keepdims=False)    \n",
    "\n",
    "    a_file = open(\"uncert_reg_l\"+s+label+\".txt\", \"w\")\n",
    "    np.savetxt(a_file, unc_reg_l)\n",
    "    a_file.close()\n",
    "\n",
    "\n",
    "data=np.loadtxt(\"dp3_data.txt\")\n",
    "#data_mag=abs(data_mag)\n",
    "obs_dataun=np.loadtxt(\"timeseries_NAO_ERA20C_1901-2010.txt\")\n",
    "data_obsmag=obs_dataun[78:110]\n",
    "data_obs=data_obsmag\n",
    "label='dp3'\n",
    "lower_tercile_dp3(a,b,time,total,data,data_obs, label)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "2df034bb",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\hss84546\\AppData\\Local\\Temp\\ipykernel_14196\\2039053944.py:75: RuntimeWarning: invalid value encountered in double_scalars\n",
      "  y_co[n-1]= E[1]/E[2]\n",
      "C:\\Users\\hss84546\\AppData\\Local\\Temp\\ipykernel_14196\\2039053944.py:76: RuntimeWarning: invalid value encountered in double_scalars\n",
      "  x_co[n-1]=F/D\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 0.1375      0.          2.          0.        ]\n",
      " [ 0.26153846  0.15384615 13.          0.        ]\n",
      " [ 0.35714286  0.28571429  7.          0.        ]\n",
      " [ 0.42857143  0.57142857  7.          0.        ]\n",
      " [ 0.525       1.          2.          0.        ]]\n"
     ]
    }
   ],
   "source": [
    "def upper_tercile_dp3(a,b,time,total,data,data_obs,label):\n",
    "\n",
    "    PCs_data_full = data\n",
    "    PCs_full=PCs_data_full.reshape(40, 40)\n",
    "    PCs=PCs_full[:,a:b]\n",
    "    PCs_data = PCs.reshape(40*time)\n",
    "    PCs_data_copy = copy.copy(PCs_data)\n",
    "    val= np.percentile(PCs_data, 67, axis=None, out=None, overwrite_input=False, interpolation='linear', keepdims=False)    \n",
    "\n",
    "    for m in range(total):\n",
    "        if PCs_data[m] > val :\n",
    "            PCs_data_copy[m]=1\n",
    "        else :\n",
    "            PCs_data_copy[m]=0\n",
    "\n",
    "        #print(ens_1o)\n",
    "\n",
    "    PCs_binary=PCs_data_copy.reshape(40, time)\n",
    "\n",
    "    #probabilities of warm winter for each year\n",
    "    binary_matrix=PCs_binary\n",
    "    #sum columns to give number of ensembles that say warm event for that year\n",
    "    column_sums = binary_matrix.sum(axis=0)\n",
    "\n",
    "    probabilities= column_sums/40\n",
    "\n",
    "    prob= np.zeros(time*3).reshape((time,3))\n",
    "    prob[:,1]= probabilities\n",
    "    prob[:,0]=np.arange(time)\n",
    "    prob_copy = copy.copy(prob[:])\n",
    "\n",
    "    bins=np.arange(0,1.1,0.1)\n",
    "    #to catagorize probabilities into correct bins\n",
    "    for n in range(1,11):\n",
    "        for p in range(time):\n",
    "            if bins[n-1]<=prob[p,1]< bins[n] :\n",
    "                prob_copy[p,2]=n\n",
    "    \n",
    "    #observed data:\n",
    "    #sort observed data and set lower tercile(warm events) to 1\n",
    "    obsv= np.zeros(time*3).reshape((time,3))\n",
    "    obsv[:,1]= data_obs\n",
    "    obsv[:,0]=np.arange(time)\n",
    "\n",
    "    val_o= np.percentile(data_obs, 67 , axis=None, out=None, overwrite_input=False, interpolation='linear', keepdims=False)    \n",
    "\n",
    "    for n in range(1,time+1):\n",
    "            if obsv[n-1,1] > val_o :\n",
    "                obsv[n-1,1]=1\n",
    "            else :      obsv[n-1,1]=0\n",
    "    #obsv_re_sort is an array of years and corresponding 1's and 0's for cold winters over UK\n",
    "    obsv_bins= np.zeros(time*4).reshape((time,4))\n",
    "\n",
    "    obsv_bins[:,0]=obsv[:,0]\n",
    "    obsv_bins[:,1]=obsv[:,1]\n",
    "    obsv_bins[:,2]=prob_copy[:,2]\n",
    "    obsv_bins[:,3]=prob_copy[:,1]\n",
    "    obsv_bins_sort= obsv_bins[obsv_bins[:, 2].argsort()]\n",
    "    #obsv_bins_sort is array with corresponding bin of year from forecast data\n",
    "\n",
    "    y_co=np.zeros(10)\n",
    "    x_co=np.zeros(10)\n",
    "    C=np.zeros(40).reshape((10,4))\n",
    "\n",
    "    for n in range(1,11):\n",
    "        A=np.zeros(time*4).reshape((time,4))\n",
    "        for p in range(1,time+1):\n",
    "            if n == obsv_bins_sort[p-1,2]:\n",
    "                A[p-1,:]=obsv_bins_sort[p-1,:]\n",
    "        B=A.sum(axis=0)#so B[1] is number of observed occurances (1's) in bin\n",
    "        D=B[2]/n #no. of years in bin\n",
    "        F=B[3] # to avg x coord\n",
    "        E=[n,B[1],D,0]\n",
    "        C[n-1,:]=E\n",
    "        y_co[n-1]= E[1]/E[2]\n",
    "        x_co[n-1]=F/D\n",
    "\n",
    "    coords=np.zeros(40).reshape((10,4))\n",
    "    coords[:,0]=x_co\n",
    "    coords[:,1]=y_co\n",
    "    coords[:,2]=C[:,2]\n",
    "    coords_new=coords[~np.isnan(coords).any(axis=1)]\n",
    "    print(coords_new)\n",
    "    size=coords_new[:,2]\n",
    "    sh=coords_new[:,0].size\n",
    "    coords_new[:,3]=sh\n",
    "\n",
    "    a_file = open(\"NAO_\"+s+label+\"_coords_upper_tercile.txt\", \"w\")\n",
    "    for row in coords_new:\n",
    "        np.savetxt(a_file, row)\n",
    "    a_file.close()\n",
    "    \n",
    "    rego= np.zeros(1000)\n",
    "\n",
    "    for i in range(1000):\n",
    "        \n",
    "        random_indices = np.random.choice(time, size=time, replace=True)\n",
    "        rand_sample = obsv_bins_sort[random_indices,:]\n",
    "        #sum columns to give number of ensembles that say warm event for that year\n",
    "        y_co_unc=np.zeros(10)\n",
    "        x_co_unc=np.zeros(10)\n",
    "        C_unc=np.zeros(40).reshape((10,4))\n",
    "\n",
    "        for n in range(1,11):\n",
    "            A_unc=np.zeros(4*time).reshape((time,4))\n",
    "            for p in range(1,time+1):\n",
    "                if n == rand_sample[p-1,2]:\n",
    "                    A_unc[p-1,:]=rand_sample[p-1,:]\n",
    "\n",
    "            B_unc=A_unc.sum(axis=0)#so B[1] is number of observed occurances (1's) in bin\n",
    "            D_unc=B_unc[2]/n #no. of years in bin\n",
    "            F_unc=B_unc[3] #avg x coord\n",
    "            E_unc=[n,B_unc[1],D_unc,0]\n",
    "            C_unc[n-1,:]=E_unc\n",
    "            with warnings.catch_warnings():\n",
    "                warnings.simplefilter('ignore', RuntimeWarning)\n",
    "                y_co_unc[n-1]= E_unc[1]/E_unc[2]\n",
    "                x_co_unc[n-1]=F_unc/D_unc\n",
    "\n",
    "        coords_unc=np.zeros(30).reshape((10,3))\n",
    "        coords_unc[:,0]=x_co_unc\n",
    "        coords_unc[:,1]=y_co_unc\n",
    "        coords_unc[:,2]=C_unc[:,2]\n",
    "        coords_new_unc=coords_unc[~np.isnan(coords_unc).any(axis=1)]\n",
    "\n",
    "        size_unc=coords_new_unc[:,2]\n",
    "        sh_unc=coords_new_unc[:,0].size\n",
    "\n",
    "        p_unc=coords_new_unc[:,0].reshape(-1,1)\n",
    "        q_unc=coords_new_unc[:,1]\n",
    "        model_unc= LinearRegression()\n",
    "        model_unc.fit(p_unc,q_unc,size_unc[0:sh_unc])\n",
    "        p_new_unc = coords_new_unc[:,0]\n",
    "        q_new_unc = model_unc.predict(p_new_unc[:, np.newaxis])\n",
    "        rego[i]=model_unc.coef_\n",
    "\n",
    "    unc_reg_l= np.zeros(2)\n",
    "    unc_reg_l[0]= np.percentile(rego, 12, axis=None, out=None, overwrite_input=False, interpolation='linear', keepdims=False)    \n",
    "    unc_reg_l[1]= np.percentile(rego, 87, axis=None, out=None, overwrite_input=False, interpolation='linear', keepdims=False)    \n",
    "\n",
    "    a_file = open(\"uncert_reg_u\"+s+label+\".txt\", \"w\")\n",
    "    np.savetxt(a_file, unc_reg_l)\n",
    "    a_file.close()\n",
    "\n",
    "\n",
    "data=np.loadtxt(\"dp3_data.txt\")\n",
    "#data_mag=abs(data_mag)\n",
    "obs_dataun=np.loadtxt(\"timeseries_NAO_ERA20C_1901-2010.txt\")\n",
    "#obs_dataun=np.loadtxt(\"NAO_EOF_ERA20C_1901-2010.txt\")\n",
    "data_obsmag=obs_dataun[78:110]\n",
    "data_obs=data_obsmag\n",
    "label='dp3'\n",
    "upper_tercile_dp3(a,b,time,total,data,data_obs, label)\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "96578179",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[2.47397322]\n",
      "[1.6591642]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmEAAAGACAYAAADlBbXxAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAA9hAAAPYQGoP6dpAACDtUlEQVR4nO2dd5gb1dWH3+PejXuva2zTsTG9md7BlFBCSICE0IJDCRBSCEm+hBY6phNKIKH3EjCdmGq6jcHYYNzt7c3b935/nBGrlaVdlZFG0p73eebRaObOnTta6ezvnnvuueKcwzAMwzAMw8gsnYJugGEYhmEYRkfERJhhGIZhGEYAmAgzDMMwDMMIABNhhmEYhmEYAWAizDAMwzAMIwBMhBmGYRiGYQSAiTDDMAzDMIwAMBFmGIZhGIYRACbCDMMwDMMwAsBEmJETiMjJIuLCtloRWSsir4vIJSIyNMo1l3llB4cduzeinvDt0Mw+Vf6RzN8pgbpnRtTdJCLrRORREdnM5+cQETleRN4WkfXec6wUkZdE5Bd+3ivsntuKyPMislxEakSkRETeFZGfxHn9aBG5XkTeFJEy7zM6OUbZbiLyFxH5TkTqReR7EblcRHpGKTtJRP4V1q6lInKtiAyKUnaiiDzh3b9KROaKyPQ4299PRH4vIm9435kqEflCRC4WkR5RyncVkT+JyDIRqRORr0TknCjlthCRW7zPstr7XGa20Y7jReRT72++2vtM+8TzDIaRKCbCjFzjFGBnYD/gbOBT4GJgkYjsG2cdNV4dkdv//G5sB8aPv1MsfufVvRdwpXePeSIyKsV6w7kc+A+wCPgFcBDwB2AdcISP9wlnE2AF+nwHAz8FlgH/EpE/xHH9JOBEoB54oZ2y/wEuBO7w7nUXcD7wcHghERkCvAfsCvzRKzsHOA14RUQ6RZR9G5gMnAocC/QA3hCRKXG0fyxwLvAx8EvgcOAx4DLgORGRiPK3AJd47TkAeBK4QUR+F1FuBjALKAFebasBInIi+tl8iP7N/wycDDwRR/sNI3Gcc7bZlvUbaggdMCPKubHAcqACGBZ2/DLvmsFhx+4FqoJ+Hh8+j15Bt8Gvv1MCdc/06j4m4vip3vHf+/F5AT2BWuC+GOc7ZfgzfQ9YHke5TmH7M7zP5OQo5Xbyzp0fcfwS7/h+Ycd+4R3bJ0bZaWHHrkIF4LiwY/2AQuDhONrfG+gd5fhvvHvtFnZsC6AZuCSi7B3ABmBgjM/lGK+umVHu0xlYDbwUcfzH3jUHZfLvblvH2MwTZuQ8zrnlwAVAX+D0dN1HRM4Wkbe84alqb6jkIhHpGlbmeu9cvyjXP+wNn4WXPy5smKTKG+6aFnHdvd65rUTkZRGpxOvRi8h+IvK0N1RWKyJLROR2CRuCDavnCBH53Bu6+VZEfi3ekG1EORGRs7whmRoRKRWRx0RkYiqfX1t/JxGZISLPeENwtSLyiYgcG2fV73mv47y6QsPQ0712lwJLE3i23kB3YE2M52gOq+sbEXkpsoyI9BGRchGZ473vJCJ/EJGvvfuWeX+LX8fxfEVAY3uFQu2Kg12910hv2XPe69Fhxxq81/KIsmXea23YsSOB15xz34e1qQL1Ih0mIl3aapRzrto5Vx3l1Afe65iwY7MAAe6JKHsPKqIPDKs33s9lJ2BElDofBarQ5zMMXzERZuQLLwBNwB7xFBaRLhFb5zguKwD+DZwEHArcjQ7p3B5W5p9AL3QoJvx+m6DDWA845xq8Y79Dhz6+9MqfhAqUt0Vk84h7dwOeAV7z6vlTWJveBc4E9gf+AuwI/C9C7B2I/jMsBo4DLgJOAH4W5TlvB64HXkH/2Z2Feh7eEZFhbXw+8bDR30lE9gLmocNxZ3jP9ynwsMSIaYpgkvdaGHH8CWAJ8COvXojj2ZxzRd51Z4nI+SIyNcpQGM45B9wE7Ccim0ac/inqBZrjvb8I9cz+BzgE/Rvc7T1zKzzB1kVEhojIWehQ25XtfQgJ0M17rYs4Hnq/ddixp1Dv5TWisVV9RGQP4LfAs865RV6be6Lfxc+j3O9zVBglK+L39l4Xhh3bEih0zq2Ncq/Q+UQJXdPqGbzf61dJ1mkYbRO0K8422+LZaGOYK6zMWuDLsPeXEX040kXZ/pdgezoBXVDh1AgMCDv3ETAvovyZ3n229N6PQb0MN0aU64N6YB6O0uZT2mmTeG0a65U/POzcB+g/024R9yrC0xPesVhDVaPRYZ4r0/B3WoTGAXWJKPcsOjzUyXs/06v7WO85ewK7A994f4OtI/7uf46oL+5nA7YHvg/7flR47TkJkLByfb1z10fUuRD1CoU/yydxfrduC7tvHXBmEr+XtoYjj/DO/STieGhY9+uI4yOAd2j9e3kE6B5WZqR3/LdR7neCd27nJJ5ja+9v80TE8ZeBr2JcUwfcHuNcW8ORv/PODY9y7qXIz8U22/zYzBNm5BMbeStiUIP+kw3fft5u5SLTvCGzYtSb0wDcj8aSTA4reg+wS0Qw8inAh865Bd77A1AhcX+4Rw4d3nkTFRyRPB6lTUNF5DYRWYEKkQZUPABs5pXpjf5Tfso5Vx+61jlXhYqDcA5F/xE9ENGutcBnMdqVKD/8nURkEjAVeNB7H37PF1ABEBnU/bD3nBuAt9DP/xjnXKQXJvLzivvZnHMfoh62A4G/o97GfdC/9zMhz5hzrhL9e5/sfc6IyN7A5sDNYff+ANhGdJbeARJluDqMv6PfyUNQz+rNIvKbNsonyouop+9Kbzh7E89T+nf0e/3D8J2IDACeRr16J6IezLOA3dDPIXKI0REb59XZOcILHfX/kIiMR4dIV6CxaVHra+teSRLr2lTqNIyotDlGbxi5gvcPcBDwRRzFm51z8xOsfyw68+tr4NforLVaYAd0yCl8av+DwD9Qr9Al3tDi9ug/rxChYb0PY7Ux4v0Gp/E14W3qhHoERgJ/RZ+9GvXSvRfWpgGo8FkX5T6Rx4a1URbg2xjH4yLK3yn0OfzD26IRGd92MTos2wQUOedWxLguMqYroWdzOgz1krchmpLhMVTMHURLTNVNwK9QkXKHt78SFS8hLkf/Nj9Bh0abROQt4OLI76LT2Lnl3tsXPL13uYjc55yLHHJNGOdcvYgcBPwL/f7gte136AzIVWHFLwa2RYPtQ5/n2yLyFfo3OBG4DyhFRcpGaSuAgd5riff6KrBn2Pn70N/KD4jIOOB1tGOxj3OuhNYUe+0i4rre6HBrZPl4KPZeB7Hxd2RgknUaRpuYCDPyhUNQj8gbaap/FhqwfZQLCzwWkW0jCzrnSkXkaeCnoqkFTkEF23/CihV5r8fQ4rlqi2i98C2BbdAhp/vC2jQpolzoH2S0eK7hEe+LvLK7s3HMEDGOJULk3yn0OVxO7DQAX0e8/zZOER35maX0bM65YhG5HvWYbYknwpxzS0TkReBs7/Vw4E/OuaawaxuBa4FrvfjAfVHP00siMsY5t6GNW3+ACreJbBz3lhTOuSXAzqJpPQaiExf6Azeg3sUQ2wKrwgRYiFDnYUuvvhoRWQJsFeV2W6He55DIPR0dxg1RFF7YE2BvoIJ5pnNuZZQ6vwCOF5HhrnVcWOj+C6Jc0x6hjsFWaJxmqD1dUG/tf6JdZBipYCLMyHk8L9U/0Blct7dTPFlC/9B/+EftDUmdFqP8PWjs0sGo9+NJ51xZ2PmX0F5+gXNuo2HGZNvk0WrmoXOuWkTmA7NE5DehIUnRBJSRCWqfQ4OuRznnHkmyXVGJ9ndyzn0tIt8A2zjnIvM7+U1cz+ZNaOjnnCuOcjqUFHZ1xPEbUK/SfaiH7s5Y9Xvfg8c8AXQ9MJ6wf/pR2Av1jKbkhYzRllV4ni8R+T/UI3Z3WJHVwD4iMsorG2Jn7zVcID0JnOuJyhVenX2Bo4BnPCGKcy5SVP+A9x15AxXqM8M7PBE8DfwfOrEkfNLCyajg+2+se7TB+6j39GRa50s7Bo2ftFxhhu+YCDNyjS29nmkXYCjq1TgF/cd3ZIzhGj9iOeaiOZD+IyJXoUkoz0SH+qLxMvoP6hbU29Rq2rtzbpmIXAr8TTQ9wn9Rj9UwdIiz2jn3J9rmK9SDcYUnCEuAw9DkpZFcCjyPel5uQP/JXYhOvQ8NF+GcmycidwD3iMgM1CtSjcZm7QZ84Zy7tZ12QWJ/p9OBF0VTPdyLioKBqOCZ7pz7URz3a5cEnq0/sExEHkVnUa5A/wnPRIeiFxHxD9k5N1dEvkQF0wPOufXh50XkWdQ7Mx/1Zo1DE5N+j04swGtbBer5WocOw/4InUl5dTxDkSJyjLcbmok4Q0SqvDY+FlbuIjQWbjn6nTsW9faeFCG25qBDjnNF5Arvs9iSlsS1D4aV/Qc6ceF577tdh4reHuhkifbaPhQdghyBxmgOldYrLKwMecWccwtF5G7gzyLShHrm9keTvP4hfPhSRHqhnSHQyRkAe4qmcal2zr3o1dnkfS7/EpHbUc/Xpmj+s7nOuWSEnWG0TdAzA2yzLZ6Nlll34bPG1qG95kuAIVGuucor2yfs2L0kmawV9Rp9iva0V3r1H0js2VZ/884tJ0aCT3Sm2muod6gWjTV7lLDkmG21GRUqL6P/vEvQWWtjvPteFlF2Fjr9vg79538x6sEpiVLvKWhcWRUaAL8E9fJs5/ffybtua9T7sA4Vu2vQ2KHTw8rMJEqy1ih1XUbErNhEng2NKboAHW783vu71KDeqisJSwQaUe+fvPvuGOXc+WgajsKwz/8uWic2PQUVhoXoxINS73P7SVvPG3EfF2uLKHep99y13n1eBHaPUec0VHSu8MovRT19Y6KULUA9YuWowH0FFdLxtH1mW+2P8n3u6v2tv/c+06+Bc6LUO76NOpdFKX8COlGjzvse3kCYDbHNNj83cc4mfBj5iReXtbNzLun1CvMZb9jtUzTmZ/+Am5PzeEO+zjm3fdBtMQwjN7DhSCPvEJHdgF3QIYib2yneYfCGb+aivfvhaLD3ZugQm5EEXqqJLVEv6XZYVnXDMBLARJiRj4SGdO5Ap90bSl80bmcIOtz1MXCwc+6VQFuV20xH45iK0eSwTwXbHMMwcgkbjjQMwzAMwwgAy5hvGIZhGIYRACbCDMMwDMMwAsBEmGEYhmEYRgCYCDMMwzAMwwgAE2GGYRiGYRgBYCLMMAzDMAwjAEyEGYZhGIZhBICJMMMwDMMwjAAwEWYYhmEYhhEAJsIMwzAMwzACwESYYRiGYRhGAJgIMwzDMAzDCAATYYZhGIZhGAFgIswwDMMwDCMATIQZhmEYhmEEgIkwwzAMwzCMADARZhiGYRiGEQAmwgzDMAzDMALARJhhGIZhGEYAmAgzDMMwDMMIABNhhmEYhmEYAWAizDAMwzAMIwBMhBmGYRiGYQSAiTDDMAzDMIwAMBFmGIZhGIYRACbCDMMwDMMwAsBEmGEYhmEYRgCYCDPyDhHZXEQuE5HxAbbhDRF5I+KYE5HLgmmRYeQGInKy91sZH3RbUkVEDg76Nx9pd0RkpndsZmCNMn6gS9ANMIw0sDnwJ+ANYFlAbTgroPsahpE9HAycDVwWYBt2BlYGeH+jDUyEGTmJiPRyzm3I1ns6575Md3sMwwieHLBF76W7PUby2HBkB0ZE7hWRZVGOXyYiLuKYE5GbReR0EVksInUi8qWIHB9RLjSUsJ+I3CMiJSJSLSLPisjEKPfaV0ReFZEKEdkgIvNEZJ9o7RGR6SLymIiUAktjPNPJwKPe29e965x33Jd7ikgnETlHRD4VkRoRKROR90Tk8LDrNxqOjNHe4SJyu4isFJF6EflORP4kItZBMowwRORUEflMRGo9u/KkiGwWdv4Q7ze7fdixo71jz0fU9bmIPB72XkTkrLDfdKn3u58Ycd0bIrJARPYQkXdEZAPwzxjtvRf1goXsZ2gb79c9RWQTEblGRL71bPJ6EXlBRKaGXR9XGISIzBCRZ7zPtlZEPhGRY9u7zkgNE2FGIhwOzAYuBY4Bvgf+IyLHRCl7N9AM/Bg4F9gBeENENgkVEJGfAC8DFcDPgGOBEuClSFHk8QSwBPgRcEaMNj4P/M7bPxt1xe/sHffrnvcCNwAfAscBxwPPAONjtCkqIjIc+AA4APgLcBD6uV0C3JlIXYaRz4jIJehvYyFwFPBrYGvgXRHZ1Cv2JtAA7Bt26b5ADbCniHT16hoKbAm8ElbuduB679gsNJxgC+AdERkW0ZwRwAPAv9HhxltiNPuvwGPe/s5h2xo/7ikifYH/AacD9wCHoTZqsVc+bkRkL2AesIlXxxHAp8DD4R1YIw0452zroBsqJpZFOX6ZfjVaHXPABmBY2LHOwCLgm7BjJ3tln4i4fhfv+O+9972AYuCZiHKd0B//+5HtAf4c53Md45WfGXE85XsCu3vH/6+dNrwBvBHlM7ws7P1tQCUwNqLcBV7ZzYP+jthmW6a3MBsy3nu/iWd7no8oNwaoBR4MO/Y28GrY+2+Aq4AmYA/v2I+9+jf13u/kvT8/ov7R3n2vDDv2hld27zif5eZIW+rXPYE/esf3bacNkXZnZqR9RO34x0CXiGufBVYDnYL+XuTrZp4wIxFedc6tC71xzjUBDwOTRGR0RNkHw984595BPWd7eYd2AQYC94lIl9CGCqL/AtuLSO+IOh8nNfy450He65wU2wJwKPA6sDqiPS965/f04R6GkevsDPREO40/4JxbAbwGhHuwXwV2FZGeIjIOmAQ8hHay9vPK7Assd859470/FBUlD0T8DtcCn6GiJZxS59xrKT6TH/c8CFjsnHuFFBCRScBUPJsd0Z4XUK/alFTuYcTG4k6MRFjbxrFBtJ6BE6vsIG8/5G5/LEq5EAOB6rD3a2IVjBM/7jkE7VVHe75k2nMYOoQSjcE+3MMwcp2QzYj2+19Ni7gCHdr7E7AbMA4oAj7xju+Leo/2ofVQ5DBAgHVE59uI96naIb/uOQRY7lNbAP7hbdEwW5QmTIR1bGqB7lGOx/rBDW/jWHGcZZd4+0Xe6zlArNk7kQbKRS0VP37csxAdhh1O6sa4CPgc+H2M86tTrN8w8oGQbYkW5zSSlt81wPtAFSq4xqPeeycirwIXeEH7Y2ktworQ3/nuQF2Ue0QeS9UO+XXPQnT40o+2AFyOxsBG42sf7mNEwURYx2YZMFREhoWGGUWkGxooHo19Isp2RgPTlzrnIvPQnEjYUJ6I7IL2TO/yDs0DytC4p5v9eZwfCBmwnhHH/bjni2jg/JnoBIVUeA4Nsl3qnCtNsS7DyFfeRYPrf0LLzGe8EIi9CfNsO+caROQt1Ds2Bvitd+ptoBH4P1TQvBpW/3NeuVHOuUd8bnud19aezrkan+/5IvAXEdk7leFR59zXIvINsI1z7nftXmD4iomwjs3D6Ky8h0TkaqAHOvuxc4zyRcBrIvJXdMjuLDSW4PgoZWeIyF2o0RwD/A1YhTeTyDlXJSLnoPFZA1FDuh51sW8DDHHOnZnkcy3wXn8pIpWox+8751xxqvd0zr0tIv8C/uDNYHoONbTTgA3OuZsSaOel6D+Ld0TkRrS32QPtwR8MnBFF3BpGh8I5V+bZnL+LyP3Af9Ahyj+hv+0/R1zyKnCNt/+KV0eNiLwD7A987pxbH1b/PBG5A7hHRGYAb6H2bQQ6rPmFc+7WJJv/hfd6sYi8iIYyfO7TPa9HO8FPi8gV6Ezrnmgs6XPOudcTaOfpwIsi8hIae7cKDc3YDJjunPtRAnUZiRD0zADbgt3Q4M5P0Bk5S2nJ7uwiyjl0ps+Z6JBiPTqj5scR5U72yu4H3A+UenU/D0yKcv89UCFT7NW50nt/TFiZy7w6ByfwXL9G4yoavWtP9uueaCD/uaiBrUO9a+8Ah4aVeYN2Zkd6xwaj6S6+9dpSDMxHe+y9g/5+2GZbpjciZkeGHf85GrQe+s09RZQZxGjqCocGrYcf/513/JoY9z0FDVOo8mzWEuA+YLuwMm8ACxJ4lm5oupn1aMqeVs+V6j3RmaPXo5Oe6tFwiueAKWFl2p0dGfa5PezVUY+GW7wKnB70dyKfN/E+fMNoE9HkrXOcc79qp9zJaM6a7Z1z8zPRNsMwDMPIRSxFhWEYhmEYRgCYCDMMwzAMwwgAG440DMMwDMMIAPOEGYZhGIZhBICJMMMwDMMwjAAwEWYYhmEYhhEAHS5Zq4gIutRFZdBtMQyjXfoCq10eB6+aTTKMnMF3e9ThRBhq7CwLuWHkDqPRDN75itkkw8gdfLVHHVGEVQKsWLGCfv36Bd2WnOO4447j4YcfDroZ/jB/PhQXw9ChiV23bh0cd5zuP/AAjPZjDd12KC6G7t1h++2hR4/03y8LqKioYMyYMZD/HiKzSUmSV/bID9asgY8+guHDoXOs1efSxNdfw+mn6/7118O222b2/mmmoqSEMUcdBT7bo44owgDo16+fGbwk6Nq1a358btXV0NAAw4ZB796JXXv//fq6ww4wdar/bYuksRE6dYJttklcMBo5g9mkxMkbe+QHzsFXX0H//rplkuZmuOEG3T/oINhtt8zePxPU1aWlWgvMNzomZWWwYQP06pXYdY2N8PTTuq+9ovSzfr32bEeMyMz9DMPIPUpLobAQBgzI/L2fegq+/FI7tL/+debvn8OYCDM6JoWF0KULiCR23Ztv6tDgwIGw557paVs4NTXaxoKCzA8vGIaRO6xdq9797t0ze9+yMpgzR/fPOAMGD87s/XMcE2FGx6O2VkVY376JX/vEE/p6xBHQtau/7YpGYSGMHWuGzTCM2NTUwKpVmR+GBBVg5eUwaRL86EeZv3+OYyLM6HiUlWlMWKKxYCtXwvvvq2dq1qx0tKw15eXaxgkTEvfYGYbRcSgshMrK5DqWqbBggQ5FAvz2tzq6YCREoCJMRPYQkWdFZLWIOBGZFcc1e4rIRyJSKyLfisgZGWiqkYs4p16v5ubWx4uKVNR0SvDrH/KC7bwzjBrlTxtj0dysYnHixMwb1g6M2SQj52hqguXLNb41k521pia48kq1s4cemnezITNF0J6w3sBnwK/iKSwiE4AXgLeBacDfgRtF5Oi0tdDITd59F847D044Ac46C15+WY1FQ4OmmEhU2NTXw7PP6n4mAvKLi2HQINAUDUbmMJtkZJamJli2TD3tyeQALS6GkpLMB+Q/+SQsWgR9+sA552T23nlEoL5D59yLwIsAEp+CPwNY7pw713u/SERmAL8BHk9HG40c5JNP4G9/g4oKFTKLF8PVV7ekeaiqSnym4euv6+yjIUPSP/26oUGnQ2+1VeaDbLOIIJLkm00yMkp1NVx1lYY5dO4M++8PZ5+d2LDe6tX6msmhwNLSlmD8M89UO2skRdCesETZGXg54thLwAwRiRolLSLdRaRfaEOXHTDymeef16G8qVOZL8KJK1Ywv6REYxeKi7W3mehMw9BQ5KxZ6Td2hYUwcqSmpeigOOe48MILg25GPJhNMpLn2Wfh1VeZ39TEiZ9/zvz774e3347/+spK9exvsknamhiVm27Se0+eDEd3DKfvomXL0lJvromw4cC6iGPrUI9erOljlwDlYZstD5LvrFqlAe0iXPfFFzyydCnXrVunPcZVq9R9ngjLlmkW6k6ddFZkOtmwQe8zcWKHTUnhnOOcc87hzjvvDLop8WA2yUie9esBuG7ZMh5ZtYrrVqzQmNVErk9mklEqfP45PPOM7neQYPyFS5dyyEUXpaXuXPz0IscoJMbxEJcD14a974sZvfxm6lT47DNobua8rbYC5zivd28YN05FzpAhidUX8oLtumv6vVNFRTrVu4O690MCbE5oqCM3MJtkJEdBAXTpwnlDh0JVFeeNG6ezoeOhoQFWrMjsxJ3GRrjiCt0//HDYeuvM3TsgFi5dyl5nnEFxeXla6s81T9hatOcZzlCgESiOdoFzrs45VxHayP916IzDD1fBtXAhMyoreXDYMGZMnAh77aVBsInk96qt1eFNSL/bvaxMvXTxGuE8I1yAiUiuCDGzSUbyHHAA/PSnzBg5kgf32osZl1wC220X37WFhWozMpkb7PHHNca2b98OEYwfEmCFpaVsU1CQlnvkmifsXeCwiGP7A/Odcw0BtMfIRiZM0N7a88+rwRgzBg48UINJq6sTq+vVVzVf1/DhmpoiXTQ16USCbbZJfLg0D4gUYHfffTdHH300Z599dtBNaw+zSUbydOkCp54KP/2phiHEmzbHOQ2t6No1c2ELxcVw6626f/bZwSyPlEHCBdj0qVN5/P/+jwnHHOP7fQIVYSLSB5gUdmiCiGwLlDjnlovI5cAo59xPvfO3Ab8SkWuBO9Gg2J8DJ2Sw2UYuMG6cpqYIUV6usV2JGo7QUOSRR6bX2BUXa1b80aPTd48sJZoAO+WUU6ioqMh4W8wmGYGQaFxVWZnGg2UyIP/GG3Vm+WabqT3MYyIF2Nw5c+jS2JiWewU9HDkD+MTbQOMkPgH+4r0fAYwNFXbOfQccDMwEPgX+CMx2ztlUcKNtyso011ciKR+WLNHYss6d0xuQX1+v26RJ0K1b+u6ThcQSYAFiNsnIftat05iwHj0yc79PPtGRBRG4+OK8njQUTYANTOOQb9B5wt6gJYg12vmToxx7E5ievlYZeYdzsGZN4gYr5AXbc8/0rt24fr16wIYNS989spAsFGBmk4zsp6ZGE7tmKhassVEz44Om6Nlyy8zcNwAyLcAgeE+YYaSf6mr1hCUyi6imJjMB+dXVGtcxYULiyyjlMNkowAwjJ8j0OpGPPKKjAv37ayxYnhKEAAMTYUZHoKxMRVXPnvFf8/LLKpBGj4btt09Pu5zTlBTjx8PAgem5RxZiAswwkqSpSb1gPXtmZp3IoiK4/Xbd/9WvMp8UNkMEJcDARJjREVi/PvFYq8e9kJ6jjkqfhyo0vXz8+PTUn4WYADOMFCgpUWGUKTF0ww3aGd1ii/Qnqg6IIAUYmAgz8p2aGjVaiaR9+Oor+PJLHSY8LDL7gE80NemQQkEB9OqVnntkGSbADCNFVq9WD3oiuQ6TZf58ePHFlmD8PAyXCBdg06ZMybgAAxNhRr5TVqZZ8hNZ1iPkBdtrr/TlwikshKFDYdSo9NSfZZgAM4wUqaqCtWszk5+rsVEXFgeNid188/TfM8NEesBeueWWjAswMBFm5DtFRdqDizd+oroaXnpJ99MVkF9Xp56wgoLM9GgDxgSYYfjAunWZWyfyoYfg22912PPMM9N/vwwT9BBkOCbCjPylrk7jwRKZRfTii+o5Gz8epqcp60BhoWbx7wApKUyAGYYPNDRoQH4mZkSuXw933KH7s2dndlmkDJBNAgxMhBn5THm5uvDj7Tk615Ib7Kij0jP7qKpKE8ZOmJCZ2U0BYgLMMHyiqEiXXcuEWLj+eu2Ibr01HHpo+u+XQbJNgIGJMCOfKfbWT443u/PChbrWZLducMgh/rfHOW3T+PF5O9U7hAkww/AJ59QL1qVL+jPVf/ihpufp1CnvgvGzUYCBiTAjX2ls1CDWRGZFhgLy99svPT3OkhIVX+PG+V93FmECzDB8pLxchwjTHZDf0NCSGf9HP4IpU9J7vwySrQIMTIQZ+Up5uaaAiFeEVVZqDxB0KNJvGhs1qHbSpMSSxuYYJsAMw2fWrs3MOpH//jcsW6aJo884I733yiDZLMDARJiRr5SWQnOzuvDj4fnnNZC/oEBjIfymqAiGD4eRI/2vO0swAWYYPlNbq0OR/fql9z5r18Jdd+n+7NmZWxIpzWS7AAMTYUY+0tysC3bHmwQ1PCD/6KP9D5ivrdU2FRTELwpzDBNghpEGCguhoiL9Iuy66zSx9bbbpiceNgByQYCBiTAjH6mo0C3e3txnn2lOnB494OCD/W9PKCXFkCH+150FmAAzjDTQ3AwrVqR/ncj33oNXX9Wg/4svzotZ27kiwMBEmJGPlJZCfX3860WGAvIPOCCxQP54qKxUI5qnKSlMgBlGmigu1jCGdAbk19e3ZMY/9ljYdNP03StD5JIAAxNhRr7hnMY3xBv8XlamvUDwPyC/uVlnRE6YkHcJD8EEmGGklTVr0r9O5AMPwPLlMGgQnH56+u6TIXJNgIGJMCPfqKpSYRXvUORzz2lvcOpU/9dHKynRmUZjx/pbbxbgnGP27NkmwAwjHVRVqQhLZz7BNWvg7rt1/9xz/R8FyDC5KMDARJiRb5SWaiB8PNO505khv7FRs04XFKR/anmGCXnAbr75ZhNghpEO1q/XlDbpFEbXXKMzwqdPhwMPTN99MkCuCjAwEWbkG+vWxR8LNn++uuJ799Z4MD9Zvx5GjMi7lBQ2BGkYaaahQQPy0ynA5s2DN97Ii2D8XBZgYCLMyCeqq3UIMN6hyJAX7MAD419fMh5qatSoFRSkf5mRDGICzDAyQFGRhlSkayiyrg6uvlr3TzhB7VSOkusCDEyEGflEWZkOAcaTH6y4GF57TfePPtrfdhQVaRzY4MH+1hsgJsAMIwM4B6tXa+ctXR24++/XBLBDhsBpp6XnHhkgHwQYmAgz8onCQk2GGo9r/dlnoakJttwSJk/2rw3l5SoC8yglhQkww8gQ5eU6uztdaSlWroR779X9887zdwQgg+SLAAMTYUa+UFurHqh4hiKbm+HJJ3XfTy9Yc7N64yZOzJtlP0yAGUYGWbdOZ2unazJPKBh/++1hv/3Sc480k08CDEyEGflCWZlO646nZ/f++7BqlQa++mmIios1386YMf7VGSAmwAwjg6R7nci33oK339bRgosuyklPfb4JMDARZuQLxcVqVDrF8ZUOZcg/9FD/epwNDdrDnDQJunf3p84AMQFmGBkmnetE1tbCP/6h+yeeqOESOUY+CjAwEWbkAw0NGkcRz5TuwkLtDYK/GfILCzUdxfDh/tUZECbADCPDNDerF6xHj/g6koly330a8D9sGPz85/7Xn2byVYCBiTAjHwgNRcYjwp5+WgPyp03T2C0/2LBBDefEiTmfksIEmGEEQEmJxrSmIy3FihUqwkCD8eOZPZ5F5LMAAxNhRj5QUqJTu7t0abtcU1NLQL6fXrBQSopBg/yrMwBMgBlGQKxZo96weBNNx4tzmhOsvh523BH22cff+tNMvgswMBFm5DpNTToUGU9A/jvv6Oyj/v1h7739uX9ZmXrgcjDGIhwTYIYRENXVOlSYDi/Ym2+q3cvBYPyOIMDARJiR65SXazBrPCkhQgH5hx3mT/B8U5Peu6Agpxe/NQFmGAGSrnUiw4PxTzoJxo3zt/400lEEGJgIM3Kd0lJdLLtr17bLrVmj66UBHHmkP/cuLtas+KNH+1NfAJgAM4wAaWzU9WvT0Yn75z91lGD4cDj1VP/rTxMdSYCBiTAjl2luVnEVT6DpU09pfMT22/vTI6yv162gwP84jgxhAswwAia0TqTfIuP77+Ff/9L93/wGevb0t/400dEEGJgIM3KZykodjmyvF9nYqLMiwb+A/PXr1QOWoykpTIAZRsA4p0mjO3Vqf1JRovVefbWm7tllF9hzT//qTiPhAmzalCkdQoCBiTAjlykr0wSp7SVcfest7XEOHAgzZ6Z+3+pqHf6cMCE9OX3SjAkww8gCKip0otDAgf7W+9pr8N57aqMuvDAngvEjPWCv3HJLhxBgYCLMyFWc06HIeDLeP/GEvh5+ePuxY/Hct6gIxo/333hmABNghpElrF3r/zqRNTVw7bW6/7Of5cQSah1xCDIcE2FGblJdrZ6w9mZFrlypvUIRmDUr9fuG4jdyaKZRCBNghpEl1NWpbYpnVnci3HWXetdGjoSTT/a37jTQ0QUYZIEIE5GzROQ7EakVkY9EZPd2yp8oIp+JyAYRWSMi94hIbmfJNBKnrEx7fe0F5YeSs+60U+qzGJuaNA6toCC+vGRZhAmw+DGbZKSddKwT+d138MADuv+b3/jrYUsDJsCUQEWYiBwHXA/8DZgGvA28KCJjY5TfDbgfuBvYAvgRsD1wVybaa2QR69e3P7TY0ADPPKP7Rx+d+j0LC2HoUBg1KvW6MogJsPgxm2SkneZmXUqoe3f/YkpDwfhNTbD77rDHHv7UmyZMgLUQtCfsfOBu59xdzrlFzrlzgRXAmTHK7wQsc87d6Jz7zjn3P+B2YEZmmmtkBTU1GpfVniv/9dc1j9iQIbDbbqnds65ODVxBQepxZRnEBFjCmE0y0ktpqdqvAQP8q3PuXPjgAxV2v/mNf/WmARNgrQlMhIlIN2A74OWIUy8Du8S47B1gtIgcLMow4Bjg+Tbu011E+oU2wOdBeCPjlJXpotntDQmGMuTPmpX6FPDCQg1yHTo0tXoyiAmwxDCbZGSE1av9XSeyuhquu073Tz45qz31JsA2JkhP2GCgM7Au4vg6IGryJefcO8CJwMNAPbAWKAPOaeM+lwDlYdvKVBptZAFFRRpo39bU62XL4KOP1N1/xBGp3a+qSnuYOZSSwjnH7NmzTYAlhtkkI71UV+usbj/XibzrLu0kjh4NP/2pf/X6jAmw6GTDfxQX8V6iHNMTIpsDNwJ/QXusBwITgNvaqP9yoH/YlrtrzBg6pXv9+vYDWkNpKXbdNbWEqs7p8kTjx6dngd00EPKA3XzzzSbAksNskpEeCgu1U+fXxJ6lS+Hf/9b9Cy/0Z03cNGACLDY+pulNmCKgiY17mEPZuCca4hJgnnPuau/95yJSDbwtIn9wzq2JvMA5VwfUhd5LDiSuM9qgrEyN2IgRscvU1sLz3mhQqgH5paU5lZLChiBTwmySkT4aG3U5od69/Umg6hxceaXGqs6cqR3OLMQEWNsk7AkTkWUicmms2ULx4pyrBz4C9os4tR8aZxGNXkBzxLGmUNNSaY+RI5SU6GvnzrHLvPqqLmc0fDjsvHPy92psVMG36aY5sfZaRxRgftkjMJtkpJnQOpF+edRfegk+/li9Xxdc4E+dPmMCrH2SGY68BjgC+FZE5orI8SKSrA/0WuAXInKqiGwmItcBY/Fc+SJyuYjcH1b+WeAoETlTRCaKyK7oUMAHzrnVSbbByBUaGzXLdHtrRYaGImfNalustUdRkQq5kSOTryNDdEQB5uGnPQKzSUY6cE4D8v1aJ7KqqiUY/+c/b3tkICBMgMVHwiLMOXeTc247NP7hS9TgrBGRm0VkeoJ1PQycC1wKfArsARzsnPveKzICNYCh8veiU8h/BSwAHgW+BnxaldnIasrLNcFhWyJsyRL47DMVX6kE5NfWqpt/4kR/F9dNAx1YgPlqj7z6zCYZ/lNRoR1Iv9JS3H67xqqOHQs/+Yk/dfqICbD4STow3zn3mXPu18Ao4M/AL4APvczRp0qcgQ7OuVucc+Odc92dc9s5594KO3eyc25mRPmbnHNbOOd6OedGOud+4pxblexzGDlEaakKo7ZEUcgLtscemh8sWQoL1cBleUqKjizAwvHLHnl1mU0y/GXdOs016EdYwzffwCOP6P5FF/mX6sIn8lKA1dXp/4Q0kHQXX0S6AkcCp6AxE++hWaNHotmm9wV+7EMbDUPz6qxZ0/asopoaeOEF3U8lIL+yUo3lhAn+BNCmCRNgLZg9MrKWujrNkO/HEkXhwfj77KPLsWUReSfAQrPja2vTNjkrYRHmufhPAU5AA1D/BZznnPsqrMzLwFvRazCMJKio0G3gwNhlXn5ZYyVGjYIddkjuPs3NGvy/5ZY6KzJLMQGmmD0ysp7QOpGprl0L2sn89FNdF/K881Kvz0fyToDV1mo6pP799f9BmtYLTsYT9iEwF13G4ynnXEOUMl8CD6XSMMNoRVmZ5ghry/UeypB/5JHJJ1UtKVGhNzblyXZpwwRYK8weGdlLczOsXOnPOpGVlXDDDbr/i1+klv/QZ/JKgIV7vyZM0NnxvXurkE4DyYiwiWFBqlFxzlWjvVPDSB3ndCiyrXiKr76CL7/UeLHDD0/uPo2NuhzS5ptrTzMLMQG2EWaPjOwltE7koEGp13XbbdpJHD8eTjwx9fp8Iq8EWMj7tckm6v0aMSLtq6QkI8KGishw59z74QdFZEegyTk335+mGYZHVZV6wtqKqQgF5O+1V9tDlm2xfr3+6LI0JYUJsKiYPTKylzVrNH4r1eD5r7+GRx/V/Ysugq5dU2+bD+SNAHNOxXJ9PRQUwKRJ0KtXRm6djMSbA4yJcnyUd84w/KWsTHsosbxT1dXw3//qfrIB+bW1GoRfUJBabrE0YQIsJmaPjOxkwwbNDZaqKGlu1mD85mbYb7/k4119Jm8EWG0tLF+u/19mzICttsqYAIPkPGGbAx9HOf6Jd84w/GXdurZ7kv/9rxq8ceNgu+2Su0dhoeYEGzw4uevTiAmwNjF7ZGQn69erFz/V+NLnnoPPP1dhkCXB+HkhwJxTu9/YqHFfBQUZFV8hkhFhdcAw4NuI4yOAxpRbZBjhbNigQZJ9+0Y/71xLQP5RRyWXUqK8XH98WZiSwgRYu5g9MrKPxkZNS5HqOpEVFXDjjbr/y19mRd7CvBBgNTUqkgcNgsmTdZJDQLY/meHIucDlIvLDpy4imwB/984Zhn+UlakQi9VDWbgQFi9WT9mhhyZef3Oz3mPixNhCLyBMgMWF2SMj+ygu1iD6VNeJvPXWFvt0/PF+tCwlcl6ANTer+CotVe/XDjtoHHCAne9kPGEXoDl3vheRT7xj2wLrgJN8apdhKIWFOuMx1o8k5AXbd9/kYi+Ki7U3NCZaWFFwmACLG7NHRvaxerXarFSWPFu0CB57TPcvvjjw5dPCBdi0KVNyT4Bt2KD/TwYNgilTYNiwrBj5SPiv6pxbJSJbAycC2wA1wD3Af2Lk6DGM5AgtFRHLQ1VZqQlaQYciE6WhQe+x1VaaxydLMAEWP2aPjKzDj3Uim5vhiis03OKgg5KPdfWJnPaANTfr/5GmJhVfBQVZlYIoKWnt5d25w+e2GEZrysp05mOslBEvvKAiqqAAttkm8foLC9UVnUVJD02AJY7ZIyOrWLdOZ9ylEr/1zDMaatG7N/z61/61LQlyWoCFvF+DB6sAGzo0K7xf4SQlwkRkMjATGEpEXJlz7i+pN8sw0KFCiJ4sL9WA/A0btN6JE7MmJYUJsOQwe2RkDfX1qa8TWVYGN92k+6efHuiM7ZwVYKHYL+dg6lS181nk/QonmbUjTwNuBYqAtYALO+0AM3pG6jQ0qEu/T5/o5z/7DL79Vn9YhxySeP1FRZqQL0tSUpgASw6zR0ZWUVios61TWSdyzhytY9IkOPZY/9qWIDkrwKqr1b4PGaLeryFDss77FU4ynrA/AL93zl3pd2MM4wfKyzXmK9ZQYcgLtv/+sYVaLMrK9JoJE1Jqol+YAEsJs0dGduCcesG6dUt+qZsFC+Cpp3Q/wGD8nBRgTU0qggE220y9X1kU6xuLZP7CA4BH/W6IYbSipESNWjQjVFYGr76q+4lmyG9qUoG3zTaJi7c0YAIsZcweGdlBaJ3IZJdNa2rSzPjOqXd/2jR/2xcnOSnAqqo0fGXoUM37lQX51OIlGbn+KLC/3w0xjB9oatI113r3jn7+uec09mLKFF1sOxGKi9U9nQUpKUyA+YLZIyM7CK0Tmaz35cknNS1Fnz4we7a/bYuTnBNgof8VVVX6v2D77XNKgEFynrAlwF9FZCfgC6DVNHDn3I1+NMzowFRU6FBktHgt51oW6z766MTG+uvrWxZoTXVB3RQxAeYbZo+M4El1ncjSUrjlFt0/80zNZZVhck6Ahbxfw4ap92vIkKBblBTJiLBfAlXAnt4WjgPM6BmpUVqqgfldu2587qOPdLHV3r3hgAMSq3f9eg2YDTglhQkwXzF7ZARPYaGKgmQ97DfdpJ3PyZMTD7HwgZxKxNrUpGlAOneGLbbQ2N6AO9WpkEyy1uyIZjbyE+d0VmSsZYpCAfkHHhh7uDIa1dUq6iZMSD5o1gdMgPmL2SMjcJqatGPYq1dys/A+/1zzgkEgwfg55QELeb9GjNBlh7JkdnsqJP3fSES6icgUEQl2LQUjv6isbJm9GElJCbz+uu4n0lt0TgNmx49PPmjWB0yApQ+zR0ZgpLJOZFOTZsYHOPzw5JJOp0DOCLDGRo392rABttxSVxDIAwEGSYgwEeklIncDG4CFwFjv+I0i8luf22d0NEpLNQt+tMR6zzyjP8Ytt1S3fbyUlWmsxrhxvjUzUUyApQezR0bgrFql3vVkPFiPPw6LF+vSbOec43/b2iBnBFhlpX7GgwbpgtuTJ+f08GMkyXjCLkfXaJsJ1IYdfwU4zoc2GR2Zdeuizy5qbtbZQ5CYF6ypSX/EBQWJDV/6iAmwtGL2yAiOigq1Wcl4wYqLW4Lxzz47tbUmEyQnBFhjo052qKnR9X1nzAhkwkK6ScZ1Pws4zjn3noiEZ6f+EijwpVVGx6SqSj1h0Rbs/uAD7Q316QP77Rd/nYWFOmV51Cj/2pkAJsDSzizMHhlBsX69ioRkZubddJPavKlT4cgj/W9bDHJCgFVU6P+CkSPV8xVgGEm6SUaEDQHWRznem9ZLhhhGYpSV6Zh/tN5OKCD/kEPiXwOsrk49YQUF0WdaphkTYBnB7JERDPX1GpAfrdPYHp9+qvkOReC3v83Y+rVZL8AaG9Wz2K0bbL21hpAEYLszSTLDkR8C4Yv1hQzdacC7KbfI6LgUFuoPLnKGUWEhvPWW7icyFFlYqFPGA0jeZwIsY5g9MoKhqEg9NomKmMZGzYwPMGuWxrhmgKwXYOXlOvw4bBjsuKOunZnnAgyS84RdAvxXRDb3rv+1iGwB7MzGeXoMIz5qatSoRetVPv20erS23VbXA4uHqiqNLQsgJYUJsIxi9sjIPKF1Irt2Tdy+PPoofPONirezz05P+yLIagHW2KhpiXr0aPF+BbRmZhAk/N/JOfcOsCvQC1iKLhmyDtjZOfeRv80zOgxlZZrLKzI/WFNT4gH5zmnQ6/jxyQXMpoAJsMxi9sgIhNJS9bQnGkxfVAS33ab7Z5+dEfuU1QKsrExjfYcP15mPBQUdSoBBcp4wnHNfAD/zuS1GR6aoSHuUkb3Kd97RGIH+/WHvveOrq7Q0kJQUJsCCweyRkXHWrlUPTqLrRN5wg3Y2N98cjjgiPW0LI2sFWLj3a9ttYezYDie+QiT81CIytq3zzrnlyTfH6JDU1+sso2gJWkMB+YcdFp/Ba2zUocjp06FnT3/b2QYmwILB7JGRcWpq1HuTqJj56CN48cWMBeNnrQArK9NYutGjNet9hkcrso1kpOcy2p51lJlpHkb+UFamwmnEiNbH16yBefN0P94p3EVF6trOYEoKE2CBsgyzR0YmWb9ecw+ObVP/tyY8GP/II9UTlkayUoA1NOioRs+e6v0aM6bDer/CSeYTmBbxvqt37Hzg9ym3yOh4lJToa2TP8KmnNL5r++3jG1qsrdUYsokTM/bjds4xe/ZsE2DBYfbIyBxNTRqQ37t3YutEPvQQfPttRoLxs1KAlZWpcB01SvN+Bd2eLCKZBbw/i3J4voisBi4Enki5VUbHIRQbEJnNvrFRZ0UCHHVUfHUVFmowfoZSUpgHLHjMHhkZJbRO5LBh8V+zfj3ccYfun3NOWgVI1gmwkPerVy+YNk2HIDOUEy1X8NNdsBjY3sf6jI5AebnGB0Qatbfe0qHFgQNh5sz266msVDf3hAmJ9VCTxARY1mP2yPCf1av1NRFP+/XXaxLqrbbSRbrTRNYJsNLSlmHbTTeFfv2Ca0sWk0xgfuQnKcAI4DLgGx/aZHQkSkvVxR9p1J7wHBiHH95+wj7ntHe6xRYZcXObAMsezB4ZGaOyUr32iQSSf/ABvPyyzvq++OK05SzMKgFWX6/er969dYKUeb/aJBlPWBkbB8IKsAI4PtUGGR2I5mYNvo/MDbZyJbz3nu7PmtV+PcXF6jHLQEoKE2BZRxlmj4xMsG5dYutENjTAVVfp/tFH6xqRaSCrBFhJiabgGDPGvF9xkows3wvYO2ybCWwOFDjnEl4mRETOEpHvRKRWRD4Skd3bKd9dRP4mIt+LSJ2ILBWRU5N4DiNoKip0i8ySH0rOutNO2otqi8ZGdfUXFMS/pmSSmADLSny1R2A2yYhCfb12DhNZJ/Lf/4ZlyzSh61lnpaVZWSPA6ut1woKIer+23dYEWJwkE5j/pl83F5HjgOuBs4B5wOnAiyKyeRv5fR4BhgE/B5YAQ/E3ts3IFGVl+uMNz//V0ADPPKP78WTIX79eU1uMHJmWJoYwAZad+GmPwGySEYOiIrVX8dqZtWvhrrt0/9e/Tm6R73bICgEWCgWpqdGRiEmT0vKs+UwyMWFxRxY6555pp8j5wN3OOe/byrkicgBwJromXOS9D0TXg5vonPPyGrAs3vYYWYRzLRmTw3n9dY0TGzIEdm/TAaEpKUTUC5bGmAMTYNmLz/YIzCYZkTinXrCuXeO3M9ddp8Jkm23g4IN9b1JWCLC6Oh2i7ddPvV+jRmV8nd58IJne2lNoDEbkFLTIY442EiWKSDdgO+CKiFMvA7vEuOxwYD5wkYicBFQDzwB/dM7VxLhPdyA81brJ9GygqkrFVmSvKRSQf8QR7c9AKizU2ZCDB6enjZgAywGewgd7BGaTjBiUlamtiTcg/7334NVX0xaMH7gAC/d+jR+vsV/RVjsx4iKZb8f+wKfAQcAmQH9v/2PgAOdcJ29rr8swGDWK6yKOrwOGx7hmIrAbsCVwJHAucAwwp437XAKUh20r22mXkQnKyvRHHL600LJlMH++Gq32AvLLyzWgf+LEtKWkMAGWE/hlj8BskhGNNWs0TCKemNP6+pZg/OOO08SkPhK4AKut1divzp1hu+3U02cCLCWS8YRdD5zhnPtf2LGXRGQDcAewWYL1RZvZFGsZkk7euROdc+UAInI+8JiInB2j53k5cG3Y+76Y0QueVatUdL36qrr6oeV1l1106aFYNDeriNt667TFH5gAyxmux197BGaTjBA1NZobLF6h88ADsHw5DBoEp5/ua1MCFWDO6Sz02lodfdh0040TbBtJkYwIK0B7b5GUA+MTqKcIaGLjHuZQNu6JhlgDrAoZO49FqJEcTZS8QM65OqAu9F4ykMjTaIc1a+DGG2HRIn3fs6cKq8+85OcDBmi8QawFu0tK1MiNGZOW5pkAyyn8skdgNsmIpLBQ84PFY2vWrIG779b9c8/11UMUqACrrdUJUP37w5Zb6kQoi/3yjWQ+yQ+B60Xkh9WWRWQ4cA3wQbyVOOfqgY+A/SJO7Qe8E+OyecBIEQn/dk8GmrGeZG7Q2AhXXAEff6xZ8idMUK9XY6Mmbe3RA77/Hh57THtfkTQ0qFGYNCm2SEsBE2A5hy/2CMwmGRE0NalXq2fP+EIerr1WO4/Tp8OBB/rWjMAEmHM6K7SoSMM+dtzRgu/TQDKf5qloz/B7EVkiIkuA5WiW6p8nWNe1wC9E5FQR2UxErgPGArcBiMjlInJ/WPl/A8XAPSKyuYjsAVwN/DNWEKyRZXz0Ebz7rgqw8CSt332nr6FA+/feaxmeDKewUHtibQ1XJokJsJzET3sEZpOMECUlug0Y0H7ZefN0ZnfnzhqM75N3MzABVlurArRbN5gxQ0M/bPgxLSSTJ2yJiGyN9g6nom73L4FXnIvmumizrodFZBBwKWo0FwAHO+e+94qMQA1gqHyViOwH3ITOSCpGc/T8IdHnMALi5Zc1ziJ8rciKCjV2IpprpkcP7X199FHrYYANG7QXNnGi7ykpTIDlJn7aI68+s0mGsmqVvrY3S7uuDq6+WvdPOEFT5vhAIAIs5P2qr9fRhkmTNl7RxPCVpBIKesbtZRF5C6hLxtiF1XULcEuMcydHOfYVGw8XGLnCokWabydsPch1X33FMKB00CAGhGZLduvWslhuiKIiNQo+p6QwAZbb+GmPvPrMJnV0Kis1B1Y8aSn+9S/12g8ZAqed5svtAxFgtbX6zAMH6qzH4cPTNvPcaCHh4UgR6SQifxSRVUAVMME7/lcRScb9b3Qk6upav29spI8ntm4PP+5c69iDsjINdJ0wwdfmmADLbcweGWlh/XpdA7G9IbhVq+Cee3T/3HN9GbLLuABrbtbnLSnRWY877qghHybAMkIyMWF/AE4GLgLqw45/AfzChzYZ+UpDgw5D1od9bVatordzrOvcmX232EKPOadlx3qjPk1Nmhds4kRfZxyZAMsLzB4Z/tLQoLmw4kl/c8012rHcfnvYf/+Ub51xAVZTo8/aq5fGfm25ZevcjUbaSUaE/RT4pXPuQXQ6d4jP0ZgMw4hOebkmL+zVS+PA4IeA/GFTpzJj4EA9VlKi06GnT9f3xcXq6vcxJYUJsLzB7JHhL6F1ItsTP2+/DW+9pfGpF12UsucoowKsuVmHHktL1SbvsIN5vwIimZiwUegitZF0ArpGOW4YSmmperd23hneeENFWWlpS0B+aDmMqio45BCNSaiv162gQOPEfMAEWF5h9sjwj3jXiaythX/8Q/dPPDHlMImMCrANG3SW+aBBMGWKjk6Y+AqMZETYQmB34PuI4z8CPkm5RUZ+0tSkgfZ9++pyHp07w7//rec22UR7ZQ0N2vs85BDdQGMVRo70LSWFCbC8w+yR4R9lZWpz2gvIv+8+jQcbNgx+kdqod8YEWCj2yzkVXwUF8S3FZKSVZETYn4F/icgotLd5lIhMQYcFDvWzcUYeUVGhM44GD9Ze5hFHwK236rmddtIkgGPH6hBkSHBVV2vZggJfEgSaAMtLzB4Z/hHqDLYlTlauVBEGcN55KaVwyJgAC3m/hgzR4cehQ837lSUkkyfsWRE5DvgdumbaX9DFcg9zzs31uX1GvlBWpsYtlJripZfUpT9uHPzlLxsbhFC+milTdMp0ipgAy0/MHhm+UVurAqtfv9hlnNMFuuvrdRbhPvskfbuMCLBw79fUqTq5ybxfWUVCIkxEugC/R7NB75meJhl5h3O6rlqox+gcPPGE7h95ZPQeWXm5Dk2OH+/D7U2A5SNmjwxfKSxUj/3YsbHLvPkmvPOOJnC98MKkvUkZEWDV1dqRHTpUvV9Dhpj3KwtJaIzHOdcIXAj4m67cyG8qK1vyfAEsXAhff62B9odGGTEKpaQoKEg5745zjtmzZ5sAy0PMHhm+0dTUkqohllCprdWUFAAnnZR0BzHtAqypCdauVbu7+eaaPsOGH7OWZAJtXgFm+twOI58pK9NcOiE3eMgLts8+0QNgQ723UaNSum3IA3bzzTebAMtfzB4ZqVNSonanrYD8f/5TPfrDh8OppyZ1m7QLsKoqHVLt31/TTkydCt27+1e/4TvJBOa/CFwuIlsCHwHV4Sedc8/40TAjj1i7tsUQVFZqPBjA0UdvXLauDhobdXmirslnGLAhyA6D2SMjdVav1jAJEfjyy5Y4qoEDYbPNVHz9619a9oILkkpomlYB1tSkbRZR79fEib6l9DHSSzIizJvSxvlRzjlsaMAIp7pac4GFsk+/8IIKrYkTdX2ySAoLYfRo9YQliQmwDoXZIyM1qqpUZH37rc56XLFCO4Kgs7KHD4clS3Ri0S67wMyZCd8irQKsqko9ecOGaeyXz2vrGuklmdmRqecKMDoOZWU6PXrQoNYB+UcfvXGMQlWV9t4mTkw6JYUJsI6F2SMjZdatg1df1ez3TU0awB7ydNXVweLFKsI6d1YvWIKxVWkTYE1N2vbOndX7NWGCeb9ykLgMmIiUiMhgb/+fIhLHolqGgXq2unZVw/XZZ7B0qQ5NHnxw63LO6fJEEybAgAFJ3coEWMfA7JHhGw0NuvzQm2+qnRo7tvVQY+fOmpQV1DtfWJhQ9WkTYKHYr8GDNfZryhQTYDlKvL3IbkAoecrPAEs0YrRPba0ardBQZMgLdsABGy+OW1qqwaRtTQ9vAxNgHQqzR4Y/FBWpCKurix4C8fXXush1r1463DdvnnYY4yAtAqyxUePXNmzQxba3286GH3OceIcj3wWeEpGPAAFuFJGaaAWdc8lNGzHyj7IyjQnbZBPdf+UVPX7UUa3LNTZqz2769KSyT5sA63CYPTJSxzn45htYtCj6Yt2VlXoeYOutNcXO0qUaPzZyZJtVp0WAVVZqZ3X4cI39GjQotfqMrCBeEfYT4DygAA127Y/1Po32KCrS2K5OneD55zXL9JQpsMUWG5cbPjyplBQmwDokZo+M1Ckvh++/VzEWkUV+fkkJPd57jy2dU9s0YoQOXYY6lm3guwBrbNTYr27dYKutdJWRFGaOG9lFXCLMObcO+C2AiHwHnOScK05nw4wcp75eDUefPq0D8o86qnVga22tBphOnKhZqBPABFjHxOyR4Qtr1+qyPl276msY73zxBbPr6qgHum29tR5satIOZRt2yncBVlGh3q+RI9X75cMSbkZ2kfDMIufcBDN4RruUl+sQY58+8NFH2uPs1QsOPLB1ucJCjQNLMCWFCTADzB4ZSVJbqwH3o0fDmDEqdEI0NHBGVRUAhePHt6zaUVqqQ4AjRkSt0lcB1tio7aur06HQ7bYzAZan2PRuIz2UlOhr587w+OO6f+CBrZchqqzUmUgTJiQ07dsEmGEYKVFY2LI+7c47q7e+tlbPffUV3erroXdvRoW8YA0NGgy/885RF8D2VYCVl2vw/bBhukh4iomrjewmmWSthtE2obXLevdWMfb663o8PCDfOT23xRbRg2JjYALMMIyUaG7WhKw9e2rnb/p0mD8fPvmkJfge1APVubN6o1auhE03hV133ag63wRYY6PazR499N7jxiUcomHkHvYXNvynvFxjGYYOhQcfVOOyxRa6jlmI4mJ1r48bF3e1JsAMw0iZ0DqRoRCI7t3h5JNVcD36qHYQQ7kKv/uuJRnqz362UYfRNwFWVqZ2c9QonbzU1hqWRl5hIszwn9LSliDW8ID8EI2N6trffPOorv1omAAzDMMX1qxRoRU+xNe3L4wfryESXbrAjBnqKRs5ErbfXm1VxJCgLwIs3Pu17bYaH2verw5FXH9tEenXfinFOVeRfHOMnKe5WY1cz57wwQcaXNqnD+y/f0uZwkINbo0R4BqJCTAjHLNHRtJUVWm8VaSnqaoKbrxR9087DX7+8zar8UWAlZXpiMHo0TrUad6vDkm8krsMzccTD7ZgbkemslLd6gMGtATkH3xwy1IgoeDXgoK4enwmwIwolGH2yEiGwkLN8xWZZf6OO3SIcswYOOmkNqsIF2DTpkxJXIA1NGj6np49Ydo0vWdn+5p2VOIVYXuF7Y8HrgDuRTNXA+yMLh9yiV8NM3KUsjLNEVZRoQviQuuhyMJCnQ0Zx1IbJsCMGJg9MhKnsRGWL1fPfDhLlsDDD+v+RRe1uQZjyh6wsjLtqI4apXm//FhH0shp4k3W+mZoX0QuBc53zv0nrMgzIvIF8EvgPn+baOQMzulQZI8e8PTTGhe2zTY6xRpUmPXsqYlZ20lJYQLMiIXZIyMpiopUBA0f3nLMObjiCrVVe++tKShikJIAa2homTE+bZoOQZr3yyC5PGE7A/OjHJ8P7JBac4ycprpajVyvXvDUU3rs6KP1tblZA/YLCjZevDsCE2BGApg9MtrHOY1P7dy5dRjECy/Ap59qx/H882NenpIAKy1VATZmjOb9GjfOBJjxA8mIsBXAGVGOn+6dMzoqZWVQU6NGbe1adbXvs4+eKynRbNNjxrRZhQkwI0HMHhntU16ucVih1BOgw4I33KD7v/hFaw9ZGEkLsPp6zUcGmots222hX9xzSowOQjJzYc8DHheRA4D3vGM7oYvpHu1Xw4wcJLTIbCgg/9BDNQdPQ4MG5G+xhb6PgQkwIwnMHhnts26dJl0NT4lz223aORw3Dk48MeplSQuwkhIdGRg7Vmc+tuP9Nzouyawd+QIwGXgGGAgMAp4GJnvnjI5ITY0mYN2wAebN02NHHqmvhYXay2wjJYUJMCMZzB4Z7RLKeB/uhfr6a03MChqMH2VZoKQEWF2der86ddL1Hrfd1gSY0SZJZYVzzq0AfudzW4xcprRUBdjLL2v814wZmvywpkYNUkFBzDgIE2BGKpg9MtqksLAlHxeofbrySn3dbz+N04ogYQEWWoatpkY9a5Mmmfgy4iKpBbxFZHcReUBE3hGRUd6xk0RkN3+bZ+QMhYVqiJ5+Wt+HAvILC9UlP2hQ1MtMgBmpYvbIiEloncgePbQzCPDcc/D55zpT+9xzN7okYQEW8n517qyxX9tsYwLMiJuERZiIHA28BNQA04FQkE9frDfaMamr0+nfCxao6Bo4EGbO1ED9Pn00L1iUlBQmwIxUMXtktElonchQNvqKipbM+L/8JQwb1qp4QgLMOa07lPtwxx114lGnpHwbRgclmW/LH4AznHOnAQ1hx99BjaDR0QglIHzBC8E5/HDtFZaXa06wyOSImAAzfMPskRGbNWvUGxZKwHrLLWqvJk6EE05oVTQhAVZbq96vLl009mvrraPaOcNoj2RiwqYAb0U5XgFsklJrjNykuFhnH73/vr6fNUt7iEOGRE1JYQLM8BGzR0Z0qqtbrxO5aFHLzO2LLmqVLyxuAeac2rvaWvV+bbqpJmA1jCRJxhO2BpgU5fhuwLeJViYiZ4nIdyJSKyIficjucV63q4g0isinid7T8JHGRhVgb7yhBmqnnWDoUM2RU1Cw0RIgJsAMn/HVHoHZpLyhsFAnC/Xpo96wK65QG3XAATpxyCNuARbyfnXrptdvvbUJMCNlkhFhtwM3iMiO6CK6I0XkROAfwC2JVCQixwHXA38DpgFvAy+KyNh2rusP3A+8mnDrDX8pK9Oe4Usv6fujj1bjN3LkRskPTYAZacA3ewRmk/KGxkb4/vsWkfT007Bwob4/77wfisUlwEKxX0VFOoy544669qPFfhk+kPBwpHPuKs/gvA70QIcC6oB/OOduTrC684G7nXN3ee/P9ZIunknbi+/eDvwbaAJmJXhPw09KS3UYsrRUF+WePl17jAUFrYyUCTAjHfhsj8BsUn4Qvk5kWRnc7H0VTj9d7RRxCrDaWvX0DxwIW22luQ7bWffWMBIhKSnvnPs9MBhdm20nYIhz7o+J1CEi3YDtgJcjTr0M7NLGdaeg2bD/HOd9uotIv9CGzpoy/KC5WQNfX3lF3x9xhBq8cePUaHmYADPSiR/2CMwm5Q2hdSI7ddK4rzlzdJLQpElw7LFAHALMOfXoFxfrdTvsoN59E2CGzyTsCRORnwGPOeeqib5wbrwMBjoD6yKOrwOiLuIlIpsCVwC7O+caJb4fxCXAn1JopxGL8nINdv30UzV4e++tWanHj/+hiAkwI534aI/AbFJ+UFHRsk7kggXw1FN6/OKLoUuX9gVYbS2sX68dyW22UW+aiS8jTSTjCfsHsF5EHhKRQ0Ukqaz7YbiI9xLlGCLSGXX3/8k5tziB+i8H+odto5NspxFJWVlLWopddtGEiJMm/RCHYQLMyAB+2yMwm5TbhNaJ7NZNM+M7BwcfDNOmtS3AmptVfIV7v2z40UgzyRisEcCBwAnAQ0CNiDwKPOCceyeBeorQ+InIHuZQNu6JgrrsZwDTRCQU69EJEBFpBPZ3zr0WeZFzrg6NEQEtnEATjZg4B8uWweuv6/t99tFZkaNGeadNgBkZwS97BGaTcp9Q9vp+/eDJJ9VT37s3zJ7dtgCrqVEBNmgQTJmiSVzt72JkgGQC8xuB54DnRKQXcCTwY+B1EVnpnCuIs556EfkI2A94MuzUfugCvJFUAFtFHDsL2Bs4BvguoQcxUqOyEp5/Xl+HDVO3/aRJ0LWrCTAjY/hlj7y6zCblOqF1Inv31sSsAGeeycLy8ugCrLlZr2lqUvE1caIuZ2QYGSIl171zboOIvAQMAMYBmyVYxbXAv0RkPvAu8EtgLHAbgIhcDoxyzv3UOdcMLAi/WETWA7XOuQUYmSV8KHLvvXV9yKFDTYAZgeGDPQKzSblLczOsXAndu6sAq6iAyZNZOH16dAG2YYMKMPN+GQGSlAgL63GeCOwLrAD+A/wokXqccw+LyCDgUnRYYQFwsHPue6/ICNQAGtnGO++oq79TJ423mDgRJ2ICzMg4ftkjMJuU05SWqqhas0bzggHfnngie519dmsB1rcvrF2rIRVTpmg6nR49Am680VFJZnbkf4DDgA3Ao8DMJGIvfsA5dwsxkio6505u59rLgMuSvbeRJNXV8J//6P7228OMGbhNNjEBZmQcv+0RmE3KWVavhoYGuOYaAEr33JOdrr++tQDr0kVjxoYMgcmTNY7VvF9GgCTjCXPAccBLXjyG0dEIzw02axZuzBhmz55tAswIArNHhnYM16yBt96CxYtp6t2bXT/9lMLychVgN93EwJoaLTt1qnq/uncPts2GQYIiTES6ou74b8zgdWAefFDjKYYNwx17LOdcfLEJMCPjmD0yfqCwUOPB7r4bgN81N7OovJxpU6Yw9+qrGVhRoV6vkPfLMLKEhESYc65BRLYkSs4co4PQ3PzDUKQ78kjO+cc/mHPrrSbAjIxj9sgAdJ3I5cvVLlVV8VnnzvyjpobpU6Yw97LLGCgCm28OEyaY98vIOpJJ1no/8HO/G2LkBhOKi+Hrr3FduvC7sjITYEbQmD3q6BQXw7x5P4RI/LKpiW0nT2buJZcwcPRoTbo6daoJMCMrSSYmrBvwCxHZD10mpDr8pHPufD8aZmQnBy5fDsDHI0ZwxUMPmQAzgsbsUUdn+XJqb76ZHsCdQGNBAXP//GcGTpum3q9u3YJuoWHEJBkRtiXwsbc/OeKcDQvkMyUl7Ll2LQC/WbHCBJiRDZg96sg0NLDm6qsZsXo1xcBD48Yx9+abGThjBgweHHTrDKNdksmYv1c6GmJkP+6f/6RXczNfA2+CCTAjcMwedWy6rl1L7//+F4Bbhwzh0YcfVg+Yeb+MHCGZmDAARGSSiBwgIj2995ZsJY9xzc2s/NvfALgDuPuf/zQBZmQNZo86Hgs/+YSjPvuMfsCCHj046623GLjjjibAjJwiYREmIoNE5FVgMfACOkUc4C4RucbPxhnZgXOOq485hjFlZdQC2914owkwIyswe9QxWfjZZ/xujz040TmagTFPPcXAqVODbpZhJEwyMWHXAQ3o0h2Lwo4/7J27wId2pZ3jjjuOrl27Bt2MrMc5x4IFC/jDsmUAPNWlCw/NnctDc+cG2zAj72loaIinWF7YIzCbFC+VJSV88P77vN+oqeFeGDOGO+bMgTlzAm6Zkc/EaY8SRpxLLHZVRNYCBzjnPhORSmAb59y3IjIB+MI51ycdDfULEekHlJeXl9OvX7+gm5PVhBbj/tecOawGegMX77ILV86bF3TTjA5ARUUF/fv3B+jvnKuIVibX7RGYTYqbxkYWzp3LXiecwE/Ky7kWKOvalU3WroWBA4NunZHnxGOPkiGZmLDe6DptkQwG6lJrjpEthATYnDlzOAn9ozNlCovM2BnZhdmjjkBpKQsfeYS9jj+eLuXl/NUL+bt/881NgBk5TTIi7C3gp2HvnYh0Ai4EXvelVUaghAswAf46aJCe+OUvbbFbI9swe5TPNDTAkiUqwM46i8KKCv7Zty+9nYPtt+fVMWOCbqFhpEQyIuxC4HQReRFNlHgVsADYA7jYx7YZAdBKgInw7CmnMKC4GHr0gFNPDbp5hhGJ2aN8paQEPvqIhc89x16//S2F5eX8cvRoDqyshE6d4LbbcNYpNHKcZPKEfSkiWwNnAk3ocMATwBzn3Bqf22dkkEgBdvc113DIM8/oycMPh002CRWEb7+F11+HTz7R3uqYMbD33rD99jZF3MgYZo/ykIYGWLZMPWDffstef/kLhWVl7DBlCrfU1GiZn/0Mpk8PtJmG4QfJzI7EObcW+JPPbTECZCMBdvvtnDJypK7JBnDWWaGCcN998OijUFYGvXtD587w1Vfw2mtqGC++GIYODexZjI6F2aMcoq5OPVzNzdCvH/Tt2/p8cTEsXgxr1rCwvJy9/vhHCktLmT51Kq/vsQed77hDO4NXXhlI8w3DbxIWYSJyIFDlnPuf9/5s4DTgS+Bs51ypv0000s1GAuzuuzllzz3hr3/VXunkybDHHlp49Wq4917o3x+22KJ1jNiGDfDuu/D3v8Pf/qYCzTDSiNmjHGHFCvWcv/SSdt6c0wW1d90V9t1XF9hevhyWLoWGBhbW1bHXBRf8IMBeuewyeoVyE15yCQwZEujjGIZfJBMTdjXQD0BEtgKuRZMkTvT2jRwiqgA79lg1hq++qoVOPlnFVkmJirDevWH48I2D9Hv1gk03hY8+grffzvizGB0Ss0fZjHPwwgtwzjlw++1qQ/r21U5cczM89RScey5cdBF8/DH06sXC2lr2OvvsHwTY3DlzGHDXXVBTA5tvDmecEfRTGYZvJDMcOQHtZQIcDTzrnPudiExHjZ+RI0QVYKecAgsWwHvvae+1Z0/4qTf5bN48qK+HESNiV9qjB3TpAv/9LxxwgM2mNNKN2aNs5uWX4YYb1A5Ees5791Z7sXo1vPIKOMfCXXbZSIANXLRIz3fqBJddpsOYhpEnJOMJqwd6efv7Ai97+yV4PVIj+4kpwIqL4fvv4c03teD++8PIkbq/apX2bDt3/qGe+YWFnPjaa8wvLGypfJNNdGihvj5zD2R0VMweZSvFxXDnnWozxo1rLcBqaqCwEEpLYcAAGD2aha+9xl6nn95agPXsCVddpdcccgjsZeu1G/lFMp6w/wHXisg8YAfgOO/4ZGClXw0z0kdMAdbcrLMei4tbRNjPf96mN+u6L77gkaVLAXhw770z0XzDCMfsUbby1luwdq3Ge4VoaoLKSqioUHHWuzeIsLC5mb0WL6awsbFFgPXvD/fcox26AQPgV7+yxKxG3pGMJ+xXQCNwDHCmc26Vd/wg4L9+NcxIDzEFGMCaNTo08MEH0Nio8V277dZy8dixKsi8NdsAzttqK44tKOC8rbZqKVdaCgUFlqrCyARmj7KVl1/W4PuQ5zzc+9W1q8aQirCwooK93n5bBVjv3sy99FIVYGvWwF136bUnnQSbbaZDkoaRRySTJ2w5cGiU4+f50iIjbbQpwOrqYMkSNY5PPaXHjjyydfzFrruqUV21SocXgBlDhrT2gG3YoB41iwczMoDZoyyloQHWr9cg/KYm9XxVVgIwv76e6xYt4rxJk+jZubMKsPp6pvfvz9xx4xjY3Kx1XHut2qWtt4bDDrMZkUZeklSeMBHpDBwJbAY44CvgKedcY5sXGoHRpgADWLkSioq097lqlfZSTzihVfwX/fvD6NFqYFeu1Fix8J5pZaUmWdx999YeNMNII2aPshTnVEStXw+1tT9M2rluwQIeWbWKsoYGPiwtVQG2ySbM3XFHBoZiS+fN05QWnTvrSh2jR+v1hpFnJJMnbEvgGWAY8LV3eDJQKCKHO+e+8LF9hg8455g9e3ZsAVZZqSkpBgyAOXP02MyZavgiGT5chwYeeAC+/FI9Y506tRjZvfeGCy4wg2lkBLNHWUpTk86SXrNGZ1P36fPDqfMmTaKsoYF3S0oobWhQAbbrrgysq1O70aMHXH21Fj7mGJg0Se2OYeQhyXjC7kLXZtsulAhRRAYA9wJ3ADv71jojZdr1gDkH330H1dWajuKtt/T4kUeq5ysaP/qRerrefBM+/1w9Y6NHw5576tCBxW0YmcPsUTbhnHq+Fi/WkIVFi7SjFkbPzp35sLSU0oYGpvXvrwKsWzcVbFtvDXPnqqd9yBA46ihdfSOWLTKMHCcZEbYNMCM8E7VzrlREfg986FvLjJRpV4CBzoRcvlwN3oMPag92s81gl100PiwWI0bA8cfrZhjBYfYoW6itVY/6t9/qMOL++2vi5nXrfsgt+EMQfmgIMiTAqqo0hnTSJF2pA2D2bJ3cM3q0xZcaeUsyLouvUdd/JEOBJak1x/CLuARYU5MaTefU2IUC8vffHwYNynibDSMJzB4FjXOaiuL99+HrrzVP4LBhmk7iiCN0os7atboWZDQBVlmp3rOdd9Ykz3V1MGMG7LCDhkgMHhz0ExpG2ojLEyYi4UkPfwfcKCKXAe95x3YCLgUu9rV1RlLEJcBA3f+hmI1331VD2q8f7LefGlLDyELMHmURNTU6q/q779RzPmZM63CEXXeF5mYWPvgge33wAYVNTUzv25e506YxsKoKysv1ut13V4/XzTerF+2ii9Q7VlDQtke+uRneeEPFX12dzsbcbjvYcksLizBygniHI8vQWUchBHgk7FjIV/wsEDadzsg0cQuw0NBBz54aQPv443p8r73UGEbEcRhGFlGG2aNgCXm/Fi/WkIahQ9WWRCLCwlGj2GvJkhYBNmUKA2tq1Mbstpt6vMaNgx//WK/58Y/Vk1ZTo6+x7v/f/8Jnn8Gf/qS5Czt1UlH273/rEkmnn66hFYaRxcQrwmytiBwgbgEGGgdWVKQJWNeu1SnhAPvua/l4jGzH7FGQhLxfy5ZpBy6UxDkKC5cuZa8zzqCwvFwz4V97LQOdU7HUu3dLHsLbb9fUOEOHwmmn6ULf48e3mlXZikcfhTvu0M7khAmtO42VlRqLdumluoUnkjaMLCMuEeacezPdDTFSIyEBVlGhwwcDB2rv8amn1Chuu61mybehSCOLMXsUEM5p+MLixSqShg1rMxXNDwIsfC3IaLMcV66E++7T/fPPbxl+DK1ZG8mCBXDvvep569VrY699377qAfv6a034esst0b10hpEFJJusdRPg57QkR/wS+Kdzrty/phnxkpAAc05nL23YoAGvjY0bB+T36hX9WsPIQsweZYANG1q8X926ten9ggQEmHOaE6y+HnbcEfbZRz30Q4bEXidy7lztSG6xhabIiYaIesi+/Rbee88W/jayloQjF0VkBrAUOA8YCAwGzgeWish0f5tntEdCAgx07bYVK9TtD/D222r0BgzQgNZYMRiGkYWYPUozzul6su+/rzGkgwapQPJDgIHmGpw3T4c1L7xQj9XUaFxq5yjhfOXlmsswntnbIQ/Za6+1X9YwAiKZ6SPXoRmqxzvnjnLOHQlMAJ4DrvexbUY7JCzAGhvVkIq0DCM88YS+HnywuvFtKNLILcwepYvqavU0ffiheqrGjGl3JYyEBFhtLVxzje6fdJLGgFVUqB2KFZdaWqoiLSJWbH5hISe+9hrzQ8sehejVS0WkYWQpyQxHzgBOC1+XzTnXKCJXAfN9a5nRJgkLMFBjtHZtS6zFypWamgJ0uaEBA2IHwhpGdmL2yG+amzX266uv1PPUTuxXiIQEGMA99+h9hg/X9SFB7zd1auwYrpB3zLlWh6/74gseWboUgAf33rvlhHPqZTOMLCUZT1gFMDbK8TFAZaKVichZIvKdiNSKyEcisnsbZY8SkbkiUigiFSLyrogckOg9c52kBFhNjXrBevduMUqhWLCddlIP2PDhlpnayDV8tUfQwW1SyPs1f756zseOTY8AW74c7r9f9y+4QEVXba3Gm7W1TuSQIRrLWlra6vB5W23FsQUFnBc5E7KqCiZPbrf9hhEUyYiwh4G7ReQ4ERkjIqNF5Hh0Dbf/JFKRiByHDhn8DZgGvA28KCLRjCrAHsBc4GBgO+B14FkRmZbEc+QkSQkwUKNXWtoS7NrQAM88o/uHHabxEzYUaeQevtkj6MA2qblZPePvvaczp0NiJ45OWcICzDm46iq1QbvsAjNn6vGyMr1vW3aoRw848EAdtmxu/uHwjCFDeHDvvZkRPoxZWanl99mn3WcwjKBIxk/7G3QG0v1h1zcAtwK/TbCu84G7nXN3ee/P9XqRZwKXRBZ2zp0bceh3InIEcBjwSYL3zjmSFmBlZWpYBw1qMapvvKHTzAcPhm220Xw9/fq1VYthZCN+2iPoiDapqgq++Qa+/15jqMaMidsjnrAAA3j9dRV7XbtqML6ILqHW0BDfOpH77APPP6+pMmJRV6fPs8ceuii4YWQpCXvCnHP1zrlfAwOAbdHe4kDn3HnOubp46xGRbmjP8eWIUy8Du8RZRyegL1DSRpnuItIvtHnlc46kBZhzKsDq6lrHe4UC8o84QoNuR4ywZT6MnMMvewQd0CY1N+tM6fffV8EydGjrjlo7hAuwaVOmxCfAampagvF/9jMVfKCxYJtsEl+i6BEjdFmjYcPU27VuXYtXrL5evf7ffKPrT15wgdk1I6tJOmLRObcB+CKFew9GlxRZF3F8HdBGUEArLgB6o0uWxOIS4E8Jty6LSFqAgS6Mu2JFa+P2/fc646lTJx2KFNGgfMPIUXywR9CRbFJlpeb9Wr5cvV/xeKDCSMoDBnD33SqaRo6Ek09u3Z5ttml7nchwpk2Dyy+HhQu1o/nVV3q8c2cYNQqOPx5mzYJ42mQYAZIN00ZcxHuJcmwjROQE4DLgCOfc+jaKXg5cG/a+L7AywTYGRkoCrKFBg/E7d26dVfrJJ/V1l11aAvXNWBlGiPy1SU1NujzQ4sUqfIYP12D4BEhagC1bBg88oPsXXNAS8F9drUIwlLswXgoKdLvzTn2ekLd/q60s4bSRMwQpwoqAJjbuYQ5l455oK7zg2buBHznnXmmrrDck8cOwhOTQ7L+UBBi0pKQYNarlWF0dPPus7h99tMaDTJ1q07gNI99tUkWFDtOtWKGdr7Gx5hrEJmkBFgrGb2yE3XeHPfdsOVdWpm3pm+So7NChiQs4w8gSAvvP65yrF5GPgP2AJ8NO7Qc8Hes6r7f5T+AE59zz6W1lcKQswDZsgM8+g0WL4KWXtNfbtasG44dy/+y4ow5XxloexDA6EHlrk5qadObj4sXqdRo2LGHvF6QgwABeeQU++EDve8EFLccbvfRusdaJNIw8J2j3x7XAv0RkPvAu8Es0589tACJyOTDKOfdT7/0J6CyoXwPviUiox1qTT+vEpSzA6uvh+uvV41VTo7FfXbtqb/TLL7XM0KEq0MaNs9QUhtFCftmkkPdr+XL1NIUC4RMkJQFWXa0LaYMG448e3XIulDYnnmWIDCMPCVSEOeceFpFBwKXACGABcLBz7nuvyAhaJ2I8HW3zHG8LcR9wctobnAFSFmB1dfDXv2rcV79+6uYPDTVWVKhBDE0Jv+sumD07qV6xYeQjeWOTmppaZgnW1OiMwniD3iNISYCB2pnCQg2L+NnPWo47px77zTePvk6kYXQAgvaE4Zy7BbglxrmTI97PzECTAiNlAQY6++i557R3GTnde9kyfR0+XLNIf/MN/PvfsPPO6hEzDCP3bVJ5uQ49rlypHbFwz1OCpCzAli5VGwOaEyw8+35lZdvrRBpGB8ASqGQJvgiwNWt0KaLu3Td27zc2as8YYMIE7SmPGaM91Bdf9OUZDMMIkMZGzQn4/vs6KWfEiJRCDVIWYKFg/KYmDcTfbbfW58vL1TsWa51Iw+gAmAjLAnwRYACvvqpCbJNNNk5QuGqVpqwITQWvr9f9wYM1aLaszI9HMQwjCMrK4OOP4dNPdWhv9Oikhx/BBwEGOiHoo4+0UxgejA+6TmTXrioUDaMDYyIsYHwTYKBLeXTqFD1Hznff6ev48RoT1tysPdChQ3WG5KefJvsIhmEERWMjfPuter/WrNFZhilOtEkqE34kVVVw3XW6f+qpG89+jGedSMPoAAQeE9aR8VWAVVXpEETv3htnvi4rg9JSGoEFm2zCto2NGqzfrVvroH3DMHKH0lKN/Vq1SsXM4MEpV+mLBwzgjjuguFhDHk46qfW5RNaJNIw8x0RYQPgqwECD7kPiKto54Ang6eXLebBfPw2QDR+usPXVDCM3aGzUpceWLNFhvVGjfEm27JsAW7IEHn5Y9y+6aOPZ1xUVvolGw8h17D9vAPguwIqL1SiPGqVTvsNpaNAM2cA3gwdz3qRJ2hMNBcPW1akAs4SthpH9lJTA/PmaiLlbN/UmZZMAcw6uuEJtzF576czrSCortd2WGscwzBOWaXwXYM3NGhPS1KRrQX7zTWuP2MqV+r5PH36/225avr6+ZS3JtWs1l9j06ak/nGEY6aGhocX7VV/vm/cLfBRgAC+8oPGlPXpsHIwP2kns2VOz9huGYZ6wTOK7AAMVUatXa4D9tGmamqKwMHTDjQPyGxpUgHXtqvuVlXDQQa3z9xiGkT2EvF+ff66/3WwVYJWVcMMNuv+LX2g+wkhKS/V4sutEGkaeYSIsQ6RFgNXXq+crJKr69FFB1dioQ5RlZZqLp1OnlsV6Gxu1J9rYCF9/DVttBfvvn/LzGYbhMw0NGnj//vtQVKTiq18/36r3VYAB3HabCsZx4+DEEzc+39ioHUNbJ9IwfsCGIzNAWgQYaKxXaAZSiD331KWJnn9es1WDGr3u3XUo0jntjX7/PWyxBfz2t7Zum2FkG8XFKsDWrNF4TZ89R74LsK+/hkcf1f2LLoqeo6ysTG2N2RvD+AETYWkmbQKsslJFVv/+rWc2isDBB2uP+Te/0WPdu2vcWGOjlt10UzjiCDj6aB3GNAwjO6iv1xCCpUs1ztPHoccQvguw5ma48kp93W8/2HHHjcs4p53DqVNtnUjDCMNEWBpJmwALxXpVV7cMM4Yjoj3opiadhXTyyTqVvbxchx+PPdaSJBpGtlFcrL/btWvVW9Snj++38F2AgXrdP/9cwxzOPTd6GVsn0jCiYiIsTaRNgIEa6xUrYhs05+Dxx3X/+ONb4sTWrtUZlCbADCP7mD9fhczo0WnxFqVFgFVUwI036v5pp8We9VheDpMnR1/NwzA6MCbC0kBaBVhTkw5VhJYdisbnn2uZ7t11aBI0o37fvibADCNb6d07bakb0iLAAG69VWNMJ0yAE06IXqauTodUo82WNIwOjs2O9Jm0CjDQ4Yo1a9qO5XriCX3df/+W2VRVVWrgU1jU1zCMNNK7d1qqTZsAW7QIHntM9y++OLZtKS3V7PgDBqR+T8PIM0yE+UjaBVhdnXq4evaMHaxbXg5z5+r+0Ufra3OzvtqsJMPoUKRNgIWC8Z2DAw6AGTNil2to0Bnctk6kYWyEiTCfSLsAA00rUVzc9hJDzz+vM6wmT9YUFKAB/L17W0/UMDoQ4QJs2pQp/gkwgGeegQUL1K7ECsYHjRnr188C8g0jBibCfCAjAqyiQmdEDhgQe7Ht8ID8o45q6XlWVupwQGipIsMw8ppID9grt9zinwArK4ObbtL9X/6ybYFVXq4zuG2dSMOIiomwFMmIAHNO83xt2NB2xuyPP1ZvWc+ecOCBLdc2Nlo+MMPoIKRtCDLEnDkqrgoK4LjjYpcLrRNptscwYmIiLAUyIsBA14JcsaJ9Yxbygh14YEuOoQ0bdFq4zYo0jLwn7QJswQJ46ind/+1v204kG1on0sellgwj3zARliQZE2CNjRqML9L2ItslJfDaa7ofCsgHHYocNMjy8xhGnpN2AdbUBFddpd71gw+GadNil7V1Ig0jLkyEJUHGBBhoOop169oPbH32WTV8m2+uS4OEqK9PW+4hwzCyg7QLMFAP2JdfajD+7Nltly0r0wlEgwf72wbDyDNMhCVIRgVYbS0sWaJerLbc/s3N8OSTun/UUa2v79HDhiINI4/JiAArLdVYMIAzz2xbXIXWiRw71taJNIx2MBGWABkVYKBB9qWlbaekAPjwQ1i5UnuoBxzQcryyUgVYGtagMwwjeDIiwABuvllnaE+eDMcc03bZqiq1OZaWwjDaxURYnGRcgJWXa0qKQYPaT3IYCsg/+ODWSxnV1MCIEZYk0TDykIwJsC++gKef1v2LL27bKw9qu0aNsjhUw4gDE2FxkHEBFkpJUVfXvherqAjefFP3wwPy6+s1N48NRRpG3pExAdbUBFdcofuHHQbbbNN2+bo6HYK0dSINIy5MhLVDxgUYwPr1OrwYT1Dr00+rodx6a5g0qeV4ZaVODbfp4YaRV2RMgIF62b/+Gvr2hXPOab+8rRNpGAlhIqwNAhFgDQ2akqJTp7ZTUoCKr1DOnnAvGGh+sBEjYmfXNwwj58ioACspgVtu0f2zzmo/NrW5WT3wY8aY3TGMOLFfSgwCEWAAq1fD2rXxecHefVdTWPTrB/vs03K8sVGNoPVGDSNvyKgAA12aqKpKU96Ez7qORWidSEtLYRhxYyIsCoEJsA0bNCVF377tB78CPPGEvh56aGuvWVWV1pFOA20YRsbIuAD79FPNPQgajB9PqomKCvWC2Rq1hhE3JsIiCEyAASxbpjOL4vFgrV0L//uf7kf2UquqNDA2HiFnGEZWk3EB1tgIV16p+7NmwVZbtX9NTY12BC0xtGEkhImwMAIVYKWlmhcsnpQUoAH5zc2w3XYwfnzL8aYmfR00KC3NNAwjc2RcgAE89hh884160n/1q/iuKSlRAWbed8NICBNhHoEKsOZmTUlRXx9fYtXGxpaA/EgvWHW11mGpKQwjpwlEgBUVwa236v7ZZ8dnR0LrRI4aldamGUY+YiKMgAUY6NqQK1fGn2H6f/+DwkIdttxrr9bnKith6FDNEWYYRk4SiAADuOEG7chtvjkccUR815SVqS0y77thJEyHF2GBC7BQSoquXeMPaA0F5B92WGux5Zx61Wx2kmHkLIEJsI8/hhdf1HCIeIPxw9eJtBhUw0iYDi3CAhdgoB6w9evjF06rVmlqCog+FNmrlw1FGkaOEpgACw/GP/JI2GKL+K4LhT8MHZq+thlGHhO4CBORs0TkOxGpFZGPRGT3dsrv6ZWrFZFvReSMZO4biABzDhYvhpdegvff12D8pUs1mDWeXidoLJhzsNNOMHp063OVlSrmwtePNAwjIYKySYu++y7zAmzDBvWAXXlliy06++z4ry8rg5EjbZ1Iw0iSQP3HInIccD1wFjAPOB14UUQ2d84tj1J+AvACcCfwE2BX4BYRKXTOPZ7IvS+88ELuvPPOzAmwxkbNPv3CC5pComtXzWi/zz7tr8cWoqEBnnlG96MlT2xosB6pYaRAkDbpkPPOo7isLHMC7Ntv4b77dFb2l1/qsX33jX+ps/p6TQo9YkT62mgYeU7QnrDzgbudc3c55xY5584FVgBnxih/BrDcOXeuV/4u4J/AbxK9cUYFGMAbb2gsV+/eGvQ6ZAh88ol6xOJJSQG6UHdxsQbA7rFH63MbNqgHzIYiDSMVArNJGRVgjY3w0EMa3lBSorGkffroChzffx9fHaWlasdsZQ7DSJrARJiIdAO2A16OOPUysEuMy3aOUv4lYIaIdI1xn+4i0i+0AX1D5zIaA/beez8Ezc8vLOTEd99lfufOOjxZXR1fHaGA/FmzNg6CraxUAda7t5+tNowOQ9A2aZtNN81cDNiKFRqP2rUrrF5NM/Dl5Mnqpf/qq/avb26GujoNibB1Ig0jaYL89QwGOgPrIo6vA4bHuGZ4jPJdvPqicQlQHratBJgzZ05mg/DDvF3Xffopj6xaxXXr18d//fLl8MEHWs+sWRufr63VYYF4vWqGYUQSqE166pprMiPAQjQ3/zAMeRvwt0TsUWidyHjT6hiGEZVs6MK4iPcS5Vh75aMdD3E50D9sGw3wk5/8JLFWpspOO2nw/bp1nDd6NMcOHcp5/frpkiDxeK+efFJfd9114xiM2lpNb2FDkYbhB4HYpIHxxmL5wZgxajeqq2no2pVPRo7kvAEDdM3ZzTZr/3pbJ9IwfCHIwPwioImNe5hD2bhnGWJtjPKNQHG0C5xzdUBd6L0E5Snac0/tdT7+ODOKinhw9GjYdFPN9dUedXVtB+RXVemspr59Nz5nGEa8dBybVFQES5YA0HXECO4cMEBtyGGHac6vtqipUfFl60QaRsoEJsKcc/Ui8hGwH/Bk2Kn9gKdjXPYuEKla9gfmO+ca/G+lj3TpAj/7mbrw161Tb9bkyfElOHz9dV3Ye9gw9YRFsmEDTJpksRmGkQIdyiZdf716wrbaCs45RwP1J06ML19haammpcik584w8pSgUxxfC/xLROajxuyXwFg0RAERuRwY5Zz7qVf+NuBXInItOiV8Z+DnwAmZbnjCOKezjvr1gy23TOzax72Z7rNmbZxPrKFBhZzNUDIMP8h/m/Tee/DKK9ppu+QS7QzGS2OjxpKNGmXxp4bhA4GKMOfcwyIyCLgUGAEsAA52zoXmSI9ADWCo/HcicjBwHXA2sBqYnWg+nkAoLtbg+kQDWZcu1VQWnTtHX8utslKFXSYDeg0jT8l7m1RfD1ddpfvHHpuYAANNzrrJJrY0mmH4RNCeMJxztwC3xDh3cpRjbwLT09wsf2lqUjHV1JR4NvtQQP7uu0dPxBpaty3ejPuGYbRJXtukBx/UzuCgQXBGgon9Q+tEbrqprRNpGD5hQUSZYM0a3RINZK2theee0/1oAfmNjTokMHBg6m00DCO/WbsW7rpL93/9a03OmgjV1TqT2wLyDcM3TISlm7o69YL17Jl473HuXJ35OHKkpriIpKpKDamlpjAMoz2uuUbt0fTpcNBBiV8fWifSEkIbhm+YCEs3y5drPFgy3qpQQP6RR0af+VhVpb3SrlETcxuGYSjz5uks686d4aKLEg+qt3UiDSMtmAhLJxUVukjugAGJp49YvBgWLFCjefjhG59vbtYYDQuQNQyjLerq4Oqrdf/44zWdTaKUlamtsdAHw/CVDhtdWVFRkd4bOKciqqhI11erqkrs+ocf1tfdd9fEiJHXh9536qRiL0M0NDSk/7MzDDLwG80yKuJdQzZR7r9f14kcOBB+/OPEbVFzsy7yPX584temGbNHRqZI1/dMnGtrNY78Q0TGA98F3Q7DMOJmgnNuWdCNSBdmkwwjp/DVHnVET1iJ9zoaqAyyIUnQF13s19qeWaztwRBqe0l7BXMcs0nBYG0Phlxte1rsUUcUYSEqnXM55ccOW2PO2p5BrO3BENg6r8GRy38ja3sGsbZnnnTZIwvMNwzDMAzDCAATYYZhGIZhGAHQEUVYHfBn7zXXsLYHg7U9GHK57YmQy89pbQ8Ga3vmSUu7O9zsSMMwDMMwjGygI3rCDMMwDMMwAsdEmGEYhmEYRgCYCDMMwzAMwwgAE2GGYRiGYRgBkJciTETOEpHvRKRWRD4Skd3bKb+nV65WRL4VkTMy1dYobYm77SJylIjMFZFCEakQkXdF5IBMtjeiPQl97mHX7SoijSLyaZqb2FYbEv3OdBeRv4nI9yJSJyJLReTUTLU3oi2Jtv1EEflMRDaIyBoRuUdEBmWqvV4b9hCRZ0VktYg4EZkVxzVZ8ztNFLNJwZCrNsnsUWbtkdeOYGyScy6vNuA4oB74BbAZcD1QBYyNUX4CUO2V28y7rh44Ogfafj1wEbA9sCnwd+/6adne9rDr+gNLgZeAT3PhO+Nd8zTwHrAvMB7YAdgl29sO7AY0AbO97/5uwALgyQy3+yDg/4CjAAfMaqd81vxOM/A3yppnNZuUeZtk9ijz9shrSyA2KaMPmaEP8n3g1ohji4DLY5S/ElgUcew24N1sb3uMOhYCl+ZK24GHgL8ClwVh8JL8zhwIlAEDg2hvim3/DbA04tg5wIoAnyEeg5c1v9MM/I2y5lnNJgUiwsweBWiPvDZkzCbl1XCkiHQDtgNejjj1MrBLjMt2jlL+JWCGiHT1t4WxSbLtkXV0QhcZzeiCx8m2XUROAQrQBHiBkGTbDwfmAxeJyCoRWSwi/xCRnmls6kYk2fZ3gNEicrAow4BjgOfT11JfyIrfaaKYTTKblAhmj3LGHoFPv9N8W8B7MNAZWBdxfB0wPMY1w2OU7+LVt8bPBrZBMm2P5AKgN/CIj+2Kh4TbLiKbAlcAuzvnGiW4xZqT+dwnom7zWuBIr45bgIFAJuMwEm67c+4dETkReBjogX7Pn0F7n9lMtvxOE8VsktmkRDB7lBv2CHz6neaVJyyMyGUAJMqx9spHO54JEm27FhI5AXWfH+ecW5+GdsVDXG0Xkc7Av4E/OecWZ6JhcZDI597JO3eic+4D59wLwPnAyZnufXrE3XYR2Ry4EfgL2ms9EI1tuC2dDfSJbPqdJorZpGDIVZtk9ig3SPl3mm+esCI0yC9SdQ9lY8UaYm2M8o1Asa+ta5tk2g6AiBwH3A38yDn3Snqa1yaJtr0vMAOYJiI3e8c6ASIijcD+zrnX0tXYCJL53NcAq5xz5WHHFqE/wNHAN343MgbJtP0SYJ5z7mrv/eciUg28LSJ/cM5lq0cpW36niWI2yWxSIpg9yg17BD79TvPKE+acqwc+AvaLOLUfOvYcjXejlN8fmO+ca/C3hbFJsu2h3ua9wI+dc4GMoyfR9gpgK2DbsO024Gtv//20NDQKSX7u84CRItIn7NhkoBlY6XsjY5Bk23uh7QynyXsNbEw4DrLid5ooZpPMJiWC2SMgN+wR+PU7DXIGQppmNYSmyJ6KThu9Dp0iO847fzlwf1j50DTTa73ypxL8dPB4234C0ACchSry0NY/29se5frLCD5FRbyfex9gBfAosDmwB7AYuDMH2n6y9505E40l2RX4EHg/w+3uQ8s/Owec5+2PjdHurPmdZuBvlDXPajYp0BQVZo8y2/ZAbFJGHzKDH+ZZwDKgDlXme4Sduxd4I6L8nsDHXvnvgDNyoe3AG96XJXK7N9vbHuXaQAxeCt+ZqcBcYINnAK8BeuZI289B0wZsAFYDDwCjMtzmmW19d7P9d5qBv1HWPKvZpOxvt9kjX9odiE0SryLDMAzDMAwjg+RVTJhhGIZhGEauYCLMMAzDMAwjAEyEGYZhGIZhBICJMMMwDMMwjAAwEWYYhmEYhhEAJsIMwzAMwzACwESYYRiGYRhGAJgIy3FEuUNESkTEici2QbcpVxGR8X58hiJyr4g81U6ZN0Tk+rD3y0Tk3LD3TkRmpdIOwwgCs0n+YTYp/8m3Bbw7IgeiSz/MBL5FF1HNSURkGXC9c+76gJuSCY5Cl+uIxQigFNQQo9mYpznnPk17ywwjNcwm5SZmkwLARFjuUwCscc7FXFC3PUREgM7OuUb/mpVdiEhXl0ULPTvnSto5vzZTbTEMnzGbFAdmkwyw4cicRkTuBW4Cxnqu4mXe8e4icqOIrBeRWhH5n4hsH3bdTK/8ASIyH133andvGOEiEflWRGpE5DMROSbinluIyPMiUiEilSLytogUeOe2F5G5IlIkIuUi8qaITI+4/jIRWS4idSKyWkRu9I6/AYwDrvPaFnM9Le/8mSLyotfO70TkR2HnQy78Yz0Xey3wExHpJCKXishK7/6fisiBUW4xVUTe8T67hSIyM6zuziJyt3fPGhH5WkR+HaOdf/L+BhUicruIdAs718r1H+MZZ3lvv/NeP/GOvyEie4hIg4gMj7juGhF5K1a9hpFOzCaZTTKblCBBLPBpm28LjvYH/ogu2DocGOIdvwFYBRwEbI4uPFoCDPTOz0QXJv0M2A/tuQ4C/gYsAg5AV7Q/GagF9vSuGwUUA48DM4DJwCnAFO/83sBP0BXlNwPuAtYCfb3zxwDlXrvGAjsAp3nnBnrP8UfvWYa38dwOHeL4hdeGvwKNwGbe+fFeme9QF/sEYCRwnnf/44EpwJXoqvebRly3Ajjae4Y7gQpgkFemK/BnYHuv3hOBauDYsPbdC1QCDwFbAIcA64G/hZV5Ax3mCL1fBpwb8YyzvP3tvff7eJ9N6O/4NXBh2DVdgHXAKUF/N23rmBtmk8wmmU1K7DcTdANsS/EPCOcCy8Le9/Z+xD8OO9YVNYAXeu9nej+gIyKuqwF2jqj/LuDf3v7f0RiPrnG2rbNnLA713p/v/UijXh/5o2+jXgfcGnHsPeAWbz9kuH4dUWYV8LuIYx8AcyKuuzjsfBfPAF7URnvmAI+Fvb8X/cfQK+zYGZ4R7OS9T8Tghdq1bcR9LwK+DHt/hHeP3kF/L23ruJvZpB+OmU0ym9TuZsOR+UcBauDmhQ44jTv4AO1FhTM/bH9zoAcwV0SqQhvwU69OgG2Bt12MOAYRGSoit4nIYhEpR3t4fdAeJsCjQE/gWxG5U0SOFJFk4xLfjfI+5vOJSD+05zkvosy8KNf9ULfTmJT54WVE5AwRmS8ihd5ndBotzxjiM+fchog6+wBj2nqoBLkXmCQiO3nvTwUecc5V+3gPw0gVs0lRns9skgEWmJ+PiPcaGb8gUY6F/zBCgvwQtHcWTp33WtPOve8FhqA94e+9694FugE451aIyBR0uGFf4BbgQhHZM5YRTZC2ni9WmWifS8y6ReRY4DrgAvTZKoELgR2TbGPSOOfWi8izwCki8i1wMOpRMIxswmxSC2aTjFaYJyz/WIK6/ncLHRCRrmi8xKI2rvsSNVBjnXNLIrYVXpnP0WDZrjHq2B240Tn3gnNuoVff4PACzrka59wzzrnZ6I9zZ2Ar73Q9OlwQDztFef9VrMLOuQpgNWGfi8cubPy5/FC31yveLqzu3YF3nHO3OOc+cc4toaVXHs42ItIzos4qYGWsNrZBvfca7bO5C40nOR1Y6pyL7FUbRtCYTYqC2SQDzBOWdzjnqkXkVuBqESkBlqPj9L2Au9u4rlJE/oHOBOoE/A/ohxqEKufcfcDNwDnAQyJyOera3wn4wDn3NWpsTxKd3dQPuJqwnqqInIz+aN8HNgAneee/94osA/YQkYeAOudcW/mFfuTd539oIOoOwM/b+XiuBv4sIkuBT9EA3m2968M5W0S+QQ3hecAA4J/euSXAT0XkADTI9iQ0SPW7iDq6AXeLyP+hM6z+DNzsnGtup43RWI9+TgeKyEqg1jlX7p17Cf07/AG4NIm6DSOtmE1qE7NJHZ2gg9JsS20jIgjWO9YDuBEoRGcS/Q/YPuz8TNQFvUnEdQLMRntY9egP7b/AHmFltkZ/ZNVogOtbwETv3DTgQ++ei9GZR8vwgjuBWWiwajnaA3sX2Ces7p3Q2VG1+tWM+cwOOAt42Su7DDg+7Px4ogeNdkKNwkrv+T4FDoxy3QleO+vQ3vjeYWW6A/cAZWjiwluAy4FPw8rcCzyFGrkidHjgTqB7WJk3iDMI1nv/C/SfVxPwRsRz/QWdiTUi6O+jbbaZTTKbZDYp/k28D8wwcgYvX8+Rzrmngm5LNiAidwLDnHOHB90Ww+iImE1qjdmk+LHhSMPIUUSkPzrscCI6FdwwDCMwzCYljokww8hdnkbjTm53zs0NujGGYXR4zCYliA1HGoZhGIZhBIClqDAMwzAMwwgAE2GGYRiGYRgBYCLMMAzDMAwjAEyEGYZhGIZhBICJMMMwDMMwjAAwEWYYhmEYhhEAJsIMwzAMwzACwESYYRiGYRhGAJgIMwzDMAzDCID/B045INRGMouNAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 600x400 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "#create reliability plot\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from sklearn.linear_model import LinearRegression\n",
    "\n",
    "fig, axs = plt.subplots(1,2, figsize=(6,4), dpi=100)\n",
    "axs[0].set_aspect(1)\n",
    "axs[1].set_aspect(1)\n",
    "fig.tight_layout()\n",
    "fig.suptitle('DJF average DePreSys3 '+s)\n",
    "\n",
    "t=np.full(50,0.3333333)\n",
    "x=np.linspace(0,1,50)\n",
    "y=np.linspace(0,1,50)\n",
    "yt=np.linspace(-0.05,1.05,50)\n",
    "for i in range(2):\n",
    "    axs[i].plot(x,t,c='black',linewidth=0.5)\n",
    "    axs[i].plot(t,yt,c='black', linewidth=0.5)\n",
    "    axs[i].plot(x,y,c='black')\n",
    "    axs[i].set_xlim([0, 1])\n",
    "    axs[i].set_ylim([-0.05, 1.05])\n",
    "    axs[i].set_xlabel('forecast probability')\n",
    "    axs[i].set_ylabel('observed frequency')\n",
    "    \n",
    "#upper tercile data unadjusted\n",
    "upperun=np.loadtxt(\"NAO_\"+s+\"dp3_coords_upper_tercile.txt\")\n",
    "wun=int(upperun.size/4)\n",
    "NAO_upun = upperun.reshape(wun,4)\n",
    "x_uun = NAO_upun[:,0]\n",
    "y_uun = NAO_upun[:,1]\n",
    "size_uun= NAO_upun[:,2]*8\n",
    "\n",
    "pun=x_uun.reshape(-1,1)\n",
    "modeluun= LinearRegression()\n",
    "modeluun.fit(pun,y_uun,size_uun)\n",
    "p_newun = x_uun\n",
    "q_newun = modeluun.predict(p_newun[:, np.newaxis])\n",
    "q_regun=modeluun.coef_\n",
    "print(q_regun)\n",
    "q_intun=modeluun.intercept_\n",
    "axs[0].plot(x,q_regun*x+q_intun, c='red')\n",
    "\n",
    "uncert_regun =np.loadtxt(\"uncert_reg_u\"+s+\"dp3.txt\")\n",
    "\n",
    "m1un=uncert_regun[0]\n",
    "m2un=uncert_regun[1]\n",
    "c1un=0.3333333-0.3333333*m1un\n",
    "c2un=0.3333333-0.3333333*m2un\n",
    "\n",
    "axs[0].fill_between(x, m1un*x+c1un, m2un*x+c2un,color='red', alpha=0.2)\n",
    "axs[0].scatter(x_uun, y_uun, c='red', s=size_uun, alpha=0.6)\n",
    "axs[0].scatter(x_uun, y_uun, c='black', s=0.5)\n",
    "\n",
    "\n",
    "#lower tercile data unadjusted\n",
    "lowerun=np.loadtxt(\"NAO_\"+s+\"dp3_coords_lower_tercile.txt\")\n",
    "wlun=int(lowerun.size/4)\n",
    "NAO_loun = lowerun.reshape(wlun,4)\n",
    "x_lun = NAO_loun[:,0]\n",
    "y_lun = NAO_loun[:,1]\n",
    "size_lun= NAO_loun[:,2]*8\n",
    "\n",
    "plun=x_lun.reshape(-1,1)\n",
    "modellun= LinearRegression()\n",
    "modellun.fit(plun,y_lun,size_lun)\n",
    "p_newlun = x_lun\n",
    "q_newlun = modellun.predict(p_newlun[:, np.newaxis])\n",
    "q_reglun=modellun.coef_\n",
    "print(q_reglun)\n",
    "q_intlun=modellun.intercept_\n",
    "axs[1].plot(x,q_reglun*x+q_intlun, c='red')\n",
    "\n",
    "uncert_reglun =np.loadtxt(\"uncert_reg_l\"+s+\"dp3.txt\")\n",
    "\n",
    "m1lun=uncert_reglun[0]\n",
    "m2lun=uncert_reglun[1]\n",
    "c1lun=0.3333333-0.3333333*m1lun\n",
    "c2lun=0.3333333-0.3333333*m2lun\n",
    "\n",
    "axs[1].fill_between(x, m1lun*x+c1lun, m2lun*x+c2lun, color='red', alpha=0.2)\n",
    "axs[1].scatter(x_lun, y_lun, c='red', s=size_lun, alpha=0.6)\n",
    "axs[1].scatter(x_lun, y_lun, c='black', s=0.5)\n",
    "\n",
    "\n",
    "axs[0].set_title('upper tercile')\n",
    "axs[1].set_title('lower tercile')\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}