Month to month correlations.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Sea Ice Autocorrelation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<xarray.DataArray 'goddard_merged_seaice_conc_monthly' (time: 469, ygrid: 332, xgrid: 316)>\n",
       "[49203728 values with dtype=float32]\n",
       "Coordinates:\n",
       "  * ygrid      (ygrid) float32 4337500.0 4312500.0 ... -3912500.0 -3937500.0\n",
       "  * xgrid      (xgrid) float32 -3937500.0 -3912500.0 ... 3912500.0 3937500.0\n",
       "    latitude   (ygrid, xgrid) float64 ...\n",
       "    longitude  (ygrid, xgrid) float64 ...\n",
       "  * time       (time) datetime64[ns] 1978-11-01 1978-12-01 ... 2017-12-01\n",
       "Attributes:\n",
       "    valid_range:    [  0 100]\n",
       "    long_name:      Goddard Edited Climate Data Record of Passive Microwave M...\n",
       "    standard_name:  sea_ice_area_fraction\n",
       "    units:          1\n",
       "    flag_values:    [-5 -4 -3 -2 -1]\n",
       "    flag_meanings:  pole_hole lakes coastal land_mask missing_data\n",
       "    datum:          +ellps=urn:ogc:def:crs:EPSG::4326\n",
       "    grid_mapping:   projection\n",
       "    reference:      http://nsidc.org/data/nsidc-0051.html http://nsidc.org/da..."
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#Sea ice autocorrelation\n",
    "\n",
    "import xarray as xr\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt \n",
    "import scipy.stats\n",
    "\n",
    "#Setting data\n",
    "data=xr.open_dataset(\"/g/data/v45/sl0603/sic_combined.nc\")\n",
    "cice_monthly=data.goddard_merged_seaice_conc_monthly\n",
    "cice_monthly"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Correlation coefficient\n",
    "The problem here is to calculate the correlation for each point in space using 2 arrays as input: month and month+1.\n",
    "Here is the correlation function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def corr2(ar,axis):\n",
    "    '''Correlation coefficient between 2 timeseries concatenated together\n",
    "    ar: 2D array, time is the right dimension. Left dimension must be size 2. \n",
    "        ar[0,:] is the first timeseries\n",
    "        ar[1,:] is the second timeseries\n",
    "    stats: 1D array containing the correlation coefficient and the p-value.'''\n",
    "    # Remove any degenerated dimensions\n",
    "    ar1=ar[0,...].squeeze()\n",
    "    ar2=ar[1,...].squeeze()\n",
    "    stats = xr.DataArray(np.asarray(scipy.stats.pearsonr(ar1,ar2)),dims='c')\n",
    "    stats = stats.expand_dims(dim='z',axis=1)\n",
    "    return stats"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This function can be called for all space point with something like:\n",
    "```python\n",
    "out = in.groupby('z').reduce(corr2, dim='time')\n",
    "```\n",
    "Where `z` is the spatial dimension. `in` is a 3D array: the dimension for the 2 different timeseries, space dimension and time. Then out is a 2D array of the stats for each spatial point.\n",
    "\n",
    "The problem here is to create an array with 1 dimension for space and the 2 months timeseries concatenated together as needed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<xarray.DataArray 'goddard_merged_seaice_conc_monthly' (time: 469, ygrid: 332, xgrid: 316)>\n",
       "[49203728 values with dtype=float32]\n",
       "Coordinates:\n",
       "  * ygrid      (ygrid) float32 4337500.0 4312500.0 ... -3912500.0 -3937500.0\n",
       "  * xgrid      (xgrid) float32 -3937500.0 -3912500.0 ... 3912500.0 3937500.0\n",
       "    latitude   (ygrid, xgrid) float64 ...\n",
       "    longitude  (ygrid, xgrid) float64 ...\n",
       "  * time       (time) datetime64[ns] 1978-11-01 1978-12-01 ... 2017-12-01\n",
       "    month      (time) int64 11 12 1 2 3 4 5 6 7 8 9 ... 2 3 4 5 6 7 8 9 10 11 12\n",
       "Attributes:\n",
       "    valid_range:    [  0 100]\n",
       "    long_name:      Goddard Edited Climate Data Record of Passive Microwave M...\n",
       "    standard_name:  sea_ice_area_fraction\n",
       "    units:          1\n",
       "    flag_values:    [-5 -4 -3 -2 -1]\n",
       "    flag_meanings:  pole_hole lakes coastal land_mask missing_data\n",
       "    datum:          +ellps=urn:ogc:def:crs:EPSG::4326\n",
       "    grid_mapping:   projection\n",
       "    reference:      http://nsidc.org/data/nsidc-0051.html http://nsidc.org/da..."
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Create a new coordinate with the month number.\n",
    "months=cice_monthly['time.month']\n",
    "cice_monthly = cice_monthly.assign_coords(month=('time',months))\n",
    "cice_monthly"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<xarray.DataArray 'goddard_merged_seaice_conc_monthly' (time: 469, ygrid: 332, xgrid: 316)>\n",
       "array([[[0., 0., ..., 0., 0.],\n",
       "        [0., 0., ..., 0., 0.],\n",
       "        ...,\n",
       "        [0., 0., ..., 0., 0.],\n",
       "        [0., 0., ..., 0., 0.]],\n",
       "\n",
       "       [[0., 0., ..., 0., 0.],\n",
       "        [0., 0., ..., 0., 0.],\n",
       "        ...,\n",
       "        [0., 0., ..., 0., 0.],\n",
       "        [0., 0., ..., 0., 0.]],\n",
       "\n",
       "       ...,\n",
       "\n",
       "       [[0., 0., ..., 0., 0.],\n",
       "        [0., 0., ..., 0., 0.],\n",
       "        ...,\n",
       "        [0., 0., ..., 0., 0.],\n",
       "        [0., 0., ..., 0., 0.]],\n",
       "\n",
       "       [[0., 0., ..., 0., 0.],\n",
       "        [0., 0., ..., 0., 0.],\n",
       "        ...,\n",
       "        [0., 0., ..., 0., 0.],\n",
       "        [0., 0., ..., 0., 0.]]], dtype=float32)\n",
       "Coordinates:\n",
       "  * ygrid      (ygrid) float32 4337500.0 4312500.0 ... -3912500.0 -3937500.0\n",
       "  * xgrid      (xgrid) float32 -3937500.0 -3912500.0 ... 3912500.0 3937500.0\n",
       "    latitude   (ygrid, xgrid) float64 -39.36 -39.49 -39.62 ... -41.72 -41.58\n",
       "    longitude  (ygrid, xgrid) float64 -42.23 -42.05 -41.87 ... 135.4 135.2 135.0\n",
       "  * time       (time) datetime64[ns] 1979-01-01 1980-01-01 ... 2017-12-01\n",
       "    month      (time) int64 1 1 1 1 1 1 1 1 1 1 ... 12 12 12 12 12 12 12 12 12\n",
       "Attributes:\n",
       "    valid_range:    [  0 100]\n",
       "    long_name:      Goddard Edited Climate Data Record of Passive Microwave M...\n",
       "    standard_name:  sea_ice_area_fraction\n",
       "    units:          1\n",
       "    flag_values:    [-5 -4 -3 -2 -1]\n",
       "    flag_meanings:  pole_hole lakes coastal land_mask missing_data\n",
       "    datum:          +ellps=urn:ogc:def:crs:EPSG::4326\n",
       "    grid_mapping:   projection\n",
       "    reference:      http://nsidc.org/data/nsidc-0051.html http://nsidc.org/da..."
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Then we can sort along the month coordinate to have all the data for January, then February etc.\n",
    "cice_monthly = cice_monthly.sortby('month')\n",
    "cice_monthly"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<xarray.DataArray 'goddard_merged_seaice_conc_monthly' (time: 198, ygrid: 332, xgrid: 316)>\n",
       "array([[[0., 0., ..., 0., 0.],\n",
       "        [0., 0., ..., 0., 0.],\n",
       "        ...,\n",
       "        [0., 0., ..., 0., 0.],\n",
       "        [0., 0., ..., 0., 0.]],\n",
       "\n",
       "       [[0., 0., ..., 0., 0.],\n",
       "        [0., 0., ..., 0., 0.],\n",
       "        ...,\n",
       "        [0., 0., ..., 0., 0.],\n",
       "        [0., 0., ..., 0., 0.]],\n",
       "\n",
       "       ...,\n",
       "\n",
       "       [[0., 0., ..., 0., 0.],\n",
       "        [0., 0., ..., 0., 0.],\n",
       "        ...,\n",
       "        [0., 0., ..., 0., 0.],\n",
       "        [0., 0., ..., 0., 0.]],\n",
       "\n",
       "       [[0., 0., ..., 0., 0.],\n",
       "        [0., 0., ..., 0., 0.],\n",
       "        ...,\n",
       "        [0., 0., ..., 0., 0.],\n",
       "        [0., 0., ..., 0., 0.]]], dtype=float32)\n",
       "Coordinates:\n",
       "  * time       (time) datetime64[ns] 2000-02-01 2001-02-01 ... 2017-12-01\n",
       "  * ygrid      (ygrid) float32 4337500.0 4312500.0 ... -3912500.0 -3937500.0\n",
       "  * xgrid      (xgrid) float32 -3937500.0 -3912500.0 ... 3912500.0 3937500.0\n",
       "    latitude   (ygrid, xgrid) float64 -39.36 -39.49 -39.62 ... -41.72 -41.58\n",
       "    longitude  (ygrid, xgrid) float64 -42.23 -42.05 -41.87 ... 135.4 135.2 135.0\n",
       "    month      (time) int64 2 2 2 2 2 2 2 2 2 2 ... 12 12 12 12 12 12 12 12 12\n",
       "Attributes:\n",
       "    valid_range:    [  0 100]\n",
       "    long_name:      Goddard Edited Climate Data Record of Passive Microwave M...\n",
       "    standard_name:  sea_ice_area_fraction\n",
       "    units:          1\n",
       "    flag_values:    [-5 -4 -3 -2 -1]\n",
       "    flag_meanings:  pole_hole lakes coastal land_mask missing_data\n",
       "    datum:          +ellps=urn:ogc:def:crs:EPSG::4326\n",
       "    grid_mapping:   projection\n",
       "    reference:      http://nsidc.org/data/nsidc-0051.html http://nsidc.org/da..."
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Since there is one year missing for January we are going to remove January for now and treat January \n",
    "# separately afterwards\n",
    "# We are going also going to remove the extra year for November and December\n",
    "#drop_year = cice_monthly.where(cice_monthly['time.year']>1978,drop=True)\n",
    "drop_year = cice_monthly.where(cice_monthly['time.year']>1999,drop=True)  #shorter timeseries for testing\n",
    "no_JAN = drop_year.where(cice_monthly.month>1,drop=True)\n",
    "no_JAN"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "18.0"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Check we likely have same number of values for each month:\n",
    "no_JAN.shape[0]/11"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<xarray.DataArray 'goddard_merged_seaice_conc_monthly' (time: 198, z: 104912)>\n",
       "array([[0., 0., 0., ..., 0., 0., 0.],\n",
       "       [0., 0., 0., ..., 0., 0., 0.],\n",
       "       [0., 0., 0., ..., 0., 0., 0.],\n",
       "       ...,\n",
       "       [0., 0., 0., ..., 0., 0., 0.],\n",
       "       [0., 0., 0., ..., 0., 0., 0.],\n",
       "       [0., 0., 0., ..., 0., 0., 0.]], dtype=float32)\n",
       "Coordinates:\n",
       "  * time       (time) datetime64[ns] 2000-02-01 2001-02-01 ... 2017-12-01\n",
       "    latitude   (z) float64 -39.36 -39.51 -39.65 -39.79 ... -41.86 -41.72 -41.58\n",
       "    longitude  (z) float64 -42.23 -42.4 -42.56 -42.73 ... 134.6 134.8 135.0\n",
       "    month      (time) int64 2 2 2 2 2 2 2 2 2 2 ... 12 12 12 12 12 12 12 12 12\n",
       "  * z          (z) MultiIndex\n",
       "  - xgrid      (z) float64 -3.938e+06 -3.938e+06 ... -3.938e+06 -3.938e+06\n",
       "  - ygrid      (z) float64 4.338e+06 4.312e+06 4.288e+06 ... 3.638e+06 3.612e+06\n",
       "Attributes:\n",
       "    valid_range:    [  0 100]\n",
       "    long_name:      Goddard Edited Climate Data Record of Passive Microwave M...\n",
       "    standard_name:  sea_ice_area_fraction\n",
       "    units:          1\n",
       "    flag_values:    [-5 -4 -3 -2 -1]\n",
       "    flag_meanings:  pole_hole lakes coastal land_mask missing_data\n",
       "    datum:          +ellps=urn:ogc:def:crs:EPSG::4326\n",
       "    grid_mapping:   projection\n",
       "    reference:      http://nsidc.org/data/nsidc-0051.html http://nsidc.org/da..."
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Stack the spatial dimensions together to create one spatial dimension.\n",
    "no_JAN = no_JAN.stack(z=('xgrid','ygrid'))\n",
    "no_JAN"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<xarray.DataArray (c: 11, time: 18, z: 104912)>\n",
       "array([[[0., 0., ..., 0., 0.],\n",
       "        [0., 0., ..., 0., 0.],\n",
       "        ...,\n",
       "        [0., 0., ..., 0., 0.],\n",
       "        [0., 0., ..., 0., 0.]],\n",
       "\n",
       "       [[0., 0., ..., 0., 0.],\n",
       "        [0., 0., ..., 0., 0.],\n",
       "        ...,\n",
       "        [0., 0., ..., 0., 0.],\n",
       "        [0., 0., ..., 0., 0.]],\n",
       "\n",
       "       ...,\n",
       "\n",
       "       [[0., 0., ..., 0., 0.],\n",
       "        [0., 0., ..., 0., 0.],\n",
       "        ...,\n",
       "        [0., 0., ..., 0., 0.],\n",
       "        [0., 0., ..., 0., 0.]],\n",
       "\n",
       "       [[0., 0., ..., 0., 0.],\n",
       "        [0., 0., ..., 0., 0.],\n",
       "        ...,\n",
       "        [0., 0., ..., 0., 0.],\n",
       "        [0., 0., ..., 0., 0.]]], dtype=float32)\n",
       "Coordinates:\n",
       "  * c        (c) int64 2 3 4 5 6 7 8 9 10 11 12\n",
       "  * time     (time) int64 2000 2001 2002 2003 2004 ... 2013 2014 2015 2016 2017\n",
       "  * z        (z) MultiIndex\n",
       "  - xgrid    (z) float64 -3.938e+06 -3.938e+06 ... -3.938e+06 -3.938e+06\n",
       "  - ygrid    (z) float64 4.338e+06 4.312e+06 4.288e+06 ... 3.638e+06 3.612e+06"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# We are going to loop through the array along the time axis 1 month at a time. Considering there are only 11 months, \n",
    "# this is not an onerous loop.\n",
    "ntime_permonth = no_JAN.shape[0]//11  # // is the integer division in Python 3.\n",
    "\n",
    "# Reshape the array to have dimensions of sizes (nmonths, nyears, nspace)\n",
    "# To get the unique values for month we use set() which is a collection of unique elements in Python. Xarray doesn't\n",
    "# like sets for coordinates so then we convert it to list.\n",
    "val_for_per = xr.DataArray(np.reshape(no_JAN.values,(11,ntime_permonth,-1)),\n",
    "                           dims=('c','time','z'),\n",
    "                           coords={'c':list(set(no_JAN.month.values)),\n",
    "                                   'time':no_JAN.time.dt.year[0:ntime_permonth],\n",
    "                                   'z':no_JAN.z})\n",
    "val_for_per"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12}"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# To see what set does. Not needed for the code.\n",
    "set(no_JAN.month.values)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create an array to save all the correlations and p-values for all months:\n",
    "stats = xr.DataArray(np.zeros((11,2,no_JAN.shape[-1])),dims=('month','stats','z'),\n",
    "                     coords={'stats':['r','p-value'],\n",
    "                             'month':list(range(2,13)),\n",
    "                             'z':no_JAN.z})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Dealing with month: 2\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/g/data3/hh5/public/apps/miniconda3/envs/analysis3-19.04/lib/python3.6/site-packages/scipy/stats/stats.py:3038: RuntimeWarning: invalid value encountered in float_scalars\n",
      "  r = r_num / r_den\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Dealing with month: 3\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/g/data3/hh5/public/apps/miniconda3/envs/analysis3-19.04/lib/python3.6/site-packages/scipy/stats/stats.py:3038: RuntimeWarning: invalid value encountered in float_scalars\n",
      "  r = r_num / r_den\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Dealing with month: 4\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/g/data3/hh5/public/apps/miniconda3/envs/analysis3-19.04/lib/python3.6/site-packages/scipy/stats/stats.py:3038: RuntimeWarning: invalid value encountered in float_scalars\n",
      "  r = r_num / r_den\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Dealing with month: 5\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/g/data3/hh5/public/apps/miniconda3/envs/analysis3-19.04/lib/python3.6/site-packages/scipy/stats/stats.py:3038: RuntimeWarning: invalid value encountered in float_scalars\n",
      "  r = r_num / r_den\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Dealing with month: 6\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/g/data3/hh5/public/apps/miniconda3/envs/analysis3-19.04/lib/python3.6/site-packages/scipy/stats/stats.py:3038: RuntimeWarning: invalid value encountered in float_scalars\n",
      "  r = r_num / r_den\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Dealing with month: 7\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/g/data3/hh5/public/apps/miniconda3/envs/analysis3-19.04/lib/python3.6/site-packages/scipy/stats/stats.py:3038: RuntimeWarning: invalid value encountered in float_scalars\n",
      "  r = r_num / r_den\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Dealing with month: 8\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/g/data3/hh5/public/apps/miniconda3/envs/analysis3-19.04/lib/python3.6/site-packages/scipy/stats/stats.py:3038: RuntimeWarning: invalid value encountered in float_scalars\n",
      "  r = r_num / r_den\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Dealing with month: 9\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/g/data3/hh5/public/apps/miniconda3/envs/analysis3-19.04/lib/python3.6/site-packages/scipy/stats/stats.py:3038: RuntimeWarning: invalid value encountered in float_scalars\n",
      "  r = r_num / r_den\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Dealing with month: 10\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/g/data3/hh5/public/apps/miniconda3/envs/analysis3-19.04/lib/python3.6/site-packages/scipy/stats/stats.py:3038: RuntimeWarning: invalid value encountered in float_scalars\n",
      "  r = r_num / r_den\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Dealing with month: 11\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/g/data3/hh5/public/apps/miniconda3/envs/analysis3-19.04/lib/python3.6/site-packages/scipy/stats/stats.py:3038: RuntimeWarning: invalid value encountered in float_scalars\n",
      "  r = r_num / r_den\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Dealing with month: 12\n",
      "CPU times: user 8min 35s, sys: 5.64 s, total: 8min 41s\n",
      "Wall time: 8min 33s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "# Since we correlate the months mn_start and mn_start+1 we need to stop the loop at the index before last.\n",
    "for mn_start in val_for_per.c:  \n",
    "    print(\"Dealing with month:\",mn_start.values)\n",
    "    \n",
    "    if mn_start == 12:\n",
    "        break\n",
    "    # Need to drop the 'c' coordinate after selection else it conflicts with the rest\n",
    "    tmp = val_for_per.sel(c=slice(mn_start,mn_start+1)).drop('c')\n",
    "    # Now the data is in the right shape to calculate the stats:\n",
    "    out = tmp.groupby('z').reduce(corr2,dim='time')\n",
    "    out.rename(c='stats')\n",
    "    stats[mn_start-2,...] = out  # -2 as the index starts at 0 but the months start at 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.QuadMesh at 0x7f9144987048>"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaUAAAEWCAYAAADGjIh1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzsvXmcHEd58P99unt6ZnZmd/bSrlbn6vIh21jYsmwDxsYGH1wm4ZKdBANOHPPiJCQkL/CGBAI4L/DjxcQhWDjBYMAgGxODAR8YHzgYC1vyrXt1r7T3NTuzc/R01++P7lmP1itppV1pdlb1/XzqMz3VVdXVvbP91FP11POIUgqNRqPRaKYDRrk7oNFoNBpNES2UNBqNRjNt0EJJo9FoNNMGLZQ0Go1GM23QQkmj0Wg00wYtlDQajUYzbdBCSaOZYYjIh0Xkd+Xuh0ZzLGihpCkrIvJ5EfnhUZS/RETaj2efKgkRaRURJSLWCb5uk4j8WEQOiMiQiDwlIuefyD5oZiZaKGk0R8mJFgDTlDjwLHAuUA/cCfxKROJl7ZWm4tFCSXNCEJFPich+ERkWka0icpmIXAn8H+CDIpISkReDsh8Rkc1B2Z0i8pdBfgx4EJgTlE+JyBwRWSUi60UkKSJdIvL1Ke77hwNN4BYR6Qc+P4E6T4jIl0Tk90E/fyEiDSJyV9DPZ0WktaT8G4K8oeDzDWPa+mLQh2ER+bWINAannww+B4PrXFhS72siMiAiu0Tkqql4FkWUUjuVUl9XSnUopVyl1O2ADZw6ldfRnIQopXTS6bgm/BfVPmBO8L0VWBIcfx744Zjy7wCWAAJcDIwA5wTnLgHax5R/Gviz4DgOXHCIfiwABg+Trj1EvQ8DBeCvAAuITuCenwDagvtIAJuAbcBbgza+D3w3KFsPDAB/Fpy7JvjeUNLWDuAUIBp8/3LJs1SANaa/DvAXgAl8DDgAyCH6+svDPJNfTvBvvALIAoly/950quykpyE0JwIXCAPLRaRHKbX7cIWVUr8q+fpbEfk1cBHw3CGqOMBSEWlUSvUC6w7R7l6g9mg7H3BAKfXvwXFhgnW+q5TaASAiDwLLlVK/Cb7/BPhiUO4dwHal1A+C7z8Wkb8G3gV8r6StbUHde4B3H+Hae5RS/xmUvxP4FtAMdI4tqJR65wTvZ1xEpAb4AfAvSqmhybSl0ejpO81xRynVBnwCXyvqFpG1IjLnUOVF5CoRWSci/SIyCLwdaDxUeeB6fC1iSzD1NamX7CHYdwx1ukqOM+N8L66/zAH2jKm7B5hb8r1UmIyU1D0Uo+WVUiPB4ZSv94hIFPgFsE4p9X+nun3NyYcWSpoTglLqR0qpNwEL8aebvlI8VVpORMLAT4GvAc1KqVrgAfypvNeUD9rerpS6BmgK2r03WH86CBFZULIWNV76k8PdwtHd8VFxAP+5lLIA2D+BupPul4g8eJhn8uBh6oWBnwX9/MvJ9kOjAS2UNCcAETlVRC4NXmJZfC3BDU53Aa0iUvwt2vhTfT1AIVigv7ykuS6gQUQSJe3/qYjMUkp5+OsglLQ/ilJqr1Iqfph011HcU9EUu3WidQ7DA8ApInKtiFgi8kFgOf5az5HoATxg8bFeXCl11WGeybgGEiISAu7F/1t+KHj2Gs2k0UJJcyIIA18GevGnlZrwre4AfhJ89onIc0qpYeCvgXvwF/uvBe4vNqSU2gL8GNgpIoPBNOCVwEYRSQH/BqxWSmWP8z3Nx59im4g2c1iUUn3AO4FPAn3A/wbeGayPHanuCHAz8FTwPC6YbH8myBvw+3w5r1r+pUTkohN0fc0MRZTSQf40mqNFRD4L9Cilvl3uvmg0MwktlDQajUYzbdDTdxqNRnMSIyJ3iEi3iLxyiPMiIreKSJuIvCQi55ScuzLYDN8mIp+eiv5ooaTRaDQnN9/DX5c9FFcBy4J0A3AbgIiYwH8E55cD14jI8sl2RgsljUajOYlRSj0J9B+myNXA95XPOqBWRFqAVUCb8l1O5YG1QdlJoT06TAGNjY2qtbW13N3QaDQVwIYNG3qVUrMm08Z8iaosE7PC7yW/EX8rRpHble+rcKLM5eDN4+1B3nj5k/YUr4XSFNDa2sr69evL3Q2NRlMBiMhY7x1HTRaP99IyobLfZk9WKbVyEpeTcfLUYfInhRZKGo1GU2EIYI4nEsZj8gbW7fj78orMw/dCYh8if1LoNSWNRqOpMASwDZlQmgLuBz4UWOFdAAwppTrw42ktE5FFImIDqynZ6H6saE1Jo9FoKgxfU5oSgYOI/Bg/JEyj+FGdPweEAJRSa/DdYL0dPxTLCPCR4FxBRG4CHsYPkXKHUmrjZPujhZJGo9FUGnIU03dHIHBmfLjzCvj4Ic49gC+0pgwtlDQajabCmEpNabqhhZJGo9FUGEdl6FBhaKGk0Wg0FYdoTUmj0Wg00wMBQlooaTQajWY6IFNo6DDd0EJJo9FoKhA9fafRaDSaaYE2dNBoNBrNtEGbhGs0Go1m2iDCVLkQmnZooaTRaDQViJ6+02g0Gs20QK8paTQajWbaIHrzrEaj0WimE1pT0mg0Gs20wN88OzOlkhZKGo1GU2EUg/zNRLRQ0mg0mgpDGzpoNBqNZloxU6fvjHJdWEQiIvKMiLwoIhtF5F+C/HoReUREtgefdSV1PiMibSKyVUSuKMk/V0ReDs7dKuL/tUQkLCJ3B/l/EJHWkjrXBdfYLiLXleQvCspuD+raJ+J5aDQazUQRAUNkQqnSKJtQAnLApUqps4EVwJUicgHwaeBRpdQy4NHgOyKyHFgNnAFcCXxLRMygrduAG4BlQboyyL8eGFBKLQVuAb4StFWPH4f+fGAV8LkS4fcV4Jbg+gNBGxqNRjONEMScWDpiSyJXBgP9NhH59Djn/0FEXgjSKyLiBu9QRGR3oBC8ICLrp+LOyiaUlE8q+BoKkgKuBu4M8u8E3hMcXw2sVUrllFK7gDZglYi0ADVKqaeDWPLfH1On2Na9wGWBFnUF8IhSql8pNQA8gi8UBbg0KDv2+hqNRjMtEAHTNieUDt+OmMB/AFcBy4FrAgVgFKXU/6eUWqGUWgF8BvitUqq/pMhbgvMrp+LeyqkpISKmiLwAdOMLiT8AzUqpDoDgsykoPhfYV1K9PcibGxyPzT+ojlKqAAwBDYdpqwEYDMqObUuj0WimB8JUaUqrgDal1E6lVB5Yiz+YPxTXAD+eorsYl7IKJaWUG0jfefhaz5mHKT7e01WHyT+WOodr6+DOiNwgIutFZH1PT894RTQVQudQms6hNCM/+SovvPcKHnvdKu5tPuM15fqGR8rQO41mHEQwzImlI3CoAfo4l5Qq/KWRn5ZkK+DXIrJBRG6YxB2NUlahVEQpNQg8gX/DXcGUHMFnd1CsHZhfUm0ecCDInzdO/kF1RMQCEkD/YdrqBWqDsmPbGtvn25VSK5VSK2fNmnWUd6wpN99MnMo3E6dyd9PyQ5bJ/89advUOs6t3mI0dSUALJs30QQxjQgloLA6gg1QqPCY8EAfeBTw1ZurujUqpc/Cn/z4uIm+e7H2V0/pulojUBsdR4K3AFuB+oGgNdx3w8+D4fmB1YFG3CN+g4Zlgim9YRC4I1oQ+NKZOsa33AY8F604PA5eLSF1g4HA58HBw7vGg7Njra2YA980+g/tbDlbIX77gzXRc90ds/v4jZAayADQvrn1NXcfz/1fbbvrAuJqURnOiEOFoNKXe4gA6SLeXNHWoAfp4rGbM1J1S6kDw2Q3chz8dOCnKuU+pBbgzWGgzgHuUUr8UkaeBe0TkemAv8H4ApdRGEbkH2AQUgI8rpdygrY8B3wOiwINBAvgO8AMRacPXkFYHbfWLyBeBZ4NyXyiR/p8C1orIl4DngzY0M4D/ufCNALhKMeR4APTnXaIjBWZ1j1Bvm2RcP982hMev/yo1825j4aXLqf/brxMyBOdHXzpk+zdKKzcnNwHQUF3Fzr+9lpHuIc6861fH98Y0JyUTsaybAM8Cy4KB/n78d+S1r7mWSAK4GPjTkrwYYCilhoPjy4EvTLZD4isHmsmwcuVKtX79lFhDao4TN0orf3KBP1X+4HOd1IVMoqYwN+qPy0wR8p7CFMi4iqgpRCIWYhrEmqo480MXUb3qInLbnkcsG1XIY9Y1EbnyhtH2gVGhZN/3VXqe38ZI9xDJ9mFizTGq59ax+JYfnehb10wzRGTDZC3VTovF1e1nvm5CZS9+5unDXk9E3g58AzCBO5RSN4vIjQBKqTVBmQ8DVyqlVpfUW4yvHYGv4PxIKXXzMdzOQWiPDpoZTVFYAKQ6UoRiNo22Sd5T1NsWpgjhan9/dG1dhPadgwBU1YQBiDfFcLIFOp/xhU2opgo3lcJumcfA8y9iv/HgNaZftJ7LxR9ZiRmycNJZ2tcdoGZeNQDmZ2/noSWvZ8Vzv+PztWdwXl2EqGlwbc+m4/wUNDMOAZki33dKqQeAB8bkrRnz/Xv4s1GleTuBs6ekEyVooaSZsXwzcero8Xl1ERILE/Ru6ycRMhhyPDYP54gYBs5gFgM4sL2fC+qjnHp2E3PPX0y4Lg7A8N4uhvYMsPvh5wjXhMklc4RiIWaffxq1+zewJXE2n/rkRdz37T+weTjH5lufotG2WLmsjuqWOEvf9XoAXnrzW7AiFr87ZRVzIq/+63V/7W9o+vt/O6HPRlPpCIY5LezUphwtlDQzkj19/r7sM2vCXPbHp5JPO3Rs6MS0DNws5D1FKHDBYgCOUtTbJnHLoOW8Vtp/v4PvP7mXbzz8T4Rrd+Dmt2CYBqGaGNULLMK1cbJ9ScQwWLbzIfYDK5c3sn5TL67y16S69w/TvCBB14bthGI2ruMykswF04RCT87lLW9oAXyz9IKraPjdHUTf/ddlemqaikGmbE1p2qGFkmZG0d7vCyPLEP7Xcz+g9767aPvF8/Ru62fIccl7yjduMA2W14QwLYOtg1l2jzg0hy3ynmLLvS/xs50DAHziii9y68AzLH79OgDEMJFwBJXL4uzbxpYvfZm9T7Xj5l1SuQLnnNqAl3cJ14TxAqOJ5L5hAAqZAv15l7wHi2IhmppizL/4NOKL5sNDtxI+YxV5wN30BObyS074s9NUDqKFkkZTPkrXhcby2b5Xxs0veIrOH3+X5O4O9m3uJWoaNMZtBtIOUVPhKkW42mYkmQN8zWZBlUXLvGrufqGLNWr36HWHjDh1c0/FrW7GyCYxkp24w4NgmAzuGRq9Zt5TpLvSBwkk1/FIdaRwci6pgkdPziVuGaQKHk34U4NWxCa5u4P0T59k2Xv9bR6Z+2/VGpPmsOjpO41mGuMBDQ9/g32PrsfLu1S11PPod56hJ+fiKoiaQjRX4NlgHxLAtu39rFG7AV/wPdWXYU3vZt4ypu1/SiznC4MbwYN6J0Ohex9D635H/+Y9WFGLqsYoe3YM8FRfBsgAcE5thFlhk2jMpms4P2rZt3vE4RPXr6CQcQCwor5A6tnURbprhOcf/T4Xf2QTnuMgj6+jekEzTjrLnM/edtyfoaZyEBHMkBZKGs0J4XCaUSn1tskDS1by5veeRi6Zoy/vb1sLxUIcWLeL/rxHxvW3PNiGsHhhgmcHsqNaUFEgAQcdl+YV+xK3TcxCltQv76TzmU1YEZvh9iEKmQLhmjBxyzioznODvvB71zzf9NwUmDunmmcHsiz+2u0Ye14g3/YSwzv2kO7oI1oXId01ginw2++uxzaEC/50xWhftDGE5iAERGtKGs2JJxEyaA5bzAqbtDTHMUwhFAvRsKweuyZK39YeBnYO4OU9RvpGcPMekboIiXnVRE0h40LU9KfLitNyML4QGo+Dy0UI3/Cv1AZOWkqF5xq1mz8at05J2b4Mt/z871CGhTd3OeGQjZtOkU+OkEvmCdfYJPpNTMvAjofoenE/APWnt5JPpif8zDQnBxPwa1eRaKGkmbaYAjWWyam1ES756gdws3n2/Ho9uaEcwx0p7GQOz/UwbZP9G3upToQxbYPqljjRxirmJ8JEU3kyrmLzsDNhQTRRjqW90NwlvDhkcEZDHWZ1A0M79uM6BayIRSFTIBQ2idRFsCIWmd4MPa90MOv1p1D7gU9Q2L8ZGe5FZdI4ezYTfc/fTun9aCoImVispEpECyXNtCRqCs1hi4teP5uG0xp4cc2v6draR1VNmKqGKIZpYMdCGKZBbihHxDQIRSwaTmug/pQ5JJbMpW9LHz9/at+UC6NjoTit95/OGXx4VgTJj+CFqkgsmUvn06/4GmA8RLgmjGmbhBNhwokwbt6l/fHnmG+uIbz0dbiZNG5fB2KaDNz2acyITc1HJu3ZRVNhiJ6+02hODKNrOJbBm1e20LxiDt0vd+DmXc79yCoyfUMM7R3EMIRv/3zbQXXf6np4rsIMGeSH0/xgmgikUl7fUoPhZLEG9qFMi9ib3on75At4riI7kCWXzBGpixCKWMRbahDTYLh9kB0/f5pQ5Flcx6X1qlWowLpPc5IiaEMHjeZ48mDr2fw8MK8+ry7CwpZqTn3/+XQ9s4Vn1neyNBHGsC2ygxnCNWFu+8nmgwTOjdLKb7rT0J2GTX58q+kmkIra0sr2B2G4F7Fs3Lp5xOfOouGMRWz/2QZCMZvahQnsmigAw+2DDLUPU8gUMG0D0zbJ/OhJckM5Wi87jUhDDdlffwereT7eyDDe8CCZHdvZ+cAG8uk8Lee1Eq6N4+ULbPvZCwetq2kqF9EeHTSa40cxRtHVCxMseusi5rzxdYhp0L1hCzt/sxtXKfYO58ms+T21jVXMXtH0mhdrpbxoSy36vvHwPyGde5l9w99hZIZI7e/lwLrd5NN5xBTcwJqwkCngOi6GKUTqImQHsrh5l033PI+Xd2k5dxOnfuTdqHSSka5eute/qkHufWI7rZeeihmxR/M+H11KZ7ZwUJ80FYbePKvRTC1jg+Wdcc05RBpqiM+fjYQjDLftZqQ7Sbrgbzr1MRnsHeHeH21kzV0nvs9TRVEwfeKKL/Jvj3wODAsv1kDLpW+g87m9iGmQTzuYIQMxBdM2MEzByRZw8y5VDYEW1ZFiJFOg84VuBj//Q5rPmoUVsTFsE9M2cLod3LzLhu8+Q2fWZchxmR2xDhJIAC9f83bm3H4v4Ifc0FQAek1Jozk+OJ5ixPFYeOMncNu34vZ1MrJ3H3sfe4XhjhSmvBpeIlXwWNefmREj+6JgsmYvQCkPGewg1Ho6Kz7xXnqfeRGA/m0HRgWUk3I4kHGw+zK8/rJWGs9YwFzbomt9G7lknsxAlo4NneSSOZycOxozqhg3KlXwcJQKvKOb1IVM0sG61NMP7eS95XoQmmNEilFlZxxaKGlOCMUpqy/9m//6a7v6HwE4u7mKpt6NjKx7CMOO8Hcf+Ba39q8jsfEmPE8RrglT3ZGivXcEdwbG/ur4yVrqTluI8jwil16DfXYdc86+mN67/4vbfrIZ8Kc1k4HfvucGs6xwPUY6++jd2sOs5c1EGwrkkp1k0w7pTAHbEOyQSVz5saGGHBfbEEzlBzjMe4oDWYe8p7ANwVWKro9/kMZ/X0t7fwrjW/+gPUhMc/zIszNTKOkgf1OADvJ3eG6rPY0Xh3wPB9eumgOAGTKJz4n7Qmd+M6n9PWQHM6Mv4q/98KN4+QJdz2xiaM8A//XgjhmhIY3lRmnl5n9/P6FYFOV5fOrPf/gaA44ih9yUO4Y3N1ZhG37QQlcp+vPe6LEvhCDtekQMwTaEZMEjJL5wilsGqz91mRZKx5GpCPJ3VnO9+tnqyydUdumtd0/6eicSrSlpjhu31Z42enx2IsKpC2qI1EXYur6DuXOqyQxkmbuqheoFzbT/bhuZEr90A5v3cPPNj45+n4kCqcg//tVP+OodHxr33JHue+z5G6WVRMQi77hkXA9ThEWxEOFqO/CGYeOk8+ztSpP3GI20WxRYGVfx3X/9DTfW/gORhgSxaz47NTepmWJm7vSd1pSmAK0pvZYbpZWzE5HR77PCJrvSDrbhrxGteP8ZZAfSDHekuPPxPQD80bJ67itxkqqZHEUt6hPXr+Ab33kBgP/7rdUUsjn+6e/uO8gZbZHTq/2Iu2e1Jjj/M+8mP5ii7mNfPnGdPgmYEk1pdoO6/0/fPqGyi//fDytKUyqbqBWR+SLyuIhsFpGNIvI3QX69iDwiItuDz7qSOp8RkTYR2SoiV5TknysiLwfnbhXxo7eJSFhE7g7y/yAirSV1rguusV1ErivJXxSU3R7UfdWWVnNYnn/P5dzdtJy7m5YD8OJQls3DORZU2ww5vnlz1PSnlVIdQzhphzsf3zP6ctQCaWpZo3azRu3mG995YfS5fuZ/rT1IIJWWA9g8nKPeNnhyez8b73wMgFf+5B0nuOeaIyEIYhoTSpVGOXtcAD6plDoduAD4uIgsBz4NPKqUWgY8GnwnOLcaOAO4EviWiJhBW7cBNwDLgnRlkH89MKCUWgrcAnwlaKse+BxwPrAK+FyJ8PsKcEtw/YGgDc0R2PE3qwE45Q3zRsNDRE0hbhlk8sFCu4AbLLav+e+tvP5nvz7IQaoWSMeHsc/4UM+5mP9UX4a6kMnvnmpn7Wd+hpN2XmPCrykzAoZhTCgdsSmRK4OBfpuIfHqc85eIyJCIvBCkf55o3WOhbEJJKdWhlHouOB4GNgNzgauBO4NidwLvCY6vBtYqpXJKqV1AG7BKRFqAGqXU08qfi/z+mDrFtu4FLgu0qCuAR5RS/UqpAeAR4Mrg3KVB2bHX1xyCh085B4CahfWAL4zsYBEdIBa1gnULv3wxnIRm+hK3BFOEIcflrl+1lbs7mnGYCk0pGNj/B3AVsBy4JlAAxvI/SqkVQfrCUdY9KqaFoUMwrfZ64A9As1KqA3zBJSJNQbG5wLqSau1BnhMcj80v1tkXtFUQkSGgoTR/TJ0GYFApVRinrbF9vgFfO2PBggVHdb8zidtqT2NxUxVtD25jS2c6sPLyfdcVvXxHG6t4bs8QZ9aEcZUiWoFTCicDpWtLs8IWGddfA4xbBne0rOCjHf66lN5gW35EBCM0Ja/vVUCbUmpn0O5a/MH8puNc95CU/e0gInHgp8AnlFLJwxUdJ08dJv9Y6hyurYMzlbpdKbVSKbVy1qxZ4xWZ8XQOpTEFXjqQwsu7LIqFAFgSC7EkFuL06jBL4yHWbuhgjdrNK8kcryRztDRE9VTdNObts+MsvGg+rz+riXrbZH7UIm4ZbLv6qnJ3TVNEjmpNqVFE1pekG0paOtQAfSwXisiLIvKgiJxxlHWPirJqSiISwhdIdyml/jvI7hKRlkBLagG6g/x2YH5J9XnAgSB/3jj5pXXaRcQCEkB/kH/JmDpPAL1ArYhYgbZU2pamhOIaw7t3b+CxpecRqYtw6iVLeMOiFsxQiPW3PkI+5eAWvIPWNG6UVq7Y9lz5Oq45IhnXw8u7iCnMbaqipzdD3vNIdaR48Y2XcM/L3XpQUW6Ozs1Q72Gs7yYyEH8OWKiUSonI24Gf4a/dT3gQfzSU0/pOgO8Am5VSXy85dT9QtIa7Dvh5Sf7qwKJuEf5DeSaY6hsWkQuCNj80pk6xrfcBjwXrTg8Dl4tIXWDgcDnwcHDu8aDs2OtrAjqH0jie/9v77bLzGHI87JjNV7/xez77Nz8l/kd/wU+39PGOfS/x7o5XDqqrX2bTl+Lf5vGeEb79820M7BykkCmMTsVuPZBiU+DJXVN+DNOYUDoChxrsj6KUSiqlUsHxA0BIRBonUvdYKKem9Ebgz4CXReSFIO//AF8G7hGR64G9wPsBlFIbReQe/PnKAvBxpZQb1PsY8D0gCjwYJPCF3g9EpA1fQ1odtNUvIl8Eng3KfUEp1R8cfwpYKyJfAp4P2tAElFphzdq3jv3BSwvgQ29eQNfWPv6q5TItfCqcC+qjZPLuqMWkH1red1GkKT8iU7Z59llgWTDQ34//jrx2zLVmA11KKSUiq/CVmT5g8Eh1j4WyCSWl1O8YX/0DuOwQdW4Gbh4nfz1w5jj5WQKhNs65O4A7xsnfib+ApxmH9L9+HIDE0rm88LM/kPcU8xP+hkvlKhoWJljTubGcXdRMguIU67r+DAAfXNFMLpkn5nrMjlps2Hu4ZV/NCUMEw5786zswALsJf/bIBO4IFIAbg/Nr8GeOPiYiBSADrA5mlcatO9k+TQvrO830p2iZdXNyE7l/+yTRU86kv+1BEiETwzbJDGRx0nmijdoyq9IpjfkUrgmTTzso14/rVG+b7B5xyto/jc9UuRkKpuQeGJO3puT4m8A3J1p3spTd+k5TWey+9mrCtXHa/vMuPFcxN2rR3Z9h+/5hfr5niEvWP13uLmqmgOL06/ef3EvtwgROzqVrOE8iZFBvm4evrDnuiAiGaU4oVRpaKGkmxBq1m6+v/UtO/9DbiDbPItWdJt4UA9BrDTOc6pY4Gdf3NJ4qeLhKse+z2tFJudFuhjQaAM9l3yPrCEUscskc8VlVLJ5TzebhnDZumGEU/563/uDlUQ/iPTn38JU0JwaZuUJJrylpJsTQdz6Lk/YXv6sXNI+amvZu7eOnW/q0QJrhuIogFpPQHLYY2tXJzdKq/+5lY+aGrtBCSXNE3E1PABCKRUm199Bw4SqamueT3/YC2dt+BVv6yttBzXHHj1ILec/3aWhFQ+Xu0kmNGFNjfTcdmZmiVjNl5J+6B4D4hW9FDAM3m6f9/ocZeOxBBjbtYOGly/VoeQZT/Nv6YdQZXVd64O5NmIfa0KE5IYhhTChVGjNT1GqmhPzTPwVAZdJ42TT54TQAPZu6eHrty2RcxYrz5zCnnJ3UnBAyriLjeqPxsPZlHFzlbxXQg5IyIIIYlWdZNxG0UNK8hjs3+D4WV+f8uEi7fvhd+tv6aFzeTO2y+dQvHeKsughr/nsrNzz0WDm7qjnOFPcrZVzFgOPSaFskQgZL4zFSBY/+vMvL17yds348pVtVNBNBCyXNTKetZ5iQ8eqcjJfsI717D71be3DSDkYoRH54BM8pEKmLlrGnmhNBaSiLRMigN19cWzKoa03QaJssCfYsdQ+laUrEytPRkxKBCpyamwhaKGnIPnS7f3DuNTie4k8b+sg9+2v2/3Y92YE0VY1VNF68gIbzVmCfdh7/suRfcK14AAAgAElEQVR9h29QU/GUCqR628RVirMTEeJxm3CNjRW1MG0TJ+VQd8pcGga2Q2JF2fp70iEgFbgxdiLMTFGrmTCdQ+nR4wUyxNzqEGqgk+Fd+8j0pcgl84RrwoTr4ljNCyjUzeOb+37Fvyb9OF6lLy/NzOPaVXO47rqzmR2zAegczDJreSN1i+sI14QRUxje20V79VK2dCV12PQThQhY9sRShaE1pZOYvuERQoYgkRgSCqFeeAQ11EfXK9vpXL+TcI1N8zmtGCGLcG01zv4ddN1zD4NtnbScfwr/+I+XcfPNj5b7NjTHgVLjhbEDj3kvdPOWh+6gsHU93U/8jtT+HiLf+Qxzr3gXQ4vfdEL7ebIiep+SZiZROpod/tLH6B0YpuXCMxncto9wXZyBbe1kB7LY8RBiGBTSWXpf2kGkoYeWt11CuHY9DW97B7u/8xon65oZyFgBde+mHu5d8C4Avr72L3GdAk46i8qOjMbZ0hxnBG3ooJl5OJ5i890vkFiYoHfT41hRi5HeEQY7UkQiFmIKsdnDAOx5YgeFbIF52/bR+Lol5Lc9j5POlvkONCeaUg/iAJ3rXsauiWFXV+GldViLE4dooaSpfMbO9/922XnELJNCtoDnKjzXI5fMA+DmPZSrSHcO4uZdaubVkB3IsuXBHSxK5kgs6sWKVN58tWbqWH1uC3ZNDC9foG/jLpx0lrmXfqjc3TppmKnTdzPzrjRHZOs7r2R/poDjeQCYIQPX8ShkCtghEytq4bkefdsH6N8+QLI9ief6ZU+97h3MfvN55JKZct6CpkwUp/PsWAgjZNHzSjuDe4bY//R2Hmw9u7ydO1kQQxs6aCqbG6WVz/7z2+h5cQ992wfo60mTCBlU1YTJpxwMU8glczieRzwRwQyZpLtHMG2Dzp4R9o74Yc8TIYMXvvFT7nx8D2vUbtbcVe4705SL7z+5l5vmVWPHQxSyBXY910nT3Gr+8La3cP4jj5e7ezMbbRJ+fBCRO0SkW0ReKcmrF5FHRGR78FlXcu4zItImIltF5IqS/HNF5OXg3K0iIkF+WETuDvL/ICKtJXWuC66xXUSuK8lfFJTdHtStvKHGIfjSFx5hx9P7Ge4bGY1/5KQccsncaJlIxMKKWFhRf03JC2IlJQsu+zIOu0ccdr3cU65b0EwzvvmjjeSGcuSSOaKmQT7tR6V9+JRzytyzmU6weXYiqcIod4+/B1w5Ju/TwKNKqWXAo8F3RGQ5sBo4I6jzLREpDhVuA24AlgWp2Ob1wIBSailwC/CVoK164HPA+cAq4HMlwu8rwC3B9QeCNiqeM2vCnFkTJuN6vDiUY3+mgCmCpxRuwWNXe5Ku4TwHkrnRGCyhiIWTctiXKRA1DU6J2yyJhXiyd0T7O9MA8MEVzQx3pAjXhInPqsIwBTsWomZedbm7NrMpWt9NJB2pKZErg4F+m4h8epzzfyIiLwXp9yJydsm53YFC8IKIrJ+KWyurUFJKPQn0j8m+GrgzOL4TeE9J/lqlVE4ptQtoA1aJSAtQo5R6WimlgO+PqVNs617gskCLugJ4RCnVr5QaAB4BrgzOXRqUHXv9iuXmqmUMOR5DjkfGVUQMwVW+y5hwtY1d5Ych8KOKgmEKYgiu49KVdch7CgPf1Ux/3tMCSTNKIVvAtE0G9wwRrgkTqY3guYpcMs8v55xV7u7NYHyHrBNJh23FH9j/B3AVsBy4JlAAStkFXKyUeh3wReD2MeffopRaoZRaORV3Nh3XlJqVUh0ASqkOEWkK8ucC60rKtQd5TnA8Nr9YZ1/QVkFEhoCG0vwxdRqAQaVUYZy2Kp4FVRa2IaQKvsFCvW0SrYvg5l2aMwUG8gVsQzBDr/6Q45ZBa1UI2/Dj6PzVC3oRSfMqkdoIdtxGDCHVnaa6JU6yPVmREU8rjqmZmlsFtCmldgKIyFr8wfymYgGl1O9Lyq8D5k3FhQ9FJf1yxoveog6Tfyx1DtfWwZ0RuUFE1ovI+p6e6bXG0jc8cpD5tykcFPumKGAAcskcI70Zko6Lq8AUIRQPEYqHiNRGiFnmaHlXQe6V32vXQppRlKfIDGQpZAsoV5EdyBKuCQPQlXX0b+V4IQZi2RNKQGPxXRWkG0paOtQA/VBcDzxY8l0BvxaRDWPaPWamo6bUJSItgZbUAnQH+e3A/JJy84ADQf68cfJL67SLiAUk8KcL24FLxtR5AugFakXECrSl0rYOQil1O4Eau3Llymmxjb34Arg5uemg/He8exkAoZjN5l9uZ2ljFQBDwzm2HkgBfrhr8IWXm3cxTN9E3FOK/ryLKULG9fjKB289Ifeimd4UN9H+6JkDXL0wQTrja9lDwzky+4dZdnojzw7ozdXHDeFoNKXew0ytHc1A/C34QqnUl9QblVIHghmtR0RkS7Asc8xMR03pfqBoDXcd8POS/NWBRd0ifIOGZ4KpvmERuSBYE/rQmDrFtt4HPBasOz0MXC4idYGBw+XAw8G5x4OyY69fsRQyBX//UXUVp121hPpldSQW1pCoDo8GbStGFO3JuWR6M6S70mQHsrhKEbcMFlRZzAqb1IVmphmq5ugo1YCcnIttCKGw/9tIFTw6to9dKtZMJYIgpjmhdAQONdg/+HoirwP+C7haKdVXzFdKHQg+u4H78KcDJ0VZNSUR+TG+xtIoIu34FnFfBu4RkeuBvcD7AZRSG0XkHvy5zgLwcaWUGzT1MXxLvii+allUL78D/EBE2vA1pNVBW/0i8kXg2aDcF5RSxf+iTwFrReRLwPNBGxXFP9b4Ico7h9IoBbO/dx9hU6jtfoXcxnXkOjtp+9kfqG6JU5/OB1qQrxGlCh77DwyTCJlkXA9ThKXza0gsTDDf9dixvoPNw+W+Q025KXU39EBnivctn0VVYxSrO409kCXvKU6vDrMv4+jotMeDqfN99yywLBjo78d/R1570KVEFgD/DfyZUmpbSX4MMJRSw8Hx5cAXJtuhsgolpdQ1hzh12SHK3wzcPE7+euDMcfKzBEJtnHN3AK/xKBos+E1a2peD8f7xJVDO857CjTXg9PUysG0fseYYnqcIGQZRU2EKRE1r1ALPVb4GX1x7Mm2DQsajZVEt9OjwBJqDBZNp+3uUYk0x3LxHXV2EeFeaetvgqT7t+WPqmRrfd4EB2E34s0cmcEegANwYnF8D/DO+Edi3gi2ghWA6sBm4L8izgB8ppR6abJ+m45qS5jigFHixBvo27qL7xXbcvIubd4nURbCyBdy8r3T2DGQBj2iJZUQumaOQKfiB3fT0naaEomBy8x5W1MCKWlQ1RhFDqF2Y8IXS43voGx6hobqq3N2dOYggVmhKmlJKPQA8MCZvTcnxnwN/Pk69ncCU+5XSQukkQQTM4S5S+/sY2DlIfsTBVYqM62tJrmLUVHzI8aftTIH+vIedK8CmXuqX1XPqB1bBU/v0lIxmlFKN6c+vWuLvV6qLUD23jlAsylf/7GIc0IJpqpHpaBIweWbmXWkACBlCyHhV4zGyw5x604exYyGqZ8doWJhgVmOU/ZkC+zMF+vMuyYKLo3yjh7ynyAROWJMjDiO9I7jZPJ/+3xcDOuqsZnwidREAlOdh11RhlniT15FppwrxhdJEUoVReT3WTJiG6ioaqqtGhdPArOUMnvkOlv/P4yhXUcgUyKccTBHSrkdXrjDqE89VvkWebfhGEEOOR/f+Yf5wy+Pse7KNT33yojLfnWY6UdSa/+vBHVTPraOQKTDSncTNF0ju7qCu66XRslowTQ1KjAmlSqPyeqw5akqnTJoSMW6ZdRYbdg6yYecgqVyBuCWERGgOW8EmWTW6b8k2BFP8jbMZV2FaBtG6CMr116C0tqQpUhRMX/3G75l15jzEEAa3tWNXV5F9/knckqi0WjBNEkFrSprKpqg1FenNFzAFmhckcJXvcmhu1KK1KkRI5KApvHrbD1kRtwzcgsfmZw7wwB0b+It3LivjHWmmI0XB9K9ffpxwTZThjhS7HniOdHsHDcmdZAuKbEHheEoLpkkh/kLxRFKFoYXSSU4o4m+KXVBtY4qQ9149ZwY/6JhlEjWFeNwmXG2TcX1N6j9/uV0bO2heQ/E38f/WrB8NHpna34uz+RmaR/bSnD1AIVDFtWA6NhSgTGtCqdKovB5rJs1X8jv5p8gSBhyXcI3NWe87g/Z1+8jsS+IqNernrhjUD8AOmYRrbCJ1ERKDWRIh9CZazWH52PtPp397P6Zt0vPKfsR8kuptW8j0DOLlC8QuWkX6ouuO3JDmtYhU5NTcRNBC6STli9kdPPWmi6ielyA2dxbxpj4G2pOj2hEwOmUHYEUtInUR7JjNrNlxrKilN9FqxqV0ndGO+ZZ30boIg22dDO3qxsu75NMOzkiGWi2Ujh0tlDQzjWW/egjzB5/HioZJdY+MuhlKux5zIiHqbRNDBNM2qF2YoKoxSj7lYEUtQrGJb9w7FmMIPS1Y2dz4x6dimL7H+aqGKMpTjPRmyA5kGen1PTwk25Oc86Uyd7Rimbma0sy8K80R6RxK43iQ/ZPPE/ngP3D+P1/Dhe9cxor5fsTQfRmHl4dydGUdPFfh5l3MkMnslYuJ1EXIBC+WIwmc0vPvPa0BgI+8dREfvWLxaP5fvHMZH1zRPJW3pykza/57Kwc2dJI6kBoNb5HuTrN3/zCvDGVZP5Dh5fZhvZl2EsxUk3CtKZ3EFOOwmYPtGM0LiDRUM9SfDQwefIuHVEGRyhWoSubwXA/DtojWRRjYOci75tXwi/Yk8KrwKd3dX+TvP76KVMcQPZt6+cBZTaS70qMbLN/ZUs3AzkHseIjV57YQDfJL29NUFqW/AcM2GenN4OZdckN5+vMuadcjVfDozBYO247mCFSgwJkIWiid5FQ/cAudW/aQ2u8HKjztra20/WwbGRfSroerfJPwkd4MyfZhYi1JEq3NJPcNM9I3wntPazhICI2nOaW7h8mnHRILE+SGcpi2SbgmPKodFbIFqlvihGvChGJhPNdj9blZ1m7o0O6MKhgzZAZB/2zEFBzPI1nwGHK80U3ammNEpsYh63TkkEJJRH7BIYI9ASil3n1ceqQ5IWw480IAYs0xf+4/FiIxr4Z9T7UHfu/8TbQevkcHwxRM2/8nmP32K+h8fg+JeIJsSSC3j7x1Ed/9zS4Arl01h0KmwJxVc8gOZChkCkTrIlS3xFGuwskWqFtcS9fLPVS3xIN1qjB1py/EDIUY6U7xR8kc9+m4PBWL8jwy6TxmR4pwTZhdaSfwQq8F0lRQiVNzE+FwmtLXgs8/BmYDPwy+XwPsPo590pwAitNus3tHWBILMXtWFaGIxdKrlrHtrpcBiJoGEUOImoJhm9QvbaB6QRPKdalf2kg+maFmXg03La3zo9TmXT5540rENMgOpOl+uYcdD+8k1hyj5ZwWGs5chJcv4DoFnKRvuddy/il4rkdVUx12QwPZji72Pv4i3/65H7ZFa0mVSyhmY6ccxDBId6X98ClB0MghxztyA5rDIEcTebaiOKRQUkr9FkBEvqiUenPJqV+IyKTC3WrKS+kUm+90VWFFLEzbxHUKgZshRd6DiCHYhmCYgud6uNk8XQ//Brs6hnI9PNcj2lhNIZMHwLAtjFCIWEuYRS31AFQ11eK5HvnkCMr1ENPA8zwMw8DNF6hqqiMUj9H/4hbywyPkkjktjGYAyvVGtWsv2HDtKn+wY4rQn3eP0ILmkBTdDM1AJrKmNEtEFgexMwgiFM46vt3SHE9KF6JPids0V9vULkxQPa92tIyroBhSyTaE3FCerhc7CddVUz2/ifT+XqJNdRSyeezqKpTroTwPMQzCdXEAQrEoRsjCCFmowHBCDAMxTHL9Q7iOP4oOVVex+1e/R0wh2lTHHQ/vrMwoi5qDKGQLiCnY8RDd3WmcYNquKIz0wGMyzFyT8IkIpb8FnhCRncH3VuAvj1uPNCccwzbJDGSJNuZpf3p/4EbI15RMEUKGgacU+ZTjazfZPCM9A9S0thCO2JghC8O2UK5H/+Y9eE4BMQ3yyRFqT12ERGMMb9mGmAZWJEy4sR7PdTFDIQzbQuwInuthmiZuNlfux6GZIvIph7zjC6CM+6pxQ2tViN0jTjm7NiNQxsy0UzviXSmlHhKRZcBpQdYWpZR+c1Q4a9RuHj7lHMyQSfWcOGIKqQPDbGlP4oyGrfDLOp4fIp3+DG2/2oiYQrgmzL6n9rHgooW0XHgmVlUU5bnkBlNk+4awYlFUOkvfi1tH4+kUBdVI9wCeUxjNG3zseaqaEoRi0YNi72gqm54gkGR+MEtPzsUU3/GvtrybAmawm6FD3pWIXBp8/jHwDmBJkN4R5M1YRORKEdkqIm0i8uly9+d4ccW256hdWIOYgmH41nWu4qCXRnGEW/zMpx1GejOM9GaI1EbY/cRuJNjwJJZNIZvHzRcopDOYIQsrFhk9X5zic7N5XKdAan8PfZvbg2keEzGN0bKaymddf4Y5NeFg35siJML8qEVWC6WpYYq8hB/pfSc+twbnXxKRcyZa91g43Bvg4uDzXeOkd07FxacjImIC/wFcBSwHrhGR5eXt1fFj5QOPYsdCWFGLfNrhQNZhKJhysQ1hVtgaja0EEK6xcQuvWk61XtJKclcHuf4hnGQSJ5lGed6oQQOAUSKUvLyvITnJNCO9vqm4FbEopDOkO/pI7urgxj8+VcdpmiGIaeAqhRM4+k0VfC38VndXubtW4UxN5NkJvu+uApYF6QbgtqOoe9QczvrucyJiAA8qpe6Z7IUqiFVAW4lhx1rgamBTWXt1HFnx04d5+lJ/DHJpay2/2zsEQMQwSBU86m0D2zCI2ibKVdQuqCEUC2HaJk46R3J3B6n9PSjPI9meZKQ3gxkyqJmfJtZUTbi2GuV6FLI5sn3D/vqRbRJvqcZJ+5tprVgUz3FIdQyRS+bL+Tg0U0isqYqeoSx1IZO4ZdCZLejpuyliivYpTeR9dzXwfaWUAtaJSK2ItODbF0z5u/Kwd6WU8oCbJnOBCmQusK/ke3uQdxAicoOIrBeR9T09PSesc8eLCx/7Lfs39uKkHN44t4a3nNLApRfMYXbEJFEdJha1sOMhInUR6hbXUt0SJ94UIxQLM9KbIT+cpX97PyO9GVIHUjhpByeVZ7h9kP5tB0ju7WZoVy8SmPQpV6E8DysawghZeI5DIZPHSeXJp3yhpLWlyqVoWXfXuv081ZfBNoRXkjmKOraOozQFTFxTaiy+q4J0Q0krE3nfHarMhN6VR8tEzDceEZG/B+4G0sVMpdRM3Wo/3iTsa4Z2SqnbgdsBVq5cOSOGfu/r2gjA+rdfNpr3uktbqVsy25+GcQp4ed9qqrgupDwPwxR6t/SRHcjiuYp0wSXdmaK/K40pkGiKAf5eFStqUdUYpWZeDWbICrQmCyedRbn+K0t56iDvEJrKY+yA4pVkjvPqItTbM9M1zolGIXjjvqrGpVcptfIQ5ybyvjtUmQm9K4+WiQiljwafHx9z4cXjlJ0JtAPzS77PAw6UqS9lIZwIA74QidRGCdVU4eUL2DVVuPmCb/JdW01ucJhs3zCZgSyFbAEvWHyyDf+3WlyL6u9K4yqFKUJtYxW5oRy5mhyhiL92pTzPF1CGP61nhkxGevVIulIpFUjvbKkm6bijroUaFiZw8y67r72ahl88MqoxaW/hR4vCmxp3TRN53x2qjD2BukfNREzCF032IhXGs8CyYJPwfmA1cG15u3RiOevHD9B20wcAMCNhQrEI4bnV5JNpjJC/Bym5qwMnncWMhAjX+EIsXOMLpirXo5At4OY93LyLW/BGBdRIvx9Tx7RNrKjlG1m4FmbIwq6JYZgGoViI4Y5UuW5fM0W8Z3EdYgoNhBhJ5phzVhOJhQmU69H1co+ewpskUzQ9M5H33f3ATcGa0fnAkFKqQ0R6JlD3qDmiUDqE+fcQ8LJSqnuyHZhuKKUKInIT8DBgAncopTaWuVsnlA3vehu1i/3YR8p18fIF9jz8LMpTRBvi2NUxsn1+LHQzEqJmXmLUE3ghU/DXi0pM9orHhggSOHbNDmQJ5f1AgaZtUsjmR03CrWiISCDotJfwysW0DayIhR2zqWqIYoYMxDBofN0SfvTjjZxV7g5WMAqYCnuRQ73vROTG4Pwa4AHg7UAbMAJ85HB1J9uniUzfXQ9cCDwefL8EWAecIiJfUEr9YLKdmG4opR7A/0OctAzu7MN1XAxDGGjrwYpaZAeyDOwcJBQLke4awXVcauZVk5hX40+75T3sWIhcMo8TxMrxXIVpG5iAk3MxXHDzLkYQN6lI0R2RmCZVsxsI16b427+Ic8t/PleGu9ccK6VTd/GWONG6CPm0g2EIkboIVU119L60g9Or/UGHnrY7dtQUeVsf730XCKPiseLg5ZvD1p0sExFKHnC6UqoLQESa8e3UzweeBGacUDrZySVf67BjuCOFGTJG142KZHoz2DEbwxQK2QKFLIgphCIWI+lXTbuL9RzPIxQYfSpP4aQcqlviGKaBmCaRBt8AQhpqyPb5pulaW6o8PnrFYgxDUK7CDJkYpmBFQ+QGhwnXVvOmrc+Uu4sVzVRpStORiQil1qJACugGTlFK9YuIdmA1w7hv9hnEAy0mM5gjFDaJt8RJHUgRqYtQPSeOm3eJ1kVwsgVyyRymbWDHbKyo/3MqZAoYRuBV3Hl1I62Xd3EdF89VvqWdaWIEAQRnndlAVUs9ViSMmAbpjj6siM2fX7WE/3pwRzkfiWaClGpJhiG4jv/3D9f4rqMiDQky3QNUL2hm1sAmuhvPOKiOHngcBQrcGSqUJrL76n9E5Jcicp2IXIe/6PWkiMSAwePbPU05SA1kCcVs+vMu2WyBQrbArDMa/Wm4kIGYQqy5CjNkEK2LEG2IYUZCWBELMQQraiGmEVjS+esKRuAvL1IbIVwT9g0dgvUG0/Y9i4thUMjmKGRyfv2ITaw5Vu7HoTlKPnrFYkb6MuRTecI1NoVMgZoFDQy2dZJL5kh39DH8xP18p2VFubta0SilJpQqjYloSt8CTgXehG+Xfie+3XsaeMtx7JumDBR32+eSOfrzLnHLIJ92WPDmZt93nevRs6kb5SqazmrGMA2sWHTUu7ebd5Eg/tLoRlkvmMKxDSzT/8l5rkco4nsWD9dESXf0EW2qw4rYmKHQaH+S7cN8fPWM9fI0YyjVeJyUw/CBFDXzqn03UlGLjmf3YUX9fWlWZIjen/yeC5fWse9l31ZKa0lHhwJmapjEiQiltfjrRn8HRIGvACvxjR80M4yY5W9uHB7KkQj5ivTc81qI1CcI18WJzqoj2rSH3OAwkfqEX3ZfF9Xzm1GeRz6ZJjeY8vct5b1RgVS0unt1fcHXqopYkTC5wRTWbD+sRbYvSbZvmMTCBDJDI2zOFIoC6Y0NUQA2buxhVtjEjttUNVaNGr30bPQ9erS91M3cOdW0nDub96QdrtzxfJl6XtlUoBI0ISYilM7HF0S/B6qBu4A3Hs9OacpHcV1oaCBDqqCImsrXhqI2VS3N4Lnkk2kkiKvUuPIMPKdApKGG3GAKMU3cvItyFWIKbtrfIBuyQ4Qilp8X7F+yIhahWBjTtnCdAlbEJp8cIVwXx4rYRJv8oIOZbj1LPN0YzwWU73DVwxQIGQaRugiZgSymbfgh7xujdAb7z4zAs0PzWTpe6LFyMhs6OEAGX0uKALsCn3iaGUhx1Prl2DIAXGXgOi6FTJ6Rji7s6hixlgaMkEXV7AbMuiYaL4iC52KYB/xAf1s7AYjWRRBTRvct5dPOqGVecT+S8jwK2TyRBl/rKoa2MEIWZPMM7eom3aU3WU4nxnps6Mo6JEK+kJlbF/XXkOZXYwRRZwd3+05293QMM+R4pAoeqbYBqhqjrHrosfLcRIWjFKNeMmYaE5kXeRZfKJ2Hv650jYjce1x7pSkrN1ctwxQZTcMHUox09jGweQ/J3R1UL1lIbMFcQvOWYjXNxUw0oBx/82tq/6vOacUUqhqrsKIWZjAyNoIpu/+/vTePk6SsD//fn6ru6u7pnp57Z2fvg11gOURYEUM0KCBINHhgRI2QxIRgNDHJzyQakogaEo9EjUFFVBSNBo3KF1SUAIIocgtyLHvB3jv31dM93V3VVc/vj+fpnt5hdnf2nJnd5/161aurn6rnqaNn6lOf4/l8vLT2G0WmNLafK1Ap+agwIvQrOF6MYt8IxeES8Uwcy+ygXiC98xydezPlOniO0NqZxvUcWle16ICXZIz0/Ga8jA5myVcUuUpEKVJsHfetQDpElJreMteY1uRZpdRjZr0HuFRE3nUEz8kyg1wty1iciuNHiuXpOItWNDOwcYjOMxbQuLgTAG/ZyUQrzqKSbCYUUD++nr5HnmZkcw8j23Iksp7RhvSfVyKbwDFBD5Wi9i+Uc2VUqGrzlyrFANdz8QsBbtxhYMMQiazHzfdus07wWUK9QHrbGZ3kdo4x7FfIxBxSzQnKubKeLB0vAZCa10J5eIxlF57Kzl+sp2PnGMUwYjQI+UzFJts9FPQ8pTkocabBdHLfPTZFm50wewxSfejsKAac3Jgg5QqNXRle9fCvavus/5M3MvDUDZzw6S/ihmUkKFIcGWFse38tzVA8GSN0I1w/3COgwXEd3Lir/Qyu7DERtzqnRYURpbxP7/ZR7u23ZrvZwGT/0bltKfyCnqLYnvHMb+zo+Wzdecq5MrFUjF0PvEDbSW0MPbeNBb+1muYVbZxV9Dnxxh8c/Ys4Bjk2RdL0NCXLcchzY2VcgeZto3u0L7vsEgBUIoOb60FGuqkUdTh4NcddOVc2k2Qj3FKF7KIskR8STyeolHxiqVgtvDyKFKEf4hd8KsVKTTABXDI/wx09NjHr0WZvdaxe1d5A1kRSAjQuyFDOlRHHQUURhb5xxHWoFCuIq3/D/mf7CQoB2+7fwZnveRXzPvCfR+kqjn2O50AHy3FKvhKRnJSjzlt9JgBqZBTPIv4AACAASURBVAcSlAkLOaKgQuiHe5jmgJqGVBwYR1zBL/jEkjHcuEOlWjspVDoDhOsQFALCICT0I7INcVLtDdCTt2mGjhBV4XOD2jqtgorZuEukFIlUvKbxxpMxolCZsH39m7qeFkzjg+P4+YCwEuFHime+8SAb/uUk3jOy/khd0nHFMWq9s0LJMkH1wV99QK0xmRfqeT5zIn6oWJXI42x7ivGnHmHouW2Evn4glXNlin5YM825nn54KVPgT4WqFg4cRVogqUjhl/xan/S8BsYHizy3eZiTGxM8N/biXHyWA2NfQqd+23svX2MmPzv6tzS/XTU9VOSHJLIJ/IKPl4kThQoHUGGE48YJShUqxUgHNgyXapOxsw1xerrzjAYhH4iv4N+DF47k5R7zKKWO6+g7y3GKI/KiCKmYA2Z+LariM943TFAoUylV9Bwk81bsB6H2MZW0FhWFkfnUgqhSnGivCiPHFeJpnSetb6hIb7nC1nEfy8FztSzbp0D6ty9cvsf3z9+yjnLOrwkhcQXH04+JaqRkOVfGS3t7BLBUX15UqIhn4gTlkGKoakIp05VhNJioq2VrKR06kZreMtewmpLlRdSbyq5vOhGAYqho9Rz+4Af/hNvSQZQbYvSJRwly45Rz2iw3tlu/CVcfPHEjeLy0RxRGtTduxxH8QqALAAYRQcHHjbvE03Ec16E4XGLI1/NZiqGyprtpMB3z2+u7Guk4pZ2v3b2Fv/7TM3G9GKNbu1+0X3UCNejglMiYWh3PNaZWISjp9EGJbAJxRdfHSmvNqTRcYjQIJwRSzGHjhkHylYglDTEuXLv4Rce0HBgKa76zHKcMmLdjB8hXhK3f/j6LXn0m/kie3JZuCj0jtSCFvr4CxVDRXw5rs/sX9I7jrfBMIb8K0WiEl/Hw0nGKfmjmJUX4+YAoVGQXNwLQ1JPHczx+M1qawaufG+xNIC1IxgiUoiXuknIdSmHEwPpB3nJSG2M7R0k0JUiGEX/7l+dQHBiraUiVos7+XvUNukYYiatzGkahAj9kbHe+JsCUeekICj6Fin4xSbkO87syjPYVaIq7dCbjzD9jHsu/88OjdGeObaJjNP7OCiXLPukxecsAfmtFC41LOol1LmH0+QcZ7xulOFyiob2B3I4cniMUQ0WgVC1bZKESkilV8MxbtPJDMp1pgpIObqhSjdwbHyhSHBgnVIrfjJaslnQAXLq0iaAcsn080BnejaaSN+XoQ6Xw/JBV7SlczyXZrAvsVYraROq4QmiK0ZRz5VpJe9dzamZViHCAMAgpjpRxYw7pzrSOuiv5BGWtIbkCq1/WRXG4RKYlSXGkTMuK5ppAssX9Dh2rKVmOa37/tHmkO9M62Wrvdgq7+hnrzlMcKJLIJoilYmQyHik/ZDQAV7TJLxOD0nBJFwE00XmjO8doaE8Rz8Rr/govHa+Zh3LjAbuKlX2dznHPVNrRxqESq1uTuMWKycahfTq5iqIQRqRch5a4y7MvjPBS4xNKtKTBrxAUAoJShSAfMF4qUhouEfqRSaTrTGhI6IKN1SAGN1IoozElsh6pZkiEioa2FMmWJOI6JJoSXPjc1Nkb6n1LVlBNn2N58uyMBDqIyFtF5FkRiURk7aRtHxKRzSKyQUQuqms/S0SeNts+JyJi2hMi8h3T/rCILKvrc6WIbDLLlXXty82+m0xfz7SLGXuziDwlImce6Xsx27kut47rcutY8bqTyS5qJJb06H9yEwCZzjSZBRm8dFy/ebck8ZyJ9ESgS2EoE+SQyHo1x3hDW4q4qacUS8ZwPJeg4Ne2W/bOVAJpWUOcfCViIO+TiQkpd+I3qN7SuAiugCtCcbiktaGiTg8FOttGFEZE/kSQQywVwzV57SJfB6+Uc2WjeenftxRG+OMBI915GtoaOO2KlzPvtE7i6QSpluRe0wnZYIeDRykIQjWt5VAQkVYRucs8K+8SkZYp9lksIveKyHPmuf7+um3XisguEXnSLJfs75gzpSk9A7wZ+FJ9o4isAS4HTgEWAHeLyGqlVIguwX4V8BC6JvzFwE+AdwPDSqkTRORydEbzt4lIK/BhdJkNBTwuIrcrpYbNPp9RSt0iIjeYMb4IvA5YZZaXM1H2/bin8R+/SCPQ8POv4SQb2PnDO0m1ZWomoPJomfzuPMVQ4TlCqNTEQ6tUIYnO2lANBx8fLJJqSVIe1bP/VRQQS8ZqmQIs++fPLl1N37oBntg1hh8Zsymw5hWLjBY7TlAOeXpU+4Y6Ei6ZmEOoFCqKKJnMC/G0VxcB6RBLxfCCuKmB5da03DAIifyQfLmC54h5AdFCrnF+ms7T5pFZ2E4s6dF68lJaLv9ztrjz2TmUZ1Fr5kXn39bYUBNMVks6UI5aSPgHgXuUUh8XkQ+a738/aZ8K8P8ppX4tIo3oZ+1dSql1ZvtnlFL/Pt0DzohQUko9B2CUnXouBW5RSpWBLSKyGThbRLYCWaXUg6bfN4A3ooXSpcC1pv/3gOuNFnURcJdSasj0uQu4WERuAV4DvMP0udn0/6IZ6xtKl2t8SESaRaRLKfXiEKXjlPjiVUS5IZKtTYjrMN4ziJ8rkmjSUVgpVxjy9T+LH6laWHA55+Ol45RzZRxTkRZ0gtZqDjxxtZnIFan5Q+zE2T2p15LGB4oMDpcYDkJSrkPScQgVtJ44n+yiIsNbhij0FsgUAoqhLimRjWuNNmWW6tyjahaNaoh+LKmT6IojqGjCbOd4LhkjwFItyZomddLl55BZcxpEOjCGmIeUC6xQm+huWs3OoTxenRbsOkJbY4MVRgfJUTTfXQqcZ9ZvBu5jklAyz8dusz4mIs8BC4F1HASzzae0EK0JVdlp2gKzPrm92mcHgFKqIiKjQFt9+6Q+bcCIUqqyr7EmbbNCyTC64rcBiJ/+Blo23kPxmcfIbe2m/fQG3Pg6wse69/AH+ZEiX4lgrIyb9wmVwjXzlLwRHVlXfQAWh0uUx3yG/JCWuEsx1CHhVjC9mHedu5h0ZwNsGGRl2iPlCu0Zj1R7AzsfeIHxwSLZRY269LxAU9xlSWeapqVNNHZNaC3Vulfjg0VAZ+HQKaF0ZvbAlBsB7ferL2XvuELDvAyuF6OwawB/7GG8xgZiqQROPEaUG8TtWEiYXY0I+KHCNdqV5RBREE6/gFC7iNTnML1RKXXjNPt2Vl/KlVLdIjJvXzsb98lLgYfrmt8nIlcAj6E1quF9jXHEhJKI3A3Mn2LTNUqp2/bWbYo2tY/2g+lzMGO9CBG5Cm1OZMmSJVPtckzjmhRC+V39uEmPxhOWccqppzD2j1/DG/MphhH6dk6Y8UKlapFZAE5x4nZX/Rqh0ts9R2j39J/njqI16VWppgT65gM7eM9bT2Z+V4ah3gIpzyWWilEcGGfpq5bQtgpGto3SuCBD5pl+PPN7Vf1DsaSOpquU/NqcMXEdMJk39DY9F6kaoOJ6bk07iiVjuuYVEBRKFAd1jkKvMUks5ZFobqRxSYVYIcf8cgl19hvpK5nZnI6wZ54Qy4FygJrSgFJq7d427utZfSDnJCIZ4PvAXymlcqb5i8DH0Kf8MeA/gD/e1zhHTCgppS44iG47gfqZdYuA3aZ90RTt9X12ikgMaAKGTPt5k/rcBwwAzSISM9rSVGNNdZw9MG8aNwKsXbv22AyD2Q+5ky7AWX0+jXd8BmnIEl+ymt/+0j+y46IPsKuoI75CpX0OeaM8VR+O+UqFTKDnM6WKFTqbk2QXNZIqJGkYLrG+pwBoIZVyxWpLU+Cl4zQtbWLpq1ew44FtgM6csOXebXgZ7RNKtSRZcUILoR/iei4qjPALASpUJFrSSFAhMlWAXc+taazxuIu4DuJIrVowYEpTmDlnJuihUqxQKVVqKYkSWY/Grgxu0qNhXgtuqUC8fzM0riQCFjWnZ/CuHRsoIDhM6Rr29awWkd6qC0NEuoC+vewXRwukbymlamnglVK9dft8GfjR/s5ntqUZuh243ETULUcHHDxi1McxETnH+IuuAG6r61ONrLsM+JnxCd0JvFZEWkzEyGuBO822e82+mL71Y11hovDOAUatP+nFVH0BbY0NxB2h9Pq/IejZgb/+MaKSnihZFT5+pGoL6DDxYqi1ptEgqmVuCIOIyI9qZqOq474QRmRis+3PdOao9yktff2r6DxzGW2nrKD9xFZKwyV2PTvA9jGf/oEi/T15BjYMkp6XJpFN1JLrlnNlglIFPzeOCiNiqdgeGlBDewrfREJWBVVNOzKJdcuj5ZpAqubIK4/5FAs+peES44PFmp/KW3EKUtaa1JIpAh4sB4GCMFLTWg6R+udr/bOyhnkmfxV4Tin16Unbuuq+vgkd5LZPZsSnJCJvAv4L6AB+LCJPKqUuUko9KyLfRTvIKsB7TeQdwHuAr6PLsv/ELKBvxjdNUMQQOnoPpdSQiHwMXTkX4KPVoAe0o+4WEfkX4AkzBuiovkuAzcA48EeH/eLnGPURUtXvU+G/6e8oRIq29XfR0JoiVazoNEEoSlFUF54se5jwXIGOhrj2T7Snam/bxTCqRRfFXxwQc9zz7//9x1AJaD39RHbc+SuGXxghn/dJp2K0LchQGi5RKlXw84FOAZSJk2pJEvqhKU+vBUaiOUPYp7M7VBOxArVUUMpoRq4p/lvVmso5nyiMCAo6E8fISIliNQmvhLSf2KaTtMZjRIUxovEc2ft+zJN3Ps4Z37+zdh12ntLBoVBHK9Dh48B3ReTdwHbgrQAisgD4ilLqEuBc4F3A0yLypOn3D0qpO4BPisgZaOVuK/Bn+zvgTEXf3Qrcupdt1wHXTdH+GHDqFO0lzI2aYttNwE1TtL8AnD1FuwLeu5/TP+7Y18Oiuq1nVJvbBk+6kPO+3cDYRX/NsuYEu3NlekphzZ9UTY0SmqQPcYEtuTKrUjHiNUe61ra06c4x/ilLPd6S1agwJOzfxdO3bSBUipPOW8rqP3gd49t3EAUVfv3Fn/PU7jzhUJElDXGWnthG6yptyquYTB3iumSXzMMfK+AmE7jxGOI6uMlh/LES+V79uzquQywZM+tSS0dULexYDE1AC9pE2/NkL9nl84k1t1Lp3Y7yS4xt15acx99wIX7Bp9A7zqm/ug+AuGNfPA6UQ5yCNC2UUoPA+VO070a/wKOU+iVT++NRSh1wlfLZFn1nmaPMb5rwEwR+idf/8N9wGlv46YXvpaekI7sitI9IR9bpN71AaXNeecwn8iNUMiKRTRDuHKMp7pCJOXiOy8a8v0f9n+OZv/ur30IFPvFFK1GlAmsuWkHz6kVkFnagigWaX/O7lJ99hDP+5Fy6P3E3u4oVNuZ9Nj7ezeJ1/Zy0pgPHcxDXIZYsEnm66GJYKhNPJ4mnkxOmt7z2G4mnnzniOloY1pUnAYxJltrcJcdzaVq5EFUsoMolVMVn5Pnumm+qyq1Lz+L3d/z66N28Y4RjOaODFUqWw078Zb9XW390YJykyfLQFHdpijvsKlYYrssi7ccU/eUKXncecXUYc6vn4keKdMylZUUzV5+7CL8QcNOdtg7PJz/7K/42iph31kkArLzizeC4tRBs8ZLEX/kW2s/O8zt9I2z5vw1s2J6jGEb4ETz8m15cERamYjS1JmlZ3ky6M40Tj5FozpBZuZxwbITclm7y3WO1siRu3K1le6iUKpRHfSKlTa2ZmOCKQ8oVmrsydJ3ZxXjPIKDnPo0+v4ueJ3uJJWM6KGKkzOa8z8a8z2+yazinNcUbtj4OWDPetDA+pWMRK5QsRxRtttPrJy1vIpFNwLoBdpcCE/4tNaHljwf4hYB4XVaIUhjhF3xasy14Jino8RqJV9UU/+R1KwEYfX4XbtKDXf01zSazcJhkFOGsaCJKZul8+WmUBkeJQsXY4Dhxx2Frwa+Z21L5gIb2BjpeuorIr5BespBY5xKIQsLSNm2qM4IonowRhhPnI67g+1EtcKXVc2ial6b9pDZUFBFrSDLeM8TQxt0Mrh/U59BToBhGjAYRG/O2VtbBcjij72YbVihZjigfKz1fW3/ija8FYOmJbWz5dQ/DgfY1DQchXkkLIbe3gJeJEypFvqIohhWKO3I0L22icVEzH/v0m/Cy6eNWMIGeL5TbNoTsHKZz7QlkFnaw+Qe/Iru0ldzWbjoKRVI923Eam/GHR1j6+lfR8+T/4sVd/CCkKe7Q6gkt6TgnXXY6S97/dzBqIn3TLUhUIRzu0yHhtXBwHUru1mfjcIWk65A0qYk6T+sgKFUY250nu6iRjT94hHx3ntGhUq02VjGMyFV0EIvn6Awg7V6Mjob4zN3QOYg131ksh4GFN99Ki69NOsMv+V12jAf0lnWUXikyQsgPCUcmHmC6LpNirDtPpisLQLy5eSYvY0aoDwUf687T2JWhdfUCMgs7iGWbOOHNv0VuSzeRX8EfG2f3A/chjtD1ijVIIsX8MzrpebIXZ7hEPu+TiTmc8vYzaT9jNVGiEScTIpUSqpiDWJyokCPydTBEfZJccYUoUohjfFIpB9dzaD2hlVKurKPyTOLWSlFH/1XNtNXCf/mKngaQiTmkXYdWz6HjlHYArsmuAazfcL8oRWQ1JYvl0HAdIZfUD5+Lf/wpRh+4h6e//ku2bhlh63hgqtYqk2UcSpEyk2/1P1+hb4xkW45K6fg1+7xxhY6ei6fjhEGF0uAo2UUr8LJttKQaENch8iusfPN59PzqSW0+e3YLqXktLHhZjF0P7aQr7dG8NEvHuWsZe+UfEiroHNlF6Yn7UVGIOC7jvQOUR3Thv1gqVquvVM3w4LhiJtw6uJ5Loa9Aobegw9CNL9D1HIqmSGQx1HPOqqa+TMyhMxGjI+Fy4kvmseKSl9YEEhy/Jtrpojg60XczgRVKlqNG/ZyngRWvZF66kZPHS5S/8gA9pYrRlnTdn6a4Q1ykZuJJtiRxPZex7X2MdeuJmMdLNF69lpTIeiRbkjTMa8JrbCDR3MjIY4+SXb2C2KKVdCw7CSedpbLodJaseBAVhWz6ry+TXd7FyOYeTn77y9j1i/V0vfwEEiedRdERnukt0D7cT1guEwUVxHHI7xqgnCvq8PGibgMI8sFEeXQzubZSrDDaV6gl4PUcXfLCS8cJ/YhiqDUu/YKh56pVf+PfjJboeXg3p7276Wjf1jmPNd9ZLIeBqmCKO8Jw10tRf3YGa3J/yWOfvr+WHw8wfg+dIW1+Uj/gvMYkYSlgaNPQHmMey2/VVYH0tjM6Cf2Ixq4Myy85Ey+bZue9v2Zo/Q7mn30ildwoidOX8tTfXsPqd1xI+ec/pemiyxiffworP/8qYiM7aDr1l0TFAq1/+XEGfAe16yFcR/idynqGfvVzeh/bSKFvnMauDFEYMT6ohRLoPHhhENZqZ0WhIgxCxgeKBFFE3HHoKQW1ydDZRY2Uc2VczyEs6wwd2ZhjXjK0MHKAlWmPjoTLd/7i23tc97H6ex4udD2lY3P+nhVKlqPO5CwRTSsX1tark2qhrjhdwkWFiljKIywFOl3OtlFOzSYA/cD7etvJPDRUPOYeZtUErN95Uk88vfqkVrysDt924jF6H+lmfKBIFEZkfvJIbR5Q4znnMdi+htFChcVJ6M+uwHvlSsqhwi8pRBRR5yoyW35F7//7Hr/5yq/YMVpmcVOC4nCJxq6MSSHkEwahKdKYqJWzcFwhKIRESpFMxgj9CD/SeQpbV7UQ+VFtgm619lJTPMbCeQ08sHUU0JnLOxIuHfMzZHNlbhhaPyP3eC5izXcWy2GmrbGhlgXing99n0zMIVfRb+XVAoGuQMpUPwUY79OJh2PJGKszHqNBVJus+dCQnqA7VVVWePGbd73pb7abAavndbUs44YfbOC6V5+Ol03jejEG+wv09RUIleJlbzyRxa8+g6GnNtDx2j9hV87HESGf9HhwxyjtDR67xkpcvLKFnB8RpttwYtu45RP3MBxM+H6WlCqkWpJ6LlKuTFjRAseLqmY7p5YOKp5wSbU3sHXLCClXaDKReUGpouc1hYqsPxHYMDRYZEEyRkvcZXk6TkNriv6ePCvXds3632G2Yc13FssR4vS3ncaGrzz+opRCKVdo8WK1h1sUKirFCkEhYElDvOafAFi7pp1KsUKyJUnbie1kl3fxkX/4cW2svQmrye3132fbw7EqQK/5i//lA+89m76ne+gvh/SWKyQdh9QPNzG4fpCz/vFdOOvv59RFJzPQsIAXhstcMvAzYstOobJyObHh50l7GcLUfMrPPswZy5v5f8/20xR3cAXijkOht0BxuERPwceP9G/hDpcI/YhE1jOJWh0a2hooDpdoijuMBhHpmEvkR7hxh6CgI/HimXgtI3x/OWTFgkZiqRg7to4w1JMnX4nYfv+O/V6/ZQJ19CrPHnVs+mXLjBF3hLgjqFAnbHWYKG3h6dI+DJQr5PvHdWbqXBkVRqTaU3SuaKa5K0N2YSMNrSkqxQqJbAIVKoKCLv/9n/det0/BcuWrl9bW37Aou1fBNZuoXs+/f/4Rnd3bFUIFhTCip1Rh03MDxDqX8MQ/fZbg0Z/Q4lZY0Oghp78alUjTU3aRYg4VT+A+cxfF7l627MzhOaIFmysEUUQ5V2aoGNSi5aqTbfN5n5HufC3cO99XYGRgnC0FvW8uCInCCdNdNeDhhLPmmyS78MS2UX69YZAHBou1JK7VsPHZ9iIwazl6WcKPOlZTssw4G25bz0uXNjEyMM6QHzIa6HQ4rugHYagUzvoh0p0NJLIJYqmYrtfjuRR6C+R2jhnBFRBLxgiDkFgqTsv738jg2HjtQTd4/d8CcMe1P2LhgkZuvncbF8xLs/j0eSSyHpdUIjbnfXIV/YCfrVQ1ppvufIHLz+pi41O9hEoXQ3RFuOG897MwFeOMjoUQBijlQhQRpZpYNLyRwR9/j5bXvA4cF9eLsyCb4LkxvxaEkDKZM1wT/VgfMVfN0LF9XIfl5weLBCbZbtrVc8uKGwbpSMRINSdoaGvQIeWjZdo60jSVKuwaLtJTCsnEHEaDkN+/6mXsemgbr3zwgRm8q3MLxbGbZshqSpYZ45rsGq7JrmHzaJnmpVk6VzSztKOBTEy/sadcnX5oNIhwPceUXDDpbgKdYSDfN04YhDiii9G5nv6Tzu8ew+t+FoBiqYQ/0seOe59ix71PESoYH9A+qK5VrTS0p2hc1EKi0aMprid0zhVOeec5tfUWU1U2VDAaRGz7xv/gFAYREXAcokQjKp7Az41T6d7K6IP30fvYel24T7S2Bbp4oLhCNu6SiTk0xV2jkU3UxspXIgb8ipl7pF8cSlFUq5c1GoSI4xDUCfdq8b/qeClXeOdfv4qmj36FhecsxTJ9lAK/Ek1rmWtYTcky4zw3VsZ9YCerOtOMjGj/RIuZ45Jy9UMvlorheC7xTJx4Ok48nWD3o7sY6c4z5If6Tb4Sks4HxDNx/HzAun/+CK0nLyG99qVIyzxCP2LXY924os2Cr+5oYKw7TzlXZn7cpf3ENm414eaHy4x0pJz3VW3pn/7mVhYkYwwHIQN+aOYK+bxhTQcLLziHaNszdJzUhEo0Elv/c/p/+kPmXfw6nv3kFwn9iEd+sYP+cqUuB6HO2gBQqIT0l/WYOlu7Nq/6kaoJo+pE52rJ+7iZ6NzRlKRlRTOZrkZiKY9kaxPF/mHGdo6wJNNJev0gO0bLBAX9crDiM9/e26VapkAxN01z08EKJcuMUR/59kyuTG+5wqnZBJ4jtbIISVc73/PDJTKA25UhKATE0zo8Oe44QFh7MOaCkMyYwo05jG7LEQZb8cfGKQ2OMrR5iJ5SBVeoVbP1xwMA3GSc+WsXwI82HbIAmco3NbmtPqJuqvbpUL1/y9Nxhkd0NoymuEOoINmSpO+Rp+k818Mp30N8wTL8/l00n7aGSv8uul6+mie/8gC7S9oXtCAZY76pl1QpVSiP+caMOuFQrwqkUFHLvDEVKVdIZD1aTuggNa+FhnktxLNZkr39pBd2EOTGcRzhtp88z/Off4jRz66xvqQDxWYJt1iODPWCqb8csqWghUR14myooL9cIQJWh4riSDcNrSlCX8+dCSLtf8pXIjoSLi3puDZHeS6O51ApVth+36ZaCpzq3I5iqGhPxMgubiS7qBHXi/Gxa+88IgLpjy5Yztfu3gLAZWs6+N66/j32vfrNJ+p78YMNB3XMBwaLvOvcxTS0p1hywUuIggptl76T0uM/w22ZR9/d95Db0kNmYTuL3v1n5H9xB7f9211szPtGkLksaYiT8lzy5Qq7+sZxBYb8kFKkTGYNfc+q5jvXtMFEkIIrQkcixqL2BlpXteJ4cZLNjSQWLiF26rnEE40ABPd/l2Rblg+saGP7/Vt5zVOPHNR1H88cyz4lK5Qss4qt4wGLU3EWZnTY8Yb+cXKViGUNcXpKFVKuw+LhEnGTwDOZjJEKJvwa8bRXm9dULXWRbEnqwXsLZl89yTPR5NHQliI9r/GwnHu9QLr8rC5SLUnmvWQR5ZEx3veOU7j+28/yvXX9ewjif7r2IvK7+vHHSod07CXnrWa8b5j0ihV86w/+k3Of3UJuZ46e3/Tx8r9+NV42xcbbn+Gpm/+UXBCyu87Xk3REl5MohDUzXDWqL1QKBy2QqqXs64VQdU5ZyhWyMT0ZtqEthRt3dOHAIECSaaJtz+I0ZMFxiC9ZTcuiE9h69+f2KBRomT7KakoWy5Gn3qQ10F2hKe7W0tM0xR0eGiriin6Axh7vofO0Dn69bZRWz6WjIc6omezpuEK6M41f8AkKgdac4g7pzjSL0l5tQmh2URZxHYJCmdS85v1qSVOZ2ia3XTAvTVOrFoK6MF6M1pOW4iY9PvKSFez8xfpan09/731QCWg6YSE9Dz5z0PfsalnGddfdw5++fhWjz6zjD7Y9RHTP1/jPd36B3aUKT1/7UwrhhCnOr3uYuSKUoomQ72obaEFTinRi3GrfqhDyHB3C74jgxrQ/aklDjOb2BoJCwNDmPnwrgQAAHK1JREFUYXI7x9j843W0nvAUKowY3jJCcaRs/F4RWwoBuUrEBQd15ZajIZREpBX4DrAM2Ar8vlJqeIr9tgJjQAhUlFJrD6R/PTMSZiQinxKR9SLylIjcKiLNdds+JCKbRWSDiFxU136WiDxttn1ORP/niEhCRL5j2h8WkWV1fa4UkU1mubKufbnZd5Pp65l2MWNvNud25tG4H8c7N6itLxIIxVDRU6owGoQsTMUohopTs0n8SNdf2tA/Tr63wBnnLiZfiWhckKFrkZ6UGYUKP+9THCjW3sRdzyWejtPQntJLq36bdz0Hv+BT7BvZ5zyl6fiJAO7uK5DIJnBcwXGE4uAo4/0jhCWfVEcLX/7RJv7roc/yhQ3fwjvhdIjFUWHEv3/+kUM2HTqu0Lh0ATz+Y3be+QDDQUjKFXKV0NQzimoCSWtB1IRUtZx5qKhNYA4VxI2AqgqllIlMdJgIGU+aUPJiqCgNlxgfKDKyPUfvhkH6B4o8/audPPfIboojZfrLFXpKFbaPa5NsdV6a5cCIlKJciaa1HCIfBO5RSq0C7jHf98arlVJnVAXSQfQHZi4k/C7gVKXU6cBG4EMAIrIGuBw4BbgY+IKIuKbPF4GrgFVmudi0vxsYVkqdAHwG+IQZqxX4MPBy4GzgwyLSYvp8AviMuVHDZgyA19WNf5U5pmUGCRU8OlyivxzSX67w+q5GhvyQYhixY/MwiazH733iLXS+ZD6JbILmpU21FDniCqXhEkHBpzRcolKs6PILyRjKpMyJJXUphpFtOh/b1bJsygVgQTLGqdkEZzYnObctVTvH1RmPc1r193ecvQCAKFTk+8YZ2TLM6Ja+WpXYT950hb6u4X6oBJT7BvjbK7++x7EPlKow+9JtG9n20weJzV9Cy8lLSbtOzeRW1XCqS/39LYZRbT5YIYyMZqQ1peo8pdIe2pXOTxjUhYhXq8nuGA/YWvDZMOazpRDQU6rUzT1TtQCT6qdrZdJBc5Qmz14K3GzWbwbeeKT7z4j5Tin1f3VfHwIuM+uXArcopcrAFhHZDJxtVMOsUupBABH5BvrifmL6XGv6fw+43mhRFwF3KaWGTJ+7gItF5BbgNcA7TJ+bTf8vmrG+oZRSwEMi0iwiXUqp7sN8CywHyI6iDoDY3T02pUaRWHgj7Zs34Xpx1n3rlyRbkhQHingm8AGg0KeTwIoruj1SqEgRz3iUcz5vWJTFD/QDtNVzuaMnzwXz0tzdV2CBiUyrRrf1lEIWp+IsT8dpSsZwXOEt8xqoFHU27cjMqQoKPsXhEqXhEmHwGG2nLGfglu/ij43T/Xg3seTh+ResN33+lfsfLD5/LWcvaKRSrHBffwFQNQEw2Y1TzTNYilTtLdUVnTDVc4SwPDGHyY9CIqjTuCbyFLZ6LvOTsVqoeH853CNyr2tVK+OD46RHfXpLAe/PbT4s1348coA+pXYReazu+41KqRun2bez+vxTSnWLyLy9nRLwfyKigC/VjT/d/jVmg0/pj9E2R4CFaCFVZadpC8z65PZqnx0ASqmKiIwCbfXtk/q0ASNKqcq+xpq07UVCSUSuQmtTLFmyZHpXapkWk5OkTmeuz8gr3gmvgMzj32fZhcOs+x/9P5jI6mAHFSn8QoAKIyrFiMgPSXema/6mWCrGaKlSZ7pSXDI/gyPCpUub2DFaptXTD+khX5vFWj2HRSuaCf0QFapa6XA37iKOIvB9ks0NeBk9aTQ9v41YMkFpcIy+dQN07xxj8QktXP/rGyg+/RDrvvnzWpRePQc61+mzX30SvvokAO86dzFhn64/NRHCrSZpS/qLg9aA0q5TM6tlYg75SkSg9lRp6oWRK/CSpiStJmu7U6yQC8Kaj6qaqDUK9b0BeLfNCH7IHEDuu4FJJrU9EJG7gflTbLrmAE7nXKXUbiN07hKR9Uqp+w+gf40jJpT2daFKqdvMPtcAFeBb1W5T7K/20X4wfQ5mrBc36jeBGwHWrl17bIbBzCD1D+DpPIznN6X1ymuuIPGyy3jp+4GbP0xuaw9b/m8DYkoqbH+il3QqRseadsR1GB8YN9pNREs6TqpUqWUlaEq6hJWIUrFSS1ba0JoiEWghlGxJ4sa1r8r1XLx03FRjdfELAaPbRgkD/WAeHxgH+thy13NEocJxhZZ0nI0bBhm+5APc0ZOf8rrqzXnTqRs1ef7TNx/YQUfCrYVzA3UTXquCRZvpWj13ovyE5zIahCZScUJIVU2nVaHmirAgGdNmvPEAR4RSGNW0soWpGCecNZ+utcsZ2rib3I4cCZOR3HLwHM7Js0qpvcaaiEhv1VokIl1A317G2G0++0TkVrTL5H5gWv3rOWJCaV8XCjoIAXg9cL4xl4HWTBbX7bYI2G3aF03RXt9np4jEgCZgyLSfN6nPfcAA0CwiMaMtTTXWVMexzBHaGhv0yvs+RRvw4yY9D+h3z19NIptgdNsoQanCot86gfyufkrDJUZ3jpHuTJP0Q5qBcs6nWND53VwBPwIv7hIUfBzPrfmiADw3TjwZw0t7pOY1EwUBY919jGzPMRqEtPbGKfoTE3wzGY9YMkaqvQF/pMzmvM85rSkWzmvg++sHpxQ+i1NxdhSDAy5oeNmaDm5fPwDs+WathYwWTtpvNDEx1o/03LBqeHioIBOTWr9ASa3CbNJUBq7OK+stBbV+SxpirDl/GV2vOJmw5NNx+jK89K6a/2462HIWU1NNM3QUuB24Evi4+bxt8g4ikgYcpdSYWX8t8NHp9p/MjJjvRORi4O+B31FKjddtuh34toh8GliADjh4RCkVisiYiJwDPAxcAfxXXZ8rgQfRvqmfKaWUiNwJ/GtdcMNrgQ+ZbfeafW9hzxt1O/A+43d6OTBq/Ulzn/eNbqjVbkoBK3/0H7iNzQSDA7ScvJSeB58hno4z/MKIFhYtSTpPa2Bk26g2zUXKFLZzKOd0BvKgVNG5+CJtCnQ9l+zSVkY29zKybZSR7jzFUPF8IWBLISBiIqqoqVgxhe+gI+Ey5Jt6UKYmFEwd2Tc/GWNBMnZAgmncjFkVSHuWJJ9ItlqNzKuepyta0FQj9EYDVSvW54pLZ0JqgRCZmENDa4pt3WOMBhG7S5WaoGtZvZDmy64itkC/GCw4oF/Osjf05NmjIpQ+DnxXRN4NbAfeCiAiC4CvKKUuATqBW01AdAz4tlLqp/vqvy9myqd0PZBA2x4BHlJKXa2UelZEvgusQ5v13quUqur67wG+jn6u/MQsAF8FvmmCIobQ0XsopYZE5GPAo2a/j1aDHtAC8RYR+RfgCTMGwB3AJcBmYBz4o8N94ZaZoWbeA3jnPwOw5Yo3ADrwodA3Tu8LIzUnfTbu0nZSay0jueM6hEHI8Au6BHjohziukB/2a2Pku0d1mXB/4mHhT4paSzo6M7b+rjWOaiRalWUNcbaa9EdVqkJlyJQn358GUfXF3dGT3yPsumpWS7liMjVMlL6oCiDQfqViGNU0qmrYeNWH5IrQFHfJZDySLUlKw6WatqU1S8WOYsDd1/+SP/zIl/fxy+wbqyHtBXV0ct8ppQaB86do341+VqKUegF4yYH03xczFX13wj62XQdcN0X7Y8CpU7SX2Iv0VUrdBNw0RfsLaJvn5HYFvHdf5245dljzjR/W1r/XeQqgs2sXQ8WuYoXU4z14jrBg/RDZxTrrw8j23B4TUGsmrt15yjmf1lUtqCjCi7s0zUtSDCN+PVKqCYbQ1dm8lzTESMdcNo6V2VWssDrj4TlCvhKRMymTXJFaCY2q8GiKu7zlpCx+PuCHO3P7vL6qYKr6ipKOcN7p8+jbOsqjw0WKStVMii1xl0Ap4qIFZL0W5EfKHF/tUcrCFWjsyhBLxXSUYajqTIB60vPqUzsO629m0dg0QxbLMc5lvc/ys9PPphXILMjwywd2knKFfEVRCiNkd55ENoHXEKdSrDDkB3v4YYphSDEsUXmqn1ygM2s3BSFLOtPkK4qt4755uCtObU8R+jrz+TmrWti6boAthaAm6JKOUAgj0q6wrCFeq2OUiTmsWtZELBlDnOlNMawKpvnJGJmYw8j2nK55FKqarwygFEW0e/pxUA3vrgY9VPEj6vLeGeFT8ImlYgwNFmmKO/Rq62ZNI2te2nQ4fh7LJJSCyjEqlOZO4RiL5Qjzmqce4TVPPcLZP/0ZfzO2sWZWW372gtqE3EqxUtNuqnWF+ss680RPKeSZ0RIb8zrD9rJzFuoKrUbTWJyKc6oJm25oT9XmKK08o7MW/VZN+5N2HSIwKYAi1sxLs3plCy0rmgmDiHSnDuaY7mTbqp9oXa7MM7lyTUOqZkxPuzqVU0fCrZUOaYo7e1SerZ+b5Ef6M5FNMLY7T7YhXjP9ZWK63wkZjy9/+5k5UdF3rlHVlGzlWYvlOOJtfev2+F7/cK1Gm61uTLC7GNRCpc8/q4t4Js7ZP/1Zbd/XTeNYLzef/5BYCcA/FzczcuM/ANB40dvoblzJotbMHufwhU238OerLp9W4MPukk7rExehMxHjpSe08MK2UZ43WdlTpkRI1Sc0XAhqGpUfKZan4zRkE7hxl3KuTCwVI26EarqzgUIvtPgRJzvC9vEKmZjD/QPj+zijI8+BRirOJZRSc7KA33SwQslimSb1E3mH/JCLz+gkChXLCz7FXXlaPaGcK3PuL39x0Mf41/LztfXmq/61tr5oin3/fNXltfV9BT5Uz7vqn3pJU5LB3XlSrkO759bKflQF0oSAEkBream0h+MKDe0pHFdItiRJtSQp5coEhYBKqUKiyaO3t1BLWzQbOJYF01zUgqaDFUoWywFQ/4Ab+8a1AIw+v4vvf/QuAK4aPrKZCqYqDjg/GSPpCFvHg70Kp/rv9bn8PEdH0Y0GEYsb4rieQ6q9gXKujFesUBzTvrD+0RJNcZd42qPjlHaSzQ3ku3NEfkQ5V8b1HLb3FthRrNSiCy1HDlu6wmKxvIjGK67Vn8ANHzm6x75BbeXalA5iXZmOszHvsyAZq9VJqhdaUwmoq2UZu0sV5idjeI7iN6NlljR6ZBdlicJIh7uXK7VyFfM7Gph32jwS2QQtqxcT+hUqJZ/RnWPE0x47to7Q6rk8NzYx4bia9eFY1VRmGmWFksVimU1cW9wzoelkDag612mqQIN6k14P8L53nMLOh3bzzQd28JFPvIEP//0POac1xQntOslsy/JmVr3lVcQu/lNiAy/gP/80qbYmyrknWHf/DvxI8cBgkZQrFMOJ+U4zKZCOZWGoFERWKFksltlMVdBEaL/Qn75+Fc0r2vjU5x560b6TBdX13362tt72rr/k8u8+hpeO8+gTvfSWK/xBewPxFacwSIpKyxpSrzgVXgHu/ZftEdCQnyhoe0wLhZlHoaafkHVOYYWSxXIMMVkQVIXPf/Xdz1/Me9W0BcV9z/YDutZSynW4/fvr2frNv6z1HxzTgui0/7mDG/7nMJy45cBQEM6SYJLDjRVKFssxzA1qK+HWJ6mkmg5Ic6lG6oHOclEdq0ot6a1lRlCAOjZlkhVKFsuxjrvsDNz977YH1vQ2+7HmO4vFcsxjS0XMEWygg8ViOV6wAmkuoGxIuMViOfaxAmluoBSE4bHpVLJCyWKxWOYgVlOyWCwWy6zBCiWLxWKxzAqUUsdsoIOtp2SxWCxzEKXUtJZDQURaReQuEdlkPlum2OdEEXmybsmJyF+ZbdeKyK66bZfs75gzIpRE5GMi8pQ5yf8TkQV12z4kIptFZIOIXFTXfpaIPG22fU5ExLQnROQ7pv1hEVlW1+dKczM3iciVde3Lzb6bTF/PtIsZe7M5vzOPxv2wWCyWA0VF01sOkQ8C9yilVgH3mO97nodSG5RSZyilzgDOAsaBW+t2+Ux1u1Lqjv0dcKY0pU8ppU43F/Ej4J8BRGQNcDlwCnAx8AURqc77+yJwFbDKLBeb9ncDw0qpE4DPAJ8wY7UCH0bXTzsb+HCdlP8E+katAobNGKDrsVXHv8oc02KxWGYVyqQZms5yiFwK3GzWbwbeuJ/9zweeV0ptO9gDzohQUkrl6r6m0VkzQN+AW5RSZaXUFmAzcLaIdAFZpdSDSuuj32Di5tTftO8B5xst6iLgLqXUkFJqGLgLuNhse43ZF/a80ZcC31Cah4Bmc2yLxWKZPSgd6DCd5RDpVEp1A5jPefvZ/3JgcjbE9xnL001Tmf8mM2M+JRG5TkR2AO/EaErAQmBH3W47TdtCsz65fY8+SqkKMAq07WOsNmDE7LvXsabYZrFYLLMERaSmtwDtIvJY3XJV/UgicreIPDPFcumBnJFxg/we8L91zV8EVgJnAN3Af+xvnCMWfScidwPzp9h0jVLqNqXUNcA1IvIh4H1oU5tMsb/aRzsH0edgxnoR5oe9CmDJkiVT7WKxWCxHBJ2Qddpa0IBSau1ex1Lqgr1tE5FeEelSSnUbq1HfPo7zOuDXSqneurFr6yLyZbS7Zp8cMU1JKXWBUurUKZbbJu36beAtZn0nsLhu2yJgt2lfNEX7Hn1EJAY0AUP7GGsAbZaL7WusKbZNvr4blVJrlVJrOzo69nYbLBaL5fBz9Mx3twPVILErgcnP73reziTT3ST3x5uAZ/Z3wJmKvltV9/X3gPVm/XbgchNRtxwdcPCIsWWOicg5xid0BRM3p/6mXQb8zPid7gReKyItxo75WuBOs+1esy/seaNvB64wUXjnAKNVe6rFYrHMJqJITWs5RD4OXCgim4ALzXdEZIGI1CLpRKTBbP/BpP6fNFHTTwGvBv56fwecqcmzHxeRE4EI2AZcDaCUelZEvgusAyrAe5VSoenzHuDrQAr4iVkAvgp8U0Q2ozWky81YQyLyMeBRs99HlVJDZv3vgVtE5F+AJ8wYAHcAl6ADLMaBPzrM122xWCyHjFKK6CjkvlNKDaIj6ia370Y/K6vfx9H++sn7vetAjzkjQkkp9ZZ9bLsOuG6K9seAU6doLwFv3ctYNwE3TdH+AjpMfHK7At67r3O3WCyW2cCxmtHBphmyWCyWOYiKwv3vNAexQslisVjmGkpZoWSxWCyW2YHCCiWLxWKxzBaUIgr8mT6LI4IVShaLxTLXsOY7i8ViscwmrFCyWCwWy6zA+pQsFovFMntQVlOyWCwWy6xBEVmhZLFYjiWulmUA3KC2zuRpWA4CpRRRxUbfWSyWY4B6YXS1LKt9r8cKqlmOUqjQakoWi2UOM1n4TCWMLHMH61OyWCzHDJM1oXoBdTi1pCM17nGPnadksVjmKvvTiA6X4LCa19Hk2BVKM1Lkz2KxHH32J3CsJjN30OXQo2ktcw2rKVksxzD12suR1mSqQs2a7I4CNvrOYrHMNfYlhCaHgx9O4WEF0VFA2XlKFotlDjFdrcj6geYmCo7ZkHDrU7JYjjGsoDkOMNF301kOBRF5q4g8KyKRiKzdx34Xi8gGEdksIh+sa28VkbtEZJP5bNnfMa1QsliOU6yZbS5zdIQS8AzwZuD+ve0gIi7weeB1wBrg7SKyxmz+IHCPUmoVcI/5vk+s+c5iOU6xGtUc5igFOiilngMQkX3tdjawWSn1gtn3FuBSYJ35PM/sdzNwH/D3+xrMCqXDwOOPPz4gIttm+jxmiHZgYKZPYoaZ8/fgS/t+6OyPOX/9h4EDuQdLD/Vgqjh4Z/Dk19qnuXtSRB6r+36jUurGQz2HOhYCO+q+7wRebtY7lVLdAEqpbhGZt7/BrFA6DCilOmb6HGYKEXlMKbVXW/PxwPF+D47364ejfw+UUhcfrrFE5G5g/hSbrlFK3TadIaZoUwd7PlYoWSwWy3GMUuqCQxxiJ7C47vsiYLdZ7xWRLqMldQF9+xvMBjpYLBaL5VB4FFglIstFxAMuB243224HrjTrVwL71bysULIcKofTNj1XOd7vwfF+/XCM3gMReZOI7AReAfxYRO407QtE5A4ApVQFeB9wJ/Ac8F2l1LNmiI8DF4rIJuBC833fx1TqoE1/FovFYrEcVqymZLFYLJZZgxVKFovFYpk1WKF0nCIinxKR9SLylIjcKiLNdds+ZNKFbBCRi+razxKRp822z4mZUSciCRH5jml/WESW1fW50qQY2SQiV9a1Lzf7bjJ9PdMuZuzN5tzOPBr342DYW2qV2YyILBaRe0XkOZM+5v2mfa/pYGby7+EI3gdXRJ4QkR8dj9c/q1FK2eU4XIDXAjGz/gngE2Z9DfAbIAEsB54HXLPtEbTDU4CfAK8z7X8O3GDWLwe+Y9ZbgRfMZ4tZbzHbvgtcbtZvAN5j1i8xYwtwDvDwTN+rvdw/19ybFYBn7tmamT6vaZx3F3CmWW8ENprf/JPAB037B2fL38MRvA9/A3wb+JH5flxd/2xeZvwE7DLzC/Am4Ftm/UPAh+q23Wn+8bqA9XXtbwe+VL+PWY+hZ7ZL/T5m25dMm5h9qkLxFcCd9fvU9dkAdM30PZrintXOear7NlcWdIjuhfX32fzWG2bD38MRuuZF6Dxsr6kTSsfN9c/2xZrvLAB/jH7Tg6lThiw0y84p2vfoo3R46CjQto+x2oARs+9ex5pi22xirpznXjFmpZcCDzMpHQxQTQcz038PR4LPAn8H1JdlPZ6uf1ZjMzocw0wnfYiIXANUgG9Vu02xv9pH+8H0OZixZhtz5TynREQywPeBv1JK5WTvue9m+u/hsCIirwf6lFKPi8h50+kyRducvf65gBVKxzBqP+lDjKP19cD5ytgN2HvKkJ1mfXJ7fZ+dIhIDmoAh037epD73oU0VzSISM2+HU4011XFmE3PlPF+EiMTRAulbSqkfmOa9pYOZ6b+Hw825wO+JyCVAEsiKyH9z/Fz/7Gem7Yd2mZkFuBidWr5jUvsp7OnYfYEJx+6j6OCDqmP3EtP+XvZ07H7XrLcCW9BO3Raz3mq2/S97Onb/3Kz/LnsGOjwy0/dqL/cvZu7NciYCHU6Z6fOaxnkL8A3gs5PaP8Wejv5Pzoa/hyN8L85jwqd03F3/bF1m/ATsMkM/PGxG27efNMsNdduuQUcZbcBEFJn2teiiX88D1zORESRp/qk2oyOSVtT1+WPTvhn4o7r2FWbfzaZvwrQLumDY88DTwNqZvlf7uIeXoKPXnkebRGf8nKZxzr+NNg09VffbX4L2a9wDbDKfrbPh7+EI34t6oXTcXf9sXWyaIYvFYrHMGmz0ncVisVhmDVYoWSwWi2XWYIWSxWKxWGYNVihZLBaLZdZghZLFYrFYZg1WKFksRwkR+aiIvGhCs4icV81WbbEc79iMDhbLUUBEXKXUP8/0eVgssx0rlCyWaSIiHwMGlFL/ab5fB/QCJwG/g56h7wA3KaW+JyJbgZvQZUKuF5GL0ZM1v2fWP4tOMfPro34xFsssxZrvLJbp81XgSgARcdApZHqBZcBpwJ+gyw7UU1JK/bZS6pZqg4gkgS8DbwBeydRJcy2W4xIrlCyWaaKU2goMishL0drPE2gh9L9KqUgp1QPcO6nbd6YY6iRgi1Jqk9IpVf77CJ62xTKnsOY7i+XA+Arwh2jt5ibgon3uDYW9tNv8XhbLFFhNyWI5MG5FZ1h/GbrC6C+Bt4iIIyKd7FmaYG+sB5aLyErz/e1H4kQtlrmI1ZQslgNAKeWLyL3oSqGhiHwfOB+dLXojuorr6H7GKInIVcCPRWQALdhOPcKnbrHMCWyWcIvlADABDr8G3qqU2mTaMkqpvIi0ocsPnGv8SxaL5QCxmpLFMk1EZA3wI+DWqkAy/EhEmtHF/j5mBZLFcvBYTclisVgsswYb6GCxWCyWWYMVShaLxWKZNVihZLFYLJZZgxVKFovFYpk1WKFksVgsllnD/w/RLh4ukLD9TAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "stats = stats.unstack()\n",
    "stats.sel(month=2, stats='r').plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's put everything you need in a function so it's easy for you to use to calculate the correlation for January. You just need to create the input array as per email instructions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "def months_correlations(arr):\n",
    "    # Stack the spatial dimensions together to create one spatial dimension.\n",
    "    arr = arr.stack(z=('xgrid','ygrid'))\n",
    "    # We are going to loop through the array along the time axis 1 month at a time. Considering there are only 11 months, \n",
    "    # this is not an onerous loop.\n",
    "    ntime_permonth = arr.shape[0]//11  # // is the integer division in Python 3.\n",
    "\n",
    "    # Reshape the array to have dimensions of sizes (nmonths, nyears, nspace)\n",
    "    # To get the unique values for month we use set() which is a collection of unique elements in Python. Xarray doesn't\n",
    "    # like sets for coordinates so then we convert it to list.\n",
    "    val_for_per = xr.DataArray(np.reshape(arr.values,(11,ntime_permonth,-1)),\n",
    "                               dims=('c','time','z'),\n",
    "                               coords={'c':list(set(arr.month.values)),\n",
    "                                       'time':arr.time.dt.year[0:ntime_permonth],\n",
    "                                       'z':arr.z})\n",
    "    \n",
    "    # Create an array to save all the correlations and p-values for all months:\n",
    "    stats = xr.DataArray(np.zeros((11,2,arr.shape[-1])),dims=('month','stats','z'),\n",
    "                         coords={'stats':['r','p-value'],\n",
    "                                 'month':list(range(2,13)),\n",
    "                                 'z':arr.z})\n",
    "    \n",
    "    # Since we correlate the months mn_start and mn_start+1 we need to stop the loop at the index before last.\n",
    "    for mn_start in val_for_per.c:  \n",
    "        print(\"Dealing with month:\",mn_start.values)\n",
    "    \n",
    "        if mn_start == 12:\n",
    "            break\n",
    "        # Need to drop the 'c' coordinate after selection else it conflicts with the rest\n",
    "        tmp = val_for_per.sel(c=slice(mn_start,mn_start+1)).drop('c')\n",
    "        # Now the data is in the right shape to calculate the stats:\n",
    "        out = tmp.groupby('z').reduce(corr2,dim='time')\n",
    "        out.rename(c='stats')\n",
    "        stats[mn_start-2,...] = out  # -2 as the index starts at 0 but the months start at 2\n",
    "    \n",
    "    return stats"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}