Topological colouring.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {},
   "outputs": [],
   "source": [
    "import geopandas as gpd\n",
    "import libpysal\n",
    "import networkx as nx\n",
    "import pandas as pd\n",
    "import operator\n",
    "import sys\n",
    "from collections import defaultdict\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.spatial import Voronoi, voronoi_plot_2d\n",
    "from random import random\n",
    "import numpy as np\n",
    "import seaborn as sns\n",
    "from shapely.geometry import Polygon"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "world = gpd.read_file(gpd.datasets.get_path('naturalearth_lowres'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def balanced(features, sw, balance=\"count\", min_colors=4):\n",
    "    feature_colors = {}\n",
    "    # start with minimum number of colors in pool\n",
    "    color_pool = set(range(1, min_colors + 1))\n",
    "\n",
    "    # calculate count of neighbours\n",
    "    neighbour_count = sw.cardinalities\n",
    "\n",
    "    # sort features by neighbour count - we want to handle those with more neighbours first\n",
    "    sorted_by_count = [feature_id for feature_id in sorted(neighbour_count.items(),\n",
    "                                                           key=operator.itemgetter(1),\n",
    "                                                           reverse=True)]\n",
    "    # counts for each color already assigned\n",
    "    color_counts = defaultdict(int)\n",
    "    color_areas = defaultdict(float)\n",
    "    for c in color_pool:\n",
    "        color_counts[c] = 0\n",
    "        color_areas[c] = 0\n",
    "\n",
    "    for (feature_id, n) in sorted_by_count:\n",
    "        # first work out which already assigned colors are adjacent to this feature\n",
    "        adjacent_colors = set()\n",
    "        for neighbour in sw.neighbors[feature_id]:\n",
    "            if neighbour in feature_colors:\n",
    "                adjacent_colors.add(feature_colors[neighbour])\n",
    "\n",
    "        # from the existing colors, work out which are available (ie non-adjacent)\n",
    "        available_colors = color_pool.difference(adjacent_colors)\n",
    "\n",
    "        feature_color = -1\n",
    "        if len(available_colors) == 0:\n",
    "            # no existing colors available for this feature, so add new color to pool and repeat\n",
    "            min_colors += 1\n",
    "            return balanced(features, sw, balance, min_colors)\n",
    "        else:\n",
    "            if balance == \"count\":\n",
    "                # choose least used available color\n",
    "                counts = [(c, v) for c, v in color_counts.items() if c in available_colors]\n",
    "                feature_color = sorted(counts, key=operator.itemgetter(1))[0][0]\n",
    "                color_counts[feature_color] += 1\n",
    "            elif balance == \"area\":\n",
    "                areas = [(c, v) for c, v in color_areas.items() if c in available_colors]\n",
    "                feature_color = sorted(areas, key=operator.itemgetter(1))[0][0]\n",
    "                color_areas[feature_color] += features.iloc[feature_id].geometry.area\n",
    "            elif balance == \"centroid\":\n",
    "                min_distances = {c: sys.float_info.max for c in available_colors}\n",
    "                this_feature_centroid = features.iloc[feature_id].geometry.centroid\n",
    "\n",
    "                # find features for all available colors\n",
    "                other_features = {f_id: c for (f_id, c) in feature_colors.items() if c in available_colors}\n",
    "\n",
    "                # loop through these, and calculate the minimum distance from this feature to the nearest\n",
    "                # feature with each assigned color\n",
    "                for other_feature_id, c in other_features.items():\n",
    "\n",
    "                    other_geometry = features.iloc[other_feature_id].geometry\n",
    "                    other_centroid = other_geometry.centroid\n",
    "\n",
    "                    distance = this_feature_centroid.distance(other_centroid)\n",
    "                    if distance < min_distances[c]:\n",
    "                        min_distances[c] = distance\n",
    "\n",
    "                # choose color such that minimum distance is maximised! ie we want MAXIMAL separation between\n",
    "                # features with the same color\n",
    "                feature_color = sorted(min_distances, key=min_distances.__getitem__, reverse=True)[0]\n",
    "             \n",
    "            elif balance == \"distance\":\n",
    "                min_distances = {c: sys.float_info.max for c in available_colors}\n",
    "                this_feature = features.iloc[feature_id].geometry\n",
    "\n",
    "                # find features for all available colors\n",
    "                other_features = {f_id: c for (f_id, c) in feature_colors.items() if c in available_colors}\n",
    "\n",
    "                # loop through these, and calculate the minimum distance from this feature to the nearest\n",
    "                # feature with each assigned color\n",
    "                for other_feature_id, c in other_features.items():\n",
    "\n",
    "                    other_geometry = features.iloc[other_feature_id].geometry\n",
    "\n",
    "                    distance = this_feature.distance(other_geometry)\n",
    "                    if distance < min_distances[c]:\n",
    "                        min_distances[c] = distance\n",
    "\n",
    "                # choose color such that minimum distance is maximised! ie we want MAXIMAL separation between\n",
    "                # features with the same color\n",
    "                feature_color = sorted(min_distances, key=min_distances.__getitem__, reverse=True)[0]\n",
    "\n",
    "        feature_colors[feature_id] = feature_color\n",
    "\n",
    "    return feature_colors"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [],
   "source": [
    "class Graph:\n",
    "\n",
    "    def __init__(self, sort_graph=True):\n",
    "        self.sort_graph = sort_graph\n",
    "        self.node_edge = {}\n",
    "\n",
    "    def add_edge(self, i, j):\n",
    "        ij = [i, j]\n",
    "        if self.sort_graph:\n",
    "            ij.sort()\n",
    "        (i, j) = ij\n",
    "        if i in self.node_edge:\n",
    "            self.node_edge[i].add(j)\n",
    "        else:\n",
    "            self.node_edge[i] = {j}\n",
    "\n",
    "    def make_full(self):\n",
    "        g = Graph(sort_graph=False)\n",
    "        for k in self.node_edge.keys():\n",
    "            for v in self.node_edge[k]:\n",
    "                g.add_edge(v, k)\n",
    "                g.add_edge(k, v)\n",
    "        return g\n",
    "\n",
    "    \n",
    "def compute_graph(features, create_id_graph=False, min_distance=0):\n",
    "    \"\"\" compute topology from a layer/field \"\"\"\n",
    "    s = Graph(sort_graph=False)\n",
    "    id_graph = None\n",
    "    if create_id_graph:\n",
    "        id_graph = Graph(sort_graph=True)\n",
    "\n",
    "    index = features.sindex\n",
    "\n",
    "    i = 0\n",
    "    for feature_id, f in features.iterrows():\n",
    "\n",
    "        g = f.geometry\n",
    "        if min_distance > 0:\n",
    "            g = g.buffer(min_distance, 5)\n",
    "\n",
    "        possible_matches_index = list(index.intersection(g.bounds))\n",
    "        possible_matches = features.iloc[possible_matches_index]\n",
    "        precise_matches = possible_matches.loc[possible_matches.intersects(g)]\n",
    "\n",
    "        for ix, r in precise_matches.iterrows():\n",
    "            s.add_edge(feature_id, ix)\n",
    "            s.add_edge(ix, feature_id)\n",
    "            if id_graph:\n",
    "                id_graph.add_edge(feature_id, ix)\n",
    "\n",
    "    for feature_id, f in features.iterrows():\n",
    "\n",
    "        if feature_id not in s.node_edge:\n",
    "            s.add_edge(feature_id, None)\n",
    "\n",
    "    return s, id_graph    \n",
    "\n",
    "\n",
    "def g_balanced(features, graph, balance=\"count\", min_colors=4):\n",
    "        feature_colors = {}\n",
    "        # start with minimum number of colors in pool\n",
    "        color_pool = set(range(1, min_colors + 1))\n",
    "\n",
    "        # calculate count of neighbours\n",
    "        neighbour_count = defaultdict(int)\n",
    "        for feature_id, neighbours in graph.node_edge.items():\n",
    "            neighbour_count[feature_id] += len(neighbours)\n",
    "\n",
    "        # sort features by neighbour count - we want to handle those with more neighbours first\n",
    "        sorted_by_count = [feature_id for feature_id in sorted(neighbour_count.items(),\n",
    "                                                               key=operator.itemgetter(1),\n",
    "                                                               reverse=True)]\n",
    "        # counts for each color already assigned\n",
    "        color_counts = defaultdict(int)\n",
    "        color_areas = defaultdict(float)\n",
    "        for c in color_pool:\n",
    "            color_counts[c] = 0\n",
    "            color_areas[c] = 0\n",
    "\n",
    "        for (feature_id, n) in sorted_by_count:\n",
    "\n",
    "            # first work out which already assigned colors are adjacent to this feature\n",
    "            adjacent_colors = set()\n",
    "            for neighbour in graph.node_edge[feature_id]:\n",
    "                if neighbour in feature_colors:\n",
    "                    adjacent_colors.add(feature_colors[neighbour])\n",
    "\n",
    "            # from the existing colors, work out which are available (ie non-adjacent)\n",
    "            available_colors = color_pool.difference(adjacent_colors)\n",
    "\n",
    "            feature_color = -1\n",
    "            if len(available_colors) == 0:\n",
    "                # no existing colors available for this feature, so add new color to pool and repeat\n",
    "                min_colors += 1\n",
    "                return g_balanced(features, graph, balance, min_colors)\n",
    "            else:\n",
    "                if balance == \"count\":\n",
    "                    # choose least used available color\n",
    "                    counts = [(c, v) for c, v in color_counts.items() if c in available_colors]\n",
    "                    feature_color = sorted(counts, key=operator.itemgetter(1))[0][0]\n",
    "                    color_counts[feature_color] += 1\n",
    "                elif balance == \"area\":\n",
    "                    areas = [(c, v) for c, v in color_areas.items() if c in available_colors]\n",
    "                    feature_color = sorted(areas, key=operator.itemgetter(1))[0][0]\n",
    "                    color_areas[feature_color] += features.loc[feature_id].geometry.area\n",
    "                elif balance == \"centroid\":\n",
    "                    min_distances = {c: sys.float_info.max for c in available_colors}\n",
    "                    this_feature_centroid = features.loc[feature_id].geometry.centroid\n",
    "\n",
    "                    # find features for all available colors\n",
    "                    other_features = {f_id: c for (f_id, c) in feature_colors.items() if c in available_colors}\n",
    "\n",
    "                    # loop through these, and calculate the minimum distance from this feature to the nearest\n",
    "                    # feature with each assigned color\n",
    "                    for other_feature_id, c in other_features.items():\n",
    "                        \n",
    "                        other_geometry = features.loc[other_feature_id].geometry\n",
    "                        other_centroid = other_geometry.centroid\n",
    "\n",
    "                        distance = this_feature_centroid.distance(other_centroid)\n",
    "                        if distance < min_distances[c]:\n",
    "                            min_distances[c] = distance\n",
    "\n",
    "                    # choose color such that minimum distance is maximised! ie we want MAXIMAL separation between\n",
    "                    # features with the same color\n",
    "                    feature_color = sorted(min_distances, key=min_distances.__getitem__, reverse=True)[0]\n",
    "\n",
    "            feature_colors[feature_id] = feature_color\n",
    "\n",
    "        return feature_colors"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 107,
   "metadata": {},
   "outputs": [],
   "source": [
    "def pysal_nx(gdf, strategy='largest_first', contiguity='queen'):\n",
    "    if contiguity == 'queen':\n",
    "        sw = libpysal.weights.Queen.from_dataframe(gdf)\n",
    "    elif contiguity == 'rook':\n",
    "        sw = libpysal.weights.Rook.from_dataframe(gdf)\n",
    "    color = nx.greedy_color(sw.to_networkx(), strategy=strategy)\n",
    "    return pd.Series(color)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 108,
   "metadata": {},
   "outputs": [],
   "source": [
    "def pysal_balanced(gdf, balance=\"count\", min_colors=4, contiguity='queen'):\n",
    "    if contiguity == 'queen':\n",
    "        sw = libpysal.weights.Queen.from_dataframe(gdf)\n",
    "    elif contiguity == 'rook':\n",
    "        sw = libpysal.weights.Rook.from_dataframe(gdf)\n",
    "    \n",
    "    return pd.Series(balanced(gdf, sw, balance=balance, min_colors=min_colors))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 109,
   "metadata": {},
   "outputs": [],
   "source": [
    "def qgraph_balanced(gdf, balance=\"count\", min_colors=4):\n",
    "    topology, id_graph = compute_graph(gdf)\n",
    "    colors = g_balanced(gdf, topology, balance, min_colors)\n",
    "    return pd.Series(colors)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 110,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "904 ms ± 20.2 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n"
     ]
    }
   ],
   "source": [
    "%%timeit\n",
    "qgraph_balanced(world)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 111,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/martin/anaconda3/envs/geo_dev/lib/python3.7/site-packages/libpysal/weights/weights.py:165: UserWarning: The weights matrix is not fully connected: \n",
      " There are 25 disconnected components.\n",
      " There are 21 islands with ids: 0, 19, 20, 22, 23, 45, 46, 47, 78, 89, 134, 135, 136, 137, 138, 140, 144, 147, 155, 159, 175.\n",
      "  warnings.warn(message)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "46.9 ms ± 1.18 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)\n"
     ]
    }
   ],
   "source": [
    "%%timeit\n",
    "pysal_balanced(world)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 112,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "56.4 ms ± 3.77 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)\n"
     ]
    }
   ],
   "source": [
    "%%timeit\n",
    "pysal_nx(world)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 128,
   "metadata": {},
   "outputs": [],
   "source": [
    "times = pd.DataFrame(index=['robust_count', 'pysal_count', 'pysal_lf'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 129,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "100\n",
      "robust_count:  0.49878664016723634 s;  14 colors\n",
      "pysal_count:  0.008840847015380859 s;  14 colors\n",
      "pysal_lf:  0.014633893966674805 s;  14 colors\n",
      "1000\n",
      "robust_count:  5.006178474426269 s;  17 colors\n",
      "pysal_count:  0.09947896003723145 s;  17 colors\n",
      "pysal_lf:  0.6779930114746093 s;  17 colors\n",
      "10000\n",
      "robust_count:  45.33429980278015 s;  20 colors\n",
      "pysal_count:  0.9851454257965088 s;  20 colors\n",
      "pysal_lf:  1.6736697673797607 s;  20 colors\n"
     ]
    }
   ],
   "source": [
    "for number in [100, 1000, 10000]:\n",
    "    print(number)\n",
    "    points = np.random.rand(number,2)\n",
    "    voronoi_diagram = Voronoi(points)\n",
    "    regions = pd.DataFrame()\n",
    "    regions[\"region\"] = voronoi_diagram.point_region\n",
    "    vertices = []\n",
    "    for region in regions.region:\n",
    "        vertices.append(voronoi_diagram.regions[region])\n",
    "    regions[\"vertices\"] = vertices\n",
    "    polygons = []\n",
    "    for region in regions.vertices:\n",
    "        try:\n",
    "            polygons.append(Polygon(voronoi_diagram.vertices[region]))\n",
    "        except Exception:\n",
    "            continue\n",
    "    regions_gdf = gpd.GeoDataFrame(geometry=polygons)\n",
    "    regions_gdf['geometry'] = regions_gdf.geometry.buffer(0)\n",
    "    timer = []\n",
    "    for run in range(5):\n",
    "        s = time()\n",
    "        colors = qgraph_balanced(regions_gdf)\n",
    "        e = time() - s\n",
    "        timer.append(e)\n",
    "    times.loc['robust_count', number] = np.mean(timer)\n",
    "    print('robust_count: ', np.mean(timer), 's; ', np.max(colors), 'colors')\n",
    "    timer = []\n",
    "    for run in range(5):\n",
    "        s = time()\n",
    "        colors = pysal_balanced(regions_gdf)\n",
    "        e = time() - s\n",
    "        timer.append(e)\n",
    "    times.loc['pysal_count', number] = np.mean(timer)\n",
    "    print('pysal_count: ', np.mean(timer), 's; ', np.max(colors), 'colors')\n",
    "    timer = []\n",
    "    for run in range(5):\n",
    "        s = time()\n",
    "        colors = pysal_nx(regions_gdf)\n",
    "        e = time() - s\n",
    "        timer.append(e)\n",
    "    times.loc['pysal_lf', number] = np.mean(timer)\n",
    "    print('pysal_lf: ', np.mean(timer), 's; ', np.max(colors) + 1, 'colors')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 130,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.axes._subplots.AxesSubplot at 0x139b3c710>"
      ]
     },
     "execution_count": 130,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmEAAAJBCAYAAAAOWWJ9AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nOzdd3hUZf7+8XtKEkqAQApFmhSRFpoiooJtEY1BEVDQVdeKDcW1fF123XXX5ae7FsCCa2Wtq9IsiKjYVgVXaQlNpNeQhPSQOjPn90dImIT0ZObMmXm/rssrpMzJJzmE3D7POXPbDMMwBAAAAL+ymz0AAABAKCKEAQAAmIAQBgAAYAJCGAAAgAkIYQAAACYghAEAAJiAEAYAAGACp1mfOCvrqDwenqIs0EVHRyojI9/sMdAEnENr4/xZH+fQ2ux2m9q3b+2TY5sWwjwegxBmEZwn6+McWhvnz/o4h6gO25EAAAAmIIQBAACYwLTtyKoMw1B+fo4KC/Pl8bjNHgfHpKXZZbc71b59rByOgPnrAgCA5QXMb9WsrHTZbDZ16NBRDodTNpvN7JEgyeGwKScnW1lZ6YqJ6Wz2OAAABI2A2Y4sKSlSVFS0nM4wAlgAsdlsat26rVyuErNHAQAgqARMCJMM2WwBNA4qEIoBAGh+pB4AAAATEMIAAABMQAhrguXLP9bs2Y80+Tg//PCd3n33raYP1ARbtmzS/PnPmDoDAAChJGDujgxlv/yyxewRtGfPbmVlZZo9BgAAISNgQ9gPG1P0fXKKT459dnxnnTW49qdbWLdujV544Rm53R516XKSwsLCtGPHr7Lb7Zo69be6+OJLJUkHDuzXnXfeotzcHI0efY5uu+0uHT6cohkzpmvRoo8lSa+++qIk6frrb9Jjj/1Vu3btlCRNnDhFgwcP0YcfLpEkderUWQkJE6qdp7i4WE8//Q8lJ2+Q0+nU7353sy64YJw2bdqoefOeVElJiaKiovTAA7PUtWs33XXXrbrxxls1fPhpSkk5VDHP7NmPqHXrSG3btlVHjqTrd7+7WWPGnKdXXvmXCgsL9frrr+r6629qlu8zAACoWcCGsECwf/8+LVq0TG+88apKSkr05pvvKzs7W7fccr369u0nSUpJOaQFC95RZGSk7r77Nn3//bfq0+eUao+3cWOScnNztWDBOzpyJF0vvPCsJkyYqMsuu0KSagxgkrR48XsqLCzU228vUlZWpu655w6NGXOeHnlklh599HH17z9QX321Uo888ke98sobtX5daWmpmj//Fe3atVMzZkxXQsIE3XzzbVq/fi0BDAAAPwnYEHbW4LpXq3ytW7ceioyM1Nq1a/TQQw9LkqKionTOOWO0fv1atW7dWmefPUbt27eXJJ1//m+0fv3aGkNYr169tW/fXv3+93dp1KizdOed99R7lg0b1mnChImy2+2Kjo7RW2+9r127dqhNmzbq33/gsc9/of75z9nKz8+v9VgjR54hm82mXr16Kzc3p94zAACA5sOF+bWIiIiQJBmGp9LbDUNyu12SJIfD4fV2j5zOsmf7Nwyj4u0uV9nHtmsXpTfffF+TJl2lffv26sYbf6u8vLx6zVJWGXT8+boOHNgvj8eo5iMNeTzuSs/tVf75y4WHl31dPP8XAADmIYTVw/Dhp+uTTz6UJGVnZ+u7777RsGGnSZJWr/5BeXl5Ki4u1sqVn+u0085QZGQb5ebmKisrSyUlJfrf/1ZLkr7//ls9+uifNXr02Zo58361bNlSaWmpcjgccrtr78scOnSYvvrqCxmGoaysTN11163q1KmzcnJytHXrZknSl19+oY4dO6tt23Zq1y5Ku3eXXXv23Xff1Pk11mcGAADQfAJ2OzKQ3HDDzXrqqX/ouuuuksfj0XXX3ah+/U7Vzp3b1aNHTz3wwD3Kz8/ThRdepJEjR0mSrrnmOt1yy3WKi+uoAQPKtgtHjTpL33zzla699kqFh4froosuUe/efZSXl6vZsx9Rhw4dNHny1GpnmDhxiubOfUK/+900SdK99z6gyMhI/e1vj+npp/+poqJCtW3bTn/722MVn3/27Ef0yScf6Zxzzq3za+zff6Bee+0lvfDCs7r99hnN8F0DAAC1sRne+2Z+lJGRX2k77fDhverUqYcZo6AWTqddLpeH82NhsbFtlJ5ev21vBB7On/VxDq3NbrcpOjrSJ8dmJSyAfPnl53rzzX9X+75///sd/w4DAAB8ihAWQC64YJwuuGCc2WMAAAA/4MJ8AAAAExDCAAAAapCZW+SzYxPCAAAAqvHT1lQ98e4Gnx2fEAYAAODFMAwt/3Gv/vXhZnWN882dkRIhzO9SUg5p8uREU2f44Yfv9O67b5k6AwAAgcjt8eiNz7Zp0Tc7NbJ/nKYnDvDZ5+LuyBD0yy9bzB4BAICAU1js0gsfbtKmXZlKOLOHJo7pJafDd+tVARvCSn/9QaXb/uuTY4f1G6OwU86q9WPWrVuj119/VQ6HUykpBzVgwEDFxnaUzWbT9Ol3SpJmz35Eo0aNltvt0TvvvCG73a4uXbro4YcflcPh0FNPPa5du3YqMzNTffr00SOPzK7XfLm5OXrssUe1b98ehYWFa8aMezVixOn64Yfv9PLLL8gwPOrS5SQ98MAsdegQrcmTE/Xssy+qc+cuWrdujV577SU999xLuuuuWzVgwEAlJW1QdnaWZs58QJ06ddaHHy6RJHXq1FkJCROa9s0EACAIZOYWad6iZB1MP6rrx/fT2KEn+fxzsh1Zi40bkzVz5n16553FKi4uUVRUlFau/EyGYaioqEhr1/6ss88eq5dffkFz5jyn1157S507n6R9+/Zo06ZkOZ1hevHFBXrvvaXKy8vT6tU/1Ovzvvzyv9S1aze9/fYiPfzw3/TSS/OVlZWpJ574f3rssSf1+uvvavDgIXr66X/WeazSUpdefHGBZsz4vV5++QWdfHIvXXbZFbrssisIYAAASNqXmqfZb65VenahZk6J90sAkwJ4JSzslLPqXK3ytaFDh6l7956SpPHjL9FHHy1Vp06dtWHDOqWmHtbo0WcrIiJCZ511jm6//SaNGXOuxo49X3379pMktW3bTosXv699+/bowIH9KiwsrNfn3bBhrf7yl7JVs969++jFFxfohx++U//+A9W5cxdJ0oQJV9T47PrezjjjTElSr169lZeX28DvAAAAwW3jrgzN/2CTWkU49dA1w9W9Yxu/fe6ADWGBwOFwVPzZ4zHkcDiUkDBBX3yxQqmpqbrxxlslSTNn3q8dOy7T6tXf69FHH9aNN96q1q1b65VXXtSUKVN1ySUTlJ2drfrWdDqdTtlstorX9+7dI8PwVPoYwzDkdrslSTabreLYbrer0seFh4ef8DEAAED6ZsNBvfXZr+oa21r3TBmi9m0i/Pr52Y6sRXLyBqWnp8nj8WjFik90xhmjdd55F2jt2p+VmXlEAwcOksvl0tSpExUVFaVrr71B48cn6Ndft2nNmp90/vkXKiFhgiIjI7V+/Vp5PO56fd4hQ4Zr5crPJJUFsPvum6EBAwZpy5aNSkk5JEn66KMlGj58hCSpXbso7d69S5L03Xff1nl8h8NREeAAAAg1HsPQwm926I0V2zTw5A76v2uG+z2ASayE1SomJlZ///tflJ6eptNPP0OJiZfL4XBo0KDB6tWrj6SyVaubbpqumTPvVEREhNq3b68//vERZWZm6q9//aNWrvxMTmeYBg+O16FDhzRiRN2f96abpusf//i7rr9+mhwOhx5++G/q0CFaDzzwR82adb9KS13q1KmTHnroz8c+/lbNmfOEFix4WSNHjqrz+EOHDtfs2Y+oQ4cOmjx5apO+RwAAWEmpy61Xlm3Vz7+k6dyhXXTNuFPksJuzJmUzTNqjysjIl8dz/FMfPrxXnTr1MGOUannfZVjOMAwVFBzV9Ok3at68+YqOjjFxQv9wOu1yuTwBd35Qf7GxbZSenmf2GGgkzp/1cQ4DR15BiZ5dslE7DuRoynm9NX5k90qX/1THbrcpOto3T9jKSlgDbN26Wffdd7duuOGWJgew9957W59++skJb4+JidGTTz7TpGMDAIDKUjMLNGdhkjJzi3X75YN0+qlxZo/EShhqx0qY9fF/4dbG+bM+zqH5th/I1rOLN0qS7p4Urz5d29X7sayEAQAANMJPW1P1yrKtim4boZlXDlHH9q3MHqkCIQwAAAQdwzC04n/7tPCbnerTtZ3unhSvyJZhZo9VCSEMAAAEFbfHo7c//1XfbDikkf3jdFNCf4U5HXU/0M8IYQAAIGhUV8Jtr+MOSLPwZK1+lpJySJMnJ9b6Ma+++qJeffVFSdJHHy3V5ZdfrOefn+eP8QAAsKzM3CI9/vY6bdmdpevH99Oksb0DNoBJrIQFvJUrP9OsWX+p15OwAgAQqval5mneomQVFrs0c0q8BvWKNnukOhHCarBu3Rq9/vqrcjicSkk5qAEDBio2tqNsNpumT79TkjR79iMaNWq03G6P3nnnDdntdnXp0kUPP/yoHA6Hnnrqce3atVOZmZnq06ePHnlkdoNmWLDgZW3dullPPfW4Zs68X2eeebYvvlQAACzNzBLupgjYEPa/lLVanfKzT459ZufTdUbnuvuDNm5M1r///ba6deuhhx9+SFFRUVq8+H3deusdKi4u1tq1P+v++/+g3/72Sr300gK1b99Bzz8/T/v27dHRo0fldIbpxRcXyOPx6O67b9Pq1T+oX7/+9Z7zhhtu0dq1P+vGG2/V8OGnNeVLBgAgKJWXcJ8U21ozTSjhboqADWGBYOjQYerevackafz4S/TRR0vVqVNnbdiwTqmphzV69NmKiIjQWWedo9tvv0ljxpyrsWPPV9++/SRJbdu20+LF72vfvj06cGC/CgsLTfxqAAAIHh7D0JJvd2n5j3s1qFcH3X7ZILWMsFasCdhpz+g8ol6rVb7kcBy/ndXjMeRwOJSQMEFffLFCqampuvHGWyVJM2ferx07LtPq1d/r0Ucf1o033qrWrVvrlVde1JQpU3XJJROUnZ0tk8oJAAAIKqUut179ZKt+2mp+CXdTWG9iP0pO3qD09DR5PB6tWPGJzjhjtM477wKtXfuzMjOPaODAQXK5XJo6daKioqJ07bU3aPz4BP366zatWfOTzj//QiUkTFBkZKTWr18rj8dt9pcEAICl5RWU6Il3N+inrWmacm5vXXtRP0sGMCmAV8ICQUxMrP7+978oPT1Np59+hhITL5fD4dCgQYPVq1cfSZLT6dRNN03XzJl3KiIiQu3bt9cf//iIMjMz9de//lErV34mpzNMgwfH69ChQxph7uIeAACWlZpVoLnvJykjt1i3XTZQI/t3NHukJqHAuwbr1q3Ra6+9pOeee6nibYZhqKDgqKZPv1Hz5s1XdHSMiRP6BwXe1kd5sLVx/qyPc9g8dhzI0TOLkyU1vIS7KSjwDhBbt27WfffdrRtuuKXJAey9997Wp59+csLbY2Ji9OSTzzTp2AAABJPyEu4ObSN0b4CVcDcFK2GoFSth1sf/hVsb58/6OIeNV7WEe8YVg9WmVbhfZ2AlDAAAhBSrlHA3RQCFMJsMwyObzZp3OAQznloDAOBP3iXcl4zqoSvGBm4Jd1METAgLD2+h7OwjatOmvRwOp2xB+M22IsMwdPRorpxO/y7/AgBCU1ZeseYuTNLB9KO6fnw/jR16ktkj+UzAhLD27WOVn5+jzMxUnk8rgNjtdtntTrVvH2v2KACAIFdewl1goRLupgiYEGaz2dSmTZTatIkyexR44YJSAIA/bNqVoeePlXD/wUIl3E0RMCEMAACEpm83HNSbx0q475kcrw5tW5g9kl8QwgAAgCmCoYS7KULnKwUAAAHDu4R77NAu+q1FS7ibghAGAAD8Kr+wVM8sTtaOAzmacm5vjT+je0g+KwIhDAAA+E2wlXA3BSEMAAD4hXcJ9wPThqpv19B+RgRCGAAA8Lmff0nTyx9vKSvhnjJEHTsERwl3UxDCAACAzwRCCXegIoQBAACfcHs8evuL7fpm/UGdfmqcbr40+Eq4m4IQBgAAml1hsUv/+nCzNu7K0MWjumvS2N5BWcLdFIQwAADQrLLyijVvYZIOpB/VdeP76dwgLuFuCkIYAABoNvvT8jV3YZIKil26Z0q8Bgd5CXdTEMIAAECz2LQrQ/M/2KSWIVTC3RSEMAAA0GTlJdxdYlpr5pTQKeFuCkIYAABoNI9haOl/d+mT1aFZwt0UfJcAAECjUMLdNIQwAADQYJRwNx0hDAAANAgl3M2DEAYAAOqNEu7mQwgDAAD1Qgl38yKEAQCAWhmGoRU/7dPCr3eqz0ntNGMSJdzNgRAGAABqRAm379T7PtJ//OMfeuihhyRJq1atUmJiosaNG6c5c+b4bDgAAGCewmKXnlm0Ud+sP6iLR3XX9MsGEsCaUb1C2OrVq7V06VJJUlFRkWbNmqX58+dr+fLl2rRpk7799lufDgkAAPwrK69Y/3h7nTbvztR14/tpyrl9ZOcpKJpVnSEsOztbc+bM0W233SZJSk5OVo8ePdStWzc5nU4lJiZqxYoVPh8UAAD4x/60fP39jTVKzS7UPVPide7Qk8weKSjVeU3Yn//8Z917771KSUmRJKWlpSk2Nrbi/XFxcUpNTfXdhAAAwG8o4fafWkPYwoUL1blzZ5155plasmSJJMnj8VR6RlzDMBr1DLnR0ZENfgzMERvLD6DVcQ6tjfNnfVY5h5/9uFfzFyere8c2+svNoxQT1dLskYJarSFs+fLlSk9P12WXXaacnBwVFBTo4MGDcjiOX5SXnp6uuLi4Bn/ijIx8eTxGwyeGX8XGtlF6ep7ZY6AJOIfWxvmzPiucw0ol3Cd30O2XD5JR6gr4uf3Bbrf5bOGo1hC2YMGCij8vWbJEP/30k/76179q3Lhx2rt3r7p27aply5Zp0qRJPhkOAAD4lncJ95ghZSXcTgcl3P7Q4OcJi4iI0OOPP64ZM2aouLhYY8eO1fjx430xGwAA8KH8wlI9uzhZ2w/kaPK5vXUxJdx+ZTMMw5Q9QbYjrcEKy+ioHefQ2jh/1heo59C7hPvmS/tTwl0D07YjAQBA8NlxMEfPLEqWYRi6f+pQndKNEm4zEMIAAAgha35J00uUcAcEQhgAACGAEu7AQwgDACDIeZdwn3ZqnG6hhDsgEMIAAAhiRSUu/evDzUremaGLR3XXpLG96YAMEIQwAACCVFZeseYtTNKB9KO67qJ+OncYHZCBhBAGAEAQ2p+Wr7kLk1RQ7NLdk+MV3zva7JFQBSEMAIAgs2l3huYv3aQW4Q5KuAMYIQwAgCDy36RDemPFNnWJaa2ZU+LVoW0Ls0dCDQhhAAAEgepKuFtG8Gs+kHF2AACwuFKXR68t36r/bUmlhNtCCGEAAFhYfmGpnlucrF8p4bYcQhgAABaVllWgOcdKuG+7bCAl3BZDCAMAwIIo4bY+QhgAABZTUcLdJkL3XkkJt1URwgAAsAjDMPTZT/v1/tc71Oekdrpr0mC1pYTbsghhAABYgNvj0TtfbNfXlHAHDUIYAAABrlIJ9xndNelcSriDASEMAIAAlpVXrHmLknQgjRLuYEMIAwAgQB1Iy9ccSriDFiEMAIAARAl38COEAQAQYCjhDg2EMAAAAgQl3KGFMwsAQACoXMLdWb8d148S7iBHCAMAwGTeJdyTxvbSJaN6UMIdAghhAACYKC2rQHMWJisjp5AS7hBDCAMAwCSVS7iHUcIdYghhAACYYM0vaXp52Ra1j4zQzCuHqBMl3CGHEAYAgB+Vl3Av/HqHep3UVjMmxVPCHaIIYQAA+EnVEu6bE/orPIwS7lBFCAMAwA8o4UZVhDAAAHysvIR7f1q+rr2on86jhBsihAEA4FN7UnL19zfWqKDIpXsmxyu+d4zZIyFAEMIAAPCRzbszNf+DTYoIs+uha4arRydKuHEcIQwAAB/4b9IhvfnZNnXr2EZ3TRxECTdOQAgDAKAZGYahpd/t0rJVezXw5A76882jdDSvyOyxEIAIYQAANJNSl0cLlm/Vj14l3K1ahBHCUC1CGAAAzYASbjQUIQwAgCbyLuGePmGgzhhACTfqRggDAKAJdh7M0TxKuNEIhDAAABqJEm40BSEMAIAGooQbzYEQBgBAA7g9Hr2zcru+XkcJN5qGEAYAQD15l3CPP6O7JlPCjSYghAEAUA+UcKO5EcIAAKjDgbR8zV2UpKOFlHCj+RDCAACoxebdmXp+6Ua1CHdQwo1mRQgDAKAG3yUd0hufbVPn6FaaOWUIJdxoVoQwAACqqFrCfcflg9Qygl+ZaF78jQIAwEt1JdxOh93ssRCECGEAABxDCTf8iRAGAIDKSrjnLkzWEUq44SeEMABAyKOEG2YghAEAQhol3DALIQwAEJIo4YbZCGEAgJBTqYS7X6xuvnQAJdzwO0IYACCkUMKNQEEIAwCEjEol3ONO0XnDu5o9EkIYIQwAEBIo4UagIYQBAILe5j2Zmr90o8LDKOFG4CCEAQCCGiXcCFSEMABAUKpUwt2zve6YOJgSbgQU/jYCAIIOJdywAkIYACCo5BeW6rklG/Xr/mxKuBHQCGEAgKCRll2oue8n6UhOoW6dMECjBnQyeySgRoQwAEBQ2HkwR88sTpbHQwk3rIEQBgCwvPIS7qjIcM2cMkSdo1ubPRJQJ0IYAMCyDMPQ5z/v1/tf7VCvLm01YzIl3LAOQhgAwJLcHo/+s3K7vlp3UCP6xeoWSrhhMYQwAIDlFJW49OKHm5VECTcsjBAGALAUSrgRLAhhAADLOJCer7kLKeFGcCCEAQAsgRJuBBtCGAAg4FHCjWBECAMABKyyEu7dWrZqDyXcCDr8TQYABCTvEu5z4jvr2oso4UZwIYQBAAKOdwn3FWN6KeFMSrgRfAhhAICAQgk3QgUhDAAQMCjhRighhAEAAsLabWl66WNKuBE6CGEAAFNRwo1QRQgDAJiGEm6EMkIYAMAUlUq4R3bX5PMo4UZoIYQBAPwuO79Y8xYma19ann477hSdTwk3QhAhDADgV94l3HdPiteQPpRwIzQRwgAAfkMJN3AcIQwA4BffJR/SGyu2qVN0K91LCTdACAMA+FbVEu7bLx+sVi349QPwUwAA8JlSl0cLPt2qHzdTwg1URQgDAPgEJdxA7QhhAIBmV6mEO3GARg2khBuoihAGAGhWOw/l6JlFZSXc9101VP26tzd7JCAgEcIAAM2GEm6g/ghhAIAmMwxDX/y8X+99tUMnd2mruyfFq21rSriB2hDCAABN4vEY+s/K7fpy3QFKuIEGIIQBABqNEm6g8QhhAIBGoYQbaBpCGACgwQ4eK+HOL3RpxqR4DaWEG2gwQhgAoEEo4QaaByEMAFBv3iXcMycPUXQ7SriBxqpXgde8efN0ySWXKCEhQQsWLJAkrVq1SomJiRo3bpzmzJnj0yEBAOYyDENL/7tLC5b/on7do/SHa0YQwIAmqnMl7KefftKPP/6ojz76SC6XS5dcconOPPNMzZo1S2+++aY6d+6s6dOn69tvv9XYsWP9MTMAwI9KXR79+9OtWr05VWfHd9Z1lHADzaLOn6KRI0fqjTfekNPpVEZGhtxut3Jzc9WjRw9169ZNTqdTiYmJWrFihT/mBQD40dGiUj393gat3pyqiWN66YaLTyWAAc2kXj9JYWFheuaZZ5SQkKAzzzxTaWlpio2NrXh/XFycUlNTfTYkAMD/0rILNfuNtdp5KEe3Jg5Q4uiesvEcYECzqfeF+XfffbduueUW3XbbbdqzZ0+lH0TDMBr8gxkdHdmgj4d5YmO588nqOIfWZsb527Y3U4+9tVZut6FHp4/WoN48BUVT8DOI6tQZwnbu3KmSkhL1799fLVu21Lhx47RixQo5HMcrKdLT0xUXF9egT5yRkS+Px2j4xPCr2Ng2Sk/PM3sMNAHn0NrMOH9rt6XrpY83l5VwTxuijm0j+DvUBPwMWpvdbvPZwlGd25EHDhzQn/70J5WUlKikpERffvmlpk6dqt27d2vv3r1yu91atmyZxowZ45MBAQD+YRiGPv9pn+Yv3ahucZH647WnqXN0a7PHAoJWnSthY8eOVXJysi6//HI5HA6NGzdOCQkJ6tChg2bMmKHi4mKNHTtW48eP98e8AAAfqFTCfUqsbkmkhBvwNZthGKbsCbIdaQ0so1sf59Da/HH+ikvcevGjzdqw44guGtlNU87rQwl3M+Jn0Np8uR3JM+YDQAjLzi/WvEXJ2pdKCTfgb4QwAAhR5SXceYWllHADJiCEAUAI2rInU88fK+H+wzUjKOEGTEAIA4AQ831yil5f8Qsl3IDJCGEAECIMw9AH3+3Wx6v2aEDP9rrj8sFq1YJfA4BZ+OkDgBBACTcQeAhhABDkjhaV6rnFG7Vtf7YmjumlS8/sQQckEAAIYQAQxNKzCzV3YZLSswt1S+IAnTmwk9kjATiGEAYAQWrXoVw9syhJbo+h+64aqn7d25s9EgAvhDAACEJrt6Xr5Y83q23rcN175RA6IIEARAgDgCBiGIa++Hm/3vtqh07u0lZ3T4pX29bhZo8FoBqEMAAIEh6Pof98uV1fri0r4b45cYAiKOEGAhYhDACCACXcgPUQwgDA4rxLuK/5zSm6YAQl3IAVEMIAwMIo4QasixAGABZVVsK9SeFOux66Zrh6dmpr9kgAGoAQBgAWVFHC3aGVZk6hhBuwIkIYAFiIYRj68Pvd+ugHSrgBq+MnFwAsoqyE+xet3nxYZw/urOvGU8INWBkhDAAs4GhRqZ5fslG/7MvWxHNO1qWje1LCDVgcIQwAAhwl3EBwIoQBQAD7dV+WZr+xhhJuIAgRwgAgQK3dlq6Xl21R21ZhlHADQYgQBgAB6POf9+u9L7frlO7tdftlAynhBoIQIQwAAoh3CffwU2L1hxtGKi+n0OyxAPgAIQwAAoR3Cfe407vpyvP6qEW4U3lmDwbAJwhhABAAcvKLNZcSbiCkEMIAwGSUcAOhiRAGACbauidTz1HCDYQkQhgAmOSHjSn696eUcAOhihAGAH7mXcLdv0d73TmREm4gFPFTDwB+5HJ7tGA5JdwACGEA4DeUcAPwRggDAGaHTXkAACAASURBVD8oL+FOyyrULZcO0JmDKOEGQh0hDAB8bNehXD2zKEluj6H7p1LCDaAMIQwAfGjdr+l66aPNats6nBJuAJUQwgDAR8pLuHt2bqt7JsdTwg2gEkIYADQzj8fQu19u18pjJdy3JA5QRJjD7LEABBhCGAA0o+pKuO127oAEcCJCGAA0E0q4ATQEIQwAmkFZCXey8gpLNOOKeA3tSwk3gNoRwgCgiSjhBtAYhDAAaALvEu57psQrpl1Ls0cCYBGEMABohBNLuAepVYsws8cCYCGEMABoIO8S7rMGd9L140+lhBtAgxHCAKABvEu4Lz/nZCVSwg2gkQhhAFBPR7ILNYcSbgDNhBAGAPVQXsLtclPCDaB5EMIAoA7eJdwPXj1EXWIo4QbQdIQwAKjFFz/v17vHSrjvnhyvdpRwA2gmhDAAqAYl3AB8jRAGAFVQwg3AHwhhAOAlJ79Y8xYlay8l3AB8jBAGAMccPHJUc99PooQbgF8QwgBAlHAD8D9CGICQRwk3ADMQwgCELEq4AZiJEAYgJFHCDcBshDAAIYcSbgCBgBAGIKRQwg0gUBDCAISM3Sm5mrewrIT7vquG6tQelHADMA8hDEBIWP9rul6khBtAACGEAQh6lHADCESEMABBixJuAIGMEAYgKBWXuPXSx5u1fjsl3AACEyEMQNDxLuG++sK+uvC0bmaPBAAnIIQBCCreJdx3XTFYw/rGmj0SAFSLEAYgaFDCDcBKCGEAgkJ5CXfHDq00kxJuABZACANgaYZh6KMf9ujD73dTwg3AUghhACzL5fbo35/+olWbKOEGYD2EMACWRAk3AKsjhAGwHO8S7psv7a/RgzqbPRIANBghDICl7E7J1bxFyXK5PJRwA7A0QhgAy6hUwj1tGCXcACyNEAbAEr5Ys1/vrtyunp3b6O7JQyjhBmB5hDAAAc3jMfTuV9u1cs0BDesbo1snDKSEG0BQIIQBCFjeJdy/Oa2brjqfEm4AwYMQBiAg5Rwt0TOLkrTnMCXcAIITIQxAwDl45KjmLUxSbgEl3ACCFyEMQEDZujdLzy3ZqDCnXf939XCd3JkSbgDBiRAGIGBUKuGeHK+YKEq4AQQvQhgA01HCDSAUEcIAmKpSCfegTrr+Ykq4AYQGQhgA0xQUleo5SrgBhChCGABTHMku1NxFyUrNLKCEG0BIIoQB8DtKuAGAEAbAz9b/mq4XP96stq0o4QYQ2ghhAPyGEm4AOI4QBsDnKOEGgBMRwgD4FCXcAFA9QhgAn6ko4U7J07QL++o3lHADQAVCGACfOHTkqOaWl3BPooQbAKoihAFodlv3Zun5JRvlpIQbAGpECAPQrMpLuOPat9S9U4ZQwg0ANSCEAWgWhmHo4x/26ANKuAGgXghhAJrM5fbo9U9/0Q+UcANAvRHCADRJpRLus09W4lmUcANAfRDCADSadwn3TQn9ddZgSrgBoL4IYQAapbyEu9Tl0e+vGqr+lHADQIPU66KN5557TgkJCUpISNA///lPSdKqVauUmJiocePGac6cOT4dEkBgWb89Xf94Z53CHHbNunYEAQwAGqHOELZq1Sp9//33Wrp0qT744ANt3rxZy5Yt06xZszR//nwtX75cmzZt0rfffuuPeQGYbOWa/Xpu8UadFNNaf7puhE6KaW32SABgSXWGsNjYWD300EMKDw9XWFiYevfurT179qhHjx7q1q2bnE6nEhMTtWLFCn/MC8AkHo+h/6zcrndWbtfQvjF68OrhahcZYfZYAGBZdYawvn37aujQoZKkPXv26NNPP5XNZlNs7PEKkri4OKWmpvpuSgCmKi516/mlG/XFmv268LSuunPiYEWEOcweCwAsrd4X5m/fvl3Tp0/Xgw8+KIfDoT179lS8zzCMBt+SHh0d2aCPh3liY9uYPQKaqCnnMCuvSI+9vU47DmTrlssHacI5vZtxMtQHP4PWxzlEdeoVwtauXau7775bs2bNUkJCgn766Selp6dXvD89PV1xcXEN+sQZGfnyeIyGTQu/i41to/T0PLPHQBM05RxWlHAfLdFdEwdr2Cmx/H3wM34GrY9zaG12u81nC0d1bkempKTozjvv1JNPPqmEhARJ0pAhQ7R7927t3btXbrdby5Yt05gxY3wyIABzbN2bpf/35lqVuDz6v2uGa9gpsXU/CABQb3WuhL366qsqLi7W448/XvG2qVOn6vHHH9eMGTNUXFyssWPHavz48T4dFID/rNqUogXLKeEGAF+yGYZhyp4g25HWwDK69TXkHHqXcJ/aPUp3XTGYEm6T8TNofZxDa/PldiTPmA9AUuUS7tGDOul3lHADgE8RwgCooKhUzy/dpK17s3TZ2SdrAiXcAOBzhDAgxB3JKdTchZRwA4C/EcKAEFaphPvKIerfs4PZIwFAyCCEASFqw/Yj+tdHm9SmZbgemDaMDkgA8DNCGBCCVq7Zr/98uV09OrbRPZPj6YAEABMQwoAQ4vEYeu+rHfpizX4N6xujWxMHKiKcDkgAMAMhDAgRxaVuvfTRZq3ffkQXntZVU8/vK7udOyABwCyEMCAE5Bwt0TOLkrQnJU/TLuir35zezeyRACDkEcKAILc/NU+z31hTVsJ9xWA6IAEgQBDCgCD2y94sPf/BJjnsNv3fNcN1cue2Zo8EADiGEAYEqfIS7i6xrTVj4mBKuAEgwBDCgCBTtYT7L7ecqcKjxWaPBQCoghAGBBGX26PXV/yiHzYeL+GObBVOCAOAAEQIA4IEJdwAYC2EMCAIUMINANZDCAMsbndKrp5ZlKwSSrgBwFIIYYCFeZdw308JNwBYCiEMsKgv1x7QOyt/pYQbACyKEAZYDCXcABAcCGGAhVQq4R7RVVMvoIQbAKyKEAZYRFkJd7L2pORSwg0AQYAQBljAoSNHNXdhEiXcABBECGFAgPtlb5aeW7JRTgcl3AAQTAhhQABbvemwXlu+VXHtW2rmlCGKpYQbAIIGIQwIQIZh6ONVe/TBd2Ul3HdeMVitW4SZPRYAoBkRwoAA413CfebATrrhklPldNjNHgsA0MwIYUAA8S7hnnBWT1129smUcANAkCKEAQHiSE6h5i1M1mFKuAEgJBDCgABACTcAhB5CGGAySrgBIDQRwgATUcINAKGLEAaYwOMx9P7XO/T5z/s1tE+Mpk+ghBsAQg0hDPCz4lK3Xv54i9b9mk4JNwCEMEIY4EeUcAMAyhHCAD9JyTiqOe+XlXDfecVgDaeEGwBCGiEM8INt+7L07GJKuAEAxxHCAB+jhBsAUB1CGOAjlHADAGpDCAN8gBJuAEBdCGFAM6OEGwBQH4QwoBlRwg0AqC9CGNBM9hzO1byFlHADAOqHEAY0g+Ml3GG6f+pwnRQbafZIAIAARwgDmogSbgBAYxDCgEaihBsA0BSEMKARvEu4LxjRVdMo4QYANBAhDGig3KMlmkcJNwCgiQhhQANQwg0AaC6EMKCevEu4H7x6uHp1oYQbANB4hDCgHlZvPqzXPqGEGwDQfAhhQC0Mw9CyVXu0lBJuAEAzI4QBNXC5PXpjxTZ9vzFFZw7sqN9d3F9hTkq4AQDNgxAGVKOgyKX5H2zUlj2UcAMAfIMQBlRBCTcAwB8IYYAX7xLue68cogGUcAMAfIQQBhyzYccR/etDSrgBAP5BCAN0vIS7e8c2mkkJNwDADwhhCGkew9D7X1HCDQDwP0IYQlZxqVuvfLxFaynhBgCYgBCGkJR7tETPLE7W7kO5mnpBX42jhBsA4GeEMIQc7xLuOyYO1oh+lHADAPyPEIaQsm1flp5bslEOOyXcAABzEcIQMlZvPqwFy7cqNooSbgCA+QhhCHqUcAMAAhEhDEHN5fbojc+26ftkSrgBAIGFEIagRQk3ACCQEcIQlDJyijR3YZIOZxboxkv66+x4SrgBAIGFEIagQwk3AMAKCGEIKpRwAwCsghCGoPHVugN6+4uyEu57JscrihJuAEAAI4TB8ijhBgBYESEMllaphHt4V027kBJuAIA1EMJgWVVLuH9zWleeggIAYBmEMFgSJdwAAKsjhMFyyku47XabHrh6mHp3aWf2SAAANBghDJby4+bDeu1YCfc9U4YojhJuAIBFEcJgCd4l3P26RemuSZRwAwCsjRCGgEcJNwAgGBHCENC8S7gTR/fU5edQwg0ACA6EMASsjJwizV2UpMMZlHADAIIPIQwBae/hPM1dmKQSl5sSbgBAUCKEIeBs2HFEL364WZEtnbpv6gh1pYQbABCECGEIKBUl3HFtdM8USrgBAMGLEIaA4DEMLfx6hz77ab+G9I7W9MsGqkU4fz0BAMGL33IwXUmpWy8v26K12yjhBgCEDkIYTJV7tETPLk7WLkq4AQAhhhAG05SXcOccLdEdEwdpRL84s0cCAMBvCGEwhXcJ94OUcAMAQhAhDH5XXsId066lZl5JCTcAIDQRwuA3hmFo2eq9WvrfXerXLUp3XjFYkS0p4QYAhCZCGPzCu4R71MCOuoESbgBAiCOEwecKilx64YON2kwJNwAAFQhh8CnvEu4bLjlV58R3MXskAAACAiEMPrP3cJ7mLkpSSalbM68cooGUcAMAUIEQBp9I2nFE/yov4f4tJdwAAFRFCEOzo4QbAIC6EcLQbCjhBgCg/vgNiWbhXcJ9/vCTdPWFp1DCDQBALQhhaLJKJdzn99FvTu/GU1AAAFCHej1bZn5+vi699FIdOHBAkrRq1SolJiZq3LhxmjNnjk8HRGBLyTiq2W+u0b60fN0xcZDGjexOAAMAoB7qDGFJSUmaNm2a9uzZI0kqKirSrFmzNH/+fC1fvlybNm3St99+6+s5EYC27cvS/3tzrYpK3Hrw6mEa0S/O7JEAALCMOkPY+++/r7/85S+Kiyv7BZucnKwePXqoW7ducjqdSkxM1IoVK3w+KALLj5sP66n3NqhNq3D98brT1LtLO7NHAgDAUuq8Jmz27NmVXk9LS1NsbGzF63FxcUpNTW3+yRCQKOEGAKB5NPjCfI/HU+maH8MwGnUNUHQ0T95pFbGxbSSVlXDPX5SkL37ap3OHd9XdVw1VmNNh8nSoj/JzCGvi/Fkf5xDVaXAI69Spk9LT0yteT09Pr9iqbIiMjHx5PEaDHwf/io1to/T0vGpLuLOzCsweD/VQfg5hTZw/6+McWpvdbvPZwlG97o70NmTIEO3evVt79+6V2+3WsmXLNGbMGF/MhgCRkVOkx95eq1/2ZeuGS07VxDG9uAMSAIAmavBKWEREhB5//HHNmDFDxcXFGjt2rMaPH++L2RAAdhzI1t/fXEMJNwAAzcxmGIYpe4JsRwa+pB1H9OJHm9WqhVMzpwyhhNui2AqxNs6f9XEOrc2X25E8Yz6qVV7C3eukdrrz8kGUcAMA0MwIYajEYxha9PVOrfhpn+J7R+tPN41Sfm6h2WMBABB0CGGoULWEe9qFfdUywql8swcDACAIEcIgScotKNGziyjhBgDAXwhhUErGUc1dmKTs/BLdMXEQHZAAAPgBISzE/bo/W88uTpbdbtOD04ap90l0QAIA4A+EsBD24+bDem35VsW0a6mZVw5RXFRLs0cCACBkEMJCkGEY+mT1Xi357y6d0i1Kd1HCDQCA3xHCQozL7dGbn23Td8kpGjWgo264pL/CnA1urwIAAE1ECAsh3iXcl47uqYnnnMwdkAAAmIQQFiIyc4s0d2GSUjIKdMPFp+qcIV3MHgkAgJBGCAsBew/nae6iJEq4AQAIIISwIJe884he+GCzWrd06g+/HUEJNwAAAYIQFsS+XndAb33xq7rHtdHdk+PVvg0l3AAABApCWBCqWsJ922UD1SKcUw0AQCDhN3OQKSl165VlW7RmW7rOG36Srr6wrxx2noICAIBAQwgLIt4l3Fed30fjKOEGACBgEcKCxOHMAs15fwMl3AAAWAQhLAhQwg0AgPUQwizuxy2H9donlHADAGA1hDCLooQbAABrI4RZECXcAABYHyHMYgqLXZq/lBJuAACsjhBmIZRwAwAQPAhhFlGphHvKEA08mRJuAACsjBBmAZVKuK8Zoa5xlHADAGB1hLAARwk3AADBiRAWoCjhBgAguPFbPQBRwg0AQPAjhAWY3IISPbs4WbsOUsINAEAwI4QFEO8S7tsvH6TTTqWEGwCAYEUICxCUcAMAEFoIYQGgvIQ7ul1L3TslXnHtW5k9EgAA8DFCmIkMw9DyH/dq8be7dErXdrprUjwl3AAAhAhCmElcbo/e+nyb/ptECTcAAKGIEGaCwmKX5n+wSZt3Z1LCDQBAiCKE+Vl5CfehIwX63cWnagwl3AAAhCRCmB/tS83T3IVJKi51694rKeEGACCUEcL8JHnnEb3w4Wa1bkEJNwAAIIT5xdfrD+qtz7epW1yk7pk8hBJuAABACPMlj2Fo0Tc7teJ/lHADAIDKSAQ+UlLq1iufbNWaX9J03rCTdPVvKOEGAADHEcJ8oLyEe+fBXF15Xh9dNJISbgAAUBkhrJkdzizQ3PeTlJVfrDso4QYAADUghDWj8hJum40SbgAAUDtCWDP535ZUvfrJFkq4AQBAvRDCmogSbgAA0BiEsCbwLuE+Y0BH3UgJNwAAqCdCWCMVl7j1/NKN2rQ7U5eO7qHLz+klO3dAAgCAeiKENUJhsUtzFyZpx8EcSrgBAECjEMIaqKCoVE+/n6S9h/M0fcJAjezf0eyRAACABRHCGiC/sFRPvrteB9OP6o7LB2nYKbFmjwQAACyKEFZPOUdL9NS763U4s1AzJsUrvne02SMBAAALI4TVQ1ZesZ58d70ycoo0c0q8BvTsYPZIAADA4ghhdcjIKdIT/1mvnIIS/f6qoTqlW5TZIwEAgCBACKtFWnahnnhnvQqKXbr/qqHUEAEAgGZDCKvB4cwCPfGf9SopdeuBaUPVs1Nbs0cCAABBhBBWjYPp+Xri3Q0yDEMPXj1c3eIizR4JAAAEGUJYFftS8/TkuxvkcNj04LTh6hLT2uyRAABAECKEedmdkqun39ugiHCHHpg6TB07tDJ7JAAAEKQIYcfsOJCjOQs3qHWLMD04bZhiolqaPRIAAAhihDBJ2/Zlae7CZEVFhuuBacPUoW0Ls0cCAABBLuRD2KbdGXpu8UZFt2uhB6YNU1RkhNkjAQCAEBDSIWzDjiOav3SjOnVorfunDlXb1uFmjwQAAEJEyIawtdvS9K8PN6trXKTuu2qoIluGmT0SAAAIISEZwv63JVUvf7xFJ3dpo3unDFWrFiH5bQAAACYKufTxw8YUvbZ8q/p2jdI9k+PVMiLkvgUAACAAhFQC+WbDQb2xYpsG9GyvGZPiFRHmMHskAAAQokImhH2xZr/+s3K74ntH686JgxTmJIABAADzhEQI+/R/e7Xw650afkqsbrtsoJwOu9kjAQCAEBfUIcwwDH28ao8++G63RvaP082XDiCAAQCAgBC0IcwwDC357y59snqvRg/qpBsv6S+73Wb2WAAAAJKCNIQZhqH3vtqhz3/erzFDuui68f1ktxHAAABA4Ai6EOYxDL39xa/6et1BXTCiq66+sK9sBDAAABBggiqEeTyGXl/xi75LTtH4kd015bzeBDAAABCQgiaEuT0evfbJVq3enKrE0T11+TknE8AAAEDACooQ5nJ79NLHW7TmlzRNHNNLiaN7mj0SAABArSwfwkpdHv3rw01av/2Irjyvj8af0d3skQAAAOpk6RBWUurWc0s3atOuTF3zm1N0wYiuZo8EAABQL5YNYcUlbj2zOFm/7M3S7y4+VWOGdDF7JAAAgHqzZAgrLHZp7sIk7TiYo5su7a/RgzqbPRIAAECDWC6EFRSV6un3k7QnJU/TJwzUyP4dzR4JAACgwSwVwvILS/XUuxt0ID1fd0wcpOGnxJo9EgAAQKNYJoTlHC3RU++u1+HMQs2YNFjxvWPMHgkAAKDRLBHCsvKK9eS765WRU6SZU+I1oGcHs0cCAAAWYhiG3IZbLo9bHsMtl+GW2+OW+9hLl1H+Z0+lt4U7HRodPcwnMwV8CMvIKdIT/1mvnIIS3XvlEPXr3t7skQAACEmGYchjeCrCjNuoEmIqvc0jt+GS2+ORy3DJbXjk9pS9dHmFH++XVY9R08eVHaP8mMcf66lmjvJjeAxPo77m2FYdNLpPCIawtOxCPfHOehUUu3T/VUPV+6R2Zo8EAECTeaoJIlVDTc0BxzuA1BBuTljdqelzVXOMY8GppmP4kk02OewOOW0OOWwOOeyVXzrtDjlsdjlsTjnsdoXbw2R3Rshpc5a9vdLHeT3O5pC96vtsDjnsdq/HVn7pPPayZVgLn329ARvCDmcW6In/rFdJqVsPTBuqnp3amj0SACCAlK3IVA4gnnqt0FRZNal4/4nhxmW45PF46l6hOSHgHD+Gx+aRy+XyWg1yy5Dh0+9NeeioGmaqhpOyPzsVcSyQOGyVA8jxYHL842sKOM5Kxz3xGBWvHwtQ5S+9Z7Tb7D79vjSG3e67HuqADGEH0/P1xLsbZBiGHrx6uLrFRZo9EgAEpfLrZCqHmdpXRmoLONVvH3kfo+r21IlhproVouPX8Bw7huFp9PZSfdlt9irBopoQ4xVEwpxhFW9z2hyyH/v41q1aqLTYUyWIVH8M77dVF6DsthNXg5zeocZml91ml83mu+CA5hNwIWxfap6efHeDHA6bHpw2XF1iWps9EgDUyvs6mZpChNs4HkDcXqGm5ouDy15GpDmVl19Qw+rO8UByQoCqdQvKazXID0GmcmDwDhInrtCE28PV0umoeyXHe3Wm2m2mqiszNa8GVfe45gwysbFtlJ6e1yzHQnAJqBC2OyVXT7+3QeFhDj04bZg6dmhl9kgA/MhzbGWkahCpe5up8vUrJ24zVROMar2I1+s6mXoew5fKr5Nx2MpWPezVXMfiHSzC7E61sEecuLpS7XZUNdtFtW4z1RJiTgg89oDcXgICRcCEsB0HcjRn4Qa1bhGmB6YNU2xUS7NHAizp+HUyZUEhrMhQVlFOLSs0lbeNqt2CquUupfpeCFyfO6J8fZ1MzVs/1V+cG+4dbrxXSaq7cNh7m6jGa2GqX6GpvM1U+X2d46KUkXHUp98XAOYIiBC2bV+W5i5MVlRkuB6YNkwd2vruTgSgPsq3l+q63qVeF+d6XQPjqeZamIat0FRZpanmjij/bC/VsIVzQjixq4WzReVg4n0xbnXbR7XdEVUeaqpc2Hv8YuKat5msep2M3c5KEoKTYRiS4ZY8Zf8ZnuN/ru5t1b/fVfl1wyN5va3qY2p+3eX1+MrHd7ZqJ13zJ598D0wPYZt3Z+rZxcmKbtdCD0wbpqjICLNHQjM5fp1MlYtxva5jOWGbqcrKSM3hxPsYJwaS2p+I78TVnxNu9/bD9lJNWzjVrZKUXyfjHW5ODCLVbx9FtW2lwqOu6i/wrWM1qHwFx3uVxopBBggGhudYwDgWFE4IEBWvl31cjYHDqC6QuKo8vpqAU93xaglRdR1fPv53thKbTbI7JJtDsjtks5e9LPvP6fW6/dj7nWUvneFSyzY+G8vUELZhxxHNX7pJnTq00v1Th6pt63AzxwlYHu9VE+8AUmsQqeXi3FrDSeVjOMPtKigqOhagql/xqXkLyuPT7SWbbNVfnFvNHUZ227HrZGwRNV6MW/VamNruUqrvhcDlYcfM27C5KBihqHGrLN6hw1N7CDFqCh0ur8eXvT01zKbiwqIqj6/j+NXN7OPt+kpsx0NK5cBS/dtsdodszohjQcdeS8gpCzrlIae24534ulOy2Su/Xt/HN+Hf3KB8iopNuzL0/JKN6hoXqfuuGqrIlmFmjdIgpR6XMgozlV54REcKM1XiLqnnhcBVtqMq3dVU/fPYlB/DL9fJ1HBxboQzTPLYKm0zhTvCqr8Y1+tC36rPB1Pf56w5MRhV/ThuwwZ8oWKVpdIqx/G3NXhbp+LxVVdBGrBNVG2IKj++97bRiceXj7flKylfZSlfPTkWRMrfVhIeJo/HdkJAsDnDawgMNQWK8uPbvd5fNeDU9vjyj7dXeX/lQCT+bfUb00LYW5//qp6d2+jeKUPVqoXpu6KVlAWtDKUXZiit4IjSCzOUXnBE6YVHlFmUXW0oKr9OptaLc722jSIc4XLaW9YQYqqGkyp3RNVxDYy9mjBT0/U3dQUZVlGAExmGUceKSNUVlcqrLPVelTE8ymzhUHF+YQ3HqCMk1bpV5Dq+LeTvVZaKX/a1rHh4BQXvVRabV1CoMYTU+fqJxz9hVaWWlaCGrrLw7yhqYlr66dG5ja4b108tI8wZodRdqoyizIqQlVZ4REcKyl5mVQlaLZ0tFdcyRie366GRnUYorlWMYltGK6ZldMVFx9yGDVSvbFvIqDWEGEZNgcRVaVun2sBR47ZQXQHEe1uoluNXF5r8uMpSYrMfCwzOEwNATSHG4ZQtLKKWAOGsZtuomuPX+rqzUdtOrLIAx5kWwm5KGCCnD/dZpbKgdaQoU+kFR5RW6L2ilXFC0GrlbKnYVjHq1a6H4jqNUGyrGMW2jFFsq2hFhvGEsfCfSqss5dssXnfsNPQ6ktyDTpXkFJy4ylJpxaSxdxQdm7GaO4oqX8viR1VXPGx1rHjY7SeustSwItKYrR2bV4Cq/9ZT+cfbFRfXjlUUIEiZFsLCnXZ5PE1f/i4PWmnHtgvLQ1ZawRFlF+dUClqtna0U0ypavdv1VGyn6IqgFdcqRq3DeGJYKypbZakaAKoLDPW5o6iOwNHUO4qqfbzvV1mKanunrWqAqCFweN9RVGmV5XhYqDnkOGq4I6mOVZaKEFT/bSdWWQBYSZNC2Mcff6wXXnhBLpdL119/va655prmmquSUndp2SpWYUZF0Eo7tqpVXdCKbRWjPlEnHwtZ0ce2DwlaUg2rLLVsCxUVRciVmef18bWveFR7cW+lgNGIVRbvcBJoqyxVrwupJoDYwsLqDiANvG263teq2J2Kwp8LCQAACM9JREFUjmmrjKzCaldZmnLHEACgaRodwlJTUzVnzhwtWbJE4eHhmjp1qs444wz16dOnUccrcZfqSHnIqnJB/AlBK6yVYlvGqE9UL8W2ilbcsW1DXwStukJFXVs7TXuiuEbcUVTrtTXusmtzGqCgKd+8+q6ylK962Ow1rLI0ZVvHe1XFXuX1+j6+/PoXmyVXWZzt2shewnYWAASaRoewVatWadSoUYqKipIkXXTRRVqxYoXuuuuuej3+u+2faU/uQaUXZyu9JFc5rqOV7s1pbQ9XjLO1ejlaKzYyRjH2Foq1tVC0LVytZC8LFEfdMvIyJE9aRegoPCGk1Of6mVpWbvx2x1D57ct1rHhUWTGp1yqL94pKDVs7NT0PTLv2bZSTV+z18d7X11QzY0WQYZUFAIDaNDqEpaWlKTY2tuL1uLg4JScn1/vxX6SvU0HeEUWXunVyiUvRpW7FlLoVfey/VrVcL1ZS/oc6byGuz9ZQXXcENfX9tYeo44EoMANLq9g2OspFwQAANLtGhzCPx1Npa8YwjAZt1fy13wS1sjklx7HngHGUhxmnbI5jQcXhrAhTNofXRcEV20fW2xqyothY31U2wD84h9bG+bM+ziGq0+gQ1qlTJ61Zs6bi9fT0dMXFxdX78bbogTpa292RhiRX1Td65LUOBj/gSQatj3NobZw/6+McWpvdblN0dKRvjt3YB44ePVqrV69WZmamCgsL9fnnn2vMmDHNORsAAEDQavRKWMeOHXXvvffquuuuU2lpqSZPnqz4+PjmnA0AACBoNel5whITE5WYmNhcswAAAISMwLwlDwAAIMgRwgAAAExACAMAADABIQwAAMAEhDAAAAATEMIAAABMQAgDAAAwASEMAADABIQwAAAAExDCAAAATEAIAwAAMAEhDAAAwASEMAAAABMQwgAAAExACAMAADABIQwAAMAEhDAAAAATEMIAAABM4DTrE9vtNrM+NRqIc2V9nENr4/xZH+fQunx57myGYRg+OzoAAACqxXYkAACACQhhAAAAJiCEAQAAmIAQBgAAYAJCGAAAgAkIYQAAACYghAEAAJiAEAYAAGACQhgAAIAJCGEh6LnnnlNCQoISEhL0z3/+U5L+f/v2F9JkG8Zx/DvdDGFRWFuTtIIQAok6UGIUsyTSNUcwggqjAw/ECJROapYkgWWGtJNOiqIIOuiPpjRMiGykTTKs8CShg9S1hk3LTGe6zfs9iFfy/eN78FrzYdfnZHuuPeO5b37s5tr9bPj9fpxOJ3v27MHj8cyd+/btW1wuF0VFRZw+fZpYLAbAx48fKS0tpbi4mKNHjzI5OZmQuSS7hoYG3G43IBlqTUdHBy6XC7vdTl1dHSAZaklra+vcOtrQ0ABIfloxMTFBSUkJHz58ABYvt/HxccrLy7Hb7ZSWlhIOh/97MEoklefPn6sDBw6o6elpNTMzo44cOaIePnyoCgoK1NDQkIpGo6qsrEz5fD6llFIOh0O9fv1aKaVUdXW1un37tlJKqfLycuX1epVSSl2+fFldvHgxMRNKYn6/X23btk2dPHlSTU1NSYYaMjQ0pHbs2KFCoZCamZlRhw4dUj6fTzLUiEgkovLz89Xo6KiKRqNq//796smTJ5KfBrx580aVlJSo3NxcFQgEFnXtPHv2rLpy5YpSSqkHDx6oqqqq/xyP7IQlGZPJhNvtJi0tDYPBwMaNGxkYGGD9+vVkZ2ej1+txOp20t7cTDAb5/v07W7duBcDlctHe3k40GuXly5cUFRXNq4vfZ2xsDI/HQ0VFBQB9fX2SoYY8fvyYvXv3YrFYMBgMeDwe0tPTJUONiMfjzM7OMjU1RSwWIxaLYTQaJT8NuHv3LrW1tZjNZmBx106fz4fT6QSgpKSEZ8+eEY1GFxyP/ldNVCxNOTk5c88HBgZ49OgRhw8fxmQyzdXNZjPDw8N8+vRpXt1kMjE8PMyXL18wGo3o9fp5dfH7nDlzhuPHjxMKhQD+lpVkuLQNDg5iMBioqKggFAqxc+dOcnJyJEONMBqNVFVVYbfbSU9PJz8/Xz6DGnHu3Ll5x4uZ28/v0ev1GI1GPn/+zJo1a/51PLITlqTevXtHWVkZJ06cIDs7G51ON/eaUgqdTsfs7Ow/1v98/Nlfj8Wvc+/ePTIzM7FarXO1f8tKMlya4vE43d3dnD9/njt37tDX10cgEJAMNaK/v5+mpiaePn1KZ2cnKSkpDAwMSH4a9CvXTqUUKSkLt1myE5aEent7qays5NSpUzgcDnp6eub9gDAcDmM2m7FYLPPqIyMjmM1mMjIy+PbtG/F4nNTU1Lnzxe/R1tZGOBxm3759fP36lUgkQjAYJDU1de4cyXBpW716NVarlYyMDAB2795Ne3u7ZKgRXV1dWK1WVq1aBfy4JXX9+nXJT4P+ms//yc1sNjMyMoLFYiEWizE5OcnKlSsXvL7shCWZUCjEsWPHaGxsxOFwALBlyxbev3/P4OAg8Xgcr9eLzWZj7dq1LFu2jN7eXuDHv4FsNhsGg4G8vDza2toAaGlpwWazJWxOyebGjRt4vV5aW1uprKyksLCQa9euSYYasmvXLrq6uhgfHycej9PZ2UlxcbFkqBGbNm3C7/cTiURQStHR0SHrqEYtZm4FBQW0tLQAP74s5+XlYTAYFry+TimlfuH8xBJTV1dHU1MT69atm6sdPHiQDRs2UF9fz/T0NAUFBVRXV6PT6ejv76empoaJiQlyc3Opr68nLS2NYDCI2+1mdHSUzMxMLl26xIoVKxI4s+TU3NxMT08PFy5coLu7WzLUkPv373Pz5k2i0Sjbt2+npqaGFy9eSIYacfXqVZqbmzEYDGzevJna2lpevXol+WlEYWEht27dIisra9HWzrGxMdxuN4FAgOXLl9PY2EhWVtaC45AmTAghhBAiAeR2pBBCCCFEAkgTJoQQQgiRANKECSGEEEIkgDRhQgghhBAJIE2YEEIIIUQCSBMmhBBCCJEA0oQJIYQQQiSANGFCCCGEEAnwB4QqzgQMtTTiAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 720x720 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "times.T.plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "geo_dev",
   "language": "python",
   "name": "geo_dev"
  },
  "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.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}