CMIP5_resample_CF.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "b9cf30f3-92ac-4bfa-ae39-3c70743f57cb",
   "metadata": {},
   "source": [
    "# Description"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a0064498-44a7-48fb-a14b-5eb554803df1",
   "metadata": {},
   "source": [
    "This notebook uses output CMIP5 pooled counterfactual histograms from the [CMIP5_create_CFdistribution.ipynb](./CMIP5_create_CFdistribution.ipynb) notebook and resamples them to recover a set of temperatures that is representative of the underlying climate distribution.\n",
    "\n",
    "We first explore the sensitivity of key quantiles in the distribution to number of resamples (with replacement) we take from the histogram. We then use the smallest number of samples where values have converged to the precision (~0.1$^{\\circ}$C) to resample the counterfactual pooled distributions over all geolocations, and output those in zarr files for use in our attribution estimates (probability ratio calculations) going forward."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "debe69d6-3bc7-492a-89cd-60dc1ff96b8d",
   "metadata": {
    "tags": []
   },
   "source": [
    "# Setup"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "2da5a6c3-82ab-452c-9116-2bcc3a0c8e3b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# import the packages we need to run this code\n",
    "import xarray as xr\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import utilities as util"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "cd56ca73-d0a0-4cdd-b85d-034f0f1abb14",
   "metadata": {},
   "outputs": [],
   "source": [
    "# add s3 feature to store data samples remotely\n",
    "import s3fs\n",
    "import zarr\n",
    "s3 = s3fs.S3FileSystem()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d4743b4f-b01d-4ffb-b3f4-ef4583acd741",
   "metadata": {},
   "source": [
    "## Find Stabilized Quantiles from Resampling of Natural Distributions\n",
    "\n",
    "We check the quantiles from the maximum temperature distribution across all geolocations and see how sensitive they are to our resampling technique. There is a certain of number of samples that will provide stability to recreating the underlying sample. We check the see whether the quantile estimate converges..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "8a0c0865-2485-4195-b5f9-4da1be7819da",
   "metadata": {},
   "outputs": [],
   "source": [
    "# define the example month we will be studying\n",
    "eg_month=7\n",
    "# define the name of the pooled distribution variable\n",
    "pooledvarname='CF_pooled'\n",
    "# define the path to saved the pooled distribution into\n",
    "pooled_loadpath='./data/cmip5/cf/'\n",
    "# load each month's pooled data\n",
    "mon_zarrname='pooledCF_mon'+str(eg_month).zfill(2)+'.zarr'\n",
    "ds_in=util.load_zarr_local(pooled_loadpath,mon_zarrname)\n",
    "# define and store the key name of the middle points of the bins, defined through xhist\n",
    "binkey=[i for i in ds_in.dims if 'bin' in i][0]\n",
    "\n",
    "# define the month names and month integer array\n",
    "monzarr='mon{:02d}.zarr'\n",
    "months=np.arange(1,13)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "b8c1db04-de1b-45d0-b48d-00ea6972e63c",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1000\n",
      "10000\n",
      "10400\n",
      "100000\n",
      "200000\n",
      "300000\n",
      "1000000\n",
      "10000000\n"
     ]
    }
   ],
   "source": [
    "# define the array giving the number of (re)samples to take\n",
    "nsamps=[1e3,1e4,10_400,1e5,2e5,3e5,1e6,1e7]\n",
    "# define the quantiles to check\n",
    "qcheck=[0.5,0.75,0.9,0.91,0.92,0.93,0.94,0.95,0.96,0.97,0.98,0.99]\n",
    "\n",
    "# initialize the array we will store the results in\n",
    "qvals=np.zeros((len(nsamps),len(qcheck)),dtype='float64')\n",
    "\n",
    "# group all of the values together across the range of latitudes/longitudes to resample\n",
    "array_to_resample=ds_in[pooledvarname].sum({'lat','lon'}).squeeze().load()\n",
    "\n",
    "# loop over the different sets of resampling N-values\n",
    "for i,n in enumerate(nsamps):\n",
    "    \n",
    "    # tell the user how many samples we are currently testing\n",
    "    print(int(n))\n",
    "    \n",
    "    # resample the pooled distribution with the new N\n",
    "    poolcheck=util.randomly_sample_histogram(array_to_resample.values,N=int(n),\n",
    "                                        bins=array_to_resample[binkey])\n",
    "\n",
    "    # calculate the quantiles from the resampled set of values\n",
    "    qvals[i,:]=poolcheck.quantile(qcheck)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "9462589a-3add-4cea-8af8-6e07498f030c",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAl4AAAGJCAYAAABb8j2MAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8/fFQqAAAACXBIWXMAAAsTAAALEwEAmpwYAABa3ElEQVR4nO3dd5wrV3338c9vRtKWu7vX3ZhibBwDNmAgMcRAQg+9hWKKKQbyEJ4HAiSQAqQYUoAQIIRAgBCwQwsOJQZCjbFNcwDbgI0x3TY2GHffu3u3SJr5PX+ckXaklbSzuyp7td/366WXNP1Io/LVOWdmzN0RERERkcGLRl0AERERkZ1CwUtERERkSBS8RERERIZEwUtERERkSBS8RERERIZEwUtERERkSBS8tgkzO8vMbjKziS7TZ81sn5md3odtHWVmbmanbmLZl5nZEzuMP83Mtu25Sczszmb2JTPbmz33J6wz/73M7GNmdq2ZrZjZFWb2djO79ZCKvC4ze4KZ/VGH8Q/MnuMDc+PONbNzh1i8dWXvw9PM7A4dpl3Rj/f6Jsp0mpk9eNjb7VIOzw0fkI379Q7znmtmX93kdhrvlbqZ3bHD9KtHtB9Ozcr1a8Pe9kaZ2avM7OfZa/idHvPNmdlrzOz72Xf5zWZ2iZm9y8wOG2KRNyTbD6eNuhzjRMFr+zgDOBB4TJfpTwams/m26hrgPsB/b2LZlwFrghfwnmyd29WbgTsAJxPKeV63Gc3sWcD5wMHAS4HfAV4HPAL4tpnddeClLeYJwJrgBVxEeI4XDbU0G3cU8FeE/dLud4G/Hmppgr8CRh68WPt5OoBQtjXBq09i4LUDWvfYMrN7A38L/Adwf+BZXeaLgf8B/i/wb8DjgOcAHwbuC2ybP3QyeKVRF0CaPg3cCDwb+FiH6c8Gfg6cu9kNmJkBZXdfAf53s+vpxN2vBq7u5zr77Djgy+7+uV4zmdmdgH8F/gs42d3TbNKXzeyjwDeAM83sBHevD7LAm+Xue+nz/h02d//2qMswSiP4PH0BONnMXufu3x3idkfGzCay78KtOC67f6e7/6zHfA8A7gU8wd3Pyo3/JPB3ZqZKkJ3E3XXbJjfgn4EqcEjb+COBFPibbPhhwGcINVeLwPeAlwNx23JXAB8Angf8AKgRahKOAhw4NTfvvYCPEr7sl4AfAn8HTLWtz9tup2fTTgtvp5btz2XP6ZfASrbOPwQsN88Ds/U8Lpv3BuD6rNwHFHjNysDfZGWrZvd/QwiY+fW33Hqs71+y1+mILtNPztbxxNw4B05rm29Tr3E237nAV4GHEmqtGvv4Cbl5Tu/wvK5oe84PbFvnuW3bOSR7vr/I9s8PgBcUfK+uuyxwK0INbWP/X0P4g3FYt/3SKHO2H0/PrevUbPp9gTOBeeBa4JXZ9EcA3wb2Ad8CfqOtLOt+ZrqU57Tc9AcAZ2fb3gd8HrjrOq/Tk7P13DY37k3ZuN/LjfudbNzx7Z+n3Hup/XZq0fdLj/I19sNjs335qbbpV3fZD0e1zdcsb9vr+TfZ63xl9pr9d7b/D8v24x7gKuBP25ZtbOf+hD9BC4Q/pm9n7edlGngDcDnhO+By4NVA1OF5PpHwx+p64JZ1Xpt7E2qpFrKynw3cu+0z1fX90uV98JsF9slGv9+fRfguWQK+AhwL7ALelb1m12bvuVKH1+NJhO+Sm4G9wAeBgzvsx/bvt7sTQuPN2Xa/Bvx22zz3Ar6YlWER+BnwjiLfL+N+G3kBdMvtjPBGdeBFbeNflY2/Yzb8wuyD+EjgQcAfE34MXt+23BWEL9PvAU8HHgIcQ+dQ8CTgzwlNnQ8A/h/wK+A/cvPcM/sy+BxwUnY7Jpt2GrkvXkIz9leyL6yXZ18mb822+3e5+RpfAJcDb8vm+4Psw3xGgdfsQ0Cd0EzyMEJzTA34UDZ9LivndYQv/ZOAk3qs70fA+T2m7wKS/BdIly+mTb3G2XznZq/zpcAzCaHii9nz/LVsnmOy53Ndbl/cs+01fWDbOs/NDc8Rvqx/Dvwfwo/2G7Pn9gfrvOaFls3K/CPgFMIP6FOAd2avzVz2/D3b343nMJd7756eW9ep2bw/Bv4i2+a7snFvAC4Bnpa9tt8n/JhXcsuv+5nJtu/A+3LluW027dHZ638W8Pjs9nXCD8/terxWBxP+ND07N+7bhB+iD+XGvQ64Njd8GqvBa4Lwh8kJQb1RtkOLvl96lK/xXnkooRnMyX0+2HrwupLwPn004Q/gXsL3x9cIn4X8fnxUh+38HPgHwmf7zwnBKl+eEuF75kZCN4iHEELXMvCmDs/zF4Rm3EfQI5gCJxC+gy4khKYnEQL9EnD3bJ7js/3h2f5pvl86rO8OhO+l72TzHthj2xv5fv85oVvEEwh/Cn8JXEwIq/9ACPR/nZXx/3V4Pa4ivN8fQfgczgPndNiPp+WGf53wvf7V7LV5FCGErZD94QFmgJuyff3YbHunAu9e7zt9J9xGXgDd2nZI+PL8Rtu4y4Cvd5nfsi+fVxN+BPL/8q4gfMHfqm2Zo2gLBV3W+UzCj8bBbev8QIdlTqM1eD2m0zayL70Vslq93BfAGW3z/TPhy9M6lTGb567tXwrZ+D/Pxp+QG9fyA9JjnUvAh9eZ51fAf+eGCwWvDbzG5xK+pI/NjTuMEGxelRt3OnB1h3U3XtMHtq3z3NzwX2Sv77Fty/4rodax1KncG1mWUFPwkh7raZTzoR2mXUHnH/y/zI0rEYJnDTg6N/5x2bwP2MRnxslqltuW+Qlwdtu4uez5/uM675fvAu/LHh+U7cc3Adfk5vlfWv/knEbr56nxfvq9Dusv9H5Zbx8Qao9/Cnyp2+eGjQevH9Fa0/LmbPyfd9iP7+uwnXe2rfPV2fNq/Al9Vjbf/TvMVwUOa3uen+j1euSW/yhwC7la92x/3wR8PDfu9zq9Hl3W+XuEz4QTPvOXEv6w3LrHMut9v98E7M6Ne0m2/ve0recicoEq93p8rm2+U7LxD2nbj6flhs8m/Cbl/9jE2bj/yoZPpO07WLfVm9qVt59/B+7dOMIo67x552w82bgjsiNhriR8udQIVfoHEL5w8/7X3X+13kazI27eYGY/JQSjGvB+wgf/2E08j/sTvlw+3Db+A0CFtR3x2zv6X0L4p3/4OttorLN9GxBqlQbBCM9tYwtt7DX+sbv/uDHg7tcRfpyO3HSpWz2C0F/tcjMrNW6E5rODCf/mt7rst4A/NrOXmtndsj6GW/XZxgMPfex+AvzI3S/PzfOD7P52jREb/My0MLNjCTWMH2x7vouE2ob791oeOIfVDvsPJDSvvRm4lZkdZ2azwG8AX1pnPb1s+f3i7jVCgHqQmT10C2XJ+6K39oVs7JvP57bb2I+3Y60z24b/g1Cbfu9s+BGEWrWvt+2bLxCC5Elty3+iYLnvD3za3W/JlXMvoWZnU98r7v4ewnN8JvDu7Hm8ArjUzO7SmG+D79Xz3X1PbnjN65sbX+T1/U/Cd1vHA6XMbIrw/P8TSHOvtxGaZRufhR8Tguu7zOyZZtZp2zuWgtf28wGypols+NmEH+mPAGSdMD9JqFH6G8IX+r0IR9YATLat75qC230foYr7nwjV0/cCXtRlnUUcBNzkazuv/io3Pe+mtuHGcr223VhH+3Psto0iribULnRkZrsI/Zt+sYl1b+Q1bn89ILwmm9kXnRxG+JKstd3+M5t+cB+WfSrhvfonhOaPX5jZX26xI/HNbcPVLuMge6028Zlp1/ix+zfWPufH0Pu1ghCojsxOm/Eg4Dx3/wWhufZBhNeyRAhom9Wv98sHCTUxf7vejAV12zedxncq67Vdhm+T3R8G3J61++Wb2fT2fVP0+/CgLvP+inD0+aa4+83u/kF3f6G7H0doIpwDXgObeq/29fV198bn6TYd5oXwusSEWu/21/zFwIFmFmVh8EGEps93AD83s++Z2ZO6rHdH0VGN24y7/8LM/gd4ppm9luzHy90bH6RjCNW4z3L3Zk2PmT222yrX26aZTRL6rJzm7m/Njb/bJp8GhB+Cg8yskn2YG26V3d+4hXXnt9FY50/7tI2zgeeb2RHu3umL99GEPyz501GsEGrx8lq+8Af0Gm/FjYQakZd2mf7DrS6b1bq8CHhRdrTocwg/MNcTOuYPy0Y/M+0a76NXEv7Vt6t2GJd3HuHP1IOz2zuz8V/Khq8EfpGvsRoVd0/N7C+Aj5vZ4zvMspzd93y/99HhhCCYH4bVPz43EvqHntxl+Svahtf9PszcxOr3SN6t6BxyN8XdzzKz77JaS7zV9+pGtbQomFmFECy7/bG8hfBefju5Vpg8z44Ed/fvAE/KasROJHx+zjSzu7v79/pR+P2Vgtf2dAbhn+frCLUr+Tf4dHZfa4wwszKhbX6zJgj/Ympt40/tMO8KMFVgnecROoU+hfBcGk4h/FD143QHjfDzNFr/oTdeiy9vYp1vBZ4LvM3M8qeTwMwOInSm/RWtTRZXEvqb5T26bXgjr3FRRfdFJ58jdKb9eRaQBrqsu/8QeJWZvZDV16pRq7nZ51DURj4z1Q7l+SHhB/wu7v76jW7c3feY2bcJ79PjWW1S/BIhhF3N+s2Mw3qtcPdPmNm3CJ2y22snr8zu70rov0X2w/qwARXnZFpfm6cRfvgbNVqfI3R8X3D3H9A/5wGPNrNZd5+HcBJrQkfxcze6MjM7JCvjctv4XYQmwIuzUYP4fu/lZOC9ueGnEPb5+Z1mdvd9ZvYVwlGNF+W/H7vJmpL/Nwv0jyOcgkPBS7adTxCO/vlDQs1C/txTlxG+/P7WzBLCB/QPt7Kx7Ifhf4GXm9k1hA7Dz6NzdfP3gd82s8cQAsgN7n5Fh/k+Szjq5Z1mdijhX+ujCB1MX+fuN2ylzFm5LzWzDwOnZV/+Xyf0TfgLQgf5i3uuoPM6LzOz3yccBHC2mb2T0ORwZ0KT2W2Bh7V9gf4H8Odm9mpCoPxtwlGk+fVu5DUu6vuEWsX/C1wALLv7JQWXfQuhNvUrZvYWQrjYRXiev+3unWo7Ci9rZrsJtUMfZPVUJo8n/Jv+QraeHxGOvHuemd1EdsqRxg9dH23kM/N9wg/u5whNLr9091+a2YuAs7IagTMJ++9wwuktfu7ub16nDF8i/BG5zt0bNTjnEppuDiYE/l6uJdTuPM3MLiYcVXa5u/ej5riTV7O6n/K+RahdfmPWLLZCODq14xU3+uBRZvbGrCz3Jhy1/O/u/qNs+gcJf5TONrM3EQ5kqBBqjh5HOHJxcRPb/WtCc9/ZZvYGQk3ZnxKC0WZONPtA4F8sXAXgK4Sao9sT/sAcROjzBwP4fl/HXczsfYTvsDsS/sCe5+5n91jmjwh/aj9vZv9G+H48hHC0Y+zuf5b9PryAcHTl5YTvh5cQjprsGOp2lFH37tet843ww+/AWzpMuwch1CwS/i2/lg5H19D9CMSjWHuqg6MIYWmeEPb+mVBr03503J0JXxyL2bTTs/GnsfaopsZ5vK4h1CT8iO7n8Xpo27Kntj+fLq9T4zxeVxK+pK4kdx6v3HyFjmrMzX8SIQBfT/iH7YTz0BzXYd5Jwg/nNdnr9xHCj8RmX+Nzga922M4VtB5htotw8MLN2TquaHtN29d5btv6DiSEqMuz/XNdtm9fVuD16bks4Yf4XYTAvUD4I/Et4Blt6/n97HWt58vc4bk23g+/1rb8mteKDkcAUvwzcz/CKQSWWXs0130I5yG7OZt+BeEH6z4FXq9HZutrP3XId9vL0OPz9ARCMKzl31tF3y9dytV4r3Q6svQccp/x3Pi7ZNtcIJzO4I+6lNdpO0K06H6k9TxeZ2XbuonO5/GazLb/A0IQvCl7r53G6hG2XZ9nj9fmN+lxHq9snkJHNRL+sL2OcFBK40jc6wnn63pw27z3YJPf792eJ21HQOfme2I27RbC99KHWHseyZbPQTbuOMJ7/7rsNb+a0DftUdn0OxG+By8nfFYaz3Xd85jthJtlL5KI9GBmf0eo8Xqyu//XiIsjIrJpFq7jeg7wO+7eqd+iDJCOahQp5tWEJqYPm9kDRl0YERHZP6mPl0gBHqqGnzHqcoiIyP5NTY0iIiIiQ6KmRhEREZEhUfASERERGZL9oo/XIYcc4kcdddSoiyEiIiKyrgsvvPAGdz+007T9IngdddRRXHDBBaMuhoiIiMi6soucd6SmRhEREZEhUfASERERGRIFLxEREZEhUfASERERGRIFLxEREZEhUfASERERGRIFLxEREZEhUfASERERGRIFLxEREZEhUfASERERGRIFLxEREZEh2S+u1SgyUO5w7uvhB/8NUQSlSShNhPu4kg03xk10mDbRNj17HE/0nhbr4ycistPom1/kZ+fAea8f/nYt7hzo1gtta0JffnqnQNg2nF+3wp+IyFDpW1fkG+8azXY9gdq+cBsV61bDV6QWrz3Q9agd7BUmoxKYje41GEfuUN0HK/OwsjfcL++B2mKYJrLTHfswKE+OZNMKXrKz3fhT+NHnswGDU/4TJmahvgz1ldxtOdySam7aMtTbhlumr/SYtgyejvSpA6EMtcVwGxWLuoSybrV4W6kdbA998eiedy/1ai40NYJTdr8yDyt7Vh83x7fdb4f3l8h29YqfKHiJjMS33gNkNQDHPgyO/Z3hbTuptwWz9rDWNpysN619eo/QV69CfWl7/Dh7GspSXxp1SUREBk7BSwYjTSCtQ1KDtBZCRlrLjctNS5PcfG3ztKwjm7f5uLHeeus2Oq63y/av/d5qmX/zBcN9jeISxDMwMTPc7eY1wl8zmK1Xi9ch2G0o9HVYNq2P7vmPs9IUTM6FGtyJ7L4yo2ZdEQhdJUa16ZFtWQL38I+/PXh0CigbCh7toaVX8Kn3WLbX9vMBqW377Gf9SA4+Fu7w4FGXYvga4W+UknouvK1XU9cr9LUHxnVqEOvLoZ/ddhSVQlianMtCUyNAzebC1CxM7O4wLnscl0f9LESkAwWvhevhsrM2EFrqfQg+beuQETN4yF+GU0nI8MWlcKvsGnVJREQGTsFrz1Xw3y8fdSnGU1QO/7qjcvhhbQ7Hucel1fs185Vy03PTolLbfNk6O26r0zrapu++LcwcNupXS0REdgAFr21RHW9dwkV5/dCwbkBpW8ea0NJl2TUBqdf224NPefseLSYiIjJCCl7Th8BvPHcTtSq9gs9GanXKauISERHZIRS85o6Ax/7jqEshIiIiO4CqWkRERESGRMFLREREZEgUvERERESGRMFLREREZEgUvERERESGRMFLREREZEgUvERERESGRMFLREREZEh0AlUR4EPf+Dlf/P6vmCjF7JoosWsi3M9MlNhVyT3ObjPZPDMTJaYrJSol/YcREZH1KXjJjnfhlTfxqk9csqV1VOKoNaw1A1rMrkqpQ5hbG+Aay+yqxJRiBTkRkXGk4CU73nu/esWW11FNUqqLKTcv1rZeIGCiFK0NcI3HldVx020Br1Ot3K5KiSiyvpRLRES2RsFLdrRf3rLE5y79VXP4b55wVyqliH0rdfat1FlYSZqP91Xr7FtJWGgMr9TD42pCknpfy7VST1mpV7lxX7Uv65sqx2sD3CZr5aYrMWYKciIim6HgJTva+//3ymZoOukOB/HMk26/4XW4Oyv1tBnIFlbqLFbbA1rSGtaycYvV/Lgwz0K1jvc3x7FUS1iqJdywsPV1mcGuSghgqzVs8dratspqkGsPeDO5+SbLkYKciOwYCl4yVtydfdWE+eUa88t15pdr7F2qs7c5XG9O27tc40uXXddc9rn3O3pT2zQzJssxk+WYQ2Ym+vIclmpJaxhbqbNY7R7gWmvg6izma+aqyZbL1Fo+WMi2dd38ypbXFxlZ37YuAa5RC5cLbdNtNXP5+SZKOyPIuTvLtZT5lfB+Xsje3wsrq+/1hZXwfl9YqbN3uc5SNcH7nepF9kP/+LR7snuqPJJtK3jJttEIHI1wtGepngtQrYGpGapy0/YuhR+YzbT63fbAKR563OH9f1KbYGZMV8LRksxufX1p6izWkrawFkJdCHPda+WaTavV1RC4VOtvkEud5j7uh1JkqwGt5YjUeE1tW2WbHsRQTdLm+3phZTVUza+0jqv3uYlbZKeoJenItj3w4GVmMXAB8At3f4yZHQR8BDgKuAI42d1vHnQ5ZLAazW17sxqmXoEpH5xah+t97ytV1MseekfiMe2AHkXGTBY0+hEtk9SzINYa2PJNq/uq7eOS1tCX6y9Xrff3C7CeOnuWauxZ6s+BDiIi/TSMGq+XApcBc9nwnwFnu/vrzezPsuE/HUI5pIeVetIzMO1dzk9bG5jml2vUku3x73uyHDE3WWZ2ssRsdt8YnpsqMztRapl29CG7OPbwPlQt7RBxZMxNlpmb7E81fS1JQ9Notb22rTXULVY7BLg1tXL1bfM+HIZKKWJuMoTq2clyCNiT2fu7Ma45vZQd4TrqUouM3uzk6Br8BrplM7st8Gjgb4E/ykY/Hnhg9vgM4FwUvLakWk/XBKa9HWqWugWmvcv9r3XYrIlSxOxkmbnJUsfgNNu8D4/nptZOK2/T5iPprBxH7J6O2D3dnyC3Uk9W+7hVu9fKLawk1EfY3NBLHFvzfZ0PVY33fiNgTZTiURdVRDaocPAys13uvm+D6/9H4E9o7alyuLtfA+Du15jZYRtc51ipN/tydG+KW9tU1zptubY9fjwqcdQSiroFprksMM12mKYfEtmqiVLMRCnmwF2VURdFRGSNdYOXmd0XeA8wAxxpZncHft/d/986yz0GuM7dLzSzB260YGb2AuAFAEceeeRGFx+KJHUWCgamvW3Nd42j7frdUXmzSpH1DExzuRqmToFpbrLMZFmhSUREpJciNV5vAR4OfBLA3b9rZvcvsNz9gMeZ2aOASWDOzD4AXGtmR2S1XUcA13Va2N3fDbwb4MQTTxxYp42lasIVN+5bG4o6NMnlA9P8cq3vh+lvVmR0DExzbU1y7YFprjl/WedSEhERGYJCTY3uflXbj/K6icPdXwm8EiCr8XqFuz/TzN4IPAd4fXZ/1saK3F+X/GIPJ7/r/JFt3wxmJnKdvzv0YeoVmGYndSZxERGR/UWR4HVV1tzoZlYBXkI4SnGzXg+caWbPB34OPGUL69qyrR7Z0H6EXLfANJdrksvPq+voiYiI7BxFUscLgbcCtwGuBr4AvGgjG3H3cwlHL+LuNwIP2cjyg3TgdIVjD5vJBabWfkuzXQJT4yijcT33k4iIiPTfusHL3W8AThlCWUbiVrsn+eIfPWDUxZARu3xxhW/sWWAyitgVR+yKY2ZKETNxzEwcsasUMR2pH5yIiGxN1+BlZm8DunZqd/eXDKREIkP286UVHnbBD5lf55xOEbArjpgpZWEsF8pm4jhMawa2bFw+vGXLNuabjExBTkRkh+lV43XB0EohMkKv/ekv1w1dACkwn6SF5i0iNnKhrDWwTecC3kxbYJvJTZvOhb0JnZJcRGTb6xq83P2MYRZEZBS+fvMCn75+T3P4UYfsZilNWUxSFpKEhXrKQpKyL0lY7vN1JBOHPfWEPfUE2Pp1BctmLTVw+Vq39tDWWjvXueaurP6LIiJ916up8R/d/WVm9ik6NDm6++MGWjKRNqk71dRZSVOq7iynTjVNWUk9u6VUU2c5u19JU1byy7RMC+POu3m+uf4nHn4g7zj+9l23X0udfUnCQhbK9tVT9jUCWpKyUE+y4fC423z7skBX9f4GuZo7N9cTbu5TkJuIrGdgazze1V4716WJNVazqohIz6bG92f3/zCMgsj2lfjawLLSKfR4jwCUzVf1tGWe1vWthqVOQanW56CSNxUZr77DET3nKUfGAVGJA/pzSUGqadoxsLWEtHou2GUhbjXAtc5f7/PLE/ZNwk19urrCVGRrD1rINbHu6hTiOtTO7crCXKQgt2O4Oynhz0/dw3dB3Z1a9r2QOC3jmvO0z98cB969C7PsAI8/7ECmRnRd315NjRdmD+/h7m/NTzOzlwLnDbJgO51nXxIhfORCjhcLQMvJ2iC0JgAlYX0toSc/fxaAkh3w/fSHR92K20wO99p+lSjioCjioPLWr1Xv2b5u1qjlat32dQhtLQEv15yaH9fvK4Aupc5SWueGrVfGAWT929Y2leYDW6fQlh+u7LDmVAfqaWsQSbLhtcEl/OlqhJVa2+PmNF8biNYEoNxw3Vlnen4czWGRfnrIwXPbL3jlPIdwHq+8UzuMGwuNH7DONTj5kJMFFm+t2cmHomaNja8TlNLOAUhfNWtVzJiIjEoUMRkZlezxRGRMWJRNMyajiEpkTESr4yaiiIrlp4Vxt52scP8DZ0b91LbEzJiMjck44pBiF6Toyd1ZajStdghsraFt/SbWfUna9/fzYhL64l1Hvc9rFhEZnF59vJ4OPAM42sw+mZs0C9w46IINy3fnF3nGd3/WDED97nczDgyykNMWbKxDsGkPPbmgNBG1P86WtbVBaaItHFUsbEPNS8NhZkzHxnQccWgfKgJTd5Ya4axDn7j2gNe1di4LfYt9OrJU9h9GOICkFFm4t+w+YvVxc9zqcONWzs0Xm7HDKjulzeQIjwLv9df468A1wCHAm3Lj54GLB1moYYqAG2vb8x9zyciFk+61OS1hxloDS3vomWgGpi5BqbE+i5phq2TofFOyJZEZu0oxu0oxh7P1TnKJe8cjTxd69IHr1MRa34F/tEIIsZYQ0ggwzfAS5QJL27jVQEOHENQaeFaXofN625aPjZbt5OfRny4ZF736eF0JXAncZ3jFGb5O5z5q1ML0qn3pGHry8/YMSp3WtzYo6Sgwkc5iM2ZLMbOlGCZGXRoRkeLW7QxiZk8E3gAcRqjtNcDdfW7AZRuKO0xNcMn97tLS7KXaHRERERmEIr1w/x54rLtfNujCjEIpMg6t9On8ACIiIiI9FOlddu24hi4RERGRYSpS43WBmX0E+C9gpTHS3T8+qEKJiIiIjKMiwWsOWAQelhvngIKXjI2lpavYs+cizEpEUYUomuhwax9fUX9AERHZkHWDl7s/dxgFERmV+flLufCip5Mk+za8bCOMmVWIowmiuC2sWSU3bm1wyw/H2b1FlebjtfO1LmNWUvgTEdmPFDmqcRJ4PnAXYLIx3t2fN8ByiQzFysq1fPfiF2wqdAGkaZU0rQL9uCz1ZkQ9QlylQ4jLhUOrrB0XVdrC4kSPENio9YtH8sxFRPZHRZoa3w/8AHg48FrgFGBsOtvv2/czvnfpS0ddDBmRavV6qtXrAYjjGQ4++P6k6UoWqFZab8kKqefHV0dceoCUNF0iTZdGVgKzci6QdWqmrfQIcQWbddvDYpwLl6YmXxHZfxQJXr/m7k8xs8e7+xlm9iHg84Mu2LCk6TILC98fdTFkxMxi7na3t3PwQb9VeBn3lDSttQSxNWGt0PgqyYaWWZ3mPvqrLrjXSJLapmsN+6F76OvdVNsSBOPJDuEwG5dNa1k2G7cdavzcU5JkH/X6PPX6Akmy0Hxcr89Tz4bTZAldBVYE7nD0SymVRnON3iLBq9GCcouZ3RX4FXDUwEokMnTGne742g2FLgCziDieII5Hd+r0NK3jXm0Ja0m3WrqkPdB1DndJuoK3rKv3MmyDH/LVsswPfduhxi8X5OLVwLY6brIlzLUGvvy4MB+QBaYsSNXnmwEqaQtTIWjtYzvsB5H9xVG3/31g+wavd5vZgcBfAJ8klPQvB1qqIZqePpp73+uT688oY6tcOZjJiVuNuhibEkUloEQcT49k++6Oe61DIKt2CXHr1Oo1wmFLWMzdvEqa5NdVxX20Tb6rNX4LI+rnJyL7kyJHNb4ne3gecIfBFmf44niK2dm7jLoYIvslM8Ms9NUa1b/HtU2+6we8pH1aspyNX14TBFfHLTeXSZLVcdulpimOd1GKZ4hLs5RKM5TiGUqlWeJSuA/TprFC580WGW9xvGtk2y5yVGPH2i13f23/iyMisjGjbPJdrfHLwlySD2jL2bjlLOAt50Lfcm7eamvgS5YBOoSm1cel0mzzFsczlEq7tkVfMxFZX5GmxnyP2UngMYzRUY0iIpuVr/ErMTvq4ojIfqBIU+Ob8sNm9g+Evl4iIiIisgGbaeyfZgz7eom4b4++OiIiMr6K9PG6hNXeozFwKOFEqiJjoV6vc+aZZ/LjH/+YKIqoVCqUy+WOt27Tio4vlUpEkTo3i4jsVEX6eD0m97gOXOvb4ayNIn3yzW9+kx/96EcAJEnC0tISS0uDOxN8qVTqS4jrNq1U0vUbRUS2qyLB61jg+OzxBe7+iwGWZ+huuukmPvGJTxBFEWa2ofvNLLPZ+2FsYydaWFjgvPPOG+o26/U69Xp9oOGuXyGu2/g4jhXuREQ2oWvwMrPbAWcRTgV9IWDAk8xsCXg88KzcOb72WysrK1x11VWjLsa2MK6Bste2vvOd77CysgLAwQcfzAte8AKSJKFWq1Gr1ahWq83H+Vu38estU68Pp7K4sc3FxcWBrN/M+hbieoU7EZFx06vG6+3AP7n76fmRZvZs4PxscL8PXupQvcrdSZKEJElGXZSRePjDH87ExGDPBZWmKfV6vXBQ20y4G8b+c3eq1SrVapV9+wZzjcYoigYe7nZqTa+IjE6v4HXn9tAF4O7/bmZ/B/z6wEo1RAcffDDPfe5zSdMUd+95X2Serd4PYxvt21L4hGOOOYZjjz124NtpdN6vVCoD20aaphsKapsJd2maDqz8+eexsrLSrJEchEa4K5VKzVv7cNFxG1muUetalLtTq9War8fKygrLy8stw0XGVaujvbySyHbx8pe/nJmZ7Xetxo5/Bc0sApbc/brBFGm4JiYmuP3tbz/qYoxUI3ztr0Fyq+uYmZnh4Q9/+Nj0WYqiiImJiYHW3uWbY7ca4rqNH8YfgmGEu07MrGdAA9YEJ/1BEhkPvYLXp8zsX4GXufs+ADPbBbwF+MwwCifDEc6+vXM72MvGxXFMHMdMTk4OZP2NZu9Bh7tRadRg1Wq6rLbITtMreP0J8DrgSjO7knAur9sDZwCvGkLZRGSHytcITU1NDWQbjXDX6HPXONq0cWsfV2Seosttpqm2VCo1azIbt8nJyQ2Nq1QqY1OzK7IVo/wcdA1e7l4DXmFmfwH8GuGoxp+4+2AOkxIRGaJ8uBtUzV03jcDXLZy5+5oA1WiCFJH9W5FrNS4BlwyhLCIj4Ylz83/+kKUf3ISVIqKJGKvE2EQcHk/ERBMlbKJ9XKf5wr2VYyxSzYJ01miqHfRRtCKy/egvlOx4+755DYvfuR4AJyFd6E+/m/ZQZpW1Aa05rfm4RFSJsCzorQa5jR0FJyIi25OCl+xo6WKNvV+8ciDr9mqCVxPS+T6szGgJbjYRE7UHu4lS50DXoVaOkoKciMgoFApeZnYbQsf65vzu/uVBFWqY0mpC/frBXbpFtrd9/3sN6WI4m3x80CSHvuBueC3FVxLSlQRfycJT9rg5biUhXanj1TTct0+v9fk8Vw6+kpCs9OnkqBFYpUNQ61Art6bJNauRa5mvpCNiRUSKWDd4mdkbgKcC3wca3/oOjEXwqt+wxHVv+/aoiyHbwO5HHk3pgP50svbUWwJbPpSl1QTPwlpLoKsma8dl99T7HORS8OU6yXKfLmEUW+c+b5V8eCt1bmbtFPZiBTkRGU9FaryeANzJ3Yd7hkGRIaocvZupux7ct/VZZNhkiWiyP635nji+Us9CW4dA1ynktYyrtwa5pM8n40w81Bwu9inIlSzX562t5q29ybVrkFutldOBDiKyXRT5VfgZUAbGMnhZOaJ8612jLoaMUDw3wQFPOGZb93my2LDpMtF0uS/r83ras5ata8jrMj9pn4Nc3UnrddhXpx+Nq1aOOvd5q4RmUitHobk0u7f17tvGkR8eYshrBvLlhHQ5q0XN3afLCb6cNYsvh/2HToAvwkFPuSPR1Gi6uXfdqpm9jfARXQS+Y2Znkwtf7v6SwRdv8MqHTnP4S8bispMihVkpIi5FsGvrQc7dQ1Bq9nnr0pSaC25ra+Va5+93OPBaGvrdLdT6EuR6im3d8EaPENfoL7caoBJ8OdRYtoSr5QH0JRTZITwZ3WenV9y7ILu/EPhk2zT9ZxIRIDsDdNmIy/258Le7h6DUtTZuNaytCXIdauV82LU8ieNJFiBFRNr0OnP9GQBm9lJ3f2t+mpm9dNAFE5GdycywSgyVGGa2vr5mkGtvSs2OPvVaitdb76l3Ht9y32WZoTJCs+lkqXkfTcahf+FE/j6bpxKzjVvURYamX/1vN6PIlp8DvLVt3Kkdxonsl9I04UvvfSc/vehblMplKpPTlCcnqUxNUZmcCvdTU5Qnp9eOm1h93BhfnpjEdMHxbSMf5OLZwW7L3UONVyOU9QhvHcNdY5p7dnBGa6hqBKhGqLKKzscmsr/p1cfr6cAzgKPNLN/UOAvcOOiCiQzL9875H777xc/2b4VmlCfagtvkVBbmpsPjtgAXHk/n5lsdLk1M6Md1P2FmUAp9vBju5R9FZD/Rq8br68A1wCHAm3Lj54GLB1moYVqa38sPz/9q+FccRURRhDVu2TizMJ7IwnTrPD0stzquMc/adcab3p7018riIl/7yPv7u1J3astL1JaX2NeH1ZlFXWrgcrVszYA33RLmmgGvMd/UFKVyRe8lEZER6dXH60rgSuA+wyvO8M3feANn/9s7Rl2MwooGPRqhzroFwTi37GoYjJrrjHPLtk/vHhY7ladXwBxEON1IGP7uFz/L4p5bAJg5+BCe9MrXUK9WqS0vUV1eorqU3bLHteXFlnG15eXs8WJzXH2lv2decU+pLi1SXVrsy/osilpDWjOgTWaPp9fU1q3ONx3ma4a9aUrl/pziQkRkJ+jV1PhVd/8tM5un9ZggA9zd5wZeuiFw378O0PQ0xUlJdcBU393/GadyyO1uv+X1pGmyGsiWllpDXDPMLYbxLaFu7Xy1pSXqtWofnt0qT1NW9u1jZV8/6uMgiku9a9naA1x7c2vWBNsYjku6hKyIjK9eNV6/ld0PuDvqaE3NznLCQx8RjnxK03CplzTB3UnTNBufhPGeZvOkq9PTNDfem8Ory7Yu0zK9yzxpl3XK4Bxxxztz5/s9oC/riqKYieldTEz358S8Sb0eglxWq9aplq0Z3Jq1cYu5+XI1dEuLJPU+nV0+kyZ1lhfmWV7ox9XAIS6X19ayTU1RmZgM4ayttq7RhNran261STaK4r6US0SkH4pcq/G1hOsynu/u/fmLvI3MHXIYv/N/XjzqYhTSGviSlrCYpgn0CG1pmq5OLxj0Ok/fxLJtITOUpXs4bZnets4iYXj9ZVuX2XXggTzsBX+wbfs9xaUS8cwMkzN9OLcCkNRr69TGNULcYi7ULWePF3MhLtynSX+rX5NajaVajaX5vX1ZX6kysbY2LquJK09MUKpMUCqXm/dxpUKpUqFUbr2PyxVKlTBfXC5Tzu5L2fxRXNq276HtxN1JkzppPSFNE5J6nTRJmuOSpDG8Oi5N6iSN4STNxof5kty6GvOmSW69adKcN03qJPWk+f0pO9eDn/tCJqanR7LtInX6VxCObnxb1uz4FeDL7n5Wr4XMbJIQ2Cay7XzU3f/KzE4D/g9wfTbrq9z9M5sr/s5iUYQBURwTruIksnFxqczUbJmp2a33FnB3klptba3bcq7mbU0/uQ5hrxHklpb6Xrtbr65Qr640+/INillEXCm3BLa4EehygW1NoGsJeWG+KI4ZyAm33FvDSlvISZJ8SOkQYtrnz6+jXidN0y7jk2ao8lS19zJ6D3jm84BtGrzc/b3Ae83sVsDJwCuAFxBOK9HLCvBgd18wszLwVTNrHLP/Fnf/hy2UW6Sv6rWEm3+1SLkSU5kqMTFVIi7rXFzrMbNmjc/03O4tr8/dqVdX1tayramNW2pOW52vtUm1urxEdXkZhlSz4Z5SX1np+8EVIjJeijQ1vgc4HriWUNv1ZOCi9ZbzUI+7kA2Ws5vqdmXbqS7VOfN132LPdUst4+NyxMRUKQSx6VLzcSX3eCKb1ngcppWpTMWUJ2I1PW2QZedAK0/05yRYnqbUqitZf7fFNbVs9ZUV6rUq9WqVpFptPq7XatSrKyTZfXO4WsvNky2X3fe7yXWcRXFMFJeI4ii7j4lKJeLm+Hh1nlJMnJsnirJl8vOXYqIoDvdxNr5Uym0nXjOvTnK8s5Wnpka27SJNjQcDMXALcBNwg7sX6p1rZjHhWo+/Brzd3b9hZo8EXmxmzyZcD/Ll7n5zh2VfQKhZ48gjjyyyOZFNuez8a9aELoCklrJYq7K4d3NHFVpkWUCLmZgut4WzXsEtC3hTJaJIwW0rLIqaR1juOuDAgW4rTRKSWo3aBgNbp6A3yBCXDzRrQkkuAFncFngaw42Qk1tHnA9FpdLqekv5ELW6Df0hkZ3MinYwNLPjgIcDfwjE7n7bwhsxOwD4BPAHhL5dNxBqv/4aOMLdn9dr+RNPPNEvuOCCXrOIbIqnzode8w1uuTacI2t6d4U0caqLddJ09BW0lcm4WePWqVZtYqrcNbhNTJeIS/pXLyIybGZ2obuf2GlakabGxwC/DdwfOBD4EqHJsTB3v8XMzgUeke/bZWb/Cnx6I+sS6aerfnBTM3RVJmNOec1JVCZLWV+jlOpSnZXFOitL9fB4qUY1P7zYGN/6uLpYp96HCyZXlxOqywkLN2+u31CpHHUJbmvHVSZLTEyXW4JbSdcCFBHpqyJNjY8kHJ34Vnf/ZdEVm9mhQC0LXVPAQ4E3mNkR7n5NNtvvAt/baKFFtsrdSeopl5xzdXPcne97BJXsivWhr1Hop7XrgIlNbSOpdwhu+XC2JrjVWoJbdXnrzU31Wkp9C82lUWRr+7VNd65dC8Gt9XFlsoSpuVREpKnIUY0v2uS6jwDOyPp5RcCZ7v5pM3u/md2D0NR4BfD7m1y/jClPPQSGakKtmlCvhsf1apoNdxm3klKrrZ1erybUVrLHtYRaNq79UI+7PaBw63khcSliarbC1GxlU8unqVNbXj+4dZyWjfctNpemqbO8r8byvtrmVmBQmYhbm0cbtWpFgtuUmktFZLwM7Noc7n4xcM8O4581qG3K4KWpdw40+ZBUywLRSnsIWg09raGpNST1o4luo253/EEccPhozunSTRRZCCnTmztnW6O5NISzGtWlpFmrlm8ebQS39ibUlaU6yVb3heeaS9lCc2nH5tEyE1ONPnBtfd4awW26RKms5lIR2T50UbQxkiQp9ZVOgWZtbdH6NUm5kLQSwlCtmpDWR9/hvF+iklGuxMwdMsVvPeXYURen7/LNpTMHbrK5tJbm+rethrMQ4JLVPm/LnYNbrV/NpXuqLO7ZQnNpl35tJdWmiexIJ/3uMc2uJcOm4DUEjf5Ea4NOCEq1atKsJWqZ3lKb1KhJ6lablG6Lo/D6pVSOKFViSpXV+3IlpjQRN6eVm9Oy6ROrj0vl3DKVmNJE7nElolSOiOLVH936DTcwf+63iCaniGZniGdniWZniWdnsR180ea4HDFdrjA9t/nm0mrH/mxtYa4Z3FZr5kKIS/rTXLpQY3lhk82lIjJ2Tnz00VT6c7rADev6i2Jml9DjhKfufsJASjRk7qv9iboFmrU1RR36EFXbQlLbfGNzWTBjNfSUOwSe3H25HAJPmL/D9E7BqhKC1TA7ZNdvvpmf/e7vklx/Q+enPD1NPDNDNDdLPDPbej87SzQ7Rzw7s/Z+bo54Zgabnt6xTV1RZEzuKjO5a/PNpbWVZAPBra3v22KdpK5L1IjI9tHrr/xjsvtG5/r3Z/enAIsDK9GQ/epne/n4Gy8cdTH6woxQI9RWG1ReJ/C0j+s4fSI8jkvj11/mlo98pGvoAvDFReqLi3DddZvbQBxnwW0uq01ru28JcW1hrhHeyjvz2pxmRmUy9Nma2eT5T+u1JNe/LTSPNoJbMkZN5yJSXHkyHtm2uwYvd78SwMzu5+73y036MzP7GvDaQRduGEqV4fTxiOLQ32a1CS0fdKLVwFTuHHhamtS6hKQotrELRYPmtRo3f+jDzeGJ444LFxLeu5dkYYF0fn7r1/pLEpI9e0j27Nn0KmxqqhnemvdttWwtYW5ujmimcT9LtGvn1rqVyjGlcrzp5lIRkX4q0nlll5n9lrt/FcDM7gvsGmyxhqdciYnLUWtTV1u4WROS1tQk9W4+iysRcaxOvNvR/Be/SD2ryYoPPYSjP/IfWGX1B9rTlHRxkXR+nmTvPOnCPMnevaQLC+F+foFkvve9Ly9vuZy+tER9aQmuv35zK4iiUJvWHt66hLiOTaYVBRcRka0qEryeD7zXzHYT+nztAXpe4md/csDh07zwbQ8cdTFkSNyddN8+0qwG6sbTz2hOO/CpT1sTLiyKiGdmiGdmKB9xxOa2Wa2G2rO9e0nmF0jnu93Pk8zPd7wn3WI/pTQl3bOHdM8e+MUvNrUKm5zs3FTavM83la5tMo2mp3VhYhHZ8YqcQPVC4O5mNke4tuPm20tE+iStVkluuSWEmT17SPaE+3TvnpbhJBtOm8N7odMFiMtlDnzqyQMpq1UqlA46CA46aFPLh7C42LO2bU2IazSVZve+tPYi4Bsux/IyyfJyz/5wPZl1rnXr1FTaqcl0dpZItW4isgHuHv64Jkl4nIQjpUfZ/aLItRoPB/4OuLW7P9LMjgfu4+7/NvDSyVjzJMma8PLhaW2YSrKamkZwSvbs6UvzXd7uRz2S0qGH9nWd/WJmxDO7iGd2Ub7VrTa1Dq/Vmn3WWppMCzSVNsJbx8C6oUKEvnPp3r3wy8JXH2thlcqmm0qjuTmiXbtU6yYD0/lHPoU0xZMk9BfNfvhJkzC+OT0F7zxfx3HZcmuWbU5vLJvNk+anpdBcJnffPl/HZbNyeNs6ui7bNi2/bMv03Hw9l11nHUnr69qtj+6x53+d0oGbPGJni4o0NZ4OvA94dTb8I+AjgIKX4O740tJqKLol1DKl69Q6JXv29Kfj+ibY9DTx3Bzx7t3Ec3NUjj6aLz3qCL74uVOZLE0yW55ltjLLTGWm5fFcZY6Z8gyzldnmbbq0f3Rat3I5fMls8ovG3fHFxc5NpgtZmJuf7xHi5vHFrR8M7dUqyQ03kNywhVq3mZlmrZpNbO7EsrKfywJS5x/vDj/yjR/xlh9+L/wjL9vQVrtvbEGR4HWIu59pZq8EcPe6mW39dNSyrXitlgWiLrVO7WEqF6CojeDElKXSanjavZtodyNI7c7GheFobo549wHN4Xhubk0/rouvv5jXf+aUTRUjsqgljLUHs07Dc5W5EOoqs8yWZynH2/9UEWaG7dpFtGsXHH74ptbh9Xros1akqbRLvzfq9a09EXfSbH11rtnaukRk/xBFEEWhtrtxv82D1z4zO5jsZKpmdhKhg71sM+4eftTyzXPNWqhcoLplNTQle/eQ3rKHtA+1EZsRzc42A1QITwfkAtUcUccwtbuv7fOfufwzm1429ZS91b3sre7d9Dom4omOoWy94UZt3HR5msi2f/OZlUpbr3VbXt7w0aX5fm+jep/LDtLpRz6KII47jIswy6aZQRxDZFgUr52vbZzFEVjrOjot2zJfY1qHZS0yiNrWsd6yLc8nWzaOwdrmi2Owtmlty1qUPf/cdMtet+ay+dckN7352nVbdpu1ShQJXn8EfBI4Jjt/16HAUwZaqh0uXVkhuSXXUTzXhNexyS5XAzWKFG+VSghGB4RQ1KhZWlMTdUB+/O5tcTmeJE34/BWfbw6/4sRXcMjUISxUF5ivzTNfDbeF6gJ7a3vD+Gx4vjbPUn3rndZXkhVWlla4YWlzzWeGrWkWzQez9abNVeaoxNu/07qZYVNTRFNTcPhhm1qH1+shrGVNpj6K2lrZHrr+8K/zI98pUJlt2x952X6K/OpdCjwAuBNgwA+B7f/3esQ8ScI/8TUdxbvUOuXm8ZWV4Rc4isLh/wfkapjm5rIw1Wiu69yEF02O6IJXfXDhtRc2A8/BkwfzzOOeSRwVP6NxLa2xUF1oCWoL1QX2VveyUFtoBrf56nzX4cS31nLveHOd7NvcOipRZd0m0m7DM5UZZsoz+02tW3zAAcQHHDDqoojIDlUkeJ3v7r9OCGAAmNlFwK8PrFTbRLNDcadap7bwlO7d0zo8Pz+SMtv0dEuNU7O5rkuTXaMWKpqZ2ZFHe332is82Hz/sqIdtKHQBlKMyB04eyIGTm28+W6ovdQ1me6urtWz5YNcYXqgusFjfevNZNa1y4/KN3Lh846aWN4xd5V1ratg6hrbKDHPlXJNpdpuI1dFdRMZfr4tk3wq4DTBlZvck1HYBzAHTQyjbUFSvuIKbPvghkj23dGzG23Jn3s0olZqdxgv1f2o04XXoOD7O3J1aWmMlWaGaVMMtra55vJKsUEtqVNNqc97Gcl+44gvN9T3y6EcO/TmYGdPlaabL0xzO5jqt19M6+2r71g1pvYbrvrX3ueMs1BZYqC1sutatHJVXg1h5tjWYtQ13qoGbKc9sODjL9pJ6Sj2tU0/r1NJa83Hd6yRpsjre66vT0myaty3TmObJmvGJJzg6AnEne+5dnst0eTRRpleN18OBU4HbAm/OjZ8HXjXAMg1Vcsst3Pz+968/4yZEs7NtTXadmvDW9n+y6e17igJ3p57WW0JMt9DTMfAkWVBKWx/n17OSZss1QlPa/XG/3GrXrbj7oXfv2/qGqRSV2D2xm90Tuze1vLuznCx3bSYt0mS6r7bJtJVTS2vctHwTNy3ftOl1NGvdCh6YsD/0beu3TiFkU2HHk7XLdFunt4WkbP3tyyoMybA87U5P237By93PAM4wsye5+8eGWKahiuZ6/1jZxETnI+zW9H9aPV1B1DhtQdyff9/uHr6kcmFks8GlZbn1aonyIScXoMbRE4994n7RR2kQzIyp0hRTpSkOm95cp/UkTZqBbL0m027D9XTrtcv7avv6EgJFRAalyCWDPmZmjwbuAkzmxr92kAUbltJhh3HYn/5p86i8ZtPe3G6Y3UWtbGuCy2LH4LLESnILtVqN6rVVVn7ZFnhyzVvtj4vUGKU+unOObGflqEwlrlCJKuG+/XHbuIl4orlM4/HtZm/HY455zKifyn4tjuIt17qtJCvr92vL9W1rn75QW+jzs5JRKEUlylGZkpUoRau32OLm43JU7j0tt2wcxc3hxnJxFBPpGLEdbao0NbJtF7lk0DsJfboeBLwHeDLwzQGXa2h+Ur2aPz7w4yHs3FClet1q+Nnq0WbjqmQlynGZiXiCSlRZfdwp+HQLPHGZSlTpvlw23D5vOW4NTTu1lmrcmBmTpUkmS5McMnXIptaRpAn76vtam0w7nBYkP9zP5ur9RTOg2GoIKUWllrDSDCi5QNMSeCy3XNv6uoWddQOUlYgs2rbdLET6pchRjfd19xPM7GJ3f42ZvQn4+KALNiyJJ1y+5/JRF2NdkUXNsNEIK90edwtAnQJPxwDVFoDyoakSVdSBWbalOIqZq8wxV5kbdVFERLoqErwaZ4hcNLNbAzcCRw+uSMNVibp3rjWsc81NXGYiKh5c8sPtNT9raoy6NJmVotGeaFRERES2rsiv+afN7ADgjcBFhEsHvWeQhRqmI+eO5L8e/18d+wWVrKRqbxEREembIp3r/zp7+DEz+zQw6e5jc63GSlzhmAOOGXUxREREZAfodQLVB7v7l8zsiR2m4e5j089LREREZBh61Xg9APgS8NgO05wx6mAvIiIiMgy9TqD6V9n9c4dXHBEREZHxte5JkMzs77LO9Y3hA83sbwZaKhEREZExVOTsk49091saA+5+M/CogZVIREREZEwVCV6xmU00BsxsCpjoMb+IiIiIdFDkPF4fAM42s/cROtU/DzhjoKUSERERGUNFzuP192Z2MfDQbNRfu/vnB1ssERERkfFT9Do03wbKhBqvbw+uOCIiIiLjq8hRjScD3wSeDJwMfMPMnjzogomIiIiMmyI1Xq8G7uXu1wGY2aHA/wAfHWTBRERERMZNkaMao0boytxYcDkRERERySlS4/U5M/s88OFs+KnAZwZXJBEREZHxVOSoxj/OLpT9W4AB73b3Twy8ZCIiIiJjZt3gZWa7gLPc/eNmdifgTmZWdvfa4IsnIiIiMj6K9NX6MjBhZrchdKp/LnD6IAslIiIiMo6KBC9z90XgicDb3P13geMHWywRERGR8VMoeJnZfYBTgP/OxhU98aqIiIiIZIoEr5cBrwQ+4e6XmtkdgHMGWioRERGRMVTkqMbzgPMAzCwCbnD3lwy6YCIiIiLjpsglgz5kZnPZ0Y3fB35oZn88+KKJiIiIjJciTY3Hu/te4AmEE6ceCTxrkIUSERERGUdFglfZzMqE4HVWdv4uH2ipRERERMZQkeD1LuAKYBfwZTO7PbB3kIUSERERGUdFOtf/E/BPuVFXmtmDBlckERERkfFU6HxcZvZo4C7AZG70awdSIhEREZExVeSoxncCTwX+gHCR7KcAtx9wuURERETGTpE+Xvd192cDN7v7a4D7ALcbbLFERERExk+R4LWU3S+a2a2BGnD04IokIiIiMp6K9PH6tJkdALwRuIhwKol/HWShRERERMZRkaMa/zp7+DEz+zQw6e571lvOzCaBLwMT2XY+6u5/ZWYHAR8BjiKcpuJkd795c8UXERER2X8UaWpscveVIqErswI82N3vDtwDeISZnQT8GXC2ux8LnJ0Ni4iIiIy9DQWvjfBgIRssZzcHHg+ckY0/g3BGfBEREZGxN7DgBWBmsZl9B7gO+KK7fwM43N2vAcjuD+uy7AvM7AIzu+D6668fZDFFREREhqLIebzMzJ5pZn+ZDR9pZvcusnJ3T9z9HsBtgXub2V2LFszd3+3uJ7r7iYceemjRxURERES2rSI1Xu8gnLvr6dnwPPD2jWzE3W8BzgUeAVxrZkcAZPfXbWRdIiIiIvurIsHrN939RcAyQHYEYmW9hczs0Ow0FJjZFPBQ4AfAJ4HnZLM9Bzhr48UWERER2f8UOY9XzcxiQsd4zOxQIC2w3BHAGdmyEXCmu3/azM4HzjSz5wM/J1yCSERERGTsFQle/wR8AjjMzP4WeDLw5+st5O4XA/fsMP5G4CEbLKeIiIjIfq9n8DKzCLgc+BNCWDLgCe5+2RDKJiIiIjJWegYvd0/N7E3ufh9C/ywRERER2aQineu/YGZPMjMbeGlERERExliRPl5/BOwC6ma2TGhudHefG2jJRERERMZMkYtkzw6jICIiIiLjbt3gZWb37zTe3b/c/+KIiIiIjK8iTY1/nHs8CdwbuBB48EBKJCIiIjKmijQ1PjY/bGa3A/5+YCUSERERGVNFjmpsdzVQ+GLXIiIiIhIU6eP1NrLLBRGC2j2A7w6wTCIiIiJjqUgfrwtyj+vAh939awMqj4iIiMjYKhK8DnD3t+ZHmNlL28eJiIiISG9F+ng9p8O4U/tcDhEREZGx17XGy8yeDjwDONrMPpmbNAvcOOiCiYiIiIybXk2NXweuAQ4B3pQbPw9cPMhCiQzV9T+Cs18DV34dzCAqQ5zdGo+jUutw83Gpy/y9puXXVWmdrzEtrqxdR1TK5u9Qnige9asoIiIFdA1e7n4lcCVwn+EVZwRuvgI+/+pRl0JGJanCT78EaX3UJdki6xLKSl1CXI8A2SvgFQqjBafFlbXzmY36hRQRGagip5M4CXgbcBxQAWJg39hcJHt5L/zg06MuhcgWOSQr4bY/s3iDAbJHiOsaIHPrUk2hyM50z2dBZXokmy5yVOM/A08D/hM4EXg28GuDLJTI0N3+fvCwv4bdR0JagyS7pfn7enZfzT2uhdqyJBvfMl/7tHrbOuvZumpt0+pd1tVjWvNUe/s5T6C+FG4iIoNylydu6+CFu//EzGJ3T4D3mdnXB1yu4TngSHjqB0ZdChml2SPgNr+xfzdzpUmXENceIPPTugTIXsutCZD5MNohQK4bVHPz7/fNvSIi6ysSvBbNrAJ8x8z+ntDhftdgizVEUwfAcY9ddzaRbS2KIZqC8tSoS7J57j1qAXsFyHwtYK/lcuGvEfo8HfWzFpFRKE+ObNNFgtezCOf7ejHwh8DtgCcNslAisgOZrXbEFxEZUz2Dl5nFwN+6+zOBZeA1QymViIiIyBjqeeb6rE/XoVlTo4iIiIhsQZGmxiuAr2Vnr9/XGOnubx5UoURERETGUZHg9cvsFhEuFyQiIiIim7Bu8HL31wCY2S5337fe/CIiIiLSWc8+XgBmdh8z+z5wWTZ8dzN7x8BLJiIiIjJm1g1ewD8CDwduBHD37wL3H2CZRERERMZSkeCFu1/VNioZQFlERERExlqRzvVXmdl9Ac9OK/ESsmZHERERESmuSI3XC4EXAbcBfgHcIxsWERERkQ0oclTjDcApQyiLiIiIyFgrclTjHczsU2Z2vZldZ2ZnmdkdhlE4ERERkXFSpKnxQ8CZwBHArYH/BD48yEKJiIiIjKMiwcvc/f3uXs9uHwB80AUTERERGTdFjmo8x8z+DPgPQuB6KvDfZnYQgLvfNMDyiYiIiIyNIsHrqdn977eNfx4hiKm/l4iIiEgBRY5qPHoYBREREREZd+sGLzOLgUcDR+Xnd/c3D65YIiIiIuOnSFPjp4Bl4BIgHWxxRERERMZXkeB1W3c/YeAlERERERlzRU4n8Vkze9jASyIiIiIy5orUeP0v8Akzi4AaYIC7+9xASyYiIiIyZooErzcB9wEucXedOFVERERkk4o0Nf4Y+J5Cl4iIiMjWFKnxugY418w+C6w0Rup0EiIiIiIbUyR4XZ7dKtlNRERERDahyJnrXwNgZrvcfd/giyQiIiIyntbt42Vm9zGz7wOXZcN3N7N3DLxkIiIiImOmSOf6fwQeDtwI4O7fBe4/wDKJiIiIjKUiwQt3v6ptVDKAsoiIiIiMta7By8xOyh5eZWb3BdzMKmb2CrJmRxEREREprleNV6Mf1wuBFwG3Aa4G7pENi4iIiMgGrNvU6O43uPsp7n64ux/m7s909xvXW87Mbmdm55jZZWZ2qZm9NBt/mpn9wsy+k90e1Y8nIiIiIrLd9TqdxB3M7JPdJrr749ZZdx14ubtfZGazwIVm9sVs2lvc/R82WFYRERGR/Vqv4HU94TqNm+Lu1xDOeo+7z5vZZYTmShEREZEdqVfwmnf38/qxETM7Crgn8A3gfsCLzezZwAWEWrGb+7EdERERke2sVx+vK/qxATObAT4GvMzd9wL/AhxD6KR/DV1q1czsBWZ2gZldcP311/ejKCIiIiIj1TV4ufsTt7pyMysTQtcH3f3j2XqvdffE3VPgX4F7d9n+u939RHc/8dBDD91qUURERERGrtAJVDfDzAz4N+Ayd39zbvwRudl+F/jeoMogIiIisp2se5HsLbgf8CzgEjP7TjbuVcDTzewegBOaM39/gGUQERER2TYKBS8zOwE4Kj9/o+mwG3f/KmAdJn1mA+UTERERGRvrBi8zey9wAnApkGajHegZvERERESkVZEar5Pc/fiBl0RERERkzBXpXH++mSl4iYiIiGxRkRqvMwjh61fACqHflrv7CQMtmYiIiMiYKRK83kt2dCKrfbxEREREZIOKBK+fu3vXi2WLiIiISDFFgtcPzOxDwKcITY3A+qeTEBEREZFWRYLXFCFwPSw3TqeTEBEREdmgdYOXuz93GAURERERGXdFTqD6PkINVwt3f95ASiQiIiIypoo0NX4693iScGHrXw6mOCIiIiLjq0hT48fyw2b2YeB/BlYiERERkTFV5Mz17Y4Fjux3QURERETGXZE+XvOEPl6W3f8K+NMBl0tERERk7BRpapwdRkFERERExl2RzvWY2QnAUfn5dQJVERERkY0p0tT4XuAE4FJWr9WoE6iKiIiIbFCRGq+T3P34gZdEREREZMwVOarxfDNT8BIRERHZoiI1XmcQwtevCNdsNMDd/YSBlkxERERkzBQJXu8FngVcwmofLxERERHZoCLB6+fu/smBl0RERERkzBUJXj8wsw8BnyI0NQI6nYSIiIjIRhUJXlOEwPWw3DidTkJERERkg4qcuf65wyiIiIiIyLjrGrzM7E/c/e/N7G2EGq4W7v6SgZZMREREZMz0qvG6LLu/YBgFERERERl3XYOXu38qe/gRd1/OTzOzQwZaKhEREZExVOTM9d80s5MaA2b2JODrgyuSiIiIyHgqclTjKcB7zexc4NbAwcCDB1koERERkXFU5KjGS8zsb4H3A/PA/d396oGXTERERGTMrBu8zOzfgGOAE4A7Ap8ys39297cPunAiIiIi46RIH6/vAQ9y98vd/fPAScCvD7ZYIiIiIuOnSFPjW9qG9wDPH1iJRERERMZUkabGY4HXAccDk43x7n6HAZZLREREZOwUaWp8H/AvQB14EPDvhI72IiIiIrIBRYLXlLufDZi7X+nup6HTSYiIiIhsWJHzeC2bWQT82MxeDPwCOGywxRIREREZP0VqvF4GTAMvAX4DeBbwnAGWSURERGQsFTmq8VvZwwXguYMtjoiIiMj46hq8zOyTvRZ098f1vzgiIiIi46tXjdd9gKuADwPfAGwoJRIREREZU72C162A3wGeDjwD+G/gw+5+6TAKJiIiIjJuunaud/fE3T/n7s8hXCboJ8C5ZvYHQyudiIiIyBjp2bnezCaARxNqvY4C/gn4+OCLJSIiIjJ+enWuPwO4K/BZ4DXu/r2hlUpERERkDPWq8XoWsA+4I/ASs2bfegPc3ecGXDYRERGRsdI1eLl7kZOrioiIiEhBClciIiIiQ6LgJSIiIjIkCl4iIiIiQ6LgJSIiIjIkCl4iIiIiQ6LgJSIiIjIkAwteZnY7MzvHzC4zs0vN7KXZ+IPM7Itm9uPs/sBBlUFERERkOxlkjVcdeLm7H0e41uOLzOx44M+As939WODsbFhERERk7A0seLn7Ne5+UfZ4HrgMuA3weOCMbLYzgCcMqgwiIiIi28lQ+niZ2VHAPYFvAIe7+zUQwhlw2DDKICIiIjJqAw9eZjYDfAx4mbvv3cByLzCzC8zsguuvv35wBRQREREZkoEGLzMrE0LXB93949noa83siGz6EcB1nZZ193e7+4nufuKhhx46yGKKiIiIDMUgj2o04N+Ay9z9zblJnwSekz1+DnDWoMogIiIisp2UBrju+wHPAi4xs+9k414FvB4408yeD/wceMoAyyAiIiKybQwseLn7VwHrMvkhg9quiIiIyHalM9eLiIiIDImCl4iIiMiQKHiJiIiIDMkgO9fvF67Zs8Sbv/AjIjOiyIgjwuPs1hyOjMggNsPMiLPhMN6I8/NE2TzWNk+HdTfmC8uQrbdxy62rsb3msqvrjhvzZ+tvrqs5fnVdndYtIiIiw7Hjg9fN+2r854VXj7oYIxVn4a01+K2GOmsPjVFrqFuzbC6krg2ERmy0LBtFrYGyaGgsEoa7lisXhm+9e4p7H32QQqiIiAzcjg9eqfuoizBySeok7OzX4bTHHs+p9zt61MUQEZExt+OD160PmOINT7obqYcA4u4kqZN6CGWpO0maPc7GJ83HTuKOZ8u2z9Oyrub8ucdrttO2LvdmufJlWV3v6vKdll277dbhdGdnrRbvOPennHLS7SnH6vYoIiKDs+OD10G7Kjz1XkeOuhgj4Y3Q2Ah5WcBM8gGyEUbbQlun0Jhf3rsF1rZ1Fwmd64XhNWXrEYZXyx+e/zk/vI6bF2tcN7/CZ7/3Kx5391uPereI9JSmzkK1zvxynYXlOgsrNfZmj+ez4YXlehi3Umd+uZbd11msJrhq+UX46Avvy4G7KiPZ9o4PXjuZZX2goq7nuR1/b/2fH/OW//kRAKd/7XIFr0ynEJu442mHAN0SaIsF8t6Bfm0wT1Na/xDkapSbob1nsGdNgE+2cZVvtZ6ysFJnbxaaVkNVuInI1iQj/AOi4CU72jN+80j++ZwfU0uci35+C+/72uUctKvSs5m2e61chybmfJjoGBrYUADp1HTdK4AkuZrMlgCyTi2jiIgMhoKX7GiHzk7w2BNuzce//QsAXvOp74+4RCLr21WJmZ0sMzNZYnayxMxEibnJMjMTpbXjssezkyV2TZSIdm4Ft0jTAVPlkW1bwUt2vOfe7+hm8JJV7afx6HSuuY6nGelxOo/8qUQ6rX/NuerWnDKk13nyOpz/Ln++vA7r365nEClFUQhPkyVmJ0otIWtXpUSs9CSy31Lwkh3vbrfdzVufdg++8P1rwelwvq/2H/LcucbafuzXnBx3vfOcdTi5breAsLqu9U+g21i+5zbNsKj1fGv586eJiEj/KXiJAI+/x214/D1uM+piiIjImNNJi0RERESGRMFLREREZEgUvERERESGRMFLREREZEgUvERERESGRMFLREREZEgUvERERESGRMFLREREZEgUvERERESGRMFLREREZEgUvERERESGRMFLREREZEjM3UddhnWZ2fXAlR0m7Qb2dFlso9MOAW7YVAEHq9fzGPW6N7p80fnXm2+z07uN174f/LLa9xszqH0/7M/8RpbZyr7fzDTt+8Eu36/9vt4823Xf397dD+04xd332xvw7n5NAy4Y9fPZ6PMY9bo3unzR+debb7PTe4zXvh/wstr322PfD/szP6x9v5lp2veDXb5f+30c9/3+3tT4qQFM224GWdatrnujyxedf735Njt9f9rvsH33/WaW1b7fmEGVd9if+Y0ss5V9Py7f9zA++75f+329efa7fb9fNDUOg5ld4O4njrocMnza9zuX9v3OpX2/c4163+/vNV799O5RF0BGRvt+59K+37m073euke571XiJiIiIDIlqvERERESGRMFLREREZEgUvERERESGRMGrCzM7zszeaWYfNbP/O+ryyHCZ2S4zu9DMHjPqssjwmNkDzewr2Wf/gaMujwyHmUVm9rdm9jYze86oyyPDY2a/nX3e32NmXx/GNndU8DKz95rZdWb2vbbxjzCzH5rZT8zszwDc/TJ3fyFwMqBDjvdzG9n3mT8FzhxuKWUQNrjvHVgAJoGrh11W6Z8N7vfHA7cBami/7/c2+Fv/ley3/tPAGcMo344KXsDpwCPyI8wsBt4OPBI4Hni6mR2fTXsc8FXg7OEWUwbgdAruezN7KPB94NphF1IG4nSKf+6/4u6PJATv1wy5nNJfp1N8v98JON/d/whQC8f+73Q28FufeQbw4WEUbkcFL3f/MnBT2+h7Az9x95+5exX4D8K/H9z9k+5+X+CU4ZZU+m2D+/5BwEmED+L/MbMd9TkZNxvZ9+6eZtNvBiaGWEzpsw1+5q8m7HOAZHillEHY6G+9mR0J7HH3vcMoX2kYG9nmbgNclRu+GvjNrH/HEwlfvp8ZfrFkCDrue3d/MYCZnQrckPsxlvHR7XP/RODhwAHAP4+gXDJYHfc78FbgbWb228CXR1EwGbhu+x7g+cD7hlUQBS+wDuPc3c8Fzh1uUWTIOu775gP304dXFBmybp/7jwMfH3ZhZGi67fdFwo+vjK+u3/fu/lfDLIiaUELqvV1u+LbAL0dUFhku7fudS/t+Z9J+37m2zb5X8IJvAcea2dFmVgGeBnxyxGWS4dC+37m073cm7feda9vs+x0VvMzsw8D5wJ3M7Goze76714EXA58HLgPOdPdLR1lO6T/t+51L+35n0n7fubb7vtdFskVERESGZEfVeImIiIiMkoKXiIiIyJAoeImIiIgMiYKXiIiIyJAoeImIiIgMiYKXiIiIyJAoeInIppmZm9mbcsOvMLPT+rTu083syf1Y1zrbeYqZXWZm57SNj8zsn8zse2Z2iZl9y8yOHnBZrjCzQwa5DREZLQUvEdmKFeCJ2y0smFm8gdmfD/w/d39Q2/inArcGTnD3uwG/C9zSnxKKyE6l4CUiW1EH3g38YfuE9horM1vI7h9oZueZ2Zlm9iMze72ZnWJm38xqlo7JreahZvaVbL7HZMvHZvbGrAbqYjP7/dx6zzGzDwGXdCjP07P1f8/M3pCN+0vgt4B3mtkb2xY5ArjG3VMAd7/a3W/OlvsXM7vAzC41s9fktnGFmf2dmZ2fTf91M/u8mf3UzF6YK+eXzewTZvZ9M3unma35LjazZ2avyXfM7F3Z846z17VRC7fmdReR7a006gKIyH7v7cDFZvb3G1jm7sBxwE3Az4D3uPu9zeylwB8AL8vmOwp4AHAMcI6Z/RrwbGCPu9/LzCaAr5nZF7L57w3c1d0vz2/MzG4NvAH4DeBm4Atm9gR3f62ZPRh4hbtf0FbGM4GvmtlvA2cDH3D3b2fTXu3uN2U1a2eb2QnufnE27Sp3v4+ZvQU4HbgfMAlcCrwzV87jgSuBzwFPBD6aK+9xhBq3+7l7zczeAZySreM27n7XbL4D1n2lRWRbUY2XiGyJu+8F/h14yQYW+5a7X+PuK8BPgUZwuoQQthrOdPfU3X9MCGh3Bh4GPNvMvgN8AzgYODab/5vtoStzL+Bcd78+u2bbB4H7r/O8rgbuBLwSSAkB6yHZ5JPN7CLg28BdCCGqoXHh3UuAb7j7vLtfDyzngtI33f1n7p4AHybUuuU9hBASv5U9z4cAd8hegzuY2dvM7BHA3l7PQUS2H9V4iUg//CNwEfC+3Lg62Z87MzOgkpu2knuc5oZTWr+X2i8m64ABf+Dun89PMLMHAvu6lM/WKX9HWTD8LPBZM7sWeIKZ/Qx4BXAvd7/ZzE4n1Gg15J9L+/NsPLdOz6u9vGe4+yvXPBGzuwMPB14EnAw8b6PPS0RGRzVeIrJl7n4ToWnu+bnRVxBqbQAeD5Q3seqnZEcXHkOo8fkh8Hng/5pZGcDM7mhmu9ZZzzeAB5jZIVnz4NOB83otkPXPunX2OAJOIDQNzhEC3h4zOxx45Cae173N7OhsvU8Fvto2/WzgyWZ2WLb9g8zs9tlBDJG7fwz4C+DXN7FtERkh1XiJSL+8CXhxbvhfgbPM7JuEINGtNqqXHxIC0uHAC9192czeQ2iOvCirSbseeEKvlbj7NWb2SuAcQm3SZ9z9rHW2fRjwr1k/MoBvAv+cleHbhP5WPwO+tonndT7weuBuwJeBT7SV9/tm9ueEvmgRUCPUcC0B78t1xl9TIyYi25u5t9dwi4jIoGRNoq9w98eMuCgiMgJqahQREREZEtV4iYiIiAyJarxEREREhkTBS0RERGRIFLxEREREhkTBS0RERGRIFLxEREREhkTBS0RERGRI/j/FbPWwG26WkwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 720x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# plot the results, and and ensure the values converge to a single number\n",
    "plt.figure(figsize=(10,6))\n",
    "plt.semilogx(nsamps,qvals,lw=3)\n",
    "plt.xlabel('Number of Samples')\n",
    "plt.ylabel('Maximum Temperature associated with Quantile')\n",
    "plt.title('Variation of Quantile estimate with Number of Samples',fontsize=16)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "86ad49eb-aea2-440d-95e9-6550930eda25",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[-1.1   -0.2    0.1    0.5    0.397  0.2  ]\n",
      " [ 0.6    0.1    0.1   -0.1   -0.1   -0.1  ]\n",
      " [-0.5   -0.1    0.     0.3    0.1    0.2  ]\n",
      " [ 0.     0.    -0.1    0.     0.     0.   ]\n",
      " [ 0.1    0.     0.1    0.     0.     0.1  ]\n",
      " [-0.1    0.     0.    -0.1    0.    -0.1  ]\n",
      " [ 0.1    0.     0.     0.     0.     0.   ]]\n"
     ]
    }
   ],
   "source": [
    "print(np.round(np.diff(qvals,axis=0),3))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3370cd06-aa1d-455d-80e5-a14f59dd0433",
   "metadata": {},
   "source": [
    "The quantile values converge after ~1e5 (one hundred thousand) samples to the order of our precision, $\\pm$1$^{\\circ}$C, which is higher than the precision of the query observations). This implies that the sampling method uses a large enough number to reproduce the quantities of the underlying distribution, at large scale for the whole domain. At smaller scale or a given geolocation with fewer outcomes/smaller maximum temperature range than the whole domain, we would expect this resampling to converge faster.\n",
    "\n",
    "Accordingly, this method is used going forward in our attribution analyses with ```N=100_000``` samples."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "34a31196-f7fa-403c-a8c3-7dc22b29dc35",
   "metadata": {
    "tags": []
   },
   "source": [
    "# Create Single Pooled Counterfactual Output File"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "eb11821f-d55e-4063-aeab-3bdc49cb4e20",
   "metadata": {},
   "source": [
    "We have created a set of pooled counterfactual distributions and stored them alongside each model's bias-adjusted natural distribution, which enables comparison, etc. in this notebook.\n",
    "\n",
    "Next, let's create a single file that has stores a resampled set of values from the pooled counterfactual distribution, for use throughout the remainder of our attribution analysis."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "27b30b08-9019-4da7-8eca-da4bef9f2fcd",
   "metadata": {},
   "outputs": [],
   "source": [
    "# define the number of samples we will take at each location,\n",
    "# given the sensitivity found above:\n",
    "N=100_000"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "8cf7919b-45c8-4095-8f64-e74a1916d0c3",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Completed resampling, month==1\n",
      "Completed resampling, month==2\n",
      "Completed resampling, month==3\n",
      "Completed resampling, month==4\n",
      "Completed resampling, month==5\n",
      "Completed resampling, month==6\n",
      "Completed resampling, month==7\n",
      "Completed resampling, month==8\n",
      "Completed resampling, month==9\n",
      "Completed resampling, month==10\n",
      "Completed resampling, month==11\n",
      "Completed resampling, month==12\n"
     ]
    }
   ],
   "source": [
    "# loop over the months\n",
    "for mi in months:\n",
    "    \n",
    "    # load each month's pooled data\n",
    "    mon_zarrname='pooledCF_mon'+str(mi).zfill(2)+'.zarr'\n",
    "    ds_mi=util.load_zarr_local(pooled_loadpath,mon_zarrname)\n",
    "    # grab the CF distribution we need\n",
    "    CF_pooled=ds_mi[pooledvarname].load()\n",
    "    del ds_mi\n",
    "    \n",
    "    # resample the distribution over space N times\n",
    "    CF_resample=xr.apply_ufunc(\n",
    "                util.randomly_sample_histogram,\n",
    "                CF_pooled,\n",
    "                kwargs=dict(N=N,bins=CF_pooled[binkey]),\n",
    "                input_core_dims=[\n",
    "                    [binkey],\n",
    "                ],\n",
    "                output_core_dims=[\n",
    "                    ['samples'],\n",
    "                ],\n",
    "                vectorize=True\n",
    "            )\n",
    "    \n",
    "    # store the result\n",
    "    store = s3fs.S3Map(root='s3://gridded-data-dev/CMIP5_resample/'+monzarr,\n",
    "                       s3=s3fs.S3FileSystem(), check=False)\n",
    "    encoding = {CF_resample.name: {'compressor':zarr.Blosc(cname='zstd', clevel=3)}}\n",
    "    CF_resample.to_dataset().to_zarr(store=store, encoding=encoding, consolidated=True, mode='w')\n",
    "    \n",
    "    # clean up\n",
    "    del CF_pooled, CF_resample\n",
    "    \n",
    "    # tell the user where we are\n",
    "    print('Completed resampling, month=='+str(mi))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c9325024-7a07-48f2-b7ce-1b050f53c259",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "AttributionKernel",
   "language": "python",
   "name": "att"
  },
  "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.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}