GWFA_analysis.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "!pip install -q -q -q geopandas\n",
    "!pip install -q -q -q rasterio\n",
    "!pip install -q -q -q py7zr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import geopandas as gpd\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "\n",
    "import rasterio\n",
    "from rasterio import features\n",
    "\n",
    "from sklearn.linear_model import LinearRegression\n",
    "from sklearn.neighbors import KDTree\n",
    "from sklearn.neighbors import RadiusNeighborsRegressor"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Chargement des fonctions et des données"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def rasterizeBlackWhite(sdf, size, points=True):\n",
    "    \"\"\"sdf : a geopandas DataFrame\n",
    "       size : resolution\n",
    "       points : if points is true return a pandas DataFrame, otherwise a numpy array\"\"\"\n",
    "    bounds = sdf.total_bounds\n",
    "    geometry = sdf.geometry\n",
    "    xmin = np.floor(bounds[0] / size) * size\n",
    "    ymin = np.floor(bounds[1] / size) * size\n",
    "    xmax = np.ceil(bounds[2] / size) * size\n",
    "    ymax = np.ceil(bounds[3] / size) * size\n",
    "    geometry = geometry.translate(xoff=-xmin, yoff=-ymin, zoff=0.0)\n",
    "    geometry = geometry.scale(xfact=1/size,yfact=1/size,origin=(0,0))\n",
    "    width = int((xmax - xmin)/size)\n",
    "    height = int((ymax - ymin)/size)\n",
    "    raster = features.rasterize(geometry, (height,width),fill=0,all_touched=True)\n",
    "    if points==True:\n",
    "        y, x = np.where(raster) \n",
    "        return(pd.DataFrame({'x' : x * size + xmin , 'y' : y * size + ymin}))\n",
    "    else:\n",
    "        return(np.flip(raster.reshape(raster.shape[0], raster.shape[1]),0))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def neighbors_near_point(df, list_of_radius):\n",
    "    neighbors = np.zeros([df.shape[0],list_of_radius.size],dtype=np.float)\n",
    "    tree = KDTree(df[['x', 'y']]) \n",
    "    dist, ind = tree.query(df[['x', 'y']], np.max(list_of_radius))   \n",
    "    for i, r in enumerate(list_of_radius):\n",
    "        neighbors[:,i] =list(map(lambda x : np.sum(x <= r), dist)) \n",
    "    result = pd.DataFrame(neighbors)\n",
    "    result.columns = [str(t) for t in list_of_radius]\n",
    "    return(pd.concat([df[['x', 'y']].reset_index(), result], axis=1))\n",
    "\n",
    "\n",
    "def neighbors_near_point(df, list_of_radius):\n",
    "    neighbors = np.zeros([df.shape[0],list_of_radius.size],dtype=np.float)\n",
    "    tree = KDTree(df[['x', 'y']]) \n",
    "    for i, r in enumerate(list_of_radius):\n",
    "        print(r)\n",
    "        neighbors[:,i] = tree.query_radius(df[['x', 'y']], r= r, count_only=True) \n",
    "    result = pd.DataFrame(neighbors)\n",
    "    result.columns = [str(t) for t in list_of_radius]\n",
    "    return(pd.concat([df[['x', 'y']].reset_index(), result], axis=1))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [],
   "source": [
    "def Gwfa_en_un_point(df, x_centre, y_centre, list_of_radius):\n",
    "    q=0\n",
    "    distance = (df.x - x_centre)** 2 + (df.y - y_centre)**2\n",
    "    weight = 1 / (1 - (distance / 16000**2))**2\n",
    "    temp = df[[str(t) for t in list_of_radius]]\n",
    "    temp = temp.mul(weight, axis=0)\n",
    "    temp = pd.DataFrame({'radius': list_of_radius, 'sb': np.sum(temp**(q-1),0)**(q-1)})\n",
    "    #temp.columns= ['radius', 'sb']\n",
    "    return(temp)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Calcul bocage"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [],
   "source": [
    "list_of_radius = np.array([25,50,100,200,400,800,1600])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bocage = gpd.read_file('building_bocage.gpkg')\n",
    "r_bocage = rasterizeBlackWhite(bocage, 4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "25\n",
      "50\n",
      "100\n",
      "200\n",
      "400\n",
      "800\n",
      "1600\n"
     ]
    }
   ],
   "source": [
    "neigh_bocage = neighbors_near_point(r_bocage, list_of_radius)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {},
   "outputs": [],
   "source": [
    "sb_bocage = Gwfa_en_un_point(neigh_bocage,387380,6893920, list_of_radius)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>radius</th>\n",
       "      <th>sb</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>25</th>\n",
       "      <td>25</td>\n",
       "      <td>0.000087</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>50</th>\n",
       "      <td>50</td>\n",
       "      <td>0.000191</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>100</th>\n",
       "      <td>100</td>\n",
       "      <td>0.000386</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>200</th>\n",
       "      <td>200</td>\n",
       "      <td>0.000731</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>400</th>\n",
       "      <td>400</td>\n",
       "      <td>0.001834</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>800</th>\n",
       "      <td>800</td>\n",
       "      <td>0.006837</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1600</th>\n",
       "      <td>1600</td>\n",
       "      <td>0.026580</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "      radius        sb\n",
       "25        25  0.000087\n",
       "50        50  0.000191\n",
       "100      100  0.000386\n",
       "200      200  0.000731\n",
       "400      400  0.001834\n",
       "800      800  0.006837\n",
       "1600    1600  0.026580"
      ]
     },
     "execution_count": 56,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sb_bocage"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 170,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([  1.33393017, -13.91206393])"
      ]
     },
     "execution_count": 170,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "lsb_bocage = np.log(sb_bocage)\n",
    "d = np.polyfit(lsb_bocage['radius'],lsb_bocage['sb'],1)\n",
    "f = np.poly1d(d)\n",
    "lsb_bocage['line'] = f(lsb_bocage['radius'])\n",
    "d"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 171,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'Mr(i)')"
      ]
     },
     "execution_count": 171,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAnoAAAG4CAYAAADSYkouAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOzdd5hU5dn48e/jEtAFsUECqIiaoFkUbNgbkcToC77RmJC8KSqIvURjizHGaNRgb1GxGxUbFiKisSG2IMVIEWsUGxpFFF2Wtsvz++Os/JAssjs7M+fs2e/nuuZyd2bOfe4zDMeb55znfkKMEUmSJOXPKmknIEmSpNKw0JMkScopCz1JkqScstCTJEnKKQs9SZKknLLQkyRJyqk2aSeQRWuuuWb89re/nXYaRTFv3jzat2+fdhpLlTufUu6v2LGbG6852xeybda+W3mTl883i8fheah0sQqNUe7t8mby5MmzY4ydG3wxxuhjuUfPnj1jXowdOzbtFL6i3PmUcn/Fjt3ceM3ZvpBts/bdypu8fL5ZPA7PQ6WLVWiMcm+XN8CkuIKaxku3kiRJOWWhJ0mSlFMWepIkSTlloSdJkpRTFnqSJEk5ZaEnSZKUUxZ6kiRJOWWhJ0mSlFMWepIkSTlloSdJkpRTIVk5QwAhhIHAwK5duw4dMWJE2ukURXV1NR06dEg7jaXKnU8p91fs2M2N15ztC9k2a9+tvMnL55vF4/A8VLpYhcYo93Z5069fv8kxxm0afHFFa6O15odr3ZaOa0yWLp5r3eZLXj7fLB6H56HSxXKt23TgWreSJEmtT5u0E5AkSVLDqhfWMnrKLGZ+Mo8e67RnQJ9udGjX+PLNQk+SJCmDJs6cw4E3TiBGqFlUR2XbCs56cAY3HbQtfXus3agYXrqVJEnKmOqFtRx44wTmLayjZlEdkBR78xbW1T9f26g4FnqSJEkZM3rKLFbUGCVGGD11VqPiWOhJkiRlzMxP5i0dyVtezaI6Zs6uaVQcCz1JkqSM6bFOeyrbVjT4WmXbCnp0qmxUHAs9SZKkjBnQpxshNPxaCDCgd7dGxbHQkyRJypgO7dpw00Hb0r5dxdKRvcq2FbRvV1H/fOMap9heRZIkKYP69libCaf2Z/TUWcycXUOPTpUM6N2t0UUeWOhJkiRlVvt2bRjUt3vB23vpVpIkKacs9CRJknLKQk+SJCmnLPQkSZJyykJPkiQppyz0JEmScspCT5IkKacs9CRJknLKQk+SJCmnLPQkSZJyykJPkiQppyz0JEmScirEGNPOITNCCAOBgV27dh06YsSItNMpiurqajp06JB2GkuVO59S7q/YsZsbrznbF7Jt1r5beZOXzzeLx+F5qHSxCo1R7u3ypl+/fpNjjNs0+GKM0cdyj549e8a8GDt2bNopfEW58ynl/oodu7nxmrN9Idtm7buVN3n5fLN4HJ6HSher0Bjl3i5vgElxBTWNl24lSZJyykJPkiQppyz0JEmScspCT5IkKacs9CRJknLKQk+SJCmnLPQkSZKy7NVXYf78gja10JMkScqi2lo491zo3Tv5bwHaFDklSZIkNdeUKTB4MLzwAvz4x3DEEQWFcURPkiQpKxYtgj/+EbbZBt57D+6+G0aOhC5dCgrniJ4kSVIWTJyYjOJNnw6//CVccgmss06zQjqiJ0mSlKb58+Gkk2D77WHOHHjgAbjllmYXeeCIniRJUnqeeSYZxXv9dTj4YDj/fFhzzaKFd0RPkiSp3Kqr4eijYdddYfFiePRRuPbaohZ5YKEnSZJUXo89BptvDldcAUcdBdOmQf/+JdmVhZ4kSVI5zJ0LQ4fC978P3/gGPPUUXHYZdOhQsl1a6EmSJJXa6NFQVQU33JBMvJgyBXbZpeS7tdCTJEkqlU8+SVqlDBwIa60F48fDsGGw2mpl2b2FniRJUimMHJmM4t15J5x+OkyeDH37ljUF26tIkiQV04cfJpMs7rkHttoKHnkE+vRJJRVH9CRJkoohxqTRcVVVck/euefC88+nVuSBI3qSJEnN9957cOihMGYM7LBDMuli003TzsoRPUmSpILFCNdcA716wdixyfq0Tz+diSIPHNGTJEkqzJtvJn3xnngC+vVLVrbYeOO0s/oKR/QkSZKaYsmSpNHx5pvDxIlw9dXJahcZK/LAET1JkqTGe/VVGDIEnn0W9toLhg+H9ddPO6sVckRPkiRpZWprk0bHffrAjBlw883w4IOZLvLAET1JkqSvN20aHHRQ0vB4333hyiuhS5e0s2oUR/QkSZIasmgRnHEGbL01vPMO3HVX0gS5hRR54IieJEnSf5s0CQYPTkbz/u//4NJLoVOntLNqMkf0JEmSvjR/PpxyCmy3HcyeDaNGwW23tcgiDxzRkyRJSjz7bDKj9suZtRdcAGuumXZWzeKIniRJat3mzYNjj4VddoEFC+CRR+C661p8kQcWepIkqTV7/PGk8fFll8ERRyT35H3/+2lnVTS5L/RCCOuHEB4PIbwcQngphHBeCCGknZckSUrR3Llw6KHQvz9UVMC4cXDFFbD66mlnVlS5L/SAWuDkGON3gS2B7YD90k1JkiSl5sEHoVev5PLsCSfAlCmw665pZ1USqRR6IYRvhxCGhxCmhBDqQghPruB9VfWjcTUhhFkhhDNDCBVN2VeM8YMY46T6nxcBU4Fst7GWJEnFN2cO/PrXMGBAcv/dP/8J558PlZVpZ1Yyac267QXsDYwH2jb0hhDCWsBjwAzgf4GNgQtJitPTCtlpCGEd4EfADwrZXpIktVD33ANHHgmffAJ/+AP8/vfQrl3aWZVcWoXeAzHGUQAhhJFAQ81pDgNWA/aLMX4OPBpC6AicEUI4r/45QgjPAOs1sP3jMcYhX/4SQmgHjAQuiTG+XNzDkSRJmfSf/8BRR8HIkbDllvDww7DFFmlnVTapFHoxxiWNeNtewD++LOjq3QEMA3YDHqiPtfPKAtVf7r0N+FeM8cKmZyxJklqUGGHECDjmGKiuhnPOSe7H+8Y30s6srLI8GWNT4JVln4gxvgPU1L/WFMOBL4DfFic1SZKUWe+/D/vsA7/8JfTsCS++CL/7Xasr8gBCjDHdBOov3cYYd1/u+cXAiTHGS5Z7/j3gbzHGUxsZfyfgGWA6UFf/9A0xxsuWe98hwCEAnTt33vquu+4q4Giyp7q6mg4dOqSdxlLlzqeU+yt27ObGa872hWybte9W3uTl883icXgeKl2sQmMUbbsY6TpmDBtfdRWhtpa3hgzhvf32S9qn5Fi/fv0mxxi3afDFGGOqD5L75p5s4PnFwLENPP8+cHYpc+rZs2fMi7Fjx6adwleUO59S7q/YsZsbrznbF7Jt1r5beZOXzzeLx+F5qHSxCo1RlO3efDPG/v1jhBh32y3G118vKGZLBEyKK6hpsnzp9lOgobVH1gA+K3MukiQpi5YsSRodb745jB8PV10FTzwB3/522pllQlqzbhvjFZa7Fy+EsD7QnuXu3ZMkSa3Pau++C7vtBs88A3vuCddcA927p51WpmR5RO8hYM8QwrJrkQwC5gPj0klJkiSlrrYWzj+fbQ4+GKZPh5tugocesshrQCojeiGESpKGyQDrAh1DCPvX/z4mxlgDXA0cA9wbQhgGbAScAVwUv9pyRZIktRbTp8PgwTBxIp/utBOd7r4bunZNO6vMSmXWbQihB/DWCl7eMMY4s/59VcAVwA4k9+VdB5wRY6xbwbbNzWsgMLBr165DR4wYUYpdlF3WZrw526108Zx1my95+XyzeByeh0oXq5SzbsPixXQfMYINbr2V2vbtef2YY3irb186rL76127XGmR61m0WH866LR1nu5UunrNu8yUvn28Wj8PzUOlilWzW7aRJMfbuncyo/dnPYvzoo2btL29oobNuJUlSa7ZgAZx6Kmy3HXz8Mdx/P9x+O3TunHZmLUaWZ91KkqTW6p//TO7Fe+UVOOgguPBCWGuttLNqcRzRkyRJ2TFvHhx3HOy0E9TUwMMPww03WOQVyBE9SZKUDWPHwsEHw5tvwhFHwF/+Ak62aBZH9CRJUqoq5s2Dww+H730PQoAnn4S//tUirwhSaa+SVbZXKT3bGpQunu1V8iUvn28Wj8PzUOliFRJj7eef5zsXXMCqc+bw3v7789ZBB7Fk1VVLtr88sr2K7VUyw7YGpYtne5V8ycvnm8Xj8DxUulhNivHJJzEecECMEKs32CDGf/6ztPvLMb6mvYr36EmSpPK6777kHryPP4bf/55Ju+7Kbttvn3ZWueQ9epIkqTw++ggGDYL99oMuXWDiRPjzn4lt26adWW5Z6EmSpNKKMWl0XFWVND3+859hwgTYcsu0M8s9L91KkqTSmTUrmVH797/DttsmPfF69Uo7q1bDET1JktQs1QtruWPCO3z4+QLumPAO1Qtrk1G8G25IRvEeeSRZ2eK55yzyyswRPUmSVLCJM+dw4I0TiBEO32QhV02awXW3jeWeF25ijWeehF13heuug+98J+1UWyULPUmSVJDqhbUceOME5i2sS55YsoQfjx/FyeNuJgALL72MdkcdCat4ATEtNkxehg2TS89GpaWLZ8PkfMnL55vF4/A8VLxYn9YsYtZnC1gSI2t8OIu9rr+Czi/P4O3Nt2DckCNZ/TvdWaty5TNqCz2OLH6/0mDDZBsmZ4aNSksXz4bJ+ZKXzzeLx+F5qHixzh0zI2544qh4Vr/BcX6btnFBZWU8Ya9j4wYnPRA3OHl0/MuYl0u27+ZslzfYMFmSJBVb77nvc/9tJ9F71qs8+u3tePfYQ7n7vW8CUNm2gh6dKlPOUF40lyRJTbN4Mfz5z+x1wP+w/qcfcPTAExm632nMW2vtpW8JAQb07pZikgILPUmS1BT/+lfSD+8PfyDsuy8zn5rAE1t+j8p2yUXCyrYVtG9XwU0HbUv7dl44TJt/ApIkaeUWLoSzzoK//AU6d4Z774V992VLYMLmGzN66iziBzP448AqBvTuZpGXEY7oSZKkrzd+fLJc2dlnw69+BTNmwL77Ln25fbs2DOrbnS4dV2VQ3+4WeRlioSdJkhpWUwO//S3suCNUV8NDD8GNN8Jaa6WdmRrJkluSJP23ceNgyBD497/hsMNg2DDo2DHtrNRENkxehg2TS89GpaWLZ8PkfMnL55vF4/A89PWxKmpq2Gj4cNb9+9+Z360br554Ip9tsUVJ87FhcvPYMNmGyZlho9LSxbNhcr7k5fPN4nF4HvqaWA8/HGP37jGGEONxx8VYXV2WfGyY3DzYMFmSJK3Qp5+yybBh8PDDsOmm8OyzsMMOaWelInAyhiRJrdmoUdCrF10eeQROPTXpk2eRlxsWepIktUYffww//zn86EfQuTOTr7oqaZ+y6qppZ6YistCTJKk1iRHuvBOqquCee+DMM2HiRKp79kw7M5WA9+hJktRafPABHH54crm2b1+44QbYbLO0s1IJOaInSVLexQg33ZSM4v3jH3D++fDccxZ5rYAjepIk5dk778AhhyQF3s47w/XXg5dpWw1H9CRJyqMlS+Dqq6FXL3jmGbj88mS1C4u8VsURPUmS8uaNN+Dgg5PCrn9/uPZa6NEj7ayUApdAW4ZLoJWeSw+VLp5LoOVLXj7fLB5Hrs9DdXWsd++9bHj99cQ2bXjj8MP5cO+9IYSy5OkSaOlwCTSXQMsMlx4qXTyXQMuXvHy+WTyO3J6HZsyIcfvtY4QYBwyI8b33Co9VjHwyvF3e8DVLoHmPniRJLdnixXDOObDFFvDaa3DbbfD3v8O666admTLAe/QkSWqhOrzxBhx/fLJs2U9+kky4+Na30k5LGWKhJ0lSS7NwIZx9Nludcw506gQjR8KPf5x2VsogCz1JklqSCRNg8GB46SU++sEP6HL77bD22mlnpYyy0JMkqSWYPx9OPx0uugi6dYMHH+SVykq6WOTpazgZQ5KkrHv6aejTBy64AIYOhZdegr33TjsrtQAWepIkZVV1NRx1FOy6K9TWwuOPJ6tddOyYdmZqISz0JEnKokcfhc02gyuvhGOPhWnT4HvfSzsrtTAWepIkZclnnyXLl/3gB7Dqqsll20sugfbt085MLZCFniRJWfHAA9CrF9x0E5xyCrz4Iuy0U9pZqQWz0JMkKW2zZ8MvfgH77APrrAPjx8O55yYjelIzWOhJkpSWGOHuu6GqKvnvGWfApEmwTcPr00tNFZK1cAUQQhgIDOzatevQESNGpJ1OUVRXV9OhQ4e001iq3PmUcn/Fjt3ceM3ZvpBts/bdypu8fL5ZPI6snIfazpnDdy65hM5PP83nm2zCqyedxLyNNipK7GLmWY4Y5d4ub/r16zc5xtjwvw5ijD6We/Ts2TPmxdixY9NO4SvKnU8p91fs2M2N15ztC9k2a9+tvMnL55vF40j9PLRkSYw33xzjWmvF2K5djMOGxbh4cXFiN0MxYhUao9zb5Q0wKa6gpnFlDEmSyuXdd+HQQ+Ghh5JJFtdfD5tsknZWyjHv0ZMkqdRihGuuSWbUjhsHl10GTz1lkaeSc0RPkqQSWnXWLNhjDxg7Nml4fO210MR78aRCWehJklQKdXVwxRX0PeUU+MY3khG9gw+GENLOTK2IhZ4kScX2yiswZAg89xyfbbcd64wcCeutl3ZWaoUs9CRJKpbaWrjggqQfXmUl3HIL09Zdl90t8pQSJ2NIklQMU6fC9tvD734HAwbAjBnwy196qVapstCTJKk5Fi1KRvC23jppn3L33TByJHTpknZmkpduJUkq2KRJcNBBMH16Mnp3ySXJWrVSRjiiJ0lSU82fDyefDNttB59+Cg88ALfcYpGnzHFET5Kkpnj2WRg8GF57DYYOhfPPhzXWSDsrqUGO6EmS1BjV1XDssbDLLsl9eY8+mvTGs8hThjmiJ0nSyjz+eNLs+O234aij4JxzoEOHtLOSVsoRPUmSVmTuXDjkEOjfH9q2Tdanvewyizy1GBZ6kiQ15MEHoVcvuP56OOkkePFF2HnntLOSmsRCT5KkZX3yCfzqV0nT4zXXhPHjYdgwWG21tDOTmizEGNPOITNCCAOBgV27dh06YsSItNMpiurqajpk6BJDufMp5f6KHbu58ZqzfSHbZu27lTd5+XyzeBxfl1OncePoeemltPn8c975xS94+xe/ILZtW7L9NVcxYxcjVqExyr1d3vTr129yjHGbBl+MMfpY7tGzZ8+YF2PHjk07ha8odz6l3F+xYzc3XnO2L2TbrH238iYvn28Wj6PBnD78MMb9948RYtxqqxhffLG0+8tg7GLEKjRGubfLG2BSXEFN46VbSVLrFSPceitUVSVNj889F55/Hvr0STszqShsryJJap3eew8OOyyZdLHDDnDDDbDppmlnJRWVI3qSpNYlRrjuumRG7RNPJOvTPv20RZ5yyRE9SVLr8dZb9DnhBHjhBejXD669FjbeOO2spJJxRE+SlH9LlsDll8Pmm7P6K6/A8OHw2GMWeco9Cz1JUr699hrsthsccwzsuisTb7wxWe1iFf8XqPzzWy5JyqfaWjj//GQG7Usvwc03w4MPsvCb30w7M6lsvEdPkpQ/06bB4MEwaRLsuy9ceSV06ZJ2VlLZOaInScqPRYvgzDNh663h7bfhrrvgnnss8tRqOaInScqHyZOTUbypU+H//g8uvRQ6dUo7KylVjuhJklq2BQvgd7+D7baDjz+GUaPgttss8iQc0ZMktWTPPZeM4r36KgwZAhdcAGuumXZWUmY4oidJannmzYPjjoOdd05G9B55JFntwiJP+gpH9CRJLcsTT8DQofDmm3DkkXDuubD66mlnJWWSI3qSpJbh88/hsMNgjz2SZsfjxsEVV1jkSV/DQk+SlH0PPQS9eiVr055wAkyZArvumnZWUuZZ6EmSsmvOHDjgANh7b1hjDfjnP5PVLior085MahEs9CRJ2XTffVBVBSNGwB/+kPTJ23bbtLOSWhQnY0iSsuWjj+Coo+Duu2HLLeHhh2GLLdLOSmqRHNGTJGVDjMnoXVVV0vT47LPh+ect8qRmcERPkpS+99+Hww+HBx6A7beH669PCj5JzdKkQi+EUAXsAWwLdAFWBeYArwHPAI/EGOcXO0lJUk7FCDfeCMcfD4sWwUUXwTHHQEVF2plJubDSS7ch8esQwkRgOnAGsB7wGfA2UAHsBdwDfBhCuCaEsGHpUpYk5cLMmbDnnsnSZVtsAVOnJqtdWORJRRNijF//hhBeqf/xFuCeGOMrK3hfJbAn8BNgH+CwGOOtRcy15EIIA4GBXbt2HTpixIi00ymK6upqOnTokHYaS5U7n1Lur9ixmxuvOdsXsm3Wvlt5k5fPt8HjWLKEbqNGsfE11xBD4M3DDmPWgAFJE+S0cmqh+ytm7GLEKjRGubfLm379+k2OMW7T4Isxxq99APtTXxA29kEy4rd9U7bJ0qNnz54xL8aOHZt2Cl9R7nxKub9ix25uvOZsX8i2Wftu5U1ePt//Oo7XXotxl11ihBj33DPGt98uWy5fLFgcb3/+7Xj7qIfi7c+/Hb9YsLgs+20p56FixCo0Rrm3yxtgUlxBTbPSfz7FGEfWB2m0GON7McbxTdlGkpRjdXVwwQXQuzdMm5bcl/fQQ9C9e1l2P3HmHLY75zHOHD2Dj79YyJmjZ7DdOY8xceacsuxfSovtVSRJpfXSS7DjjnDiick9eTNmwIEHQghl2X31wloOvHEC8xbWUbOoDoCaRXXMW1hX/3xtWfKQ0rDSWbchhAnAgTHGGfUTMr52dC/GaNtySRIsXswGt9wCt9ySLF92xx3w05+WrcD70ugps1jRdakYYfTUWQzqW56RRancGtNe5SVg/jI/N+kyriSpFXrhBRg8mA2nTIGf/Qwuuww6d04llZmfzFs6kre8mkV1zJxdU+aMpPJZaaEXYzxomZ8PLGk2kqSWbcECOOssGDYMvvlNpp11FpufdlqqKfVYpz2VbSsaLPYq21bQo1NlCllJ5eE9epKk4hg/HrbaCs45B379a3jpJT7Zeee0s2JAn24rvFocAgzo3a28CUll1JiGydeHEL7T2IAhhG+EEAaHEH7VvNQkSS1CTU2yssWOO8K8efDww3DDDbDWWmlnBkCHdm246aBtad+ugsq2STPmyrYVtG9XUf+8q4Eqvxrz7a4BpoQQJgMjgeeA6XGZpc5CCBsAW5OskPEj4H3gkOKnK0nKlCefhIMPhn//G444Av7yF1h99bSz+i99e6zNhFP7M3rqLOIHM/jjwCoG9O5mkafca0wfvaOBTYCngN8CzwPVIYR5IYQ5IYRa4E3gDpL1b4fGGLeIMU4oYd6SpDR9/jkcfjj065f8/uST8Ne/ZrLI+1L7dm0Y1Lc7XTquyqC+3S3y1Co06lseY3wX+D3w+xBCT5LRuy7AqsAc4FVgQozRqUuSlHcPPwyHHALvvZdcsj3rLKh0QoOURY3+50wI4RvAtsBbMcbbS5eSJCmTPv00Kexuugm++1147jnYfvu0s5L0NZoy67YOeAL4bolykSRl1ahRUFWVND/+/e/hX/+yyJNagEaP6MUYl4QQXge+VcJ8JElZ8vHHcPTRcOed0KcPjBkDW26ZdlaSGqmpffR+D5weQti8FMlIkjIixmTJsqoquO8++POfYeJEizyphWnqlKPTgHWAF0MI7wP/Ybkl0VzrVpJauFmzklYpo0bBttsmPfF69Uo7K0kFaGqh9xIwvRSJSJJSFmMy0eK442DhQrjgAvjNb6CiIu3MJBWoSYWea91KUk69/XbSMuWRR2DXXeG66+A7jV4USVJGrbTQCyGc3oR4McZ4VjPykSSV05IlMHw4nHRSMqL317/CYYfBKi6FLuVBY0b0zgDmA/OAFSwLvVQELPQkqSV44w0YMgSeegq+/3245hro0SPtrCQVUWMKvTeB7sBkkmXO7osxfl7SrCRJpVNXB5deCqedBm3bJpMtDjwQwsr+LS+ppWnMWrffBnYkmYhxFvBhCOHeEMJPQgirlTpBSVIRzZgBO+0Ev/0t9O+f/H7QQRZ5Uk416iaMGOOkGOMJMcbuwA+BD4ErgI9CCLeFEHYtZZKSpGZavBjOPjvpg/fGGzBiRNI+pVu3tDOTVEJNba9CjPEp4KkQwm+As4HjgNWAp4qcmySpGF58MRm1e/FFGDQILrsMvvnNtLOSVAZNLvRCCDsBPwP2B1YHRgJXFTkvSVJzLVwIZ50Fw4ZBp05w772w775pZyWpjBpV6IUQtiIp7gaRrHX7MMlI3t9jjDWlS0+SVJDnn4fBg5N78A44AC66CNZeO+2sJJVZY/rovQpsCDwB/BG411m3kpRRNTVw+ulw8cWw7rowZgzstVfaWUlKSWNG9L4DLAC2BrYCzgtfMzsrxuiNH5KUhqeeSvrivfFG0vR42DDo2DHtrCSlqDGF3p9KnoUkqXBffAGnnAJXXgkbbQRPPAH9+qWdlaQMWGmhF2O00JOkrHrkERg6FN59F447Lpl80b592llJyggXM5Skluizz5LLtHvuCZWV8OyzyYQLizxJy7DQk6SW5u9/h6oquPlm+N3v4F//gh12SDsrSRnU5D56kqSUzJ4NxxwDt98OvXvD6NGw1VZpZyUpwxzRk6SsixHuuisZxRs5Es48EyZOtMiTtFKO6ElSln3wARxxBNx/P/TtCzfcAJttlnZWklqI3I/ohRDGhRCmhBCmhhBGhhBsKiUp+2JM7sGrqoKHH4bzz4fnnrPIk9QkuS/0gH1ijH1ijL2Bd4AT005Ikr7WO++w+SmnwIEHJoXdlClwwgnQxoswkpqm7IVeCOHbIYTh9aNsdSGEJ1fwvqoQwuMhhJoQwqwQwpkhhIqm7i/GOLc+3ipAeyA26wAkqVSWLIHhw2GzzVhz6lS4/HIYNw569kw7M0ktVBr/POwF7A2MB9o29IYQwlrAY8AM4H+BjYELSQrT05q6wxDCGKAv8BLw24KylqRS+ve/4eCD4cknYY89mDhkCNv//OdpZyWphUvj0u0DMcb1Y4w/ISm8GnIYsBqwX4zx0Rjj1SRLsR2/7D12IYRnQggzG3hcv2ywGOPeQBdgAnBESY5KkgpRVwcXXwybbw4vvADXXQePPsqCrl3TzoH8r84AACAASURBVExSDpR9RC/GuKQRb9sL+EeM8fNlnrsDGAbsBjxQH2vnJuy3LoRwM3AncF7jM5akEnn5ZRg8GMaPhwED4OqrYd11085KUo5kdTLGpsAryz4RY3wHqKl/rVFCCGuFEL61zFM/BqYXJUNJKtTixXDuubDFFvDaa3DrrclqFxZ5koosxJje3IQQwkigU4xx9+WeXwycGGO8ZLnn3wP+FmM8tZHxNwLuIrkXMAAvA0fHGP/TwHsPAQ4B6Ny589Z33XVX0w8og6qrq+nQoUPaaSxV7nxKub9ix25uvOZsX8i2WftutRTt33iDTc87j9Vff52PdtuN1485hsVrr/1f78vL55vF4/A8VLpYhcYo93Z5069fv8kxxm0afDHGmNoDGAk82cDzi4FjG3j+feDsUufVs2fPmBdjx45NO4WvKHc+pdxfsWM3N15zti9k26x9tzJvwYIY//CHGNu0ifFb34px5MivfXtePt8sHofnodLFKjRGubfLG2BSXEFNk9WmTJ8Cazbw/BrAZ2XORZKaZ8KE5F68l16CX/86mXzRwCieJBVbVu/Re4Xl7sULIaxP0gfvlQa3kKSsmT8fTjoJdtgB5s6FBx9MVruwyJNUJlkt9B4C9gwhrL7Mc4OA+cC4dFKSpCZ4+mno0ydZuuzgg2H6dNh777SzktTKpLEyRmUIYf8Qwv7AukDnL38PIVTWv+1qYCFwbwihf/1EiTOAi+JXW65IUrZUV8PRR8Ouu0JtLTz2WLLaxRprpJ2ZpFao7LNuQwg9gLdW8PKGMcaZ9e+rAq4AdiC5L+864IwYY10JcxsIDOzatevQESNGlGo3ZZW1GUnOditdPGfdpm+tyZPZ5PzzaffRR7y/3368OWQIS1ZbraBYefl8s3gcnodKF8tZt+nI7KzbrD6cdVs6znYrXTxn3abo009jHDIkRoixZ88Yn3mm2SHz8vlm8Tg8D5UulrNu08HXzLrN6j16ktQyjB4NvXrBjTfCySfDiy/CTjulnZUkAdmdjCFJ2TZ7NvzylzBwIKyzDjz/PPzlL1DgpVpJKgULPUlqihjh7ruhqgruvBPOOAMmTYJtGr49RpLSlNWGyZKUPR9+CEceCffeC1tvncyo7d077awkaYUc0ZOklYkR/va3ZBTvwQdh2DAYP94iT1Lmlb29SpbZXqX0bGtQuni2VymNdh99RM+LLmKd559nbq9evHLSSczv3r3k+83L55vF4/A8VLpYtldJh+1VbK+SGbY1KF0826sU2ZIlMQ4fHuPqq8dYWRnjpZfGWFtbtt3n5fPN4nF4HipdLNurpAPbq0hSE7z5JuyxBxx6KPTtC9OmwTHHQEVFyXddvbCWOya8w4efL+COCe9QvbC25PuUlF8WepL0pbo6uPRS2HzzZCbtNdckEy422qgsu584cw7bnfMYZ46ewcdfLOTM0TPY7pzHmDhzTln2Lyl/LPQkCeCVV5L1aX/zG9h9d3jpJRg6FEIoy+6rF9Zy4I0TmLewjppFyUqPNYvqmLewrv55R/YkNZ2FnqTWrbY2aXS8xRbw8svJ7NrRo2H99cuaxugps1jR3LgYYfTUWWXNR1I+2EdPUus1dSoMHgyTJ8N++8Ff/wpduqSSysxP5i0dyVtezaI6Zs6uKXNGkvLAET1Jrc+iRcmKFltvDe++m6x0cc89qRV5AD3WaU9l24Yne1S2raBHp8oyZyQpDyz0JLUukyYlBd6f/gSDBiX34u2/f9pZMaBPtxXeDhgCDOjdrbwJScoFGyYvw4bJpWej0tLFs2Hy11tl4UJ63HQT6991F4vWXpvXjj+eT3bYIe20vqJmUR1vzZ4HQOdVIx8vSCq/DTuteLQv67L4PfE8VLpYNkxOhw2TbZicGTYqLV08GyZ/jaefjrFnzxghxoMPjvHTT9POaIWqFyyOd0x4O94+6qF4x4S3Y/WCxWmn1CxZ/J54HipdLBsmpwMbJktqlaqrk0bHu+6a3Jf36KNw7bWw5pppZ7ZC7du1YVDf7nTpuCqD+nanfTvnzEkqnIWepHx6/PGk8fHll8NRRyWrW/Tvn3ZWklRWFnqS8mXuXDjkkKSo+8Y34Kmn4LLLwPt4JLVCFnqS8uPBB6FXL7j+ejjxRJgyBXbZJe2sJCk1FnqSWr5PPoFf/QoGDEjuvxs/Hs47D1ZbLe3MJClVFnqSWrZ77oGqKrjjDjj99GSVi759085KkjLB6VySWqb//AeOPDIp9LbaCh55BPr0STsrScoUGyYvw4bJpWej0tLFazUNk2PkW489xrevuIKK+fOZeeCBvDtoELGiZTYUXpGs/d0tVBaPw/NQ6WLZMDkdNky2YXJm2Ki0dPFaRcPkd9+N8X/+J2l8vMMOMc6YUf4cyiRrf3cLlcXj8DxUulg2TE4HNkyW1KLFmDQ67tULnngCLr4Ynn4avvvdtDOTpEzzHj1J2fbWWzB0aNIAeffd4brrYOON085KkloER/QkZdOSJcmqFpttBhMmwNVXJ8WeRZ4kNZojepKy59VXYcgQePZZ+OEPYfhw6N497awkqcVxRE9SdtTWJo2O+/SBl16Cm26CMWMs8iSpQI7oScqGadNg8GCYNAl+9CO48kro2jXtrCSpRXNET1K6Fi2CP/0Jtt4a3n4b7rwT7r3XIk+SisARPUnpmTwZDjooGc37+c/h0kuhc+e0s5Kk3HBET1L5LVgAp5wC220Hs2fDqFEwYoRFniQVmUugLcMl0ErPpYdKF6+lLIHWcfp0Nj3vPCrffZcP9tqLfx9xBLUZ+o5mRdb+7hYqi8fheah0sVwCLR0ugeYSaJnh0kOli5f5JdCqq2M89tgYQ4ixe/cY//GPJu+zNcna391CZfE4PA+VLpZLoKWDr1kCzXv0JJXeE0/AwQcnq1wceSScey6svnraWUlS7nmPnqTSmTsXDj0U9tgDKipg3Di44gqLPEkqEws9SaUxZkyyfNl118FvfwtTpsCuu6adlSS1KhZ6koprzhw44AD4n/+Bjh3huefgggugsjLtzCSp1bHQk1Q8994LVVVw221w2mnwwgtJCxVJUiqcjCGp+f7zHzj6aLj7bthiC3joIdhyy7SzkqRWzxE9SYWLMWl03KtX0vT4z3+GCRMs8iQpIxzRk1SQth9/DPvsA6NHJ5dnb7ghuWwrScoMCz1JTRMj3HAD2x57LNTVwYUXwrHHJu1TJEmZYqEnqfFmzoShQ+Gxx6ju3Zs1R46E73wn7awkSStgoSdp5ZYsgSuvhFNOIYbA5JP+zBs7bAmftmPAwlo6tPNUIklZ5GQMSV/vtddgt93g6KOZu9W2fH/wFfy63dZ8PG8xZ46ewXbnPMbEmXPSzlKS1ICQrIUrgBDCQGBg165dh44YMSLtdIqiurqaDh06pJ3GUuXOp5T7K3bs5sZrzvYNbRvq6ljv7rvpceONLGnbljcOP4In+uzIkvrXv7Ua/Gd+8vMqIfDdrh1ZJRScvpaTtb+7hcricXgeKl2sQmOUe7u86dev3+QY4zYNvhhj9LHco2fPnjEvxo4dm3YKX1HufEq5v2LHbm685mz/X9tOmxZj374xQoz77BPj++/H259/O373Dw/FDU4eHTc4eXS87Nb7l/783T88FO+Y8Haz8tdXZe3vbqGyeByeh0oXq9AY5d4ub4BJcQU1jZduJf1/ixfDWWfBVlvBW2/B7bfD/fdDt27M/GQeNYvqGtysZlEdM2fXlDlZSdLKeAe1pMQLL8DgwTBlCvzsZ3DZZdC589KXe6zTnsq2FQ0We5VtK+jRybVsJSlrHNGTWrsFC9jw2mth222Tpczuuy8ZyVumyAMY0KcbYQX34IUAA3p3K0OykqSmsNCTWrN//hO23JINRoyAX/0KZsyAH/2owbd2aNeGmw7alvbtKqhsmzRHrmxbQft2FfXPe4FAkrLGM7PUGtXUwGmnwSWXwHrrMWXYMPqcdNJKN+vbY20mnNqf0VNnET+YwR8HVjGgdzeLPEnKKEf0pNZm7FjYfHO4+GI47DCYPp1Pt9220Zu3b9eGQX2706Xjqgzq290iT5IyzEJPai0+/xwOPxy+973kprqxY5PVLjp2TDszSVKJWOhJrcHDD8Nmm8Hw4XD88TB1Kuy+e9pZSZJKzEJPyrNPP4WDDoK99oIOHeDZZ+HCC6HSViiS1BpY6El5df/9UFUFt9wCp56a9MnbYYe0s5IklZF3UUt58/HHcPTRcOed0KcPPPhgstKFJKnVcURPyosY4Y47klG8e+9NljKbONEiT5JaMUf0pDyYNYvNTjsNnnsO+vaFG2+EXr3SzkqSlDJH9KSWLMakqKuqYq1Jk+D885NizyJPkoSFntRyvf02/PCHMHgw9O7NpOuvhxNOgDYO1EuSEhZ6UkuzZEnS6HizzZJ2KVdcAU8+yfz11ks7M0lSxoQYY9o5ZEYIYSAwsGvXrkNHjBiRdjpFUV1dTYcOHdJOY6ly51PK/RU7dmPirfb++2xy3nmsOXUqc7bZhtd++1sWdOnS7HwK2TZr3628ycvnm8Xj8DxUuliFxij3dnnTr1+/yTHGbRp8McboY7lHz549Y16MHTs27RS+otz5lHJ/xY79tfFqa2O84IIYV1stxjXWiPH662NcsqRo+RSybda+W3mTl883i8fheah0sQqNUe7t8gaYFFdQ03gzj5R1M2Yk9+E9/zwMHAhXXQXrrpt2VpKkFsB79KSsWrwYzj4bttwS3ngDRoyAUaMs8iRJjeaInpRFL76YrFH74ovw05/C5ZfDN7+ZdlaSpBbGET0pSxYuhNNOS5oef/BBssLFnXda5EmSCuKInpQRq8+YAUcemdyTd8ABcNFFsPbaaaclSWrBLPSktNXUwOmns9XFF0O3bjBmDOy1V9pZSZJywEJPStO4cXDwwfDGG3wwcCDdbr0VOnZMOytJUk5Y6Elp+OILOOWUZIWLjTaCJ57gtRDoZpEnSSoiJ2NI5fbII8nyZVddBb/5DUydCv36pZ2VJCmHLPSkcvnsMxgyBPbcEyor4Zln4OKLoX37tDOTJOWUhZ5UDn//O1RVwc03w+9+B//6F+y4Y9pZSZJyznv0pFKaPRuOOQZuvx1694YHHoCtt047K0lSK+GInlQKMcJddyWjeCNHwp/+BBMnWuRJksrKET2p2D74AI44Au6/P1nh4oYbkskXkiSVmSN6UrHEmNyDV1UFDz8M550Hzz1nkSdJSo0jelIxvPMOHHpoUuDtvDNcfz307Jl2VpKkVs4RPak5liyBq6+GXr3g6afh8suT1S4s8iRJGeCInlSgVd9/H/bYA558MvnvtdfChhumnZYkSUtZ6ElNVVcHl19O31NOgXbtkgJvyBAIIe3MJEn6Cgs9qSlefjkp6v75Tz7dYQc63XUXrLde2llJktQgCz2pMRYvhgsugDPOgA4d4NZbmd6tG7tb5EmSMszJGNLKTJkC228Pp54K++wDM2bAL37hpVpJUuZZ6CmXqhfWcseEd/jw8wXcMeEdqhfWNj3IwoVw+umwzTbw/vvJChd33w3f+lbxE5YkqQRCjDHtHDIjhDAQGNi1a9ehI0aMSDudoqiurqZDhw5pp7FUOfKpWVTHW7PnAdB51cjHC5KRtw07taeybUWjYqz+8stset55tJ85kw+//33eOPJIatdY4yvvKfaxNDdec7YvZNusfbfyJi+fbxaPo9w5lXJ/xYxdjFiFxij3dnnTr1+/yTHGbRp8McboY7lHz549Y16MHTs27RS+otT5fLFgcaw6/aG4wcmj4wYnj46X3Xr/0p+rTn8oVi9Y/PUBampiPPHEGFdZJcb11ovxwQdX+NZiH0tz4zVn+0K2zdp3K2/y8vlm8TjKnVMp91fM2MWIVWiMcm+XN8CkuIKaxku3ypXRU2axokHqGGH01Fkr3viZZ6BPHzj/fDj4YJg+HfbeuzSJSpJUBhZ6ypWZn8yjZlFdg6/VLKpj5uya/36huhqOPhp23RVqa+Gxx2D4cFjuUq0kSS2NhZ5ypcc6K74Pr7JtBT06VX71ycceg803h7/+FY45BqZNS1a5kCQpByz0lCsD+nRbYdeTEGBA727JL3PnwtCh8P3vQ9u2yTq1l1wC7duXL1lJkkrMQk+50qFdG246aFvat6tYOrJX2baC9u0q6p9vA6NHQ1UV3HADnHwyvPgi7LRTyplLklR8royh3OnbY20mnNqf0VNnET+YwR8HVjGgdzfaV8+FX/4SbrstuVw7alTSI0+SpJxyRE+51L5dGwb17U6XjqsyqG932j9wfzKKd+edyTJmkyZZ5EmScs8RPeVa2zlzYP/94Z57YOut4dFHoXfvtNOSJKksLPSUTzHCrbfS96ijkqXMhg2D44+HNn7lJUmth//XU/68+y4cdhiMGUNNr16scc89sMkmaWclSVLZWegpP2KEa6+FE06Aujq49FL+1asXu1vkSZJaKSdjKB/efBP694dDD4W+fZPGx8ccAxUNN0+WJKk1sNBTy7ZkCVx2WdIuZeJEuOaaZLWLjTZKOzNJklLnpVu1XK++CkOGwLPPwt57w9VXw/rrp52VJEmZ4YieWp7a2mQWbZ8+MGMG/O1vyWoXFnmSJH2FI3pqWaZNg4MOgsmTYb/94K9/hS5d0s5KkqRMckRPLcOiRcmKFltvnbRPufvupAmyRZ4kSSvkiJ6yb9IkGDw4Gc375S/hkktgnXXSzkqSpMxzRE/ZNX8+nHIKbLcdfPIJPPAA3HKLRZ4kSY3kiJ6y6dlnk1G8116Dgw+G88+HNddMOytJkloUR/SULfPmwbHHwi67JPflPfpostqFRZ4kSU1moafsePzxpPHxZZfBUUcl9+T17592VpIktVgWekrf3LnJ0mX9+0ObNvDUU0mx16FD2plJktSiWegpXQ8+CL16wXXXwUknwZQpyWVbSZLUbBZ6Sscnn8Cvfw0DBiT3340fn6x2sdpqaWcmSVJuWOip/O65B6qq4Pbb4fTTk1Uu+vZNOytJknLH9ioqn//8h6ozzoBx42CrreCRR5L1aiVJUkk4oqfSixFuuw2qquj03HNw7rnw/PMWeZIklZiFnkrr/fdhn32Spcs22YRJ116brHbRxsFkSZJKzUJPpRFjMpO2qirpj3fxxfD009RssEHamUmS1Go4rKLie+stGDo0KfB23z0p+DbeOO2sJElqdRzRU/EsWQKXX56sbjFhAlx9dVLsWeRJkpQKR/RUHK+9BkOGwDPPwA9/CMOHQ/fuaWclSVKr5oiemqe2Fs4/P5lBO3063HQTjBljkSdJUga0mkIvhHBlCCGmnUeuTJ8OO+6YLF22114wYwYccACEkHZmkiSJVlLohRB2AdqnnUduLFoEZ56ZND2eORPuvDNZ7aJr17QzkyRJyyh7oRdC+HYIYXgIYUoIoS6E8OQK3lcVQng8hFATQpgVQjgzhFBRwP7aAX8BTmhm6oL/v1zZH/8IP/lJMor30586iidJUgalMaLXC9gbeK3+8V9CCGsBjwER+F/gTOC3wJ8K2N/pwPUxxo8LylaJBQvg1FNhu+3g449h1KhktYtOndLOTJIkrUAas24fiDGOAgghjAQaqhQOA1YD9osxfg48GkLoCJwRQjiv/jlCCM8A6zWw/eMxxiEhhN7AdsBppTiQVuO555IZta+8AoMHw4UXwpprpp2VJElaibKP6MUYlzTibXsB//iyoKt3B0nxt9sysXaOMfZo4DGk/i07AVXAWyGEmQAhhJkhhM5FOZi8mzcPjjsOdt4Z5s+Hf/wDrr/eIk+SpBYiq5MxNgVeWfaJGOM7QE39a40SY7wqxtjtywKw/rkeXsZthLFjoXdvuOQSOOIImDYNfvCDtLOSJElNEGJMr+PIl5duY4y7L/f8YuDEGOMlyz3/HvC3GOOpBe4vxhgbnDUQQjgEOASgc+fOW991112F7CJzqqur6dChQ6PfXzFvHhsPH063Bx6gZt11efXEE5nbp09q+WR5f8WO3dx4zdm+kG3L/WfZ2uTl883icXgeKl2sQmOUe7u86dev3+QY4zYNvhhjTO0BjASebOD5xcCxDTz/PnB2qfPq2bNnzIuxY8c2/s1jxsS43noxrrJKjCecEOO8eenmk/H9FTt2c+M1Z/tCti33n2Vrk5fPN4vH4XmodLEKjVHu7fIGmBRXUNNk9dLtp0BDN4KtAXxW5lzyb84cOPBA2Htv6NgxmXxx/vlQWZl2ZpIkqRmyWui9wnL34oUQ1idpevxKg1uoMPfdB1VVcOutcNpp8MILSQsVSZLU4mW10HsI2DOEsPoyzw0C5gPj0kkpZz76CAYNgv32S1a0mDgRzjoL2rVLOzNJklQkZe+jF0KoJGmYDLAu0DGEsH/972NijDXA1cAxwL0hhGHARsAZwEXxqy1XWrTqhbWMnjKLmZ/Mo8c67RnQpxsd2pX4jyRGuP12OOYY+OILOPtsOPFE+MY3SrtfSZJUdmWfdRtC6AG8tYKXN4wxzqx/XxVwBbADyX151wFnxBjrSpjbQGBg165dh44YMaJUuwGgZlEdb82eB8CSGFmlfgmxDTu1p7Jtk1d6W6FlZyS1/fhjel5yCZ2ee47Pv/tdXjnpJGp69CjavpqaT0vfn7Nune1WSnn5fLN4HJ6HShfLWbfpyOys26w+Sj3r9osFi2PV6Q/FDU4e/V+PqtMfitULFhdtX2PHjo1xyZIYr78+xjXWiHG11WK88MIYa2uLto8m55OT/TnrtvD9aeXy8vlm8Tg8D5UulrNu00ELnHWba6OnzGJFA6kxwuips4q2r3Yffgh77pksYbbFFjB1Khx/PFQUb9RQkiRlUxpr3bZ6Mz+ZR82ihq9A1yyqY+bsmubvZMkSuOoq+p54YlLUXXklHHoorGJtL0lSa2Ghl4Ie6yT34TVU7FW2raBHp2b2r3v99WQE7+mn+bxvX9a++27YYIPmxZQkSS2OwzspGNCnG6HBhdggBBjQu1thgevq4MILkzVqp02DG29k6rBhFnmSJLVSFnop6NCuDTcdtC3t21UsnWFb2baC9u0q6p8vYKD1pZdgxx3hhBOSe/JeeilZ7WJFFaUkSco9L92mpG+PtZlwan9GT53FzNk19OhUyYDe3Zpe5C1eDMOGwZlnwhprwB13wE9/aoEnSZLK30cvy8rZR68YOrz+Opucdx6rv/EG//ne93jj6KNZvOZXlwjOWo8h+1eVLp599PIlL59vFo/D81DpYtlHLx320ctYH71mW7Agxt//PsaKihi7dInxvvtW+Nas9Riyf1Xp4tlHL1/y8vlm8Tg8D5Uuln300sHX9NFzRK8BIYS5wOtp51EkawBz005iGeXOp5T7K3bs5sZrzvaFbNsJmF3g/rRyWfu7W6gsHofnodLFKjRGodt5HkpsEGPs3NAL3qPXsDtjjIeknUQxhBCuydKxlDufUu6v2LGbG6852xeybQhhUlzRpQI1W9b+7hYqi8fheah0sQqN0YztPA+thLNuG/ZA2gkUUdaOpdz5lHJ/xY7d3HjN2T5r3xPl588ki8fheah0sQqNkcXvSS546VZSQfyXtKS0eR5aOUf0JBXqmrQTkNTqeR5aCUf0JEmScsoRPUmSpJyy0JMkScopCz1JRRdCWD+E8HgI4eUQwkshhPNCcF0+SeUTQhgXQpgSQpgaQhgZQuiYdk5psNCTVAq1wMkxxu8CWwLbAfulm5KkVmafGGOfGGNv4B3gxLQTSoOFnqSlQgjfDiEMr/9XcF0I4ckVvK+qfsSuJoQwK4RwZgih4svXY4wfxBgn1f+8CJgKrF+Wg5DUYhXrHAQQY5xb/95VgPZAq5x96soYkpbVC9gbGA+0begNIYS1gMeAGcD/AhsDF5L8w/G0Bt6/DvAj4AelSVlSjhT1HBRCGAP0BV4CfluyrDPM9iqSlgohrBJjXFL/80igU4xx9+Xe8zvgJJK1FT+vf+4k4Aygy5fP1T/fDngYGB1jvLAsByGpxSr2Oaj+tQrgXGB2jPG8kh9ExnjpVtJSX55gV2Iv4B/LnUzvAFYDdvvyifqT623AvyzyJDVGMc9By8SsA24Gfl2UJFsYCz1JTbUp8MqyT8QY3wFq6l/70nDgC1rp5RJJJbPSc1AIYa0QwreWecuPgellyzBDvEdPUlOtBXzWwPOf1r9GCGEnYAjJifVf9Z1VbogxXlauJCXl1krPQfX/vSuE0BYIwMvA0eVJL1ss9CQVoqGbe8OXz8cYn63/XZJKYWXnoDeBbcqaUUZ56VZSU30KrNnA82vQ8L+yJamYPAc1gYWepKZ6ha/ei0cIYX2SPlWvNLiFJBWP56AmsNCT1FQPAXuGEFZf5rlBwHxgXDopSWpFPAc1gffoSVoqhFBJ0qwUYF2gYwhh//rfx8QYa4CrgWOAe0MIw4CNSPpXXbR8/ypJagrPQcVnw2RJS4UQegBvreDlDWOMM+vfVwVcAexAck/MdcAZ9f2qJKkgnoOKz0JPkiQpp7xHT5IkKacs9CRJknLKQk+SJCmnLPQkSZJyykJPkiQppyz0JEmScspCT5IkKacs9CQ1W0iMCiHcl3YukqT/z0JPagFCCGeEEGL9Y0kI4dMQwsQQwtkhhC7LvbdH/fsGlCu/mHRe/xVQFUI4oKH3hBAOCSH8qFw5Lbfvzeo/k93T2H9DQgg3hRAmNeH9Z4QQZpcyp1IKIVSGED4MIeyWdi5NFUJYLYTwUQhhl7RzkZrKQk9qOeaSLPezI/Az4F6S4mpaCGHrZd73Qf37nilncvVrTO4LnB1C6N7AWw4BUin0Muos4MAmvP86YM/SpFIWRwNvxRhb3KLzMcb5wOUkf2ZSi9Im7QQkNVptjHH8Mr//I4RwFfAUcGcIYZMYY12McSEwvuEQpRVjnAGs15wYIYTV6v/Hmmsxxn838f3vAe+VKJ2SCiGsAhxJyy6UbgL+FELYPMY4Le1kpMZyRE9qwWKMnwEnARsDnzLP7gAACYNJREFU34eGL92GEH4dQngmhPD/2jn3YK+qKo5/vogJJiI+0mJMMgdzcpx8lo+KsPKB6eQDChsBsyYVHR/pOD7oyvRSLC3TUAwIc7yiSWgCSTwyX3FxBDIDZFB8ogKCCmEoqz/WPrI5/l5X7gz33tZn5szvd9bZ771/c9Zvrb32quT2nSXpkLyswpUoaYCkpyWtk/SApJ0l7ZPyrE1pDijl7SLpMklLJL0jabGkM7Pns4GDgSGZC3poevacpF9IukrSi8CbWb6jJP0ttWWlpDGSetQbF0nnSHohtfd+4OOl5xXd2/XcqZKuTu7HLiX5Cam8fdL9NsnV+nwaj39JGlyrLkk7SbpN0suS1qe8Y7Lnm7luJfUr3NGS7pb0tqSlks4p1XO4pPtSuWslzZN0eoW+9ZO0INXdIukwSSskNZXSnZTWwPo0FtdK2rbamCX6A71xK3ReVqV1M6SUZrakeyQNTunelDRVUs0/FJKGpvE5KJWxLvX9IEkflTRO0po0Zt+u037M7AWgBTijXtogaE+EohcEHZ9ZwLvAF2qk6QNMAE4DBuOWoYck7V1K90lgJHAl7mo9ArgVaE7XqbgnoFmSsnw3pjy3AicAk4HbJJ2Unp8DLASm4G7lw4EHsvyDgS+ndIMAJB0JzACWp3ovAI4HxtUajFTnTcCfgZOBfwJja+VpBc3A7qmtOQOBJ8xsSbofCVyBj8eJwCPAHXUUil8CRwEX4i7aywFroE1jgPm423w2cJOkw7Lne6X6zwK+AfwRGJe3RVJvfG5ew8f6FuAOoHtekaSBuLI2J/Xranyd/KxOG48GFpvZypI8XzcDgEnA2LICDnweGA5cnOo7KOVphN8DdwKnAALuAX4HvIz39R/AhHqKY+JR4KsN1hsE7QMziyuuuNr5BTQBK2o8fwX4bfreB1cQTqiStguurC0ERmTy8bjC+OlMdm0q64xMdnyS7Zfu9wE2AsNK9dyBKz/F/VxgfIX2PJfa360k/zswqyTrn+rev8ZYzAGmlmRjUr5+tcYojcHcOnMxHxid3W+H75/8YbrfGVgL/KiUbwqwqFpdwFPAeY2uAaBf6sPITLYt8Drw8yplKM39LcDMTD4KWAF0z2QDU/lNWd5lwLhSmWcC/wF2qdH2B4G7S7Ji3QwpyScALdn97DS+vTLZBalt3WvUOTSlGZLJirU7NpP1BDYAZzfwOxyK/0a61UsbV1zt5QqLXhB0DlTzobSfpEmSXgXew19s+wJ9S0mfs833jhUWqpkVZL3T59H4C3tiqay/AAdIamQv8AwzW5+1d3vc6jdRUtfiwgNMNuBu4Er93AY4ELco5txbIfmH5S7glKxfxwE92NT//YHtgbsr5Osr6WNVyp0HXJLczuV5qcWDxRcz2wA8Q7ZPUlIvSb+WtAwfuw24VSyv41Bgum2+N/K+Uj19cYtveU5mAt3wfldjD1yRzCnWzaRSeTOAz6W5LGgxszey+6fTZ2/qMyP7/oH1bGZrcOW4kbJWANsAuzWQNgjaBaHoBUEHR1I3YBfg1SrPe+DKwJ7ARcAX8Rf7fPwFnbO6dP/fCvJCVuTdFX/5rUz7ttZLWo9HiXbFX/L1KLe9VyrzZjYpJxuAd3Cr1Z5Vytkt1flaSV6+3xKa8T73T/eDgMfM7Pl0X+wHLPepuO9VpdzhwJ+AEcAiSc9I+lYD7ak0Z/m8jk9tHAV8HZ/7saU0e+DKzvskxfvtTLRr+pzC5nPybJJXmxNSXe+UZMW6WVMqbzw+h/m+ymrrsrx+K1Fp7dYbs2oUfWgkbRC0CyLqNgg6Pl/Bf8uPVXl+OG7h+ZqZLSyEknq2Uf2rcHfWwbi1sExFBbREeS/a6iRrwhWLMi9XKef11Jay1ax8X1gPP1KS71yzlYCZLU1BFIMkPYzve7s8S/JKVme+J2339LmqSrmrgfOB8+XBLpfi+/oWmEczt5r0J2AAMNzMRmfy8p/85ZSsVCnvDpmoaPf3gScrVPdsBVmed6cKsneBI3HLXpm2VM7biqIPFecwCNojoegFQQdG0k7ANbhL6q9VkhUb6t+3qEg6At+n9kQbNGMmbpnZ08ym1UjXqNUEM1sr6XFgXzMb2WhDzOw9SfOAk4DR2aOTS0lfw61H+xUCSTvgSvGyBqpqxoMtZuLjm7tpnwLW4YEvedsH4gEJm1nOqvRjgaRLgNOBz7DJVdlatsPnJp/7HnggRa5ctwDDtPnRNieWyloEvAT0MbMxtI5FwKdKsmLd9DSz6a0sb2vRB1hpHwwqCYJ2Syh6QdBx6CqpiKztgVvQzsb3gx1rZpWsaeBn6r0NjJF0LW7da8Jf2luMmS2SNBq4M5U/h017tvY1s+KYlYXAMZKOwS1dz9Z5YV4KzJC0EY+UfAvfIzYAuMLMFlfJ91PgXvkZg5PwCNljS23eKGkycGHau7Yaj+hs9Py+ibgrdBTwkJkVVjzMbJWkG4ArJb2LB6GcjAcCVI26TdbBSbiiaMD38KCOOQ226QOY2RpJLcAISW/ilrPLcHfpjlnSG/Bz7u6XdD3uyr0MV1g3prI2SroYuF3SjsBUXHnfGz8I+1QzW1elKY8A35TUxcyK8op105zWzVx83XwW6GtmZ33Yfm8pkkbggUrld+QheORtEHQYQtELgo5DT9w9a/hZc0uAPwA3mtnyapnM7FVJpwHX4UEKzwA/wBWptuJcYDGunDSl9j2NH2NR8GPSZn5cyRiG78eq1u6HJX0JP8Ljdtz6swyYRg13sJlNknQerqgMwaM2v4sHh+QMx4/ouBl4A/gJfpxMraCCoo4XJD2Kux2vrpBkBO6WPBt32S4BvmNmzTWKfQyP6uyDu8CfBI4zPyh5SxiM93MCrmD/Bv9zMLxIYGYvSRoA/AoPXPk3Hk07nexcQzO7KymMl6fn7wFL8aNsiv1vlZiMH3lzJB5NXZCvm5FUXjdbgy74enufFChyNP6HIAg6DDJr5JimIAiC4P8JSUfhSll/M5vVBuVNBl40s3O3uHFbgWSJngh8wszWbu32BEGjhKIXBEEQIOka3Iq4HD965yrcAnhg4W7dwvIPxY862at0VEqHQNI04HEza9rabQmC1hCu2yAIggA8cGMU7mp+Cz+S56K2UPIAzKxF0qW4+75DKXqSuuOu9eu3dluCoLWERS8IgiAIgqCTEgcmB0EQBEEQdFJC0QuCIAiCIOikhKIXBEEQBEHQSQlFLwiCIAiCoJMSil4QBEEQBEEnJRS9IAiCIAiCTsr/AOPQeg7U6Kt8AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 720x504 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "sb_bocage['exp_line'] = np.exp(lsb_bocage['line'])\n",
    "font = {'family' : 'normal',\n",
    "        'weight' : 'normal',\n",
    "        'size'   : 15}\n",
    "import matplotlib\n",
    "matplotlib.rc('font', **font)\n",
    "ax = sb_bocage.plot.scatter('radius','sb', loglog=True, figsize=(10,7),s=50)\n",
    "ax.set_xlim((20, 1700))\n",
    "ax.set_ylim((sb_bocage.sb.min()-0.00001, sb_bocage.sb.max()+0.01))\n",
    "sb_bocage.plot.line('radius', 'exp_line',color='Red',ax=ax, legend=False)\n",
    "plt.grid(which='both')\n",
    "ax.set_xlabel('Diamètre du voisinage (en m.)')\n",
    "ax.set_ylabel('Mr(i)')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Calcul openfield"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {},
   "outputs": [],
   "source": [
    "openfield = gpd.read_file('building_openfield.gpkg')\n",
    "r_openfield = rasterizeBlackWhite(openfield, 4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "25\n",
      "50\n",
      "100\n",
      "200\n",
      "400\n",
      "800\n",
      "1600\n"
     ]
    }
   ],
   "source": [
    "neigh_openfield = neighbors_near_point(r_openfield, list_of_radius)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {},
   "outputs": [],
   "source": [
    "sb_openfield = Gwfa_en_un_point(neigh_openfield,786645,6879505, list_of_radius)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 173,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([  1.22020562, -12.40344766])"
      ]
     },
     "execution_count": 173,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "lsb_openfield = np.log(sb_openfield)\n",
    "d = np.polyfit(lsb_openfield['radius'],lsb_openfield['sb'],1)\n",
    "f = np.poly1d(d)\n",
    "lsb_openfield['line'] = f(lsb_openfield['radius'])\n",
    "d"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 174,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'Mr(i)')"
      ]
     },
     "execution_count": 174,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAnoAAAG4CAYAAADSYkouAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOzdeZxd4/3A8c8jkZCMXUqsURpLCa2tQUusRWJpEbRVpVRr158tqKDWlqqltRVVIiViC2pNqFqy0KRksUaQWGKfjEySyfP749yJyZhJZrn3nnPPfN6v133JnHue53zPmevkm+d7n+eEGCOSJEnKnyXSDkCSJEmlYaInSZKUUyZ6kiRJOWWiJ0mSlFMmepIkSTlloidJkpRTndMOIIuWX375uN5666UdRlHMmjWL7t27px3GAuWOp5THK3bf7e2vPe3b0jZrn628ycv1zeJ5eB8qXV9t7aPc7fJm3LhxM2OMPZp8M8boq9Grd+/eMS9GjhyZdggLKXc8pTxesftub3/tad+Wtln7bOVNXq5vFs/D+1Dp+mprH+VuV0m+mD033v78W/HCByfG259/K34xe+7X9gHGxmZyGkf0JEmSMmjM1I859KbRxAg1c+ro1qUT5z0wkZt/sRVb9lqxRX34HT1JkqSMqa6dx6E3jWZWbR01c+qAJNmbVVtX2D6vRf2Y6EmSJGXMiPHTiRE6za/jsDH30uvjdxe8FyOMmDC9Rf1YupUkScqYqR/NYuPXx3PuI39lg5lvsWxtNZdv9xMgGdmbOrOmRf2Y6EmSJGXJjBkceMUgTvvXPbyzbA+O3PcMHvnW9xa83a1LJ3qt3K1FXZnoSZIkZcHcuXDVVXD22axdW8s13z+Iy7f8MbOXXGqh3UKA/n1Wa1GXfkdPkiQpbU89Bd/9Lpx0Emy3HeGll9j8lqvoVNWdbl06AclIXveunbj5F1vRvWvLxuoc0ZMkSUrLjBlw8slw222w9tpwzz2w114QAlsCowftzIgJ05k6s4ZeK3ejf5/VWpzkgYmeJElS+TUo01JbC2edBaedBt0W/u5d966dGbjlWm0+jImeJElSOT35JBx9NLz8Muy+O1xxBZTo0asheXKGAEIIA4ABPXv2PGLIkCFph1MU1dXVVFVVpR3GAuWOp5THK3bf7e2vPe3b0jZrn628ycv1zeJ5eB8qXV9t7aPc7dLSZeZM1r3mGlZ5/HFmr7IKrx5zDB9tu20yu6Id+vXrNy7GuEWTbzb3bLSO/PJZt6XjMyZL15/Pus2XvFzfLJ6H96HS9eWzbpsxZ06Ml14aY1VVjF27xnjWWTHOmlW07vFZt5IkSSkYNQqOOSYp0+6xB/z5zyUr0zbF5VUkSZKKbfp0OPhg6NcPZs2Ce++FESPKmuSBiZ4kSVLxzJ0Ll14K668Pw4fD734HEycuWDKl3CzdSpIkFcOoUcls2okTkzLtFVfAuuumGpIjepIkSe3RsExbU/NVmTblJA9M9CRJktomY2Xapli6lSRJaq2RI5PZtBkq0zbFET1JkqSWevddOOgg2HHHzJVpm2KiJ0mStDhz58If/wgbbAB3353JMm1TLN1KkiQtSsMy7Z57JoseZ3QErzFH9CRJkprSuEx7332ZLtM2xURPkiSpocZl2rPPTkbzBgxIO7JWs3QrSZJU74knkjLtpEnQvz9cfnlFjeA15oieJEnSu+/CgQfCTjvB7NlJmfb++ys6yQMTPUmS1JHNmQN/+EOy6PE998DgwfDyyxVZpm2KpVtJktQxPf54UqadPDlJ7C6/HL75zbSjKipH9CRJUsfyzjswcCDsvDPU1iYl2vvuy12SByZ6kiSpo5gzBy65JJlNe999X5Vp+/dPO7KSsXQrSZLyrwOUaZviiJ4kScqvhmXaOXNyXaZtiomeJEnKn8Zl2nPOyX2ZtimWbiVJUr489hgce2xSpt1rr6RMu846aUeVCkf0JElSPrzzDhxwAOyySzKiN2IE3Htvh03yAEKMMe0YMiOEMAAY0LNnzyOGDBmSdjhFUV1dTVVVVdphLFDueEp5vGL33d7+2tO+LW2z9tnKm7xc3yyeh/eh0vXV1j7a2y7Mncsaw4bR65ZbYP58ph18MG8fdBDzu3RpdZ+VqF+/fuNijFs0+WaM0VejV+/evWNejBw5Mu0QFlLueEp5vGL33d7+2tO+LW2z9tnKm7xc3yyeh/eh0vXV1j7a1e7RR2Ncf/0YIca99orxjTfa1FclA8bGZnIav6MnSZIqz9tvs9HgwfDkk8kM2hEjYM89044qc/yOniRJqhxz5sDFF8MGG7DSs8/Cuecms2lN8prkiJ4kSaoMjz6azKadMgX23psxAwfyvYMOSjuqTHNET5IkZdvbb8P++8Ouu8K8efDAA3DPPczu2TPtyDLPRE+SJGXTnDlw0UXJoscjRsB558FLL8Eee6QdWcWwdCtJkrLnkUeSMu0rr8A++8Cf/gS9eqUdVcVxRE+SJGXHtGmw336w225QVwcPPgh3322S10YmepIkKX21tXDhhbDhhklyV1+m3X33tCOraJZuJUlSqlYYMwZ+9SvLtCVgoidJktIxbRqcdBKb3nUXrLdeMpLnCF5RWbqVJEnl1ahM+8bhh8P//meSVwImepIkqXwefhg22QQGDUomXEyaxLSf/hSWWirtyHLJRE+SJJXetGnw4x/DD38IMcJDD8Hw4bD22mlHlmsmepIkqXRqa+GCC5JFjx96CH7/+2Q27Q9/mHZkHYKTMSRJUmk8/HCy6PGrr8K++yazaR3BKytH9CRJUnFZps0MEz1JklQUYc6chcu0559vmTZllm4lSVK7VNfOY+y1t/Pdi0+H6e8yb5996Xy5ZdoscERPkiS12finx/PsZtuzw/GHUDc/csTBv2fTPkcyJi6TdmjCRE+SJLVFbS2155xH735bs91r47jkB4cw5MI/8+iamzGrto5DbxrNrNp5aUfZ4ZnoSZKk1vnXv2Djjek6+Hc8td4W7HTEX/lL3wOYv+SSC3aJEUZMmJ5ikAK/oydJklpq6lQ48US45x7o3Zuhv7+e077o2eSuNXPqmDqzprzx6Wsc0ZMkSYs2e3ay0PFGG8EjjyQzaydMgF12pVuXTk026dalE71W7lbmQNWYiZ4kSWreQw8lz6Y96yzYYw+YNAlOPx26dqX/pqsRQtPNQoD+fVYrb6z6GhM9SZL0dVOnJk+z2GMPWGKJ5CkXw4bBWmst2KWqa2du/sVWdO/aacHIXrcunejetVNhu98QS5u/AUmS9JXZs+GPf0wWO15iCbjwwuR7eV27Nrn7lr1WZPSgnRkxYTpxxkTOHrAR/fusZpKXEf4WJElS4qGHkmfTvv467LcfXHrpQiN4zenetTMDt1yLUaPeYIctF7+/ysdET5Kkjm7qVDY+80z4z3+gd+9kwsUuu6QdlYrA7+hJktRRzZ4N550HG27ICuPGJWXaCRNM8nLEET1JkjqiBx+E445bUKYdvf/+9D3ggLSjUpE5oidJUkcydSrssw/suSd07pyUae+8k9pvfCPtyFQCJnqSJHUEDcq0PPooXHSRZdoOwNKtJEl598ADcPzxSZl2//2T2bRrrpl2VCoDR/QkScqrN9+EvfeG/v2TMu2jj8Idd5jkdSC5T/RCCGuGEB4PIUwKIbwcQrgkhOYe2CJJUg7Mng3nnps8m/bxx+Hii5My7c47px2ZyqwjlG7nAafGGMeGELoAjwI/Au5KNyxJkkrggQeS2bRvvAEHHJCUaddYI+2olJJURvRCCOuFEK4NIYwPIdSFEEY1s99GhdG4mhDC9BDCuSGETq05VoxxRoxxbOHPc4AJgGPWkqR8efNN2GuvpEzbpUtSpv3nP03yOri0RvS+DewBPAd0aWqHEMIKwGPARGBvYF3gUpLk9My2HDSEsBKwD7BrW9pLkpQ5X34Jl1ySzKLt1Ckp055wQpLsqcNLK9G7P8Z4L0AIYRiwchP7HAUsDfwoxvg58GgIYVlgcAjhksI2QghPA039c+XxGOPh9T+EELoCw4DLY4yTins6kiSV30rPPguHH26ZVs1KJdGLMc5vwW67Aw/XJ3QFQ4GLge2B+wt9bbe4jgrl3tuAF2OMl7Y+YkmSMuSNN+CEE9jk/vthgw3gscdgp53SjkoZlOVZtxsAkxtuiDFOA2oK77XGtcAXwG+LE5okSSn48ks455xkNu0TT/D6r34F48eb5KlZIcaYbgCF0m2McYdG2+cCJ8cYL2+0/R3glhjjoBb2vy3wNPASUFfYfGOM8YpG+x0JHAnQo0ePze+44442nE32VFdXU1VVlXYYC5Q7nlIer9h9t7e/9rRvS9usfbbyJi/XN4vnUan3oZWefZb1rrySpWfM4IN+/Xj917/mo6WXLtq5FCPOtvZR7nZ5069fv3Exxi2afDPGmOqL5Htzo5rYPhc4vont7wLnlzKm3r17x7wYOXJk2iEspNzxlPJ4xe67vf21p31b2mbts5U3ebm+WTyPirsPvf56jP37xwgxbrhhjI89Vry+GyhGX23to9zt8gYYG5vJabJcuv0EWL6J7csBn5Y5FkmSyuvLL2Hw4KRMO3JkMrP2v/+1TKtWyfKCyZNp9F28EMKaQHcafXdPkqRcuf/+5Nm0b74JAwfCH//obFq1SZZH9B4CdgshLNNg20DgS+DJdEKSJKmE3ngDBgxIFj5eaqnk8WVDh5rkqc1SGdELIXQjWTAZYHVg2RDCfoWfH4wx1gDXAMcBw0MIFwPfBAYDl8WFl1yRJKmyffllstDxRRdB587whz8kjzFz0WO1UyqzbkMIvYA3m3l7nRjj1MJ+GwFXAX1Jvpd3AzA4xljXTNv2xjUAGNCzZ88jhgwZUopDlF3WZiRV6my3cvTtrFs1lJfrm8XzyNp9aKVnnmG9q65i6RkzeH/HHXn9qKOY06NHUfouZpyl7MNZt+2T6Vm3WXw567Z0Km62Wxn7dtatGsrL9c3ieWTmPvTaazHuuedXs2mfeKJ4fbeBs24rFxU661aSpPz58ks4+2z49rfhySeTiRbjx0O/fmlHphzK8qxbSZLyI8avZtNOnQoHHZR8F2/11dOOTDnmiJ4kSaX2+uvQvz/svTd06wZPPAFDhpjkqeQc0ZMkqVRqauh1441wxx2w5JJJmfa445I/S2VgoidJUrHFCPfdByecQK/6Mu0f/wirrZZ2ZOpgUlleJatcXqX0srasQZb6dnkVNZSX65vF8yh1TEu/+y7rXXklKz3/PLPWXpsJv/oVtX37luRYLq+Svc9XGlxexeVVMiMzyxpksG+XV1FDebm+WTyPksU0a1aMZ50VY5cuMVZVxXjppTHOmVMx9yGXV6lcLGJ5FUu3kiS1R32Z9vjj4a234OCDk9m0lmmVAc66lSSprV57DfbcE/bZB6qqYORIuO02kzxlhomeJEmtVVMDZ52VLHr89NNw6aXw4ouwww5pRyYtxNKtJEktFSPcey+ccIJlWlUER/QkSWqJ+jLtvvsmZdpRoyzTKvNM9CRJWpTGZdrLLkvKtNtvn3Zk0mJZupUkqSmNy7Q/+UlSpu3ZM+3IpBZzweQGXDC59FwwuXT9uWByvuTl+mbxPFoS09LvvJMsejx6NNXrrMOrxx3HZ5ttVrLjtZULJmfv85UGF0x2weTMcMHk0vXngsn5kpfrm8XzWGRMs2bFeMYZyaLHyywT45/+FOOcOaU7Xju5YHLb2uUNLpgsSdIixAj33JOUaadNg5/+FC65xDKtKp6JniSpY3v1VTj2WHj4Ydh4Y3jySfjBD9KOSioKZ91KkjqmWbPgjDOS5O6ZZ+BPf4IXXjDJU644oidJ6lhihOHD4cQTLdMq90z0JEkdxyuv0OfUU2HMGNhkE8u0yj1Lt5Kk/Js1izmnnk7dxptQ9fLLvHDi2VQ/O9okT7lnoidJyq9Cmba29/p0ueQiHthgO/5x8VX8tOp7bH3JKMZM/TjtCKWScsHkBlwwufRcMLl0/blgcr7k5fqmeR5Lv/0237riClYcO5aZa67Nkz8/kukbfJtVlob3v0z2WSIENuy5LEuE0sVRKfchF0yuXC6Y7ILJmeGCyaXrzwWT8yUv1zeV86iujvH002NccskYl102jjvx7LjxGSPi2qcmrytuvWfBnzc866E4dPRbJQ2nUu5DLphcuVjEgsmWbiVJ+RAj3HUXbLghXHghHHggTJnCI7sM5It5TTepmVPH1Jk15Y1TKiMTPUlS5XvlFfjhD2G//WD55eGpp+CWW2DVVem1Une6denUZLNuXTrRa+VuZQ5WKh8TPUlS5Zo1CwYNShY9fu45+POfk0WPv//9Bbv033Q1QjPfwQsB+vdZrUzBSuVnoidJqjyNy7QHHQRTpsBxx0HnhZeIreramZt/sRXdu3ZaMLLXrUsnunftVNjukrLKLz/dkqTKMmVK8mzaRx+FPn1gyBDYbrtFNtmy14qMHrQzIyZMJ86YyNkDNqJ/n9VM8pR7juhJkirDrFlw+unJEy2efx6uuALGjVtskleve9fODNxyLVZddikGbrmWSZ46BD/lkqRsqy/TnngivPMO/PzncPHFsMoqaUcmZZ4jepKk7Jo8GXbdFfbfH1ZaCZ5+Gm6+2SRPaiETPUlS9lRXw2mnJd/BGzMmKdOOHQvbbpt2ZFJF8RFoDfgItNLzEWil689HoOVLXq5vq88jRno8+STr/uUvLPXhh7y32268fuSRzF1xxfRiyvDxfARaPv4/aS8fgeYj0DLDR6CVrj8fgZYvebm+rTqPSZNi3HnnGCHGTTeN8emn048p48fzEWhta5c3+Ag0SVJmNS7TXnmlZVqpSJx1K0lKR4wwbBicdFIym/bQQ+Gii5xoIRWRI3qSpPKbPBl22QUOOABWXhn+8x+46SaTPKnITPQkSeVTXQ2nnpqUaceOTcq0Y8bANtukHZmUS5ZuJUmlFyPceWdSpn333aRMe/HF8I1vpB2ZlGuO6EmSSmvSpKRMO3Ag9OjxVZnWJE8qOUf0JEmlUV3NN6+9Nplw0b07XHUVHHUUdOqUdmRSh+GIniSpuGKEO+6ADTZgraFD4Wc/g1degaOPNsmTysxET5JUPJMmwc47LyjTvnDVVXDjjZZppZSY6EmS2u+LL+CUU5LZtC+8AFdfDWPH8vm3v512ZFKH5nf0JEltV1+mPekkmD4dDjsMLrzQETwpIxzRkyS1zcSJSZn2wAOThY6feQb+9jeTPClDQvIsXAGEEAYAA3r27HnEkCFD0g6nKKqrq6mqqko7jAXKHU8pj1fsvtvbX3vat6Vt1j5beZPl69uppoa1b7mFNYYNo27ppXnz8MOZPmBAkxMtsnge3odK11db+yh3u7zp16/fuBjjFk2+GWP01ejVu3fvmBcjR45MO4SFlDueUh6v2H23t7/2tG9L26x9tvImk9d3/vwYb789xtVWixFiPOywGD/4YJFNsnge3odK11db+yh3u7wBxsZmchpLt5KkxZs4EXbaCQ46aOEybY8eaUcmaRFM9CRJzfviCzj5ZNh0U3jxxWQ27Zgx0Ldv2pFJagFn3UqSvi5G+Oc/4be//Wo27UUXOYInVRhH9CRJC3v5Zcu0Uk6Y6EmSEl98Af/3f7DZZkmZ9i9/sUwrVThLt5LU0cUIQ4cmSd706XD44cmix47gSRXPET1J6shefhl23BEOPhhWXRWefRZuuMEkT8oJEz1J6ogalmnHj0/KtKNHw/e+l3ZkkorI0q0kdST1Zdrf/hZmzLBMK+WcI3qS1FG8/DL065eUaXv2tEwrdQAmepKUd59/nozgbbopTJgAf/2rZVqpg7B0K0l5FSPcfnvyXbwZM+CXv0zKtCuvnHZkksrERE+S8uill+CYY+DJJ2HzzeHuu2HrrdOOSlKZWbqVpDz5/HM46aRkNu2ECXDNNfD88yZ5UgfliJ4k5UGMMGRIUqZ9//2kTHvBBZZppQ7ORE+SKt1LL8HRR8NTT8EWW8C998JWW6UdlaQMCDHGtGPIjBDCAGBAz549jxgyZEja4RRFdXU1VVVVaYexQLnjKeXxit13e/trT/u2tM3aZytvWnJ9O82aRa+bb2aN4cOZV1XFG4cfzow994ROncoU5eJl8XPifah0fbW1j3K3y5t+/fqNizFu0eSbMUZfjV69e/eOeTFy5Mi0Q1hIueMp5fGK3Xd7+2tP+7a0zdpnK28WeX3nz4/x1ltjXHXVGEOI8YgjYvzww7LF1hpZ/Jx4HypdX23to9zt8gYYG5vJaZyMIUkZUl07j6Gjp/He57MZOnoa1bXzFt7hf/+DHXaAn/4U1lgDnnsOrrvO7+JJapKJniRlxJipH7P1BY9x7oiJfPhFLeeOmMjWFzzGmKkfw2efwYknwne+k3wn79prkyTP7+JJWgQnY0hSBlTXzuPQm0Yzq7ZuwbaaOXUQI8NOvIjNn/sHS7z/PhxxRDKbdqWVUoxWUqUw0ZOkDBgxfjqN58at/+FUzn3kr2z9zst8vGEfVrzvPthyy3QClFSRTPQkKQOmfjQrGcEDlqmdxXa33sqvH36AL7p25/TdjmGFY47ilC2/nXKUkiqNiZ4kZUCvlbrTbckl2PW/jzNo5I30qPmU2/vsxiXbH8Kc5Vbg7FWWSTtESRXIRE+SMmCvJWay3i2nsMW0l/hvz28x8pRBDKrZEIDuAfr3WS3lCCVVImfdSlKaPvsMTjiBbt/bis0+f5ez+x/HTw6/nA/W/RbdunSie9dO3PyLreje1X+XS2o97xySlIYY4dZb4eST4YMP4Mgj6Xz++ZxStRwbTZhOnDGRswdsRP8+q5nkSWozR/QkqdwmTIAf/AAOOQTWWgtGj4ZrroGVVqJ7184M3HItVl12KQZuuZZJnqR2MdGTpHL57DM4/nj47ndh0iS4/vpk0eMtmn5EpSS1l/9UlKRSa1ym/dWv4Pe/d9FjSSVnoidJpTRhAhx9NDz9dPK4shEjHMGTVDaWbiWpFJoq0z77rEmepLJyRE+SiilG+Mc/4JRTvirTnn8+rLhi2pFJ6oBM9CSpWMaPh2OOsUwrKTMs3UpSe3366cJl2htusEwrKRMc0ZOktqov0558Mnz4IRx1VDKb1jKtpIww0ZOkthg/PplN+5//wNZbw4MPwuabpx2VJC3E0q0ktcann8JxxyVl2ilTkjLtM8+Y5EnKJEf0JKkl5s//ajatZVpJFaJViV4IYSNgJ2ArYFVgKeBj4BXgaeCRGOOXxQ5SklJlmVZShVps6TYkDgkhjAFeAgYDawCfAm8BnYDdgbuA90II14UQ1ildyJJUJo3LtH/7m2VaSRUlxBgXvUMIkwt//AdwV4xxcjP7dQN2A/YH9gKOijHeWsRYSy6EMAAY0LNnzyOGDBmSdjhFUV1dTVVVVdphLFDueEp5vGL33d7+2tO+LW2z9tkqqvnzWeWRR1j3uutY8rPPmL7XXrx52GHMW2aZsoWQl+ubxfPwPlS6vtraR7nb5U2/fv3GxRibXs8pxrjIF7AfhYSwpS+SEb/vtaZNll69e/eOeTFy5Mi0Q1hIueMp5fGK3Xd7+2tP+7a0zdpnq2hefDHGbbaJEWLs2zfGF15IJYy8XN8snof3odL11dY+yt0ub4CxsZmcZrGl2xjjsEInLRZjfCfG+Fxr2khSqj79FI49NinLvvoq3Hhj8oSL73wn7cgkqc2cdSupY5s/H265JZlN+9FH8Otfw3nnwQorpB2ZJLXbYhO9EMJo4NAY48TChIxFju7FGLcqVnCSVFIvvpjMpn32WejbFx5+2BE8SbnSkhG9l4EvG/y5VWVcScqcTz6Bs86Cv/4VVloJbroJDjkElnANeUn5sthEL8b4iwZ/PrSk0UhSKc2fD3//O5x6alKm/c1v4NxzLdNKyi2/oyepY2hYpt1mG3jkEdhss7SjkqSSasmCyX8LIXyrpR2GEJYMIRwWQvhZ+0KTpCL45BM45hjYYgt47bWkTPvvf5vkSeoQWjKiVwOMDyGMA4YBzwAvxQaPOgshrA1sTvKEjH2Ad4Ejix+uJLVQU2Xa886D5ZdPOzJJKpuWrKN3LLA+8BTwW+B5oDqEMCuE8HEIYR7wBjCU5Pm3R8QYN4sxji5h3JLUvBdegG23hcMOg299C8aNgyuvNMmT1OG06Dt6Mca3gTOAM0IIvUlG71YFlgI+BqYAo2OMNaUKVJIW65NP4Mwz4Zprktm0N98MP/uZs2kldVgtnowRQlgS2Ap4M8Z4e+lCkqRWmj8/SepOPRU+/jiZdHHuuY7gSerwWvPP3DrgCWDDEsUiSa1XX6Y9/HBYf/2kTHvFFSZ5kkQrEr0Y43zgVWCV0oUjSS30ySfJyN2WW8IbbyQTL5xNK0kLae0XV84AfhdC2KQUwUjSYs2fDzfeCL17J9/FO+YYmDIlebJFCGlHJ0mZ0toFk88EVgL+G0J4F3ifRo9E81m3kkrmhReSUbznnkvKtVdfDZtumnZUkpRZrU30XgZeKkUgktSs+tm0f/0r9OiRlGl/9jNH8CRpMVqV6PmsW0ll1Xg27bHHwjnnONFCklposYleCOF3regvxhjPa0c8kpR44YXkaRbPPw/bbQdXXWWZVpJaqSUjeoOBL4FZwOLqJBEw0ZPUdh9//NWix9/4BtxyC/z0p5ZpJakNWpLovQGsBYwjeczZ3THGz0salaSOZ/58uOkmOO20JNk77rikTLvccmlHJkkVqyXPul0P2IZkIsZ5wHshhOEhhP1DCEuXOkBJHcC4cbDNNvDLX8IGG8CLL8Lll5vkSVI7tWgdvRjj2Bjj/8UY1wJ+CLwHXAV8EEK4LYTwg1IGKSmnPv4Yfv3rZNHjqVOTMu1TT0GfPmlHJkm50OonfccYn4ox/gZYE7gGGAicUOzAJOXY/Plwww3JosfXX5+UaadMcckUSSqyVid6IYRtQwhXAm8BvwaGAX8udmCSsqm6dh5DR0/jvc9nM3T0NKpr57Wug/oy7RFHwIYbJrNrLdNKUkm0KNELIXw3hHBJCOEt4HGS0bwTgW/EGA+MMT5ZyiAlZcOYqR+z9QWPce6IiXz4RS3njpjI1hc8xpipHy++ceMy7T/+YZlWkkpssYleCGEK8BzQBzibJLnbJ8Y4NMZYU+oAJWWQcj0AACAASURBVGVDde08Dr1pNLNq66iZUwdAzZw6ZtXWFbY3M7LXuEx7/PFJmdYlUySp5FoyovctYB6wOXAJ8FoI4YPmXiWNVlJqRoyfToxNvxcjjJgw/etvjB0LffsmZdqNNkpm0/7pT5ZpJalMWrKO3jklj0JS5k39aNaCkbzGaubUMXVmgwH+jz6CM86A666DVVaBW2+Fgw92BE+SymyxiV6M0URPEr1W6k63Lp2aTPa6delEr5W7JWXav/0NTj8dPv0UTjgBBg+GZZctf8CSpNbPupXUMfXfdLVmB+RCgL3mTk/KtEceCd/+dlKmvewykzxJSpGJnqQWqeramZt/sRXdu3aiW5dOQDKSt1rdLEa+9k+W3m4bmDYtKdOOGgWbbJJuwJKkFn1HT5IA2LLXiowetDMjJkwnvvsSf68bzxY3XkL47DPLtJKUQSZ6klqle9fODOR9Pj/9tyw7eTL84Adw9dWw8cZphyZJasTSraSW++gj+NWvYOut6frBB3DbbUmZ1iRPkjIp9yN6IYQngeWBALwCHBZj/DzdqKQKU1f31Wzazz6DE09k9I478v0990w7MknSInSEEb29Yoybxhj7ANOAk9MOSKooY8bA976XjORtsgn8979w6aXUde+edmSSpMUoe6IXQlgvhHBtCGF8CKEuhDCqmf02CiE8HkKoCSFMDyGcG0Lo1NrjxRg/K/S3BNAdaGZtf0kLmTkzWSpl663h3XeTMu3IkZZpJamCpFG6/TawB8nzc7s0tUMIYQXgMWAisDewLnApSWJ6ZmsPGEJ4ENgSeBn4bZuiljqKurrk2bSDBi0o03L22c6mlaQKlEbp9v4Y45oxxv1JEq+mHAUsDfwoxvhojPEakkexnRRCWPC3TQjh6RDC1CZef2vYWYxxD2BVYDTwm5KclZQHo0cnZdqjjlqoTGuSJ0mVqeyJXoxxfgt22x14uNGkiaEkyd/2DfraLsbYq4nX4U0ctw74O3BIO09Byp0lP/ssKdN+73tJmXbIEMu0kpQDWZ2MsQEwueGGGOM0oKbwXouEEFYIIazSYNOPgZeKEqGUB3V1cO21bHXIIXDjjXDSSTB5Mhx0EM0+70ySVDFCjOnNTQghDANWjjHu0Gj7XODkGOPljba/A9wSYxzUwv6/CdxB8l3AAEwCjo0xvt/EvkcCRwL06NFj8zvuuKP1J5RB1dXVVFVVpR3GAuWOp5THK3bf7e2vte2XmTSJ3pdfzjKvvMLMTTbhjRNPpGaddUp2PLVOXq5vFs/D+1Dp+mprH+Vulzf9+vUbF2Pcosk3Y4ypvYBhwKgmts8Fjm9i+7vA+aWOq3fv3jEvRo4cmXYICyl3PKU8XrH7bm9/LW7/4Ycx/vKXMYYQY8+eMd5+exz5xBOlO57aJC/XN4vn4X2odH21tY9yt8sbYGxsJqfJaun2E5JFjhtbDvi0zLFI+VBXB9dcA717w803w29/C1OmwIEHWqaVpJzK6pMxJtPou3ghhDVJ1sGb3GQLSc17/nk4+mgYNw769YOrroKNNko7KklSiWV1RO8hYLcQwjINtg0EvgSeTCckqQLNnAlHHJHMpp0xA26/HR5/3CRPkjqIso/ohRC6kSyYDLA6sGwIYb/Czw/GGGuAa4DjgOEhhIuBbwKDgcuiz6mVFq+uDq6/Pln0+Isv4P/+D373O1hmmcW3lSTlRtln3YYQegFvNvP2OjHGqYX9NgKuAvqSfC/vBmBwTNbDK1VsA4ABPXv2PGLIkCGlOkxZZW1GkrPdStdffftlJk5MZtO++iqffOc7vHrccdT06lX0Y2fts5U3ebm+WTwP70Ol68tZt+nI7KzbrL6cdVs6znYrXX9P3313jIcfHiPEuNpqMQ4dGuP8+SU7dtY+W3mTl+ubxfPwPlS6vpx1mw4WMes2q5MxJLVUXR1cdx1bnXoqfPklnHwynHWWZVpJkomeVNGeey6ZTfvCC1R/5zuscOutTrSQJC2Q1Vm3khblww/h8MOhb1947z0YOpTxl15qkidJWoiJnlRJ6urgL39JFj2+5ZakTDt5Mgwc6KLHkqSvsXQrVYoGZVp23DFZ9HjDDdOOSpKUYWVfXiXLXF6l9FzWoPX9Lfnpp3zzuuvo+dBD1K68Mq/95jd8uMMOXxvBa088Lq+SPXm5vlk8D+9DpevL5VXS4fIqLq+SGS5r0Ir+5s2L8eqrY1x++Rg7d47xlFNi/OKLksTj8irZk5frm8Xz8D5Uur5cXiUduLyKVGGefTYp0774Iuy0E1x5pWVaSVKrORlDypIPPoDDDoNttkn+fMcd8OijJnmSpDZxRE/Kgro6Vrv7bth3X6iuhlNPhTPPBL97IklqBxM9KW2FMm3vF1+EnXdOyrQbbJB2VJKkHLB0K6WlYZn2ww95+eyz4ZFHTPIkSUVjoieVW10dXH01rL8+3HprUqadNKnJJVMkSWoPS7dSOT3zTDKb9r//tUwrSSo5F0xuwAWTS6+jLlS65CefJIse/+tfzO7Rg9ePPpoPf/CDhUbw2hurCybnS16ubxbPo6Peh8rRlwsmp8MFk10wOTM63EKlc+fGeOWVMS63XIxLLhnjaac1u+hxe2N1weR8ycv1zeJ5dLj7UBn7csHkdOCCyVIK/vOfpEw7fjzssktSpl1//bSjkiR1IE7GkIrt/ffh0ENhu+3go4/gzjvh4YdN8iRJZWeiJxXLvHlfjdoNGQKnnQaTJ8N++zmbVpKUCku3UjFYppUkZZAjelJ7NCzTfvwxDBtmmVaSlBmO6EltMW8eqw8fDrfcAjU1cPrpcMYZ0L172pFJkrSAiZ7UWk8/DUcfzbcmTIBdd4UrrnAET5KUSS6Y3IALJpdeJS9UuuTHH7Putdey6iOPMPsb3+Clww+nepddijbRwgWT1VBerm8Wz6OS70Ol7NsFkyuXCya7YHJmVORCpXPnxvjnP8e47LLJoseDBsVYXV30c3HBZDWUl+ubxfOoyPtQGfp2weTKhQsmS21UKNNSX6a98kro3TvtqCRJahFn3UpNef99+PnP4fvfh08/hbvugn/9yyRPklRRTPSkhubNSyZX9O4NQ4fCoEEwcSL86EcueixJqjiWbqV6//43HHNMUqbdbbevEj5JkiqUI3rSe+/BIYfAD36QlGmHD4eHHjLJkyRVPBM95VJ17TyGjp7Ge5/PZujoaVTXzvv6TvPmwZ//nKyB989/JgseT5oE++5rmVaSlAsmesqdMVM/ZusLHuPcERP58Itazh0xka0veIwxUz/+aqd//xu++1044QTo2xdeegl+/3vo1i29wCVJKjITPeVKde08Dr1pNLNq66iZUwdAzZw6ZtXWcehNo6l56x342c+SMu1nn31Vpv3Wt1KOXJKk4nMyhnJlxPjpNPWwl07z6/jpc/fR+dIDYG5tUqYdNMgRPElSrvkItAZ8BFrplTqe9z6fzYdf1C74eZWlodOLL7P9369j5bffYvp3t+DtE4/nyzXWaPexin0uPgJNDeXl+mbxPHwEWun68hFo6fARaD4CLTNKHc/tz78VNzzrobj2qSPiFkffEidtu32MEN9etkc8Zv+z4tDnpxbtWD4Cre3H0+Ll5fpm8Tx8BFrp+vIRaOnAR6Cpo+i/6WpccP//OGzMvZz49K10mz+XK/oO5C9996dTVXcu2nT1tEOUJKlsTPSUK1XPP8Nzd55MtymTeGq9LZh61GFc89k6dApw8y+2ontXP/KSpI7Dv/WUDzNmwMknw2230W3ttZl9513MWGtzurw3ibN7bkT/PquZ5EmSOhyXV1FlmzsX/vSnZNHjO++EM8+EiRNZar8fMXCrtVl12aUYuOVaJnmSpA7Jv/1UuZ56Co4+OlnsePfdk6dcuB6eJEkLOKKnyjNjBvz0p7D99vDFF3DPPfDAAyZ5kiQ1YqKnytG4THvWWTBxIuy9t8+mlSSpCZZuVRmefDIp0778clKmveIKWG+9tKOSJCnTHNFTtk2fDj/5CeywA1RXf1WmNcmTJGmxTPSUTXPnwmWXJWXau+6yTCtJUhtYulX2jBoFxxyTlGn32COZTesIniRJreaInrJj+nQ4+GDo1w9mzYJ774URI0zyJElqo5A8C1cAIYQBwICePXseMWTIkLTDKYrq6mqqqqrSDmOBpuIJ8+ax+vDh9Lr5ZpaYN49pBx3EtIMPZn7XriU5XrEUu+/29tee9m1pm7XPVt7k5fpm8TzKHVOl3IeK0Vdb+yh3u7zp16/fuBjjFk2+GWP01ejVu3fvmBcjR45MO4SFfC2ekSNj3GijGCHGPfaI8bXXSnu8DPfd3v7a074tbbP22cqbvFzfLJ5HuWOqlPtQMfpqax/lbpc3wNjYTE5j6VbpaFimran5qky77rppRyZJUm6Y6Kmswrx5cOmlyWza4cPhd79LZtPutZezaSVJKjJn3ap8Ro5ki1/+Et56C/bcM5lN6wieJEkl44ieSu/dd+Ggg2DHHVmithbuu88yrSRJZeCInkpn7txk1O6cc5I/n302Y/r25Qe77ZZ2ZJIkdQiO6Kk0Ro6EzTaDk09OHl/28ssweHBRlkyRJEktY6Kn4mpQpuXLL5My7f33W6aVJCkFJnoqjrlz4Y9/hA02gLvvhsGDk1G8AQPSjkySpA7L7+ip/Z54Ink27aRJSWJ3+eXwzW+mHZUkSR2eI3pqu3ffhQMPhJ12gtmzkxLtffeZ5EmSlBEmemq9OXPgD39IFj2+996vyrT9+6cdmSRJasDSrVrn8ceTMu3kyZZpJUnKOEf01DLvvAMDB8LOOycjepZpJUnKPBM9LdqcOXDJJcls2vvuSxY/tkwrSVJFsHSr5jUs0+61V1KmXWedtKOSJEkt5Iievq5xmXbEiGTShUmeJEkVJcQY044hM0IIA4ABPXv2PGLIkCFph1MU1dXVVFVVtWjfMHcuawwbRq9bboH585l28MG8fdBBzO/SJZV4sn68Yvfd3v7a074tbcv9u+xo8nJ9s3ge3odK11db+yh3u7zp16/fuBjjFk2+GWP01ejVu3fvmBcjR45s2Y6PPhrjBhvECDHutVeMb7yRbjwVcLxi993e/trTvi1ty/277Gjycn2zeB7eh0rXV1v7KHe7vAHGxmZyGku3Hd0778ABB8Auu1imlSQpZ0z0Oqo5c+Dii5PZtPffD+eem8ym3XPPtCOTJElF4qzbjuixx5LZtFOmwN57w5/+5AieJEk55IheR/L221+VaefNgwcegHvuMcmTJCmnTPQ6gsZl2vPOg5degj32SDsySZJUQpZuc26FsWPhqKOSMu0++yRl2l690g5LkiSVgYleXr39Npx0EpsOGwbrrgsPPgi77552VJIkqYws3ebNnDlw0UVJmfaBB3jzsMOSMq1JniRJHY6JXp488ghssgmcfjrsuitMnMhbP/sZLLVU2pFJkqQUmOjlwbRpsN9+sNtuMH9+Uqa9+26/iydJUgdnolfJamvhwgthww2T5O73v4f//c8yrSRJApyMUbkeeQSOPRZeeQX23TeZTbv22mlHJUmSMsQRvUrTuEz70EMwfLhJniRJ+hoTvUrRVJn2pZfghz9MOzJJkpRRlm4rwcMPJ2XaV1+1TCtJklrMEb0smzYNfvzjZNQuRsu0kiSpVUz0sqi2Fi64IFn0+KGH4PzzLdNKkqRWs3SbNQ3LtD/6EVx2mSN4kiSpTRzRy4qGZVqAf/0L7rrLJE+SJLWZiV7aGpdpL7ggWfR4t93SjkySJFU4S7dpu/pqOOOMZDTvsstgrbXSjkiSJOWEiV7afv1r6NMHdt457UgkSVLOmOilqLp2HiP+9yFT565Or9HT6L/palR19VciSZKKI8QY044hM0IIA4ABPXv2PGLIkCElPVbNnDrenDkLgPkxskQIAKyzcne6delUtONUV1dTVVVVtP7aq9zxlPJ4xe67vf21p31b2mbts5U3ebm+WTwP70Ol66utfZS7Xd7069dvXIxxi6beM9Frwvrrrx+nTJlSsv6ra+ex9QWPMau27mvvde/aidGDdqZ7kUb2Ro0axQ477FCUvoqh3PGU8njF7ru9/bWnfVvaZu2zlTd5ub5ZPA/vQ6Xrq619lLtd3oQQmk30nHWbghHjp9Ncfh0jjJgwvbwBSZKkXDLRS8HUj2ZRM+fro3mQlHSnzqwpc0SSJCmPTPRS0Gul5r+H161LJ3qt3K3MEUmSpDwy0UtB/01XozD34mtCgP59VitvQJIkKZdM9FJQ1bUzN/9iK7p37bRgZK9bl05079qpsN0lViRJUvuZUaRky14rMnrQzoyYMJ2pM2votXI3+vdZzSRPkiQVjVlFirp37czALX3kmSRJKg1Lt5IkSTlloidJkpRTJnqSJEk5ZaInSZKUUyZ6kiRJOWWiJ0mSlFMmepIkSTlloidJkpRTJnqSJEk5ZaInSZKUUyZ6kiRJOWWiJ0mSlFMmepIkSTlloidJkpRTJnqSJEk5ZaInSZKUUyZ6kiRJOWWiJ0mSlFMmepIkSTlloidJkpRTJnqSJEk5ZaInSZKUUyZ6kiRJOWWiJ0mSlFMmepIkSTlloidJkpRTJnqSJEk51WESvRDCX0IIMe04JEmSyqVDJHohhO8D3dOOQ5IkqZzKnuiFENYLIVwbQhgfQqgLIYxqZr+NQgiPhxBqQgjTQwjnhhA6teF4XYGLgP9rZ+iSJEkVpXMKx/w2sAfwHNClqR1CCCsAjwETgb2BdYFLSRLTM1t5vN8Bf4sxfhhCaGvMkiRJFSeNRO/+GOO9ACGEYcDKTexzFLA08KMY4+fAoyGEZYHBIYRLCtsIITwNrNFE+8djjIeHEPoAW9P65FCSJKnilT3RizHOb8FuuwMP1yd0BUOBi4HtgfsLfW23mH62BTYC3qwfzQshTAW2jDF+2LrIJUmSKktWJ2NsAExuuCHGOA2oKbzXIjHGv8YYV4sx9oox9ips62WSJ0mSOoIQY3orjtSXbmOMOzTaPhc4OcZ4eaPt7wC3xBgHtfF4McbY5Bf1QghHAkcC9OjRY/M77rijLYfInOrqaqqqqtIOY4Fyx1PK4xW77/b21572bWmbtc9W3uTl+mbxPLwPla6vtvZR7nZ5069fv3Exxi2afDPGmNoLGAaMamL7XOD4Jra/C5xf6rh69+4d82LkyJFph7CQcsdTyuMVu+/29tee9m1pm7XPVt7k5fpm8Ty8D5Wur7b2Ue52eQOMjc3kNFkt3X4CLN/E9uWAT8sciyRJUkXKaqI3mUbfxQshrEmy6PHkJltIkiRpIVlN9B4CdgshLNNg20DgS+DJdEKSJEmqLGVfXiWE0I1kwWSA1YFlQwj7FX5+MMZYA1wDHAcMDyFcDHwTGAxcFhdeckWSJEnNKPus2xBCL+DNZt5eJ8Y4tbDfRsBVQF+S7+XdAAyOMdaVMLYBwICePXseMWTIkFIdpqyyNiPJ2W6l689Zt/mSl+ubxfPwPlS6vpx1m47MzrrN6stZt6XjbLfS9ees23zJy/XN4nl4HypdX866TQcVOOtWkiRJ7WSiJ0mSlFMmepIkSTlloidJkpRTJnqSJEk5ZaInSZKUU2VfRy/LXEev9Fy/qnT9uY5evuTl+mbxPLwPla4v19FLh+vouY5eZrh+Ven6cx29fMnL9c3ieXgfKl1frqOXDhaxjp4jek0IIXwGvJp2HEWyHPBZ2kE0UO54Snm8Yvfd3v7a074tbVcGZrbxeFq8rP2/21ZZPA/vQ6Xrq619tLWd96HE2jHGHk29UfZn3VaIf8YYj0w7iGIIIVyXpXMpdzylPF6x+25vf+1p35a2IYSxsblSgdota//vtlUWz8P7UOn6amsf7WjnfWgxnIzRtPvTDqCIsnYu5Y6nlMcrdt/t7a897bP2OVF+fidZPA/vQ6Xrq619ZPFzkguWbiW1if+SlpQ270OL54iepLa6Lu0AJHV43ocWwxE9SZKknHJET5IkKadM9CRJknLKRE9S0YUQ1gwhPB5CmBRCeDmEcEkIIaQdl6SOI4TwZAhhfAhhQghhWAhh2bRjSoOJnqRSmAecGmPcEPgOsDXwo3RDktTB7BVj3DTG2AeYBpycdkBpMNGTtEAIYb0QwrWFfwXXhRBGNbPfRoURu5oQwvQQwrkhhE7178cYZ8QYxxb+PAeYAKxZlpOQVLGKdQ8CiDF+Vth3CaA70CFnn/pkDEkNfRvYA3gO6NLUDiGEFYDHgInA3sC6wKUk/3A8s4n9VwL2AXYtTciScqSo96AQwoPAlsDLwG9LFnWGubyKpAVCCEvEGOcX/jwMWDnGuEOjfU4HTiF5tuLnhW2nAIOBVeu3FbZ3Bf4FjIgxXlqWk5BUsYp9Dyq81wm4EJgZY7yk5CeRMZZuJS1Qf4NdjN2BhxvdTIcCSwPb128o3FxvA140yZPUEsW8BzXosw74O3BIUYKsMCZ6klprA2Byww0xxmlATeG9etcCX9BByyWSSmax96AQwgohhFUa7PJj4KWyRZghfkdPUmutAHzaxPZPCu8RQtgWOJzkxvpiYWWVG2OMV5QrSEm5tdh7UOG/d4QQugABmAQcW57wssVET1JbNPXl3lC/Pcb4n8LPklQKi7sHvQFsUdaIMsrSraTW+gRYvonty9H0v7IlqZi8B7WCiZ6k1prMwt/FI4SwJsk6VZObbCFJxeM9qBVM9CS11kPAbiGEZRpsGwh8CTyZTkiSOhDvQa3gd/QkLRBC6EayWCnA6sCyIYT9Cj8/GGOsAa4BjgOGhxAuBr5Jsn7VZY3Xr5Kk1vAeVHwumCxpgRBCL+DNZt5eJ8Y4tbDfRsBVQF+S78TcAAwurFclSW3iPaj4TPQkSZJyyu/oSZIk5ZSJniRJUk6Z6EmSJOWUiZ4kSVJOmehJkiTllImeJElSTpnoSZIk5ZSJnqR2C4l7Qwh3px2LJOkrJnpSBQghDA4hxMJrfgjhkxDCmBDC+SGEVRvt26uwX/9yxReTldd/BmwUQvh5U/uEEI4MIexTrpgaHXvjwjXZIY3jNyWEcHMIYWwr9h8cQphZyphKKYTQLYTwXghh+7Rjaa0QwtIhhA9CCN9POxaptUz0pMrxGcnjfrYBDgSGkyRX/wshbN5gvxmF/Z4uZ3CFZ0zuC5wfQliriV2OBFJJ9DLqPODQVux/A7BbaUIpi2OBN2OMFffQ+Rjjl8CVJL8zqaJ0TjsASS02L8b4XIOfHw4h/BV4CvhnCGH9GGNdjLEWeK7pLkorxjgRWKM9fYQQli78xZprMcbXW7n/O8A7JQqnpEIISwBHU9mJ0s3AOSGETWKM/0s7GKmlHNGTKliM8VPgFGBdYBdounQbQjgkhPB0COHjQtl3ZAhhi4Z91ZcSQwh7hhAmhhBqQggPhBBWDCGsV2gzq7BPn0ZtlwghnBZCeC2EUBtCeCWEcFiD90cBmwM/b1CCPrTw3tQQwqUhhLNCCO8Anzdot10I4clCLB+FEK4PISyzuOsSQvhNCOHtQrz3Az0bvd9keXtx5dQQwjmF8uMSjbb3L/S3XuHnToVS67TC9Xg5hHDwoo4VQlg+hHBDCGF6CGF2oe31Dd5fqHQbQtihvhwdQrgzhFAdQngjhPCbRsfpG0K4r9DvrBDCf0MIP2ni3HYIIUwoHHtMCGGrEMLMEMLgRvvtXfgMzC5ci0tCCEs2d80KdgRWJxmFbthXU5+bnzfaZ1QIYVgI4eDCfp+HEB4KISzyHxQhhEML1+e7hT5qCuf+3RBC9xDCTSGEzwrX7KDFxE+M8W1gDHDI4vaVssRET6p8I4F5wPcWsU8v4BZgf+BgkpGhp0II32y031rAucCZJKXWbYDrgKGF134klYChIYTQoN2VhTbXAf2Be4EbQgh7F97/DTAZeJCkrNwXeKBB+4OB7Qv7DQQIIWwLPA68VzjuCcAewE2LuhiFY14NjAB+BPwPuHFRbVphKLBKIdaGDgDGxRhfK/x8LnAGyfXYC/gPcNtiEorLgO2AE0lKtIOA2IKYrgfGk5TNRwFXhxC2avD+2oXj/xIYANwF3NQwlhDC6iS/mw9IrvW1wG3A0g0PFEI4gCRZG104r3NIPicXLibGnYBXYowfNdre8HOzJ3A3cGPjBBzYGjgG+G3heN8ttGmJvwO3Az8GAjAM+BswneRcnwduWVziWPAMsHMLjytlQ4zRly9fGX8Bg4GZi3h/BvDXwp97kSQI/ZvZdwmSZG0y8LsG228mSRjXbbDtkkJfhzTYtkdh24aFn9cD5gO/aHSc20iSn/qfxwI3NxHP1EL8SzXa/m9gZKNtOxaOvfEirsVo4KFG264vtNthUdeocA3GLuZ3MR64psHPXUm+P/l/hZ9XBGYBZzdq9yAwpbljAS8Bx7b0MwDsUDiHcxtsWxL4ELiomT5C4Xd/LfBEg+1/AGYCSzfYdkCh/8EN2r4F3NSoz8OAL4GVFhH7I8CdjbbVf25+3mj7LcCYBj+PKlzfFRpsO6EQ29KLOOahhX1+3mBb/Wf3xgbblgPmAr9uwf+Hh5L8P7LU4vb15SsrL0f0pHwIi3wzhA1DCHeHEN4H6kj+Ylsf6N1o16lx4e+O1Y9QPdHEttUL/92J5C/sOxr19TDQJ4TQku8CPx5jnN0g3m4ko353hBA6179IJpjMJSkDN3WenYDvkIwoNjS8id3b6p/Ajxuc1+7AMnx1/hsD3YA7m2jXO4TwjWb6/S9wcqHs3Pj3siiP1P8hxjgXeJUG35MMIawQQrgihPAWybWbSzIq1vAYWwKPxoW/G3lfo+P0Jhnxbfw7eQJYiuS8m7MqSSLZUP3n5u5G/T0ObFb4XdYbE2P8pMHPEwv/XZ3Fe7zBn7/2eY4xfkaSHLekr5lAJ6BHC/aVMsFET6pwIYSlgJWA95t5fxmSZGBN4CTg+yR/sY8n+Qu6oU8b/Tynie312+rbrkzyl99Hhe9tzQ4hzCaZJdqZ5C/5xWkc+wqFPv/CV8nJXKCWZNRqr8xE9QAABTdJREFUzWb66VE45geNtjf+uT2GkpzzjoWfBwLPxhinFX6u/z5g43Oq/3mFZvo9BrgH+B0wJYTwagjhwBbE09TvrOHv9eZCjH8AdiX53d/YaJ9VSZKdBQqJd3WDTSsX/vsgC/9O3ixsb+53QuFYtY221X9uPmvU380kv8OG36ts7nPZ+PPblKY+u4u7Zs2pP4eW7CtlgrNupcrXj+T/5Webeb8vyQjPLjHGyfUbQwjLFen4H5OUszYnGS1srMkEtJHG30X7tLBtMEli0dj0Zvr5sBBL41Gzxj/Xjx52abR9xUVGCcQY3yhMohgYQnia5HtvgxrsMqPBMRt+J22Vwn8/bqbfT4HjgONCMtnlFJLv9U2IyWzmViv8I2BP4JgY4zUNtjf+R/57NBqlKrStarCpPu4jgRebONybTWxr2Hb5JrbNA7YlGdlrrJjJebHUn0OTv0Mpi0z0pAoWQlgeuJikJPVYM7vVf6F+wYhKCGEbku+pjStCGE+QjMysGWP81yL2a+moCTHGWSGE54D1Y4zntjSQGGNdCOG/wN7ANQ3e+lGjXT8gGT3asH5DCKGKJCl+qwWHGkoy2eIJkuvbsEz7ElBDMvGlYewHkExIWGjkrJnzmBBCOBn4CbABX5UqW6srye+m4e9+GZKJFA2T6zHAL8LCS9vs1aivKcC7QK8Y4/W0zhRgnUbb6j83y8UYH21lf2npBXwUvz6pRMosEz2pcnQOIdTPrF2GZATt1yTfB/thjLGp0TRI1tSrBq4PIVxCMro3mOQv7XaLMU4JIVwD3F7ofzRffWdr/Rhj/TIrk4HdQgi7kYx0vbmYvzBPAR4PIcwnmSn5Bcl3xPYEzogxvtJMuwuA4SFZY/BukhmyP2wU8/wQwr3AiYXvrn1KMqOzpev33UFSCv0D8FSMsX4U7//buZ8XG+MojuPvz+xsxt6GLGysFCuyYCUrYiMloWjGhtIkbkMsRM2I/AHGYsZmmixMWdj4lVF2FJMSamwIsWKOxXke3ca9tzvurdt9fF7beeY731vfms/9nuccIuKTpHHgjKSfZBPKbrIRoGnXbXE7OE0GxQCOkE0dT9vc018i4oukOaAm6St5czZClksH6x4dJ+fc3ZE0RpZyR8jAulistSjpJDAhaRC4S4b3teQg7D0R8aPJVh4CuyQNRES5XnluJotz84w8N+uBdRFx+F8/d6ck1chGpaX/IzeSnbdmfcNBz6x/rCTLs0HOmpsHbgHXImKh2S9FxEdJe4ErZJPCa+AoGaS6ZQh4RYaT0WJ/L8gxFqULFC/zkyHjIPk+VrN9P5C0lRzhMUHe/rwFZmlRDo6IaUnHyaBygOzaPEQ2h9QbJkd03AA+AxfJcTKtmgrKv/FO0iOy7HiuwSM1six5jCzZzgP7I2KyxbKPya7ONWQJ/DmwI3JQcif2kZ/zJhmwr5NfDobLByLig6SdwFWyceUl2U17j7q5hhExVQTG08XPfwFvyFE25ftvjcyQI282k93Upfpzc57G56YXBsjz9kfRKLKd/EJg1jcU0c6YJjMz+59I2kKGsm0Rcb8L680A7yNiqOPN9UBxE30bWBUR33u9H7N2OeiZmRmSLpG3iAvk6J2z5A3ghrLc2uH6m8hRJ6uXjErpC5JmgScRMdrrvZgth0u3ZmYG2bhxmSw1fyNH8pzoRsgDiIg5SafI8n1fBT1JK8jS+liv92K2XL7RMzMzM6soD0w2MzMzqygHPTMzM7OKctAzMzMzqygHPTMzM7OKctAzMzMzqygHPTMzM7OK+g1/LHSphqqhigAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 720x504 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "sb_openfield['exp_line'] = np.exp(lsb_openfield['line'])\n",
    "font = {'family' : 'normal',\n",
    "        'weight' : 'normal',\n",
    "        'size'   : 15}\n",
    "import matplotlib\n",
    "matplotlib.rc('font', **font)\n",
    "ax = sb_openfield.plot.scatter('radius','sb', loglog=True, figsize=(10,7),s=50)\n",
    "ax.set_xlim((20, 1700))\n",
    "ax.set_ylim((sb_bocage.sb.min()-0.00001, sb_bocage.sb.max()+0.01))\n",
    "sb_openfield.plot.line('radius', 'exp_line',color='Red',ax=ax, legend=False)\n",
    "plt.grid(which='both')\n",
    "ax.set_xlabel('Diamètre du voisinage (en m.)')\n",
    "ax.set_ylabel('Mr(i)')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Rasterization des fonds à 4 m"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'openfield' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-5-3b1ef4aa4b7a>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mr_openfield\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mrasterizeBlackWhite\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mopenfield\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m4\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;31mNameError\u001b[0m: name 'openfield' is not defined"
     ]
    }
   ],
   "source": [
    "r_openfield = rasterizeBlackWhite(openfield, 4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Calcul des voisinages"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "neigh_openfield = neighbors_near_point(r_openfield, np.array([25,50,100,200,400,800,1600])) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Calcul de la GWFA en un point"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def Gwfa_en_un_point(df, x_centre, y_centre, q):\n",
    "    distance = (df.x - x_centre)** 2 + (df.y - y_centre)**2\n",
    "    weight = (1 - (distance / 16000**2))**2\n",
    "    \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = r_openfield.copy()\n",
    "list_of_radius = np.array([25,50,100,200,400,800,1600])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "neighbors = np.zeros([df.shape[0],list_of_radius.size],dtype=np.float)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "tree = KDTree(df[['x', 'y']]) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "dist, ind = tree.query(df[['x', 'y']], np.max(list_of_radius))   "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "for i, r in enumerate(list_of_radius):\n",
    "    neighbors[:,i] =list(map(lambda x : np.sum( x <= r), dist)) "
   ]
  },
  {
   "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.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}