Results_on_synthetic_linear_Erdos_Renyi_graphs.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "Oex92m5gPgTK"
   },
   "source": [
    "Code for generating the results for synthetic linear Erdos-Renyi causal graphs (Fig. 16)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "8Xqc_6q-WNW1"
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import networkx as nx\n",
    "import numpy\n",
    "import pandas\n",
    "import networkx\n",
    "from itertools import combinations, permutations\n",
    "import logging\n",
    "from causallearn.utils.cit import CIT\n",
    "from sklearn import linear_model\n",
    "import collections\n",
    "import matplotlib.pyplot as plt\n",
    "import pickle as pkl\n",
    "from datetime import datetime\n",
    "import copy\n",
    "\n",
    "from causal_discovery.utils import *\n",
    "from causal_discovery.pc_alg import PCAlgorithm\n",
    "from causal_discovery.mb_by_mb import MBbyMBAlgorithm\n",
    "from causal_discovery.sd_alg import SequentialDiscoveryAlgorithm\n",
    "from causal_discovery.ldecc import LDECCAlgorithm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "pnj3tN45cbZq"
   },
   "outputs": [],
   "source": [
    "def get_erdos_renyi_graph(n, p, seed=None):\n",
    "  nx = networkx\n",
    "  np = numpy\n",
    "\n",
    "  def relabel_graph(g, labels):\n",
    "      return nx.relabel_nodes(g, {i: labels[i] for i in range(len(labels))})\n",
    "\n",
    "  graph = nx.erdos_renyi_graph(n, p, directed=False, seed=seed)\n",
    "  max_cc = list(max(nx.connected_components(graph), key=len))\n",
    "  graph = nx.relabel_nodes(graph.subgraph(max_cc).copy(), \n",
    "                          {val: idx for idx, val in enumerate(max_cc)})\n",
    "  graph = graph.to_directed()\n",
    "\n",
    "  # assume label number is topological order and convert to DAG.\n",
    "  D = len(graph.nodes())\n",
    "  for i in range(D):\n",
    "      for j in range(D):\n",
    "        if i > j and graph.has_edge(i, j):\n",
    "          graph.remove_edge(i, j)\n",
    "\n",
    "  assert(nx.is_directed_acyclic_graph(graph))\n",
    "  return graph\n",
    "\n",
    "def get_valid_graph(n, p, min_size=1, seed=None):\n",
    "  nx = networkx\n",
    "  np = numpy\n",
    "\n",
    "  def get_successors(g, node, depth):\n",
    "    if depth < 1:\n",
    "      return set()\n",
    "\n",
    "    depth_one_succ = set(g.successors(node))\n",
    "    succ = set(depth_one_succ)\n",
    "\n",
    "    for n in depth_one_succ:\n",
    "      succ |= get_successors(g, n, depth - 1)\n",
    "\n",
    "    return succ\n",
    "  \n",
    "  def get_tmt_and_outcome_idx(g):\n",
    "    potential_tmts_to_outcomes = {}\n",
    "\n",
    "    for tmt in list(nx.topological_sort(g)):\n",
    "      if len(set(g.predecessors(tmt))) == 0:\n",
    "        continue\n",
    "      \n",
    "      succ = get_successors(g, tmt, depth=6)\n",
    "      if len(succ) == 0:\n",
    "        continue\n",
    "      \n",
    "      potential_tmts_to_outcomes[tmt] = succ\n",
    "    \n",
    "    if len(potential_tmts_to_outcomes.keys()) == 0:\n",
    "      return None\n",
    "    \n",
    "    tmt = np.random.choice(list(potential_tmts_to_outcomes.keys()))\n",
    "    outcome = np.random.choice(list(potential_tmts_to_outcomes[tmt]))\n",
    "    return tmt, outcome\n",
    "  \n",
    "  while True:\n",
    "    graph = get_erdos_renyi_graph(n, p, seed)\n",
    "    if seed is not None:\n",
    "      seed += 1\n",
    "\n",
    "    if len(graph.nodes()) < min_size:\n",
    "      continue\n",
    "\n",
    "    tmt_outcome = get_tmt_and_outcome_idx(graph)\n",
    "    if tmt_outcome is not None:\n",
    "      return graph, tmt_outcome[0], tmt_outcome[1], seed\n",
    "\n",
    "\n",
    "def get_params(n, p, min_size, seed=None, get_next_seed=False):\n",
    "  nx = networkx\n",
    "  np = numpy\n",
    "\n",
    "  def graph_to_adj(g):\n",
    "    D = len(g.nodes())\n",
    "    adj = np.zeros(shape=(D, D))\n",
    "\n",
    "    for i in range(D):\n",
    "      for j in range(D):\n",
    "        if g.has_edge(i, j):\n",
    "          adj[i, j] = 1\n",
    "    \n",
    "    return adj.T\n",
    "  \n",
    "  graph, tmt, outcome, seed = get_valid_graph(n, p, min_size, seed)\n",
    "  A = graph_to_adj(graph)\n",
    "  A *= np.random.uniform(low=.25, high=1, size=A.shape) * (2 * np.random.binomial(1, 0.5) - 1)\n",
    "  noise_vars = np.random.uniform(low=0.1, high=0.2, size=(A.shape[-1]))\n",
    "\n",
    "  B = np.linalg.inv(np.eye(A.shape[0]) - A)\n",
    "  noise_cov = np.diag(noise_vars)\n",
    "  cov_matrix = B @ noise_cov @ B.T\n",
    "\n",
    "  labels = []\n",
    "  other_count = 1\n",
    "  for i in list(graph.nodes()):\n",
    "    if i == tmt:\n",
    "      labels.append(\"X\")\n",
    "    elif i == outcome:\n",
    "      labels.append(\"Y\")\n",
    "    else:\n",
    "      labels.append(\"V%d\" % other_count)\n",
    "      other_count += 1\n",
    "\n",
    "  if get_next_seed:\n",
    "    return A, noise_vars, labels, cov_matrix, seed\n",
    "\n",
    "  return A, noise_vars, labels, cov_matrix"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "u682hu7SP3ji"
   },
   "source": [
    "Below we demonstrate the various (local) causal discovery algorithms on a randomly generated graph."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "0qSr_-irdHwp"
   },
   "outputs": [],
   "source": [
    "A, noise_vars, labels, cov_matrix = get_params(n=30, p=0.08, min_size=15)\n",
    "graph_true = adj_matrix_to_nx_graph(A, labels)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 516
    },
    "id": "NhhHfF9PdQ1w",
    "outputId": "cdfc50d8-c135-4097-a0df-636cff0ffc57"
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAApQAAAHzCAYAAACe1o1DAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nOzdd3xT5fcH8E/uTaG0pZQhG1mCVLaUL1OW7CmiIAiITFGWA7cCPxRUXCAgogXZoAgIQhkqS4ZA2VsqgqyWVaDQ0qa5vz+OobRN0qRtcjM+79crL0hyb+7TUpqT8zzPOQZN0zQQEREREWWTovcAiIiIiMi7MaAkIiIiohxhQElEREREOcKAkoiIiIhyhAElEREREeUIA0oiIiIiyhEGlERERESUIwwoiYiIiChHGFASERERUY4woCQiIiKiHGFASUREREQ5woCSiIiIiHKEASURERER5QgDSiIiIiLKEQaURERERJQjDCiJiIiIKEcYUBIRERFRjjCgJCIiIqIcYUBJRERERDnCgJKIiIiIcoQBJRERERHlCANKIiIiIsoRBpRERERElCMMKImIiIgoRxhQEhEREVGOMKAkIiIiohxhQEkeQdP0HgERERFll1HvAZD/SEkBDh0C9uwBoqOBP/8EYmKApCTAZAICAoDAQOCRR4B69YA6dYCICKBKFUDhRx8iIiKPZdA05obItf79F5g5E/j6a+DqVXnMaJQg0paAAAlAAaBsWWD4cKBfP6BwYZcPl4iIiJzEgJJcZtcu4MMPgVWrJMOYmpr911IUCUJ79QLefhuoVCn3xklEREQ5w4CScl1iIvD++8BnnwGqaj8T6SyjUYLLjz4CRoyQ1yciIiJ9MaCkXLVjB9CnD3D6NGA2u/Za9esDc+YAlSu79jpERERkH7c6UK759lugUSPgn39cH0wCsrmnZk1g7VrXX4uIiIhsY0BJueKzz4DBg6X8T07WSjrDZAKSk4FOnYBly9xzTSIiIsqMU96UY19/Dbz4on7XNxhkLeWqVUDbtvqNg4iIyF8xoKQcWbsWaN9e/8LkiiKlhvbulTqWRERE5D4MKCnb4uOBhx8Grlxxz5rJrBiNQI0aUjDdyJL9REREbsM1lJRtL78shco9IZgEZE3lvn3Ap5/qPRIiIiL/wgwlZUtUlEx1eyKjEdi/H6haVe+REBER+QdmKMlpmga89ppn99ceM0bvERAREfkPZijJadu2AY0b6z0K+1QVOHsWKFlS75EQERH5Pg/OMZGnmjbNOza9fPut3iMgIiLyD8xQklPi4oBSpXK3P7erFC0KnDsn5YSIiIjIdZihJKesXm09mFy2DLhzByhQwPa58+dLZ5uiRaXf96JFwPHj0lnn9Gnb59WtK/Uub94Ebt0Cfv8daNgw67HGxQG7dmV9HBEREeUMA0pyyp491jN+kZFAvnxAr17WzwsNBbp2BX75RQK9Pn1kF/auXUBMjO3rRUQAW7bIa/fpI7fAQOC334D69e2P1WAAoqMd/9qIiIgoezjlTU6JiLAepCmKbIK5eFEyihkNGQLMmAF07ChZToMhrbvOqlVAtWpA+fKZz4uKAmrVAipUABIT5bGQEODvv4GTJ+1vDjIaJcCdM8f5r5OIiIgcxwwlOcxkAg4dsv6c2SyBW0SEBIcZPf88cOGCBIiA460aGzUCNm1KCyYBICFBspaNGgHFi9sf786djl2HiIiIso8BJTns779lDaQts2ZJYNm/f/rHw8OBevUk4HS2q06ePMDdu5kftzxWvbr980+dkjWaRERE5DoMKMlhCQn2n4+Jkcxh797pywpZAsxZs5y/5tGjslbSYEh7TFUlQAWAwoXtn282A0lJzl+XiIiIHMeAkhxmLVOYUWQk8MADQOfOcl9VJcDcskWyhc766ivg4YeBqVOlSHnp0rIWs2xZed6RjCcDSiIiItdiQEkOc6Se49KlQHy8rJkEpN938eISaGbH7NnAG2/I7u7z54F//wUeeQT49FN5/vz5rF8jT57sXZuIiIgcw4CSHBYYmPUxSUlSX7JtWwkk+/eX+pE//pj9637yCVCkiGz2KVtWNuMULChT8I6UBXJk3ERERJR9DCjJYWXLpl/LaEtkpKyhHD1aMpSLF6ffpZ0dycnAkSNSmqhMGaBHD2mtmNV0dtGi7JRDRETkal7QkZk8Rf78Ug/SXiFyQLKGBw4Ao0ZJfUpr093h4TJ1DUgmMygI6NZN7h89Chw7Jn+vWlUe37NH1nDWrAm8+Sbw11/Ae+/ZH4eiZF38nIiIiHKOGUpySv36stEmK5GREtAdOWK9/WH37rLeculSqV1ZtGja/e7d045LTgZatADmzpWC6C+8IJtymjUDbt+2PwaDQV6biIiIXIudcsgpX3wBvPaa8/Uk9bJ6tUy7ExERkeswQ0lOefxx7wkm8+QBGjTQexRERES+jwElOaVGDSkqrnj4T46lj3fBgnqPhIiIyPd5eFhAnmj4cM/PUppMwIsv6j0KIiIi/8A1lOS0u3eBEiWA69f1Hol1iiKZ1H379B4JERGRf2CGkpyWN6/UmHSkJqUezGbg3Xf1HgUREZH/YIaSsiUlBahTR2pGpqbqPZo0qip9xH/6yXMDXiIiIl/DgJKy7eBBCSpNJr1HIgwGICwMOH5c6loSERGRe3DKm7KtRg3g/fc9JxOoaVL0nMEkERGRezFDSTliMgFPPQWsWqX/zu/XXgMmTdJ3DERERP6IASXl2N27QMeOwO+/6xdUDhoEfPON52RLiYiI/AmnvCnH8uaVDGXHjvoEdKNGyVQ3g0kiIiJ9MKCkXBEYKDur33lH6kCqqmuvZzTKNb/+Gvj8c8/v3ENEROTLOOVNuS46GujdGzhxQjbKuEKTJsD33wPly7vm9YmIiMhxzOtQrqtTB9i/X4qLBwfLVHROp6MtGc/ChYFp04CNGxlMEhEReQpmKMmlEhKABQuAKVOkCLrRKIXQHfmpUxQJRFNTgUaNpId4165AnjyuHzcRERE5jgEluYWmAdu3Axs2AHv2AH/+CVy5Yvv4UqWAevWAiAigUyegWjX3jZWIiIicw4CSdLN+/WGMHTsf7703ASaTgrx5ZYq8ShWZ2iYiIiLvwICSdKFpGsLDw3HixAmsW7cOrVu31ntIRERElE3clEO6WLx4MU6cOAEAiIyM1Hk0RERElBPMUJLbxcfH46GHHsLVq1cBAEajERcvXkSRIkV0HhkRERFlBzOU5HZvv/024uPj791PTU3F3LlzdRwRERER5QQzlORWu3fvRr169ZDxx65SpUo4ceIEDOyfSERE5HWYoSS3MZlMGDBgABQrfRL/+usvbN++XYdRERERUU4xoCS3mT59Og4dOoTU1NRMzxmNRsycOVOHUREREVFOccqb3OLChQuoVKkS7ty5Y/OYvHnz4tKlSwgLC3PjyIiIiCinmKEktxg5ciTu3r1r95jk5GQsXLjQTSMiIiKi3MIMJbnc2rVr0a5duyyPMxgMqFq1Kg4ePMjNOURERF6EGUpyuYkTJwKQdZIBAQE2g0VN03D48GFER0e7c3hERESUQ0a9B+BumgacPw9ERwPHjwO3bwN37wJ58gBBQUDlykBEBPDggwCTZLlj+vTp2Lp1K+Li4hAbG4vY2Fjs378fMTExCA4Oxu3bt9Mdv2bNGkREROg0WiIiInKWX0x537wJzJ8PrFoF7N4N/NegBaoK3F/BRtMAk0n+XqAAULcu0L498NxzQKFC7h+3L5s0aRI++OAD3LhxA0lJSYiLi7t3q1+/PgrxG05EROQ1fDqgPHQImD4dmDMHSEqSx5z5ai0ZyoAAoFcv4KWXJHtJOffaa69h5cqVOHnypN5DISIiohzyyTWUcXFAt25AjRrAd98BiYkSSDobOlvOSU6WDGfdukDbtsC5c64Ztz+Ji4tD0aJF9R4GERER5QKfCyh/+AF4+GHg55/lvmUKO6csr/Pbb0CVKsCsWc4HqJQmNjYWxYoV03sYRERElAt8JqBMSACeegro0QO4cQOw0owlV5hMspFnwADJVl675prr+DoGlERERL7DJwLKa9eA5s2BFSvkvrsyh7/9BjRqBFy86J7r+ZLY2FhOeRMREfkIrw8ob9yQYHLfPtdlJW1JTQVOnQIaN5Z1m+QYs9mMy5cvM0NJRETkI7y6DmVSkpT1OXLE/cGkhckEnDkDPP44sG0bEBqqzzg8kabJ92bPHuDo0bSanyZTElJT38eRI3Vx6BAQHg4YvfonkYiIyL95ddmgN98EJk0CzGa9RyI1LZ9/Hvj2W71Hoq9bt2RH/E8/Sc3PmzflcaPx/kLxGkwmEzQtAACQN6/syG/XDhg0CChdWpehExERUTZ5bUD5559Agwaet9N63TqgdWu9R+F+hw4BX38NfP999mp+AhKUm81Ap07AsGGS9VW8flEGERGR7/PKgDIpCaheHTh9Wr+pbmsUBShWDDh2TDrt+IOrV4Hhw4FFiyQLmRtlmiyv8+ijwNy5QNWqOX9NIiIich2vzP9MmgT8/bdnBZOAZNfi4oAxY/QeiXusWCE1P3/4Qe7nds3PAweAWrWAiRNz77WJiIgo93ldhjI5GShZMq0ftyfKlw+IjQXy59d7JK6RmChrHRcskHWRrv4JMhiA2rWB5cuBBx907bWIiIjIeV6XoVy+3LODSUCm5OfP13sUrnHrFtCmjUxxA+5Zw6ppwMGDQP36wIkTrr8eEREROcfrMpSNGwM7dnjGzm5bDAaZCj569P6dzd7v9m2gZUvZva3HcgOjUdam7tgBVKrk/usTERGRdV4VUP71F1C5st6jcNyOHZJV8wUmk5T12bhR37WrRqNsfNq9GyhRQr9xEBERURqvmvLets2545ctA+7csb/jev58WZfZvbtMrdq6ff21c9dWFGD7dufO8WSffiqtJvXeCGUyAZcuAf37e17JKCIiIn/lVQFldDQQEOD48ZGRskGmVy/rz4eGAl27Ar/8AkRFSTYx423OHDl2+XLnxmowyHh9wZEjwHvveU4Al5oKrF0rJYWIiIhIf1415f2//8lUp6MUBTh7Frh4EahbN/PzQ4YAM2YAHTsCq1dbf41Tp6TgdoUKzgdUFSoAMTHOneNpTCb5vh88qH928n4GAxAcDBw/DpQqpfdoiIiI/JvXZChTUyWocYbZLBnGiAigWrXMzz//PHDhgmQnrWneHKhYEZg9O3vZudOnZVe0N5sxA9i/37OCSUD+PZKSgFGj9B4JEREReU1AefMmcPeu8+fNmiWBZf/+6R8PDwfq1ZOA09aO8QEDJJCaPdv56wIS9MTFZe/c3Hbz5k04m4w2m4HPP3fRgHKBySTrZM+d03skRERE/s1rAkpLf2hnxcQAW7YAvXvLDmELS4A5a5b18woUAJ58EtiwAfj33+xdG5Ai4Hq7ceMGihcvjubNm+OEE4Ucf/1VsqyevCjCYABmztR7FERERP7NawLKnNSdjIwEHngA6NxZ7quqBJhbtsgaSWuefVY29Hz3XfavC3hGvcx///0XiYmJ2LJlC6pVq4b33nsPiQ5EulOnpg/CPVFqquzAT07WeyRERET+y2s25Vy9ChQpkr1zAwNlY84ffwCdOslt5Urgueds7xSOjgZKl5YNHznpI9279//h4YcVFC1aFEWLFkWxYsXu/RkcHAyDGyqf//7773j88cfv3TcYDChTpgy++eYbtG3b1uo5cXFA8eKenZ2834oVQJcueo+CiIjIP3l4/ilN/vyyazs7Gb+kJGkVOGiQBEn9+8uazB9/tH58rVrAo49K7cWcBJMAEB39G9avP44rV67AnGHw+fLlyxRkZgw8LX8vXLgwVFXN1hhiY2PT3dc0DefOnUO7du3w5JNPYsqUKSiVYav0rl3Wg8lly4C2baWo+I0b1q83f77U9SxdWto0tm8vvbgrVZJd9+XLZz3mAQMkO5yQkHVPdKNRisgzoCQiItKH1wSUefJIl5zjx7N3fmQkMHQoMHq0BDjff297feOAAWnn5ETRosDRo5sBAKmpqbh69Sri4uIQGxuLuLi4dH+PjY3FoUOHEBsbi9jYWCRlWDSqKAqKFCliM+DM+FhQUNC9c+Pi4qAoSrqA1vL3n3/+GVFRUZgwYQKGDRsG439z3Hv2SKCWMaCOjJTanb16WS/2fn9tz7g4oE8fCeJ37ZIPBI7UES1ZUoL58+ftF6W3MJmAP//M+jgiIiJyDa+Z8gakzM/8+dnPGu7fD1SvLoFNvXoS5GSUN6+UEjp2TPqGZ5eiAB06yNS6szRNw+3bt9MFm9YCUMuf165dy/QaISEh94LMK1euICYmJlOGNKNq1arh22+/Rf369dGhg5RTyvjT4WxtT4Mh7TVWrZLyTVllKFeulHOuXQOeeirrDCUgx9y44Vu904mIiLyF12QoAaBOnbTONdkRGQlMmSKdX6wFk4Ds7C5UKOebcRRF6l9mh8FgQEhICEJCQlCxYsUsj09JScGVK1dsBqD//vtvlsEkABw+fBgNGjRA27ZtER0dZXXK21Lb8+23JTg8fDj98xlrezr7ceXZZ4GmTYFHHgE++MDx827dkh3pFSo4dz0iIiLKOa8KKP/3v5xtEvnqK7nZs2iR3HLKZLKewXOFgIAAlChRAiVKlLD6fKdOnXDOwWKNBoMBd+7cwdWrto+ZNQt4801Zi/rKK2mPW2p7TpyYvbWuDzwAfPmlvPb5886fHxvLgJKIiEgPXlM2CJCMX4UK3jGtWaQIcN/Gal2dtxOdWdZMVqhQAZMmTcKlS5fw66+b7QaE2a3tmZXp04ETJ6yvzXREdmuVEhERUc54VUCpKMDw4XqPImuqKhuA8uTReyQi4y5vQDKRAQEB6N69OzZt2oRTp07htddeQ9GiRR3KAmentqc9Tz4p5ZwGDXL+XAtPqPlJRETkj7wqoASkdqSnBGq2aBoweLDeo0hj2bRjyUZWqVIFX375JS5duoQFCxagadOm6ephBgRknQVeuhSIj5c1k4DsnC9ePHs744ODgWnTZDnChQuys7tAgbR/5wIFgPs2rduUL5/z1yYiIqKc87qAsmBByYRlsySjyxmNkmkrXVrvkQiTyYTU1FQEBgaid+/e2LFjB44ePYoRI0agUKFCVs8xGKT8jz2W2p5t2zpW29OeIkXkNV57TYJUy61XLyAkRP6+YEHWr1OwoPPXJiIiopzzqrJBFhcvAlWqyM5eTxt9YKDsfHZgc7bbHDlyBGXKlEFoVlHifZo3BzZtsn9MnTpSr/Lzz4Fhw6S255Ahto+3VTYob16gfv3Mx7/5puz4btcOuHJFdufbEhgoRdA99YMGERGRL/OqXd4WJUrIBo7evfUeSWaffOJZwSQAVK1a1elz6tWTVpX2an5GRwMHDgCjRsn6VmvT3eHhUgIIkCxkUBDQrZvcP3pU6n3evQts3pz53H79pFe3tecyqlmTwSQREZFevDKgBGQ6dMkSYM0aCTr0ZjRKEPbSS3qPJHfUqeNYAfmsant27w6MHZv+saVL5c+xY4Fx43I6UlnzWa9ezl+HiMgfpKTITFp0NHDmjHSNM5lkpickRBqA1KkDlCrlHVVVyDN45ZS3xaVL0nf7yhV9g0pVlU4te/Z4XnYyu86cAcqV03sUjps3zzMz1kREnuDUKWDmTODXXyWYTEmRxzO2w9W0tGRCoUJS//npp4FnnnFscyT5L68OKAGZMm3YUNZT6hFUqqrsRt682X2FzN2lUSPpke0JGWB7goLkw4UjLRqJiPxFaqq0wP3qKwkkVdX53+eKIiXZQkOBgQOBF14AKlVyzXjJu3ndLu+MwsNl80iBAumLbLuDqkqpmg0bfC+YBKTmp6cHk0ajlC5iMElElGbXLlm/3qULsHGjPJad3+eW+r43bwKTJwOVK0tZvFu3cm+s5Bu8PkNpcfIk0Lo18O+/7ilwrapS7mbtWpl290XJybIB6r8ylh7r8GEgG/uOiIh8zt27sj79448lu+iKpICiyCbLOXOAli1z//XJO3l9htKicmXZGGLZFOOqHb/Kf9+xPn2A48d9N5gEZCr/pZfSvmZPo6rAY48xmCQiAuQ9qUYNCSY1zXUzTGazLDNq1UreIxzZwEm+z2cylPfbvBno21eylUDu1Kq07HQrWlR6Vbdvn/PX9AYJCbKs4MIFz2ttaDTKLsUaNfQeCRGRvqKjJVvo7v0EBgPQsSPwww+yS5z8l4fmnnKmaVPZrDN9uhRAB7K/vtJyXvnywBdfyNS6vwSTgJSQmDvX84JJgwF4/30Gk0RE+/YBzZrpszlV02TjT9euskyK/JdPZijvp2nAtm3SK3rVKuD2bXk8ICCtbML97n88Xz5ZlzlsGPD44/5dj+vFF4FvvvGMwFJVZZp7z57MJS+IiPzJqVNS2ufmTX03USqK1B1euNC/3yv9mc8HlPczm4GYGJkaiI6WLGZCgvSlzpsXCA4GHn5YCrpGREhpBHZfEQkJ8n35+29918soikyr/PmntHEkIvJXqalAgwaSofSUdYzffQcMGKD3KEgPfhVQUs6cPy+/vC5e1OeXl6LIEoS2bT/HlSs/4dVXX0XLli2d6lFOROQrPv0UGD1a71GkFxwsyZoyZfQeCbkbA0pyyj//AM2bA+fOuTeoVFWZ3l61CpgwoQU2/ldYTVVVNGzYEB07dkS7du1QrVo1GDjfQkQ+zrKj29rSLT0ZjbKec/16Tn37GwaU5LRLl4B27YADB3JnB31WVBUICwN++QWoXx/Yu3cv6tSpc+95g8EAg8EAs9mMokWLomPHjmjfvj1atmyJAgUKuH6ARERu1rQpsH2750x1Z7RokbRrJP/BgJKyJSVFap2NG+e6emeWll+9eknrsEKF5HFN01C9enUcOXLE6nlGoxEmkwmqqqJ+/fro2LEjBg4ciCJFiuT+IImI3OzgQaBmTb1HYZuiyJr7Xbv0Hgm5k0+WDSLXCwgA3n1XNjdZCovnVgF0g0FuBQsCy5cDCxakBZPyvAGDBw+2ObVt+u8je2pqKrZt24a33noLM2fOzJ3BERHp7Ouv3d9q2BlmM7B7N7B3r94jIXdihpJyLCUFWLZMsojbtskvuuxMw1jOK18eGDEC6NdPprqtuXr1KooXL34veLRFVVU8+uij2LhxI4KDg50fFBGRB7l5EyhWTKqTeDKjEXjuOdn1Tf6BGUrKsYAAoEcP4I8/gEOHgEGDpKOQhdGYOXtpMKSvIRkSAnToAGzYIHXVRo2yHUwCQOHChfHEE0/AaOdjuqqqqFixItauXctgkoh8wg8/SL9uT2cyAfPnA4mJeo+E3IUZSnKZS5dkSnzPHuDoUallmZgoNT/z5QMqVpR1NnXqyN+dnTJfu3Yt2rVrZ/U5VVVRrFgx/PnnnyhdunQufDVERPobMEC6l2WcnFm2DGjbFihRArhxw/q58+dL8fHSpYE2baTrW+3aUnP57FmZHcooJAR47z2gVi059oEHgLFjZf28I3buBOrVc+pLJC/lwaswyNsVLy5Zxw4dXPP6rVq1QrFixRAbG5vpudTUVIwfP57BJBH5lB07rC8pioyU9oe9eskay4xCQ+X5X34B4uKAPn3kd/SuXfJh3lbXscKFgcGDparHihUyA+Uog0ESCgwo/QOnvMlrqaqKgQMHQr2vnZGiKAgMDETDhg3xwgsv4IcfftBxhEREuScxEThxwvpzUVHSfKJ/f+vP9+wJBAVJ4AlIhrJGDaBvX9uvCQBnzsgGyWbNgLfecm68qiqzVOQfGFCSV+vXrx9S76tZpCgKVq1ahU2bNqFHjx545plnMH36dB1HSESUOw4dkh3U1pjNwJw50jbYWlva558HLlyQwBNwTw1hk0mmvMk/MKAkr/bQQw+hcePGAKSc0MKFC9GyZUsEBARgzpw5GDlyJF566SWMHTsWXC5MRN7s0iX7z8+aJYFlxixleLhMO8+ZYzsgdZW4OPdej/TDgJK83tChQwEA06ZNw9NPP33vcUVR8Pnnn2PixIkYN24chg0bli6bSUTkTbIqFRQTA2zZAvTunb5OpSXAnDXLdWOzxRt2pFPu4KYc8no9e/ZEgwYNUN7KFkWDwYA333wTRYoUwZAhQ3D16lXMnTsXefLkydUxpKQAd+7IGidVlV3swcHsZUtEuceRSZbISGDePKBzZ9n5raoSYG7ZIiXZ3M3dGVHSDzOU5PUMBoPVYPJ+AwcOxNKlS7FixQp07NgRCQkJ2b6e2Sw7I6dNk0/+jzwCBAZK3cwSJaQGZ/78cmvSBHjtNelre/Zsti9JRITAwKyPWboUiI+XNZOAlAYqXjxtM4675c2rz3XJ/ZihJL/RtWtXrF27Fp07d0aLFi2wZs0ap/p7X7sGfP+9dAT65x/JPqqq7a5At28DW7fKovSUFHmsZUvpAtS+vZxLROSo+1vQ2pKUJB9gBw2SQLJ/f+mu8+OPrh+fNY6MmXwDM5TkV5o1a4bNmzfjzJkzaNy4Mc46kDa0lOIoUUKyjf/8I49rmmMtJi3BJABs3ChTUQ8+CHz+efZaVBKRf6pZ07HjIiNlDeXo0fLhdfFifTrWqCpQt677r0v6YEBJfqd27drYtm0bkpOT0ahRIxw9etTqcZoGzJ4NVKkia5KSk3NeasOyJ+jCBQlOIyKAw4dz9ppE5B9CQ613s8koOloKkY8aBeTJY326Ozwc6NZNbsWLS41Ky/3w8PTHtm0rj3fqJPcfeSTt2Hz5bI9D0+R3HPkHtl4kv3Xx4kW0adMG58+fx+rVq1G/fv17z50/DwwcCKxdK1Pbrvpfoqry+mPHAm+8kX5nJhFRRr16ST/vrApWDB8OTJkCHDlivS7lmDHye8eajK0VT58GypWzfmy5clL83JaNG6UoOvk+BpTk1+Lj49GpUyfs3bsXy5YtQ5s2bXDgANCihaw7cteUtMEgnSt++kkyBUTO0jR544+Olg9EiYkSdAQGAgUKyHRpjRqObewgz/Xll8Crr3rH7mlFkbXnBQroPRJyBwaU5Pfu3LmDHj16YO3atXj//V/wySdt7r0Zu5OiAPXrA+vWASEh7r02eadz56S24ObNwO7dwK1b8rgl8w1IoGn5WVZV4OGHgQYNpBVfixYsbeVtLl4EypRx/5NaZvMAACAASURBVO8nZxmNQMeOwPLleo+E3IUBJRGAlJQUPPXUOKxcORqKEgqzWZ93WVWV6aHVq1lug6wzm4HffwemTgVWrpQPIs4GF0ajZN8rVpSp0eeek7JX5B26d5dAzdM39f36K/D443qPgtyFASURpG5beLiGuDgNZrO+e9UUBRgwAJg5U9dhkAfaswfo10/WxVmCwpywZCfz5gXGjwdefpnlrLzB5s2evS7RYAAqVAD++osZcH/CXd5EkDfSy5cNugeTgGSgvv0WWL9e75GQp7h7F3jnHenHfPy4PJYb2SlNk1tSkpSYadAg7fXJczVpItUnPDn4HzGCwaS/YYaS/N6aNUCHDnqPIj1FAYoVkzf30FC9R0N6OnpUyrOcOOG6agMWRqMEAZ99JlPh5Lm2bAGaNtV7FJmpqgS7e/dKySLyH/qnY4h0dOOGFC1XPOx/gtkMxMVJrUryX7t2AQ0bytShOz76m0xSiH/ECODNN91zTcqeJk0k6Pe0310AMH8+g0l/5IE/ikTuM3MmcPmyZ5bgSE0FvvsurTMP+Zfdu4HmzYGEBH129H78sZSnYVDpuSZOlB3fnjL1bTAA774L1Kql90hIDwwoyW+ZzdKX2xODSQtF4eYcf3TsGNCqlayd1LM8zBdfAB98oN/1yb7gYGDhQkBRNBgM+kb+ljaLb7+t6zBIR1xDSX4rKkr63Hq6ggWl9hzLCPmH5GTg0Udl/ayn1BrctMkz1+sRkJqailatvsLGjcNhMCjQNPfvhDEagUqVgG3b5PcV+SdmKMlvTZ3qOVNF9ly/Lh10yD9MmCAbcTwlmFRVoG9f4PZtvUdCGaWkpODZZ5/Fli2vYeTI7VBVg9vXVFqK5W/cyGDS3zGgJL90966U5bH2pr1sGXDnjv12YfPnSyapaFGgTx9g0aK0jNLp09bPad4ciIyU6cyEBOlysmKFZKPsUVXgl18c/9rIe+3fL1PMnjRvlJoqP6tvvaX3SOh+d+/eRffu3bFs2TIsWbIEX375GFavli5bRqP7xtGoEfDHH1KVgvwbA0ryS4cP267jFxkJ5MsH9Opl/fnQUKBrVwny4uIkoKxaVXbkxsTYvubQoUC5csDkyTLVPnKkBKQ7d0qwaUtqqhxDvk3TpGi5J7KsN961S++REAAkJiaia9euWLNmDZYtW4Zu3boBAFq3lg+2bdrIca6qA2k0Sk/4KVMkM8kuSwRwDSX5qZkzgSFDrD+nKMDZs7JusW7dzM8PGQLMmCF9alevll/alv9Fq1YB1aoB5ctnPu+BB2RH+f2Cg4FTpyTAbdXK/phv3GBNSl/mqXUFLYxG4OmnZRMI6ef27dvo0qULtm/fjp9//hmtrPzi0DSZRXnpJVmqkFsbD1VVPuA2agTMmSOtO4ksmKEkvxQdbXtayGyWX5YRERIcZvT888CFC7KpB3B8ejJjMAnIL/ujR6X0R1b27XPsOuSdpk1z71Sls0wm4McfJStP+rh16xbatWuHnTt3IioqymowCciH3D59ZMbkww+BUqXk8eyuGbec17y5LNPZsoXBJGXGgJL80pEj9lvXzZolgWX//ukfDw+X9ndz5uTOp/7QUFlDeeSI/eMMBll7Sb7p0iXZeJUb7RRdyWyWJSHkfvHx8WjdujUOHDiA9evXo6kD6ewHHpAC9WfOyOxJ69aynMciICDztLiqpg88S5YERo2S4vobNgBdunhmMXXSnwd/HiZynVu37D8fEyOfwnv3Bl5/Pe2N3hJgzpqVO+OYNk2mvT/80P5xiiIbhcg3zZnjWRtxbDGbgenTJUhhn2b3uXr1Klq3bo3Tp0/jt99+Q0REhFPnq6os0enYUf4NT56UWZroaGmccOeObDIMDpZNPdWrywzNo48CRYq45msi38OAkvzS3btZHxMZCcybB3TuLDu/VVUCzC1bZN1jTv3f/8nrDRsmfW/tMRgcGzN5py1brD++bBnQti1QooSsobVm/nyge3egdGnZjNG+PVC7ttQFPHvW+npeQIKHDz6QcwsVks0cH30ELFlif6znzsmSD8s0KrlWbGwsWrVqhUuXLmHjxo2oWbNmjl5PUaTXdpUqwLPP5tIgicApb/JTjvSZXboUiI+XNZOAvFEXL547U37vvw+89550lZg2LevjNY29cX3Zrl3Wl1C4quIAIMHqc88B48YB7dpJq8fFi4GePbMeb3R01sdQzp0/fx7NmjXDlStXsGnTphwHk0SuxICS/FJwcNbHJCVJfcm2bSWQ7N8fuHlTNibkxPvvy5v4mDHSi9cRZrOU6SDfc+ECcOWK9eeiooDz5zOv5bXo2RMICkr7kNOmDVCjhhQiP3HC9jXbtZP1dC++KBUPNm0CBg+W2qyTJtlfI2c0MqDMDfHx8fj000+RmJho9fmzZ8+iadOmuH37NjZv3oxHHnnEzSMkcg4DSvJL4eGO7aiNjJTjRo+WDOXixYCN3/8OefddCSbHj5cpb0dpGlC5cvavS57LXnDmqooDXbvKOuKMH45mz5ap7Hr1bJ+bmgrs2ePYdci2r7/+GqNHj0bnzp1xN8N6lpiYGDRp0gRmsxlbtmxBpUqVdBolkeMYUJJfqlPHsdZ20dHAgQOyyzFPHuvT3eHhQLduciteXDJGlvvh4WnHvfKKBJJRUVK/sl699DdHxky+5+JF+8+7ouJAtWpSNSDj/4GDB9Oet0XTgH//de56lNncuXMBAL///ju6deuG5ORkAMCJEyfQpEkT5MmTB5s3b0a5cuV0HCWR4xhQkl+qU8fxbE5kpEwBHjlivVNI9+6y3nLpUskkFS2adr9797TjOnWSP9u1k843GW/2lC4tGyfI9yQm2p9ivr/iwP1Z9ZxUHChcGLh2LfPjlscKF7Z/PisO5MyhQ4dw/PhxAIDZbMaaNWvQs2dP7N+/H02bNkVYWBg2b96MMo4UqCXyEAwoyS/VrOl4LbWvvpJd1rayNuPGyfPWbuPGpR3XvLnt4+yVYFEUxzKY5J3M5qxL8ERGSk3Bzp3lfm5UHLD3gSqrD1uOZPfJtoULF0K9r9ijpmlYvnw56tWrh2LFimHTpk0oUaKEjiMkch4DSvJL+fIBjRtnv3OEO5nNsjGIfFNgYNZT1rldceDqVetZSEsW3Fr28n5BQdm7LklGct68eUjNEJVrmobk5GRUr14dhbNKERN5INahJL81fLjt+n+eJCTEsVIu5J0KFMg6I2ipODBoUO5UHDh0SH6mLL2ZLapXlz8PH7Z/fsGC2buuN7l8WdZQ//23LEtISQHy5pUKEVWrArVqOVYtIqMdO3bg/PnzNp9fsGABChQogKlTp8LA6vHkRRhQkt/q0kWmEa312PYURqMED9l54yLvUKuWY8dFRgJDh6ZVHPj+++xXHFi+XMoEdesG/PBD2uPPPSdliv780/a5AQHSQcXX3LghReJ//VXWNF+6JI8bDLLsxGCQwN8SgCuK9LOuXx944glZjuBI5YgFCxbAaDTCZKfP5vTp05E3b1589tlnDCrJaxg0zRsafhG5xtixsvM6N/pyu8qxY9LVgnxTaqpkoZOSsj52/37JIlrW1WbcJBYeDljKFb75JvDgg1JrEgCOHk3fD37dOtlE9sYbsg6zZ08JMp99Fli40P445syRWpe+4OBBaSc5Z450ozIYnPt9YDRKa9ZixYCXXgIGDpTORtakpKSgaNGiiI+Pd+i1J06ciDfffNPxwRDpiGsoya+9+KK8mXtiEkBVJYPEYNK3qapzWcrcqDgAAE8+Ka1F/+//gLVrJUB95pmsg0nAN0pYnTsnmd6aNeX7mpQkGUhnP1xaEo2xsfIBtUwZKRFmbSf8r7/+mmUwafwvzWk0GnH79m3nBkOkI2Yoye8tWCA7Zj2JwSBr606ckKCAfNvLLwNTp6YFJ54sXz4piu4NG9qs0TQp4D5ihGQkXfE9VxSgbFkJ2Bs1Snu8d+/eWLRoEcwZolbLFHhQUBA6deqErl27ol27dggNDc39wRG5CANK8nuaJusp16zxrHIoS5ZkziqRb9q9G/jf//QeRdaMRlln+d13eo8ke65ckQ+P69a5/lqqKtnOV18FPvoIuHv3DsLCwpCSkgIgLYgsXLgwnnzySXTt2hUtWrRA3rx5XT84IhdgQEkEWYBfpYpkXvReT6koZnTubMCyZQaPnIon16hdW9bz6f3zl5XoaO/clHP+vNSC/ftv935wNBhk087//vcF3nrrFQDAgw8+iO7du6Nr166oV69eupqURN6KayiJIKVYVq2SDIyeQZyqajCbD8Js7ou7dx3YpUE+Y8QIzw4mFUXWZXpjMHnhAtCwIXD6tPtnITQN+PlnICpqIHr27IsDBw7gn3/+waRJk9CwYUMGk+QzGFAS/eexx4AVKySodLSLTm5SVeChhwyYO/cSNmz4Ca1bt8a1rCpMk8/o0UPWzXoqs1lqt3qb69eBFi0kqNRrjarZDGzdmh/AHFSvXoOlgMgnMaAkuk+7drK+Km9ex2rK5RZFkZ2+27YBffq0xW+//YajR4+icePGOHv2rPsGQroJCgI+/1zvUVhnNMoaz2ef1XskznvxRSmLpPeGJ02T4vQzZ+o7DiJX4RpKIisOH5bF+wcPZt3FJCcURbIXQ4cCkyalL2B+8uRJtG3bFklJSYiKikLNmjVdNxDyCJoGtGkD/P67Z20QCwiQ7joPP6z3SJyzfLmUR/Ik+fJJTdBy5fQeCVHuYoaSyIpq1WTn7fjxkp1xxTInRQFKlpTgYfr0zN1wKleujO3bt6NEiRJ47LHH8Ntvv+X+IMijGAzArFnS39uTTJjgfcHk1avSqtLTZpdTUqQnO1M55GsYUBLZEBAAvPMOsG+f7A4Fch5YGgxyCwqS2oPHjqW9tjXFixfHpk2b0LBhQ7Rr1w4LFizI2QDI45UuDXzzjd6jEKoqa4tfflnvkTjvtdeA+HjPC9xMJmDTJunMQ+RLOOVN5KC//pI3+m+/BW7elMxlamrWb1iqmtaBIzwcGDlS1qKFhDh+7ZSUFAwePBjff/89Pv74Y4wePZoL+33cpEnA66/rd31VlTaOW7d69mYha2JjJTDXe92kLQaD/C44fNjzMqhE2cWAkshJiYmycWf37rSbrW5q5csD9etLq7rGjWVjQ3bfQDRNw5gxYzB+/HgMGzYMX375JUuO+LgPPwTefdf911VVoGpV4NdfgQcecP/1c2riRPm+eXIZJgD444/0nXSIvBkDSqIc0jTg4kXgxg0JNlVVFt4XLw64onPaN998gxdffBFdunTBggULkC9fvty/CHmM6dOlRiXgvo06TZoAK1d6X2YSkO/Rgw9KmSBPZjQCTz/tWO90Im/AgJLIC61atQo9evTAo48+ipUrV6JQoUJ6D4lcaN8+qTpw7Jjr1gRa6q9+9JEEsN6a/I6KAtq313sUjjEagbg4oGBBvUdClHPclEPkhTp16oSNGzfixIkTaNSoEc6cOWP1OLPZjHHjxuGPP/5w8wgpN9WuLUHle+9JoJebNVItRfwjImRN38sve28wCciaT2vfn2XLgDt37Gdd588HkpMlwzl+PLB9O3D5ssw+7Nkju8YzNj1o3hyIjJRgPyEBOHdOGiQ40lHIZAJ27XLu6yPyVAwoibxUvXr1sH37diQnJ6N+/frYv39/pmNef/11jB07FiNHjtRhhJSb8uQBxo2TepCDBsmyCkvVgOywBF0NGgCLF8t6vkqVcm+8etm92/rSgMhI+Z716mX9vNBQoGtX4JdfgOrVgb59gd9+kz+7dQM2bwa+/lo25d1v6FCpKTl5smRGR44EihYFdu60X8EBkMA9OjpbXyaRx+GUN5GXi42NRYcOHXDy5EksW7YMLVu2BABMnjwZo0aNunfcvn37UKtWLb2GSbns5k3JqE2bJoWyAUtmUUNqauYoMyBAaiACEjz17i3BULVqbhuyy2maTB/fuJH5OUUBzp6V9c5162Z+fsgQYMYMoGNH6ViVkJB5l/hXXwHDhgFlykgmEpBNS5cvpz8uOFi68xw+DLRqZXu8igJ06SLZUyJvx4CSyAckJCSge/fu2LBhA2bPno3AwEB0794dlv/eRqMRAwcOxNdff63zSMkVbtwA9u6VbNfOnSn46afdKFs2HCEhBREUBISFybR5RIRUHChf3jfL1Zw5Y78DzYcfAm+/LRnIw4fTP7dzpwSKZcrY3h3epw8wd65kdXfutD+W334DSpUCqlSxf1zJksD58/aPIfIGDCiJfERKSgqGDBmC2bNnQ1EUaJqG+/97BwUFIS4uDsEZW/KQT/nnn39Qvnx5rF+/Hq3spcd80I4dQMOGtp+vWBE4eVKmp195Je3x8HDJ8k6cKAGnLbNnSw3Z4sWBa9dsHxcaKsHt77/LdLk9RmNa5pjIm+Xi0m4i0lNAQABeeeUVLFiwAMnJyZmev3PnDhYvXowBAwZk6/VNJnnTjY6WKb6kJJliDAwEihSRDFi1arLWj/Rz7b9Ixx93/icm2n8+JgbYskWm+19/PW1Ku39/+XPWLNvntmolGcrJk+0Hk4AsQwgOloxoVkwmWfPpzRuhiAAGlEQ+4/z582jdujVSbRQrVBQF06dPdyqgPHAA+P572Tl76JDsgAUkq2KZMtW0tDdmo1EKYjdqBPTrZ32tGrnW9evXAfhnQOlInc7ISGDePKBzZ1m7qKoSYG7ZIuseraldG/jhB5nmfust+6//f/8nrzdsmCxDcHTcDCjJ23GXN5EPuHHjBlq3bo3Lly/bDCjNZjP27t1rdTf4/e7elWLLDRoAtWoBU6dKVvL+pKfJJNN0KSnpNy6YTBKEzpwpXYFq15Zpwjt3cuOrJEf4c4bSkRr/S5dKZ6vnn5f77dvLFHZkpPXja9UCNmyQ1qvt26f/f5DR++9Laae335YspSMURTZMEXk7ZiiJvJzJZEKXLl1w/PhxmLPoNWc0GvHNN9/Y3JyzYgUweLBMaVvq7WWnH7LlnIMHZTrxlVckMO3Vyzc3g+hlzJgxOHfuHAoVKnTvFh0dDUVRcOrUKRQuXBiFChVC/vz5/aL3uyOdfZKSgEWLpPRS8eLy83nzJvDjj5mPrVVL2k+eOQO0bi3H2fL++1LWacwYWYvpqOBg/p8g38BNOURe7uLFi6hYsSISExNhNBphyiICDAoKQmxsLEJCQu49dvUqMHy4vNEqSu73QDYYZGq8Y0fJXpYokbuv768qVqyIv//+GwEBAdA0zea/vaIoKFiwIJYuXYpmzZq5d5BulJwMhIRkvcmlTh0pVP755zI1/f33UjbofjVryk7tc+eAFi3sr5t8910phD5+vASWznjsMZluJ/J2nPIm8nIlSpTApUuXsHLlSgwdOhQVK1a895xqZWHWnTt3sGTJknv3168HHn5Y1ogBuR9MAmntAqOipIwK6+7ljqFDh0JRFKSkpNj9IGE2m3Ht2jUULlzYjaNzvzx5ZA1vVqKjZWnGqFFyTsbp7sqVJTMJAO+8IwXf69VLuxUpknbsK69IIBkVBaxenf64evXsjyMgQJaGEPkCZiiJfNC///6LX3/9FevXr8e6detw/fp1KIpyb0q8bNmy+Oeff7B4sWwgMJtd1yM6I8v03owZMr1O2Xf16lWULFnS6q7++6mqiiFDhmCaowv7vNgLL8hu7ayylMOHA1OmAEeOZC7u/txzkrW0pV8/YM4c+fvGjYC9pG9W09kLFwI9e9o/hsgbMKAk8nFmsxmHDh3Chg0bEBUVhU2bNsFsNmPGjBsYOjQUgPuCyYymT5duLZR9/fr1w4IFC+xmKMPCwhATE+MXG3UiI2V9pLe8s5086RstL4kYUBL5mcTERMyc+S9efbWyQ2VWXMlgkAzNM8/oOw5vtnv3bvwvi3nTmTNnYtCgQW4akb6uXpXuM1kkbXWnKECNGsC+fXqPhCh3MKAk8jOxsbJm8uZNz8jiBAZKwfTy5fUeifeqXbs2Dh48mGmXv6qqqFatGqKjo62up/V0p06dQps2bQAABQoUQKFChRAWFoYCBQogNDQUoaGh9/5u+RMAFixoiUWL1GxVKHCn2bNl+pzIF7BsEJEf0TTZzZqQ4BnBJCAlhvr1k7VoCrcJZsvIkSPxvKWw4n1SU1MxY8YMrwwmASAkJARnz55NN51vMBigquq9MkhmszlT7dXu3SfBZHrNrWN1Vmgo0KOH3qMgyj389U3kRxYvBn7+2bGOIu5iMknZlBkz9B6J9+rRo8e97JyFqqro27cv6tevr9Oocq548eLo3r07jMa03IelPFJKSgpSUlIyBZOlS5fGzJmDUKuW53afUVVZ5+lIIXYib8EpbyI/cfMm8OCDnjPVnVFgIHD6tBSbJueNHj0aX3zxxb0AKygoCDExMSju5d9QR9aIAlJrs3Dhwti3bx9KlSqFbdukxqOn/awrClCwoGzG8YM9UuRHmKEk8hPz5nluMAnIJorvvtN7FN7rhRdeSLeG8oMPPvD6YBIA6tati4iICCh21kMYDAYEBARgzZo1KFWqFADpJ//yy563jMJslp3oDCbJ13jYfzUicgVNk5p7nsxslv7Hnr6RwlNVrFgRDRo0AACUKVMGw4YN03lEuePOnTuoVauW3baimqZh/vz5iIiISPf4Bx8A5cp5ztS3qkr70S5d9B4JUe5jQEnkB7ZulSk2T81OWly6BPzyi96j8F6WjTkff/wxAgICdB5Nzpw9exZvvPEGSpcuje+++w558uSxeeyHH36Ip556KtPj+fIBCxZIeSq9+2Wrqizn+OorfcdB5CpcQ0nkB3r3BpYs8fzsn6oCrVpJGzuyLTFRWgdaWgjeuCGPqSqQmnoLlSvnR5060rP6oYc8b9rXFk3T8Mcff2Dy5MlYvnw58ufPjwEDBuCll17CwoULMWbMmHSZSkVR0LNnT8ybN+/erm9rVqwAunWTD1R6vOOpKhAWBuzYwSLm5LsYUBL5gdKlgfPnMz++bBnQti1QooQEJdbMnw907y6v0aYN0L49ULu2vDGePWu7fmStWsCYMdKrOCxMjl24EPj0Uwl+bAkJkbWeemeUPM3ly9JScN484NgxWSJgMEiwYjan9WBXVQkgLa0Hg4OBxo2lXFSnToDRA4vFJSUlYdGiRZgyZQr279+PKlWqYMSIEejTpw9CQkIAAJcvX0apUqWQ8t8Xpqoq6tati02bNiFv3rxZXmPJEuDZZ93bZhSQ73dYmJTFytjikciXeMnnViLKrqtXrQeTgGwOyJdP1nVZExoKdO0q09BxcUCfPkDVqsCuXUBMjO1rhocD27fL+rVRo4COHaVk0fvvA4sW2R9vQoL91/Y3O3fK971kSeDtt6X3tCV41DTJOt+/vDA1NX0f69u3gV9/BZ58EihTBhg/XpYWeIILFy7g3XffRZkyZdC/f3+ULFkS69atw5EjRzB06NB7wSQAPPDAA+jVqxeMRiNUVUWJEiWwcuVKh4JJQGo+rlwp1QTcFVQrClCqlGQmGUySz9OIyKetW2eZ6Mt8UxRNO3dO03bvtv78kCHyGh06yH2DIe25Vas07fRp6+eNHy/nVaiQ/vEZM+TxsDDbYwI0bfFiXb5VHuXKFU3r2VO+H0aj/e+XMzdF0bR8+eTfwmx2/9dlNpu1HTt2aM8884xmNBq1/PnzayNGjNBOnjyZ5bl79+7VAGhBQUHa4cOHs3X906c1rVmz3Pt+Wrupqvw5bJimJSRka5hEXocZSiIft3ev7V2uZjMwZw4QEWE9g/L888CFC2lrGh2dKrRkyDJOo8fHSwbNXp/lgABZG+jPfv5Z2mP+8IPcz821r2azLDl44QWgRQvgzJnce217kpOTsWDBAtSrVw8NGjTA7t278dlnn+HcuXOYPHkyKjmwuLB27dqYMGECoqKiULVq1WyNo1w54PffpZB+vny5uwPcsvmnVClg0ybZgBMcnHuvT+TJGFAS+bgLF+xvypg1S4KM/v3TPx4eDtSrJwGnnYotVs2ZA1y/Dnz9tayxDAkBOnSQdXzTpgF37tg+12wGLl507nq+IiVFgvgnngCuXXN9R6M//pB/56VLXXeN2NhYjBs3DmXLlkXv3r1RsGBB/PLLLzh58iRGjBiRqcNPVt566y00adIkR2MyGORn8eRJ4LXXZI0jkP3NS5Yp9AoVgC+/lN70TZvmaIhEXocBJZGPS0qyn1mMiZHWh717p19bZgkwZ81y/ppnzgANGkjW8++/gVu3ZB3mnDnAyJH2zzWbZcyOOnnyJD755BP89NNPzg/UgyQmSn3COXPkvjs2jphM8r3u3j33i8pHR0ejb9++ePDBB/HJJ5+ga9euOHLkCNatW4cOHTrYLVTuLqVLAx99JB9g5s2TTL1lM5ii2F5reX9FpoAAoHnz21i69Br++gsYMYJZSfJPHrjfj4hykyPZxchIeUPt3Fl2fquqBJhbtgCnTjl/zbJlgVWrgNhYKddy+bJkO999V7KVAwfaP99eZk7TNOzduxfLly/Hjz/+iJMnTwIAmjdvjm7dujk/WA+QnCybnzZscH9ZG8v1Bg2SAKpfv+y/VkpKCpYvX47Jkydj+/btKFeuHCZMmID+/fujYMGCuTJeVwgMlJ/33r1lU9j+/bLsIjoaOHFCgv3kZDkuJASoUUOCzzp1gEceAUqVKo8+fRJQt+5xPPjgg3p/OUS6YEBJ5OMCA7MuwbN0qaz3ev55CSjbt5cizG+8kb1rfvSR7BCvVSttenvrVuDKFWD2bGDuXAlWrVEUWdt2P5PJhD/++APLly/H0qVLceHCBaiqeq9vtdFoRHh4ePYGqzNNA557ToJJZ5cW5Lb+/YHChaW8kDOuXLmCmTNnYvr06Th//jyaNWuG5cuXo1OnTlA9pU2Ng0JCpMxS48aOn5OQkIDExEQ0bNgQ27ZtQ9myZV03QCIPxYCSyMeFhWWd9UpKknI+gwZJINm/v9SC/PHH7F2zVi1ZR5ZxreTu3fJntWr2C2xjGQAAIABJREFUA8qwMKlNuGHDBixbtgwrVqxAfHw8jEYjTP/tUEm9L42ZmpqKKlWqZG+wOlu4UEoqeYrnnpOs3AMPZH3sgQMHMGXKFCxYsAAGgwHPPvssRowYgRo1arh+oB7i5s2bSPyvsOrFixfRuHFjbN26FeXKldN3YERuxoCSyMfVrOnYLuHISGDoUGD0aMlQfv+9/QLk9ly4IEFjcLDUQbT4r9U0zp2zfW5KCpCQsA1BQY9B0zQoinKvO4rJxheiaZpXZigvXgRefFEyyJ7QYkLT5IPEiy/a/jCRmpqKlStXYvLkydi8eTNKlSqFMWPGYNCgQShSpIh7B+wBYu4rmmo2m3Hx4kU0atQIW7duRYUKFXQcGZF76b8qmohcKiLCseMsbfxGjQLy5JEAM6PwcFkT2a2bZDKDgtLu3x/PffklUKSITOM+/TTQvDnw1lvA559LYe6sWiu2b1/03po7s4PzwN6WodQ0YPBgCbg9IZi0SE2VJRAZA8rr169j0qRJqFixIp588kmYTCYsWbIEp0+fxltvveWXwSSQPqAEJOCOjY1F48aNMz1H5MvYepHIx2maTCHfvJn1scOHA1OmSNBnrS7lmDHA2LHWzx07Fhg3Lu1+s2bAm2/KBoYCBYB//5WNOhMnSkkcW4xGCbLM5iS8/vrr+Oqrr2AwGJDVr6rQ0FCEh4ejSpUq6f6sUKECjB7YbzAqSjLBnshgAAoVkg5LMTFHMWXKFMybNw8mkwnPPPMMRowYgTp16ug9TI/w8ccf45133km3BAOQdb2FCxfG1q1bHaqxSeTtGFAS+YFWraSYs96bPhxRu7YUY7eIiopC7969cePGjUxv2hbly5fH4MGDcfz4cRw7dgzHjx/Hzf8i6ICAAFSqVClTsFmlShUE61jfpW1baYno6lqTOVGt2gQcPvwOihcvjqFDh2LIkCEoVqyY3sPyKIMHD8bs2bOtLsdQVRWFCxfGli1b8PDDD+swOiL3YUBJ5AfmzQP69tV7FFkzGIAvvshcqzI2NhZ9+/bF+vXrM52jqip69eqFuXPn3ntM0zRcvHgxXYB57NgxHDt2DBcuXLh3XJkyZaxmNYsWLQpDVlvjcyAmBqhUybOmujMyGMwIDj6KGTMO4Omnn0aePHn0HpJHatasGTZv3mzzeVVVUbBgQWzdutXrlmUQOYMBJZEfSEqSNY8ZWyF6mrx5gUuX0jqX3E/TNEydOhWvvvoqzGbzvWylqqoYP3483nrrLYeucfPmTRw/fjxTsHnq1Kl7r1mwYMFMQWaVKlVQvnz5XCmD8/rrsp7Uk7OTFvv3y8Yusq5kyZK4mEVrJ1VVERYWht27d6N8+fJuGhmRezGgJPITnh7EGI1Ssiarji2HDh3C008/jb/++uvehp3ly5fjiSeeyNH1k5OTERMTky7ItASeCQkJAIC8efPemz6/P9isXLkygoKCHLqOpklJnqtXMz+3bJlMhZcoYTv4nz9fOtuULg20aSPrMGvXlozn2bPS6tKWRo2At9+W3faBgbLbfu5c4IMPrB9vNErnl88+c+hL8zt3795Fvnz57K7vtdRLzZs3L9avX5/jtpFEnooBJZGf+Ptv4KGHPHuade9eCY6ykpiYiNGjR2PatGkAgGPHjrlsOlHTNJw/fz5ToHns2DFcunQJAGAwGFC2bNlMWc3w8PBMu5/PnAFslSjs0EFaVL74ovRBzyg0VEoNrVsHPPkksH69ZJ737wfq15c2gLYCyp49ZenDDz9I7cuEBKBiRaBkSWD8eNtff4MGwPbtDnyj/NCJEyfs/twFBASgXbt26NmzJzp27IiQkBA3jo7IvRhQEvmR11+XbJOnbc5RVWn552w/6dWrV+OXX37B1KlTdenIEh8fnynIPH78OGJiYu5lTwsXLpwuyLx6tSkmTLC+Q1pRJMt48SJQt27m54cMAWbMADp2BFavTl+/ctUq2ZlvLaAsWVKKlc+dC7z0knNfY2CgBJ9e1vDGLdasWYMOHToAwL16qZaMZN++ffHVV18hNDRU51ESuQcDSiI/kpgIVK8O/POP50x9KwpQtChw/LiUF/IFd+/exalTpzJtCDpx4gTu3HkHBsPr0DTrpYw+/FCmpatXBw4fTv/czp1AmTJyy/ihwF5A+f77UtKpbFkJWJ115Ij0rKb0vv32WwwePBiKoqBFixbo1asXnnjiCTz77LNISEjAFlvtoIh8kOcVZyMil8mXT9bgNWyo90jSmM3SlcdXgklA1lpWrVoVVatWTfe42WxGs2Z38ccfttN9s2ZJ/c7+/YFXXkl7PDwcqFdP6ng6m2Fu0kTWbFapAvz8swSe167Jms3XXwdu3bJ//r59DCit6dq1K/Lnz4+WLVumW9rQo0cPPP/88zh//jxKlSql4wiJ3Iedcoj8TP36wLvv6j0KYTDIesE2bfQeiXsoioJbt/JB02yXJIqJkT7nvXvLphiL/v3lz1mznL9uqVLS1ejHH4ElS4CWLYFJk6SU1Jo1WY3ZfiF6f1akSBE888wzmdbJdunSBQEBAfjRVv9KIh/EgJLID40bJwGKC0stZklRpGXjlCn6jUEPSUlZHxMZKTvBO3eW+6oqAeaWLcCpU85fU1EkOz1hAvDRR8DmzcCnn0o7zMaNgccft39udnu6+6uwsDC0adMGS5Ys0XsoRG7DgJLIDxkMwMyZUqZHr+t36QIsWMDNHtYsXQrExwPPPy/327eX3dzW+qs7wlKiaN269I9beqo/+qjtczVN3w8e3qpHjx7YuXMnzpw5o/dQiNyCASWRn1JVCVBGj5b7iht+G1iuMXCglK/xx+Yr+fJlfUxSErBokdSkLF5cssk3b8qUdXYcPGj9cUugaG9NpqY5NmZKr3PnzggMDMQPP/yg91CI3IIBJZEfUxTgk0+Ab745jWLFNJdmC1UVKFgQWL5csqNGP90SWKKEYxm/yEj5Ho0eLRnKxYuzP/X800/yZ7t26R9v317+3LnT9rlms0y/k3Py58+P9u3bc9qb/Iaf/konIkAKgg8cOBDbt2/Hp5/OxPHjg/DddxLImEy5cw3La3XvDnz1FVC4cO68rreKiJCC5Fl9f6OjgQMHgFGjJPC3Nt0dHp62+7p4cdl4062b3D96FDh2TP6+YQOwcqWUD1IUCSAjIoAxY6Tc0LZt9sdib0qcbOvRowd69OiBfftO4/r18tizR/5dd+8Grl+XTLTBILU+H3hAdvFHRAB16gC1agHBwXp/BURO0IjI71y8eFEbPHiwpiiKBkADoO3YsUPTNE3bsUPTevXSNKNR0xRF02TS0/mbositWzdN27RJ5y/Yg/z8s+Pfw+HD5ZzDh60/P2aM7euMGZP+2MBATZs4UdPOnNG05GRN++cfTfvwQ03Lk8f+GIKDNS011R3fGd+Smqppq1bd0RRlpWYwmO/9n1BV29/r+//P5cmjaf36adquXXp/JUSOYWFzIj+SkJCAzz77DB999BFSUlKQel918+vXryMsLOze/bg4KVHzzTdSCB2QaWtNy7zmTlEk02J5udKlgQEDgEGDpGQNpblwwbu+J02ayK5wckxysnQzmjxZ2p0qihlmc/ZWl1my+7VqAa++CvTq5Z61zkTZwYCSyA+YTCZERkbinXfewfXr1++1BbQICwvD9evXbZ5/9ar02Y6OBvbsAWJjgTt3JLgMDgaKFJFpOsutaFFXf0XerVQpCSw9nWUN54QJeo/EO+zbJ+WdLEsNcuvdVVHkQ1zTptIEwFYveCI9MaAk8mGapmHVqlV49dVXccpOAcO6deti165dbhyZfxs7Fhg/3vN6qltz4gRQubLeo/BsycnSMvPDD+W+q9qaGo1AQADw+efS153lnMiTMHlO5KP27t2Lxx57DF26dMHff/9t8zhVVREeHu7GkdGgQXqPIGuqCrRowWAyK9evA489Jh8QUlNdF0wCMv2dmAgMHSqbr+7edd21iJzFgJLIR/Xv3x/b/tu+m3GK+34GgwEPPfSQu4ZFkCnvJ57w7NJJqanA8OF6j8KzxcUBjRrJUhB3z/X9/DPQoQO7GJHnYEBJ5KMWL16MWrVqwZDFvJjJZEKlSpXcNCqyGDky90oz5TZFkY1VHTvqPRLPFR8vLSv/+su1WUlbzGZg40bgySdlyp1IbwwoiXxUlSpVsGvXLrz77rtQFAWqnarlDCjdr0kToE8fz2w9aTanFVanzMxmyTAfO6bvhwKzWdppDhum3xiILLgph8gP/Pnnn2jfvj2uXbtm9fn4+HgUKFDAzaOi69eBKlWAK1c8Z4OOqkoP8W+/1XsknmvqVM9bDhAVJa06ifTCDCWRHwgODsatW7dQs2ZNAIByXzG7ggX/n73zDovi6uLwb2YWQUCK2DD2jkYsWIgajcausUVNNHbjZ42xtxSNNUZT7NHYNTZsiSVqjL2hYldib8QaUQSVtnO+P04WWFi2sbuzwLzPcx/YnXbm7sydM+ee4qsqkwrh68u5Pp1FmRRFTvn0/fdKS+K83LzJqZScCVHkl4CoKKUlUcnOqAqlikoWJy4uDl26dEGZMmVw4sQJ/PHHH/Dz80uaAlenu5WlRQvnsHYJAlsn168HvLyUlsY5kWWgRw/n832VZeDpU2DoUKUlUcnOqAqlikoWZ8KECbhy5QpWr14NNzc3NG3aFOHh4WjTpg0AqCmDnICffgI++US5vIKCwG3TJk6Bo2KYnTuBI0ecT6EEODBo2bLkpOoqKo5G9aFUUcnCHDlyBHXr1sWUKVMwduxYvWVEhD179qB06dIoUaKEQhKq6EhM5HKVK1c69riiyG3jRqB1a8ceO7PRpAnw11/KRHWbg0YDDBjAZR9VVByNqlCqqGRRdD6T/v7+OHTokNEobxXnQJaBr74Cpk1jJc/eioskAT4+QEgIUL++fY+V2bl5E8gM6Vo9PLg0qoeH0pKoZDfUKW8VlSzKsGHD8OTJE6xcuVJVJjMJosjl+44cAYoW5c/2Og4AfPQRcO2aqkyaw88/O2eKp9S8egWsWaO0FCrZEVWhVFHJgmzbtg2LFy/Gjz/+iJIlSyotjoqF1KoFXLrEQRa6YBlboFMkfX2BceNOYv/+gjh79i/b7DyLs26dYYvx5s3A69eAsUQJq1dz8vEiRbhE47FjHEQTFQWcPs2lOFO/PHh6AtOnc57JJ0+4Es/48ablFAQOrFJRcTSqQqmiojD//ssPjalTuepFQAA/ePz92UpVoQLQuTOncjlwAHj50vj+nj59ik8//RQtW7bEp59+6pBzULE9OXMCM2cCN24Aw4cnKyzWWC11CcorVOD8knfvAoUKheHhw4do2LAhevbsiadPn9pO+CzGs2dARIThZUuW8G/VubPh5V5eQNu2wPbtQMWKQLdu7IfZrRvX4z54EFiwIG3eTz8/4H//A1xdga1bzZeVCDh1yvGlIFVUQCoqKg4nJobol1+IAgOJeOgnkiQiUUz+nLKlXlanDtH69URxcfr7lWWZ2rRpQ3ny5KGHDx8qc3IqduHNG6IVK4gaNSLy9k6+FjQaIhcX/ZbyWilWjKhDh5c0aNBqkuXk/R09epQAEACSJIm8vLxo4cKFpNVqlTtJJ2XPHsP3JcB9HRFBdOqU4eV9+/I+WrQg8vHh3yv1OnPm8DqFChneh58fLx8/Pn05Urfbt5XoKZXsjFpYS0XFgdy5wyliliwBYmL0rU3GAjBSLzt+nP3s8uQB+vfnPIZ58wIrVqzA1q1bsXnzZhQoUMAu56CiDG5ubNXq1o1Vhnv3gLAw4Nw5tlq/ecNT4zlzAgULAkFBQJUqbNns0qU/5s79FTlznsd3330HAPD390/at1arxcuXL9G3b18sXLgQv/zyC6pWrarUqTodYWHct4buUVkGVqwAxo0D3n6bXRVS0rMn8OABV7JJL4H9yZP8t1Ch9C2h1shcrJht9qWiYg6qQqmi4gBkGZg/nytsJCYm57GztkKK7sH27788VT53LjBp0hOMHTsY3bt3R9u2bW0juIpTIgjsDlG0KLtJmIL+m/+cMWMGPDw8MH78+HRfOM6fP49q1aphwIABmDx5Mnx8fGwpeqYktZKYmqVLgTFjgF69gGHDkr8PCABq1uSofWP3eoMGQEICB0jZAo2GZf7wQ9vsT4UhYneRsDD2fb12jf1nY2P5Rc7Dg0upBgUB1aoBb72lXG5ZJVAVShUVO3PrFtC9O1sU7YFWC7x4AQwalA85c67FV1+pmalV9Ll9+3bS/xMmTAAAjB8/Hp6enoiJidFbV/vf28rPP/+MdevWYdasWejcuTOE7PRkTEV0tPEZhJs3gUOHgC5dgFGjkl8Ye/Xiv0uXpr9to0ZA166cOzIy0jbyCgLLrJJxiICjR9kgsGNHsg+7iwv/zil9VXU5XXW/f+7cQKtWwMCBrGBmddSgHBUVO7JjB0+DnThh3+PoBrX4+OaoUcMLp0/b93gqmYtrqUxfEyZMwDfffIN8+fKlu41Wq0VkZCS6dOmCevXq4c6dO3aW0nmJjTW9zpIl7HbSqhV/liRWMA8d4sAqQ1SpAmzYwONDqroDGcYcmVXS59UrThVVoQJXjwoJ0Q+ITEhIG/gky/pVlCIjOcK/enWgalVg+fKs/buoCqWKip1Yv54rj8TGOq5Um1YrICoKqFePH2QqKjExMXj27Fma7ydMmIBEExembqr88OHDWJ+Nc9G4uJheZ+NGnino2ZM/N28OFCjAiqYhKlcG/vwTuH6d142Pt528gHkyqxhm3z6gbFmuOvT33/ydtWO4brvz5/naqFgRCA21jZzOhqpQqqjYgd9+4zQisuz49B1aLSuxTZty+hCV7M2N9MxjAO7du2d0KluSJHh4eGDevHkYMWKEPcTLFOTMadoXLjYWWLuW77sCBXi6++VLtmylpnJlYO9e9sdr3Nh0KjBrZVaxjOhooF8/4P33gYcPk2PmbYHOh/b2beCdd4DRo7OetVJVKFVUbMypU0CHDrYdjCxFltni0agRP7RUsi/GFEog2QppiFatWuHatWsYMGBAtq62VKxYci5PYyxZwuuNHMlWx3XrOPo+JZUqsTIZEcH354sXtpc3MZEDtlTM59IloHz55Hyg1gZMmkKr5efCzJlAYCArmFkFNShHRcWGxMYqZ5lMjVbLfkC9e/PUWjaOqcjWXL9+HRqNxuT0tg5BEEBE+PHHHzFkyBCD6xBxCqywMG43bvC1FheXHO0aEMDRrkFBnMYoMxMUxD5zpggL46nNIUM4OCP1dHeZMqxMAsAXXwClS+svv3mTMzfoaNqU+zJXLv5cvnxy5PbOnWmVVR1E2SMIxFacPMnK/atX9lMkUyPLrEwGBwP79/Nvm9kRyNjrqYqKikWMGQPMmOG4QclcFi3i8m4q2Y/evXtj5cqVZiuUgwcPxl9//QU3NzeEhoYmWSaJ2C93/nxg167kaVqNJtnqoiN1tGvevFwtZsAAttBlNq5fZ2XQHD77DJg9G7h8mQPyUtK9OwdmpEePHpzTUsft2+nnkixWLP3ZB42G89y6uponc3bm7FkOuomNNR7Jby8kCfDx4dzCqV8wMhuqQqmiYiNCQ9k3xhnvKHd3IDycSzqqZC9q1aqF48ePm1xPFEWMHj0aU6dOxYkTJ1CrVi38+OOP6N37c6xcyUrS1ausrFgToKDbrmZNVro6dABy5LDihBRAlrmE4qtXSktiHoGBbClVMc7duxxp//KlMsqkDo0GyJ+fldu8eZWTI6OoPpQqKjZi6FDr6iw7gvh44KuvlJZCRQmuX7+e7jIPDw/4+vqiT58+iI6OxtSpUwEAwcHB6N+/P8aM2Y8SJbQYNCg56XZGo11PneJ0OlWq8AM0MyCKQP365vlRKo0k8fStinFkmS3GpnKMOoLERODRI7bgZ2ZUC6WKig24cMH5p/JcXDhy0c9PaUlUHEVMTAxy/eeAp/Oj9PT0hJeXF/z8/HD69GnkMGAmjIoCBg+Ox8qVOSAIBCLbO+DqYny++IKbs1sr//iDA20yA9euZf7pU3vz889cttbZ2LCBrfeZESe1p6ioZC4WLHB+64VWCyxbprQUKo7kxX8hxJ6enujcuTN27NiBZ8+eYdKkSbh8+TIiDZRmOXeOA2p+/ZU1PHsokwBfj1otMGkSJ322VQ1re9GkCVC4sNJSGEeSgIYNVWXSFHfu6JfIdBYEAejbF3j6VGlJrENVKFVUMsjLl+xI76jk5dYiy8CcOc4XMKRiPwoVKoSrV6/i2bNnWLFiBZo3b44cOXKgTZs2EEURmzdv1lv/6FEOUHjyxHHTgETsmxkczMEvzooosu+ns7q1APybDRyotBTOz8iR5kXtOxoifp58/bXSkliHE98aKiqZA2PpOwyxeTPw+jXg7Z3+OqtXs99j6sp4+fJxWhGi5PQhlnDvHtSyjNmMMmXKpJnWzp07N95//32EpMi6rUud8vq1433KEhOBx4+BunWdL29qbGwsHj16hKtXryIg4Bhy5XLOFFySxJbJli2VlsS5+ecfHoOd1QCg1XImgKgopSWxHCefpFNRcX5OnWL/RHPfeJcs4RQqnTvzVHlqvLx4+fbtbClKybx5GauuIAisUNaoYf0+VLIGHTt2xKeffopHjx4hJqYAGjfmlxilLNiJifyy9N57HKzj42Pf4z18+BDLly/H8+fP8eLFC0RFRSEyMhLPnj3D8+fP8fLlS0RHRyMh1Y09duwFTJtW0b7CWYEs84uos7veKI0ucbkzExcHrFoFDBqktCQWQioqKhmiTh1dTRzzmigSRUQQnTpleHnfvrzfFi30v2/XjujlS6KuXXn5hx9adlyASKMh6tVLub5ScR6ePXtGGo2G5syZT8HBfG1Yej3Zo0kSUffu9j//ZcuWEQCSJIk0Gg2JokgAjLYKFSqQVitT69Ysp9J9lXJMGTVKJlmW7d9xmZj4eKK8eZX/vUw1QSAqXZoos/2capS3ikmiojinWWQkW8cEAXBz43q1gYHZu2asLHMVi9evLdtuyhRg3DigYkUu+ZWSEyfY+b9w4WRrka8vcOUKMHUqR5QfOAC0bw9s2mS5zOXLc9JlFZVmzZrh+vUWuHnT+Uwh27cDLVrYb/+xsbEoW7YsIiIiIJthlnVzc8OVK1dQvHhxPH4MlCvH/m5K+yRLElCiBPD6dSk8fnwXQUFBqF69OipWrIi3334bb7/9Nry8vJQV0knYvZurD6Vm82b+3t8//anm1auBjh2BQoU4QKt5c059Vbo0uxIVL552m3r1eKw2RHAw5y42xqlTmavikWocV0nDw4dcg/bkSVZu7txJf11R5AoSwcFArVp8wxnzDcxq3L1ruTIJAEuXclWdXr30ow0DAjjx87Rp+g+q2bO5asbcuexnlhGuXmU/nWxcmlnlP+rW7Y1du5zP6U4U+d74+29+mbIHbm5umDVrFtq2bWvW+jNmzEDx/7SG/PlZCWnSJNmupASSxC+0v/0GTJ5cE2vW3ERoaChOnz4NbQpH2IIFC6Jy5coIDAxMUjTLlStnMGVUVubkSe6z1D7Clrohde3KBpWTJ/ladXExftyxY7m8YkpSGxJSIwiZT6GE0iZSFedAlokOHCBq356nckTRsikdjYbN9G5uPGV77pzSZ+QYzp2zflpj/36iJ0/0pxpnzOD9liqV/F3z5kRxcUQVKvDnevV4HWumvHUtOlqpHlNxJho2jCMgXvEpPkNNkoiGDbPv+cuyTO+99x5pNBoCDE9zS5JEtWrVIq1Wm2b7rVtZTkFQpo80GqLVq1mWZ8+ekbu7e7rnIQgCubi46J3XuHHj7NvBTkarVvxsS92Plrohpfy9t20jun3b8HYZGas1GqLevZXqKetQo7xVcOgQT4O+9x6wdSu/vcmyZZGeiYl8G8TG8tte5cpsSfv7b7uJ7RRkJEBmyRIus9WqFX+WJK4gcugQcOMGf+flBSxcCEyfbttpakui0lWyJrdvA3/9lQOACfOKQmi1HEBhzQyAuQiCgBkzZhitcy5JEpYvXw7RQL6g1q2B339nC5USFn8iHjNq1AB27MiN4cNHG5ST1yW9ACOtVou33nrLUaI6BaGhhl0UZJlTv1Wrlrb+OgD07Ak8eMDJ7QHud3uTmMj1vTMTqkKZjXn1Cvj8c/bz0OV/s0UqBd0+jh9nH8sZM5QvbWUuZOFIkZGHyMaNwIsXPFgB7JNToAArmjqmTOHo8blz2ZXA2xvw9ORl7u7WuxeokaAqCxc6d05FgMvirVtnn30TETZu3Ij27dtDEAQI6eQCmjp1KkobyRTevDm/BBYt6vj+1I2rYWFAt27A1q3j4OZmeo5UkiR88MEH6O+MpWLsxL//cmqq9Fi6lBXLXr30v9e5Ia1YYb2/7Lx5PI5HRQG7dgG1a5u33dWrnHkh06CwhVRFIY4eJSpa1LD53x6tWjWia9eUPmvjhIeHk6urKwUGBtKIESNox44d9PLlS6PbXLqUsX6ZP58oIYGoQAGiLVuIoqKIcuZMXr5/v2m5vb0tP+6bN7bpM5XMyZs3RD4+ykzTWtJEkSgw0PbRrufOnaP33nuPAFDLli0pNDSUvLy8CNCf6q5WrRolJiaatc/Xr3mKXhCUiwBnd6VEEoRxBEh655PyvAoWLEjPnj2zbac6OVevmu4/c92QUjZjU96VKxP9+CNR69acDaRHD6LLl3nMb9zYvN80M/1MUFoAFcezdSuRi4vjlEndQOfrSxQWpvTZp8/9+/eTBl2dT5UoilS9enUaN24c7d27l16/fq23zfPnGeuXoCDez/ffs5/kwoX6yytVYj+clO3zz3mbr7/mz5Y+vLy9M186ChXbsmWLZdfM5s2sMBl7eVm9mtOyFClCNGkS0bFjRE+f8kvS6dNEffpYP+Zcvmyb83769Cn169ePRFGkcuXK0R9//JG0bO7cuXqKl0ajoctWHPjoUaK332a5lUvFpCUglAA/Sq1QiqJIR48etU15B1GIAAAgAElEQVSHZiLOnzfdb1268Lrt2iU/tx4+JDp4MP1tjCmU6Y2/9+6Z738fEaFot1kElBZAxbFs3MiDuhJO5JJE5OHBjs/OSt68edMMwCkVTI1GQ3Xq1KFvvvmGDh8+THFxcVS0aMb65dw5Ip2/f40aptfPiKO3IBA1aqRU76o4C2PHWqbstGjB2/Xvb3i5lxfRq1eseLZoQXT3LiuVzZoRNWzIL0yJiURLllh3jyxblrHzjY+Pp59++ol8fHzI29ubfvrpJ4qPj9dbJyEhgcqWLUuiKJIgCDR16lSrjyfLrFh27uzYF3f9lkDA3wT4641lbm5utGnTpox1aCbkwgXTfebmxkaCbdv48wcf8LbduqW/jaUKJcAzU0R8PFPrPnigWJdZDJQWQMVx7N6dHI2tzADHSqWPD9GVK0r3hmHatm1rVoJjQRBI97bfpk1shqa4PvuMj23u9HlGIwezWWCnigEaNrRsHLAkCtbHx7CyOmcOr1OokGXXrIsL0aBB1p/r7t27KSAggARBoL59+9KTJ0/SXXfv3r0EgAIDA9MonNbw88/KjbXc4gm4SoAfSZJEdevWpbZt2xIAGjhwIL3JRr4v166Z12em3JBSN2sUygULWCZXV9PrPn+uaLdZBJQWQMUxPHhAlCuXkm/LyU2SuAqAM45l06dPJ0ky7HtkqDVs2JCmT5edol/NbdnQOKGSAlm2zu92yhTeXjedm7KdOEH0zz/GxxddhafgYMuPXb265ed57do1+uCDDwgA1a1bl86ePWvWditXrqRbt25ZfsBUbN2q7Mt7cksgIIxy5y5Ajx49IlmWad68eeTq6kqVK1emq1evZvhcnR1Zlunvvx+a1V+m3JBSN0sVSh8fovv3ic6cMb1ujhxs2c8sQGkBVOyPLLPlwFlKqwH84Bk9WumeSSYhIYHOnDlDQ4cOJVNKpCiKpNFoaOnSpUREFBqqfH+a3+8yPXyocGerKMq9e9ZdOyVLsmvGDz/ofx8QwPudOtX49suWsY9l7tyWH9vV1Xy/36ioKBo1ahS5uLhQkSJFaMOGDQ4vSfjkCfuMO4dCSQRoqWvXm3oynj17lsqUKUMeHh60atUqh/aPvdBqtXTz5k3avn07zZgxg3r16kXvvPMO+fj4/Dd+3zarv0y5IQUE8AzRhx+y1f7x4+TPAQHJ6/36K9G0afx9vXpEn35KFB7O98H775uWo1o1RbvTYqC0ACr2Z+VKpQczw00Q2LLhaGRZpnv37lFISAiNGDGC3n33XcqZMyfpIiCNKZOSJFG+fPnoRArBZZmTjjvPwyO9/k4kV9cdNHLkSLrirD4HGSUmhujsWXZg27eP6MgRjgSLilJaMqfh7FnrryFromAB9ttNTGSrj7XHThUPlwatVktLly6l/PnzU86cOembb75JE0TnKHQFIpS+51M2SWKrWEqio6OpW7duBIB69OhBMTExivSXpcTFxdGlS5coJCSEJk6cSJ06daJKlSqRm5tb0ljt4eFBQUFB1KVLF5o6dSpt2bKFmjePMut3MeWGNH58+rKNH5+83ujR3OfPn/M0+uPHPENUrZppGVxciAYOdERv2g4oLYCKfXn6lKe6nVHZ0U19JyTYtw+io6Np//799O2331Lbtm2pYMGCSYNOkSJFqEOHDjRz5kw6cuQIvX79mgICAsiQMikIAgUHB9NDAya+hQuds49Tt7Zt51Lu3LkJAL3zzjv0yy+/UFRmVrauXyeaNYu95suUMT7nWrw4R0l8/z0/KbIpx49bf/1YEwVbpQo/UI8c4Sk8a49tLH3KsWPHqFq1agSAOnXqRPfu3XNMZxpgwwbl7/P0xtvy5XkaNzXLly8nd3d3KleuHJ0/f15v2f79+6lRo0Z0//59B/VgMtHR0XTq1ClauXIljR07ltq0aUNlypTRe/H38/OjOnXqUJ8+feiHH36gXbt20d27dw1WNpoxwzncvsxtGQ1GczRQWgAV+zJtmvPfQL/9ZrvzTUxMpIsXL9LixYupT58+FBgYmBRk4+npSfXr16cxY8bQli1b6EE64XN9+vQxWIqtb9++FGdoNCYuZejhoXxfptcEgacsZZkoNjaWNmzYQE2bNiVRFMnd3Z26d+9OBw8edPjUoFUkJLCDWsOGfHKiaL4/h66uKEBUuzbRunWGn7BZmBMnrL+OLI2CrVyZ6N9/iU6e5EjwjFzDhoIT7t+/T5988gkBoKpVq9Lhw4dt0kcJCQnUrVs3WrdunUX3REICUf78zvtyKQhEc+calj08PJwCAwPJzc2Nfv75Z4qLi6OxY8cmBSCuWLHCJn1riKdPn9KhQ4do4cKFNGTIEGrSpAkVLlxYb/wtVKgQNWrUiAYPHkwLFiyggwcPGg2wMsSRI8r/Bpa0v/+2U4fbCSgtgIr9SEwkeust5W8KY02SWC+wlkePHtHWrVtp7Nix1KBBA8qVKxfprIkVK1akTz/9lH755Re6cOGC2QmKlyxZkjSI6fwlFy1aZHK70aOdW3lfuDCtzPfv36cpU6ZQyZIlCQCVLFmSJk+erIg1wiSyTLR4MZG/f/LFk9GLDyDy8+Psw5nJ+z0DmJOPz1gzNwpWp0yGhdkmiXpsbPI5vH79miZPnkzu7u6UL18+Wrx4sdn3tzlcvXo1aQxo166d2YrLxo3K3+fGmiCwIT89Hfn169fUv39/AkA+Pj5JyqSLiwuNGTMmQ32qczXavXs3/fTTT9S3b1+qW7cu5cmTh1KOt6VLl6ZWrVrR6NGjafny5RQaGmqzWRRZ5hdrZ1X4dU0UOXgtswGlBVCxHzorQmZo16+bPp/Xr1/TkSNH6Pvvv6eOHTtS0aJFkwaiAgUKUJs2bWjatGm0b98+kxVujHHp0iXSKaV58+alY8eOmbXdq1dExYo5n++URsNVGgzMACUhyzIdPHiQunfvTu7u7iSKIjVt2pQ2bNhAsSmf5Epx926yRdJeLTiYy2lkcR4/zlg3mRMFW6kSK5PnzlkXhJO2RVNgYCVq3bo1tWzZknLnzk2SJFH37t3pn3/+sXkf7dq1K2lskSSJfH19aePGjSa3s6bQgBLtwAHD8suyTCtXrqQcOXJQanefVq1amdV3CQkJdO3aNdq6dStNmzaNunXrRtWrVydPT8+k/bm6ulLFihXpo48+ogkTJtD69evpwoULDkljNGuW8yuUAAf0ZDYEIiIDFRlVsgDNmwN79jh/HW1JAoYNA777Lvk7IsL169dx4sQJhIaGIjQ0FOfPn0diYiLc3NwQFBSE4OBg1KxZEzVr1kThwoXTrcVrKVqtFv7+/ihWrBh+++03+Pv7m73tkSNA3bo8JDgLbm7A5ctAiRLmrf/y5Uts2LABS5cuxfHjx5E7d2506dIFPXv2ROXKlS069tatW7FmzRosW7YMHh4elgtPxEV2Bw/mora2KDafHhoNF2P+9ltgyBDARteTM5I/P/DkifXbnzsHVKzI3VWzJnDyZPKyMmWAo0e5+7p35xrKKbl5M+13xiEUL34fQUHDsW/fPkRGRkIQBKR8dOXLlw/Fixc32IoUKQIXFxeLzu/nn3/GgAEDko6hO17Hjh0xf/58+Pn5pdkmPBwoX96iwyiCRgO0bQts2KD/fVRUFPr164d16RRPL1q0KO7cuZP0OTY2FteuXUN4eDiuXLmC8PBwhIeH49q1a4j/rwB1rly5EBAQgICAAJQvXz7p/+LFi0OSJHudolFevAD8/YHYWEUObxa+vsDDh4Crq9KSWIiS2qyK/dBqjSdjTd0sKa2WLx/nlFu7ln08tFrLE7umbm+/nUA7duygr7/+mpo0aUK+vr5Jb7Nly5albt260fz58yksLMwmCYdNERkZSQlWRgsNGeJcU9+6hNLWEB4eTqNGjaL8+fMTAKpSpQrNmTPHrDrAsixT+fLlCQA1btzY8t9NlolGjFCm03r0sH+0mIOJjo6mP//8kyZMmEDe3hcIkK3uHmNRsN27G5eje3fLjqXRyFS58p8kiiKVLVuWdu7cSVqtliIiIujw4cO0cuVK+uabb6hHjx5Ur149KlKkiF5xAlEUqUiRIlSvXj3q0aMHTZgwgVasWEGHDh2i+/fvGwze0KUd0u1D1yRJIj8/P9q6dWuabb7+2rB10p5lKz082Fvjn384r+/Zs0QffWS6TyVJP2r+2LFjVKhQIZNZLgYMGEAffPABlSxZUq+P8+XLR/Xq1aN+/frR7Nmz6c8//6SIiAib+GRrtcZnV6xhwADntSSLIl9LmREoLYCKfQgPt+witqS0GkC0Zw+Xslq5kisQZFSh5IoOLuTn50fNmzeniRMn0u7duykyMlLRfrSG16+JatVSfsASRX642GIwjo+Pp99//53atGlDGo2GcuTIQR999BHt3r07Xd+1sLCwpAeOIAj08ccfG3x4G0SWifr1U67zBIFzv2Riv8r79+/TunXr6LPPPqOqVasmKQDu7l2Jaz1br1A6uuXM2ZN++OGHdIPiUhMfH083btygP//8kxYtWkTjxo2jTp06UXBwcNLLka7lyJGDSpcuTY0bN6a+ffvSt99+S8HBwUn+g6mb7vtOnTrpvVg1aWJ4KtWeZSt37yaKjCT63/+I3nuPSOfq3amT6T49fpyDGCdOnEiiKJpVIUwURapRowYNHTqUFi1aREeOHDHr5dK865X7YNw4TjPl56cfa+fiwu7TbdoQTZ5M9Mcf7FZhDZGRbBhxphd/gJ8ZJUqYTpHlrEBpAVTsw+rVll3IlpRWA/QHTmtKTxlqv/12L3NEGZvBixfsR6aUUimK/ECyhzH30aNHNHPmzCTrY+HChenLL7+kmzf1EycPGDAgTbT8Z599Zt5vPHq08qO7IBD17Gl+Rm0FSUxMpHPnztG8efOoc+fOVKRIkaQ+L1myJHXr1o0WLlxI06bdJ0GQM4UPWcp26tRTm/bXq1ev6PLly7R9+3aaM2cODRs2jNq2bUuVK1cmb29vMqVY6ZQrPz8/2rZtG8ly+r6i9ipb2awZf/fxx/rr7t7NxzOmLAkC0YwZr8nf39/oOaZuAQEBJIoiTZw40SZBUPHxRCEhRHXrJstmqjywKCaPq5LE733791t+m27frvx1beh3OXo0w92qGFBaABX7MHQov9FZcjFbW1rNVgrlL78o2mU25/lzopo1Ha9UCgJR27b6UbH2QJZlCg0Npb59+5KXlxcBoPfee49WrlxJz549S4q4T90mT55sfMfr1ik/sqdsVvoMHDx4UC8Bvi2JiYmhv/76iyZOnEiNGzdO6n+NRpNkQdq0aZNeztQtW5zPImOqiaJM771nly40iq4/zW0dOgw1eh72KFu5aBHRy5dpx5ePP+Z133kn/f1pNEQdOkSTh4cHpbS6GmsajYbGjx9P48ePJ0EQqEGDBummXjNFTAzRhAlEefOyPBkZI3UKeOnSXDvdEj23a1flZ5KSr3WiYcOs6k6nAUoLoGIf2ra1PJLN2tJqtlAoXVyIvvhC0S6zCzExyTn67P0wlyT+zUePdrz736tXr2jVqlVUv359AqBXscJQ+/nnnw3v6NEjNtM4kwnNzY3oxg2z+0KWZfr2229JEASqUKGCTfo3IiKC1q9fT4MHD6agoKAkXzdvb29q1qwZTZkyhQ4cOECvXr0yuP2VK3yPOVO3mts2b7ZJF5pNVFSU0WtXp4DpWokSJWjUqINGz8EeZSuPHTNc9rV8ed5nnz7G96m7NMPDw2natGlUtWrVpHMzNP0tiiJ99NFHRES0b98+8vf3p7x589KuXbss6t8DB9hX1Nbjoe7arlnT/GQNz58TlS2rfFliSeJ69Zl1qlsHlBZAxT40bWrdhW1NaTVbKZTDhyvZY/Zl+3Z+G7fX27Ao8u+jRCnL1Ny8eZOKFy+e5oGU+oEcEhKiv6Ess4OUs5gMdM2cvEv/8ebNG+rcubPeeVqawkqr1dKFCxdo/vz59Mknn1CxYsWS9le8eHHq2rUr/fzzz3Tx4kWzfFITEjjVj9IPTWta/vyOfzmaNm1autetTtGqXr06zZ49O8kCbI6Lka3LVl69yn6EqdcvUID3OWaMcXkKFUp77hERETR//nxq0KBB0ktLykCdAJ0GTESPHz+mpk2bEgAaPXq0yaC7mBiiQYP42Pa8xTUarsik8z01xT//EBUurNz9odHwi0UmDBdIA5QWQMU+mFN43lCzprSarRTKwYMV7TK7ExnJwcOCYLsBVRS570aN4ihPZ+Du3bsmp9AEQSCNRkN//fVX8obOWrNO1+bPN3reDx48oGrVqqWx7uzZs8fodjExMbRv3z6aNGkSNW3aNMmHT6PRUPXq1WnIkCEUEhJi9fTi9OmZ0zIpCGy5czQDBw40eM0WLVqUpkyZQrdu3UqzzeLFps/H1mUrr14l2rkz7TY6hdKUG7Kfn/F+eP78Of3666/04YcfJs04uLm56flAa7Va+u6770ij0dA777xDd+7cMbivBw+IypVzvMtF8+bmWf3u3WMrsqOVSkliX/untnURVgwoLYCKfWjZ0roL3NLSaoDtFMrRoxXtModx9y5P7+umrywdxHTrFy5M9N13zjcYTZw40WT6EV1zc3OjU6dOsXXy7bedW/MpWDBdk8fp06cpf/78ac5bkiQaP3683roPHjygkJAQGjJkCFWrVi0pcMnb25uaNm1KkyZNon379lFMTEyGf4vwcMt9qZ2hSRIrIErk1B8wYECa61T3uxYqVIjGjRtHl1LVgl+61PZjq6mylRmd8s6Xz/w+efPmDW3fvp1WrVplcPnx48epaNGi5OPjQ5tT+SjcucNT3EpYAEWRJxfMuZWePNEPDrJn0w1zLVpwaqisApQWQMU+dOpkvRXM3NJqumYLhVKj4VQZ2Ym4OM7l2bMnT3mkfXvX6n12cWGLRb9+PIXujBlttFptmhq8ppokSXTf2a2Tuvb772nOecOGDeTq6mpQiRYEgWrUqEELFiygLl266LkCFCtWjD755BOaP3++RaVBLaFHj8w51S2KnH9RCfbs2WPyegVA5cuXp2+//Zbu3LlDa9fadmw1p2zlwoWGg3L+c3M0GpQDsJJnSyIjI6ldu3YEgAYNGkRv3ryhR4+4epiS16Ak8YydOS8nWi3XOndzs5/MGg1Rrlycci8TJJCwCCgtgIp9mDrVeoXSnNJqKZutorwN5ArOFmi1Wjp79iwdPHiSFi06T19+GU7DhoWTRjOY6tZdSSNHnqEff9xP69dvpg0bNlBoaKjSIqfL/v379R6+Go0mTeognaLl5+dHhQsXpuLFi9PLVq2cX/ORJHZo+w+tVktfffVV0vkYU0JEUaSgoCAaPHgwbdiwwS7lAlPz7FnaadLM0ARBVjRALz4+3qycjClb48aTzTo3W5atbNqU99Wxo/73O3eaThsEENWubfu+k2WZ5s2bR66urlSpUjUqX/6NU9zWosgvV+Zy82aytdJW8uv207w5+21mRTRQyZIEBVlfcjEsDDh/nqvPiSKwZEnadQICksuMFSgAuLsDH37In69c4TJk1sicHVm+fDl69+5tcNmhQ9xSUqBAATx8+NABkllOQkICChUqBD8/PxQuXBhvvfUW/P394e/vj4IFCyb9ny9fPmg0/w0/z57xRWTPsoq2QKsF/vwTuHkTrwoUQNeuXbFlyxYAABEZ3fTgwYOoU6eOI6RMYvlyICHBoYfMMIIgg+g48ua9CKCfw46r1Wrx4MED3LlzB3fu3IG7uztiYmLM2tbLyws9e1bBnj2m1zU1tpYpA+zdy/9/8QVQurT+8pRlK3ft4tK6CxYAXl7AjRtAp05As2bAJ58Aspy+HC4uQI0aZp2eRQiCgAEDBqBWrVpo2HA/nj3LYfuDWIEs8/3Qvj3QooXp9UuUAA4c4BKi8+YBGzfyPoz1aXoIAvf3J58AAwYA1apZvo9Mg9IarYp9+PffjL1NGSutBhClcgvTY/x4y4+XO3fWM/+bS0REhMEyb4aaKIo0duxYpUW2LVu22MYEsG0bO6mlzP6sa76+HBlw5EiG/TRf/vAD5c6d26zfS9eGDh3q0C6VZZ5qtKWVx95NkogCA2Xq1280AaBhw4bZzA1AV6rxyJEjtHr1apo8eTL17t2b3n//fSpZsmSa+8+QVT31fQhwon6dr2vx4hkfWy0tW+nhQfTTT3xpx8ayVdOc0osAR6bbiwsXiCTJuSoxiSJn2rAmmvrRI84lWru2vouCi4t+S2nN9PRkK+f33/NsQXZAICKyv9qqogSFCgH//KO0FKYRRaBRI37jzq4MHToUc+bMgdaEWdnV1RX3799H3rx5HSSZA/jyS2D69IxbKPPnBy5dYjNQ06b6y379FWjVCqhcmc081uLighft2qHIzp2Ijo4GAGg0GiSakN3V1RWxsbHWH9dCrl0DypY1f/3Nm7nL/P2BqCjD66xeDXTsyOPKlClAcDD/7+oK3L8P/PYb/4zPnlkuryCw5WbXLiB3bmDOnDkYMmQIWrdujdWrV8Pd3d3o9rIs4/Hjx0kWxtu3byf9f+fOHdy9exfx8fFJ6+fJkwfFihVLasWLF0/6v2jRohgwYADWrFlj8HcVRRGFChXCqlWrULdu3aTvO3cGNmywfmbI0YSHA+XK2X6/iYn8W16+7HyTDpLEv9PKldbvQ5b5/goL479v3gDx8XwfuLtznwYFASVL8nWdrVBao1WxHyNGOF9Kv/Ta0qVK95ayPH782GQycEmS6PPPP1daVNvTuLHtors7dOB9/u9/yd+1acPf2ao2eOXKRET077//0m+//UbDhw+nypUrJ1mt0gvOcSRr1lh2SpbUmwZ4/599xuX/6tcnGjmSy41eumRdVHmePETR0frn8Pvvv5O7uztVr16dHjx4QI8ePaITJ07Q2rVradq0adS3b19q0qQJlS1bllxdXfX6O3fu3FS1alVq164dDR8+nGbPnk3btm2jixcvmpUX9Ouvv05jtdT9voMGDTIYgT9vnnMnKUjZfH3NSqtqFeYGKCnZrlyxz7lnd6C0ACr248YN5W9cc5ruYZXdGTdunNFgAI1GQxEREUqLaVtkmZ9utryg1qzh8NeiRdmX4tEjLnBsq/1rNBxRkYo9e/YkVRsx1O7du+ewbh0+3DLFzpJ60+nto18/Xqd+fWvGAJkeP35CoaGhtH79epo+fTr179+fateuTZIkpQl68vHxocqVK1ObNm1o6NChNGvWLPrtt9/o/PnzFGWDPCzLli2j1MpkkSJF6MCBA+lu8/w5kaur8uOpqSZJnPTcXtSq5dyGDI0m6+c8VgooLYCKfWnUyLlvbknK/PVLbUVkZCR5enqSIWVEkiTqq3uqZyVevrT9ReXry2GUf/3FdcEjI4neesu2x7h5k4g4qnXv3r307rvvJv1Ohn6/nDlz2iSvpLm8+67lp5TRetMffmj9sbkVTeqvXLlyUWBgILVu3Zp69epF/v7+5O7uTgsXLqTnz5/bvf902Qp0iuzAgQPN+v1693b+ZAWCwFk57MGFC8qfnznNw8O83JQqlgGlBVCxL7//rvzNa6pdu6Z0LzkPEydOTLeO7m17PQWU5PFj+1xUupwqRESffGLz/cuXLtH27dupevXqZEyR1LWCBQs6tFvz5rX8tKypNy1JRO7ubJW6coXo0CHrq6F8+eVhOnPmDEVGRupVYyHi+tpNmjQhjUZDS5YssXv/3bt3jwBQ4cKFaf/+/WZvFxam/HhqrEkSuynYi/79nV+hBlipXrjQfv2QXVGDcrI4sgzUrQuEhjqfg7QoAp9+CixcqLQkzkN0dDQKFy6MqBSRERqNBp07d8aKFSsUlMxO/PMPR3bYg2PHAD8/y6JTzOSzOnUw98gRs9fXaDQoVaoUJEmCRqOBJElp/rf0s7Fl06d/gfh4V4vPa/9+oEIFoGDB5PFixgxgxAhOYXPjRvK6NWsCJ04kf96xA/j4Y8DMbDtp2LIFaNMm/eUJCQkYNGgQFi1ahC+++AKTJk2CYMeoh6NHj6JSpUrw9PS0aLt33+V+cbbxVsfevcD779tn3/nyAU+fmreupYFgoaFAsWJp1/n5Z6B/f8vkFASgcePsHQhqD1SFMhtw8ybw9tuAA4NMTSKKPJCEhwO5ciktjXMxc+ZMjBo1CrpbUxAEhIeHo6wdFCPF+fdfwF4R6/v3A3nyABUr2nzXYatWodu0abhy5QpEUYRsIkGdt7c3evfuDa1WC61Wi8TERIP/G1tmzv+6zxER1wG4WXxeXboAq1ZxTtnNmzkqNiKCo1nr1dNf192d89G6u3Pw/JgxwL17QIMGHPlqKevWAR99ZHwdIkq6Pzp16oSlS5fCzc3y87QnV68CgYEc+etMSBLnQrTXe+mDB8Bbb5m/fosWwPbtnJtxwYK0y728gIcPgd27gXbtgNu3+VocMUJ/vcePgTt3LJc3d24efrJdJLY9UdI8quI4Zs9Wfpohddu7V+lecU5evXpFefLkIYB9uDroIpezIq9e2e8C27+f6OJF++z77l0iIjpy5IhJ/0kAVLJkSYd2q6endadlab3plK1GDV53yBDrjr1li/nnFxISQm5ublSnTh166mzF7Ilo5kzlx9eUTRS5drc1ORjNxVL3KksDwW7fTr4ubdUcGCeXLRCVU2VVHMnAgWxhkCSlJeE3wgED7Dftktlxd3fH+PHjAQBEhC+//FJhieyIu7tlZg1nwN09aZq+du3aOHToEP766y9UrVoVAOcpTE1SVSAHkTOnddvFxgJr1/JUZIECQK9ewMuXQEiI6W1Pn+YcjGXKWHdsE6km9Wjfvj3279+Pq1ev4p133sH169etO6idGDKEK9E4+GdPF1kGli4FfH3td4ywMMvOV5bZWlqtGs+gpaZnT7Z6/vGH7WRMTViY/fadHVEVymyCKAJbt3K5RCUHOVEEPvgAmDVLORkyA3369IGrqytKlSqFwMBApcWxLzVrOsebjrlUrcoXcgoaNK1CiSAAACAASURBVGiA0NBQ7NixAxUqVACgr1g6WqEMCLB+2yVLeIwYORJo3pynos2Zwta9sKb0s7QESz06goODceLECWg0GgQHB+OIBT6t9kaSuN98fJzj0h4+3LySgxnh8mXLSxMuXcrb9Oql/31AAA8LK1bo77NuXX7BiY/n4w0bluZWNBuNhvehYjtUhTIb4ePDDtllyigzyIkiO0KvX+88b+5ORUwMcPgw8NNPcP3iCzzu0AEX69Rhc8f48cCmTewsRKS0pOmSmJiIx48fgyyRMTMVtzVSBFkQBDRv3hznzp3Dxo0bUapUqaRlkoNvuBo1WFRrSFlvOkeOtPWmW7Tgl9NevXiWoWlTLnYUEgJcvw4sXmz5Mb29gSJFLN+uRIkSOHbsGAIDA/H+++9j7dq1lu/EThQvDuzbB3h4KKdUCgL/TjNm2P9Y0dGWK5Q3bwKHDrHvbspngk7BXLo0+bsdO4DPP+frr3174MwZ4PvvuUa3NQgC8OqVdduqpIPSc+4qjicykqhmTcdXdfj4Y4P5oLM3p0+zs1CpUsk/iCgaLxLr5UXUsCGXF3KyjPCff/45ASBXV1cqVaoUNWvWjAYOHEgzZ86kTZs20ZkzZ9LmEbRl0nFHtF9/NasvEhMTafXq1VS0aFFq2bKlHXo7fTJarcRYvemyZYk2bGCfttevuV25QjR9unU56gWBL+eMEBcXR926dSMANHny5DRph5Tk/HkiPz9l0un072+/ijipqVvXOhm7dOHt27Xjz5JE9PAh0cGDpredPZu3rVzZ8uO6uBANHeqYvskuQGkBVJQhPp5owgS+ee2Z+Fyj4SSyixdzURQV4ifw8uVEVasmd5KlHatL9uflxZnhnSSZ5+LFiwnQD0hxcXFJk1vT09OTKlSoQMHBwdTlww/5PBz9tLWmubgQWRgEIsuywxWca9eU7ypLxoixYzN+zrIs0zfffEMAqGfPnhTnRG+v168TVarkmJd4jYbbtGmOHXMbNrROXlsEgllTVdXFhWj0aMf1T3YASgugoiznznFlDEGw7WCn03caNlQj6fT44w+iAgX0OymjTfdGMGiQ4uUfIiMjSaPRUGql0ljz8PAgeehQ5y7ppHtSd+2qaP+aiywTlSiReWpLHz5su3NftWoV5ciRgxo0aOCQqjrmEh9PNHlyssJnj34UBKKKFdkq6mhat7Ze7vnziRISeGjcsoUoKoooZ07T29Wsycf+3/8sP6aLC9H48Y7vp6yM6kOZzalUiX1RZs8GSpbk7zLi36jbNiiIo0X37AEKF864nJmeqCh2DGrWDHjyhL+z1OEoPbRa/jt/PkddHTxom/1aga+vL5o2bWqWz6AgCMifPz/Cw8Mh9O+ffB7OSmIip0vIBAgC8NlnSkthGkHgAIzatW23zy5duuDPP//E2bNnUatWLdyxJkmhHXBxAb74gsdbXWpUW/mSiyLvf8IE9oFVIo6vRAnr/XatDQTr1o3/pkywby4JCSyzig1RWqNVcR5kmWjfPq7JqzOeiaJxw5FGk2wFcXMj+vRTojNnlD4TJ+PYMX71doQFTneM4cOJEhMVOd1169YRYNwqKUkS+fv7061bt5I3bNjQea2UokgUGJip/DYiI4lcXZXvOmNNEIgWLLDP+f/9999UokQJypcvH4WGhtrnIFYiy+wj2KEDX/KiaLk1WWflLFyY6LvvLPbEsDmrV2fsWjh3Ltnfs0YN/WWdOhGFhBD16EFUvz5R27ZEa9bwukuXWn/MS5eU66+siFopR8UgL18CZ8/y225YGHDyJBAZCcTFsVXBzY3LbAUHszUyKIjfiq3Nf5dl2bMHaN2a81zYyiJpDoLANctWrbLebGAFt2/fxooVKzBx4kSkN7RIkoTcuXPj2LFjepHQCA0F3nmHx3pn5I8/OKQ5E9G7N7BypfOWAXR3Bx49sl+1rKdPn6J169Y4d+4cfv31V7Rt29Y+B8oADx8Cy5ZxYadTp5LLEIpickocIh4+iPi70qX5VunYEWjSxPrUObbk6lWgXDnrt//sM54pu3w5bV7KmjWByZO5LKifH1sXL1/mflu40Lohw82NE2s4Q1qnrIKqUKqo2Iu9e3mKW5Ydq0zqEASgbVtgwwa7jpqPHj3Chg0bsHbtWpw4cQLu7u7IkycPIiIi0pQklCQJ3t7eOHr0KMoZevqMHg3MnKlMf6WHJPHcWsocJpmE69f54exsZQABvjwnTAC+/tq+x3nz5g26d++OjRs3YubMmRg6dKhda4BnBCIuXxkWxn9jY1l5ypmTle7AQHZTsiQJvKOQZZbx9WulJTGP4GDg+HGlpchiKGofVVHJqpw8yT4Atgq8ycic4qef2vz0nj9/TkuWLKGGDRuSKIrk4uJCLVu2pDVr1lB0dDRt376dgLTT3D4+PnTx4sX0d/zmDVHp0s4z9S2KRPnzcxhqJuX775XvxtRNkjgYMD7eMX2g1WppzJgxBIAGDBhACQkJ6a67Zs0aunnzpmMEy2K0bKlMeiRrrr+vv1a6t7IeUFoAFZUsx6tXRMWLO49SBLADUoZP6xWtX7+e2rRpQzly5CBBEKh+/fq0aNEievbsmd66cXFx5OXlRSmVSS8vLzp79qzpA4WGct85S4jyH39kuO+UJDGRo2Gd6UEvSewz52gWLVpEkiRR8+bN6eXLl2mWz5s3jwBQ8+bNHS9cFmDnTuWvLXOaKBLdv690b2U9oLQAKipZjqFDlbdMpmyCwBmnnzyx+FTi4+Npx44d1KVLF/L09CQAVL16dfrhhx8oIiLC6LZ9+vQhjUZDoiiSh4cHnTp1yvwDr1njHArlnDkW95kzcvUqB+g4Q5cKAtE33yjXF3v27CEvLy+qXLmy3jW8fft2EgSBAJAgCHT79m3lhMyEnD59mjp37kLAXQJkxa+z9JpGwymOVGwPlBZARSVLcfiwczy1UzdJ4pBSM9BqtXTw4EHq27cv+fn5EQAKCAigiRMn0vXr183uiv379xMAypkzJx07dszyvly82PYJUi1p06dbLrMTs3278kZzQeCIXUdVb0mPixcvUuHChemtt96ic+fO0ZkzZ8jNzS1JoZQkiUarWa9NkpiYSFu2bKG6desSACpWrBh98MEhEkXnVSgBoj//VLrnsiZQWgAVlSyDVutc/n+G2q5dBkWXZZnCwsJo+PDhVKhQIQJARYoUodGjR9O5c+esqvSi1Wpp8ODBdOTIEev7dN06+5dzStl0U+3z51svsxOzZo2WBEFLSliQBIGoVSvH+U2a4sGDBxQUFETu7u7k7e1NkiRRSp9fb29vevPmjdJiOiXR0dE0e/ZsKlmyJAGg2rVr08aNGykxMZGePeMJEWd8r9ZoiKpUUf6FJquiRnmrqNiKP/8EGjdWWor0kSSgUSNOf/MfV69exdq1a7F27Vpcu3YNefPmRceOHdGpUye88847EJ0hH8mZM0CXLsDff/NzwV6IImfhX7UKePdd+x1HQUaMGIHvv78OSdoCQHRoLvlPPuE0Lw7MYmWSR48eoUyZMoiOjja4fMWKFeimy55tIU+fJqddu3mTE3UnJACuroCHByd0r1YNqFIF8PTMyFk4jvv372POnDlYtGgRYmJi0KFDBwwdOhQ1atTQWy8khFMaORsaDafDS52WSMVGKK3RqqhkGVq1cq7Ih3TMRP8cPkwzZsygqlWrEgDy8vKi7t270+7du41GvypKXBzRl1+azrRvTdPtb/BgxUtX2pOZM2cSAJo9ezadOUNUvrz9rUgaDSc7mD3b+axCCQkJ1KhRozSWSV0TRZGqVq1q0T5Pnybq1YvI31//8tJokpOXp/z83y1JJUtyXemUef6didDQUPr4449JkiTy9vamkSNH0j0TNXXbt3euyRpBIJoyxUEdlk1RFUoVFVtw755zzvGkalpBoG8BcnV1pfbt29OmTZsy17ReWBh71OuezLZQJBs3Jjp0SOkzsyurV68mADRmzJik7+LiiL76yj46uq7Vrk1044aCJ54OsixT3759k3wmjTVTwWRv3hAtX05UtWqyEm3t5SgIRM2aEe3YobwCnpiYSBs3bqRatWoRACpZsiTNmTOHoqOjzdr+yROe+naG+ETdVLezvi9nFVSFUkXFFkye7Fyv40ZarLs7RWXivIpExAr8V18R+fnxeQmC6Se5rsYdQOTlRTRiBJEFQUaZld27d5NGo6Hu3bsb9IUNC2Od2hY6uihqCeBygPPnK68UpceMGTNMKpIAkvotPQ4cICpSRHfutrlFdb/BO+9k/PKMjY2lrl270rZt28zeJioqin788UcqVqwYAaC6devSli1bKNGKUq6HDhHlyKHsu7ZGw5Vv79yxWHwVC1EVShUVW9C8uW1GzSVLiGJjOetz6mW6qNOWLTN+nKyiSMXFER0/TjR3LlHPnkTlyqVVLEWRqFQpoi5diH78kSPxM5NVNgOcPn2aPD09qVmzZhRvIhrmxg3Wsb290+rfxpqLS7JO//bbdwloRjdv3nbMCVpBTEwMiaJolnUSAOXIkYP+/fffVPsg+uwzfQXQHoqQqyvRTz9Zp5jLskw9evQgnXVRa2Int2/fpmHDhlGuXLlIo9HQJ598QqdPn7b8wKnYuVN/it/RymTu3ETh4Rk+DRUzUBVKFRVbkCePbUbAXLn4VTosTF8xevttVjSXLrXNcdatU7rH7EtCAj/1nSWkWAFu3LhB+fLlo+rVq1OMBb6hr1+zEjBpEnsXFChg+BJycyOqVYtoyBCiVauI7t5lZc3b25vGjh1rxzPLOIcOHaLRo0dT1apVSRRF0lkjDSmUgiDQd999l7TtpUtERYs6VkGqXZsolU5rkh9++EHvPHbs2GFwvWPHjlGHDh1IFEXy9fWlsWPHmswxayl79/L14shJHI2GqGBBzsGq4hhUhVJFJaP8849tR8L332eTxIQJySPj2bP8xPbyyvj+XVyIRo5UutdU7MijR4+oZMmSVLp0aXpiRUL71Lx4QRQRwVbMu3fZPy69GdDBgwdT3rx5KTY2NsPHdQRRUVG0Y8cOGjFiBFWqVCnJcplSwcydOzclJiZSaChbcB3t3SJJRGXK8FBjDjt37tSzwEqSRO+9917S8oSEBFq/fj3VrFmTAFDp0qVp3rx5Fr14WMrZs0QBAY6b/m7cmK9ZFcehKpQqKhnl999tPxrOm8fWtapV2T+TiKhhQ9vtv149RbtMxX68fPmSgoKCqECBAnRLgbDh8PBwAkBr1qxx+LFtwfPnz+n333+noUOHUoUKFZKUsp9+OkQeHsq5Sms0HA3++LFx+cPDw8nDw8PglP7h/zI8FClShABQ/fr1adu2bSanw22FPZM16PrIw4M9h6xInauSQdQ8lCoqGWXpUqB3b9vu090duHAhOTfiL78AAwbYbv9lygBXr9pufypOQXx8PD744AMcP34chw4dQuXKlRWRo0GDBkhISMDhw4cVOb4tefbsGZYtO4ipU9vi5UvBobk7UyNJQPnywIkTPESkJjIyEkFBQbh//z60qQQVBAGSJEEQBHTq1AlDhw5V7PoICwMGDgRCQzk3ZGJixvYnSYAsA61bA7Nn85Cp4nicIGuxikom580bQBBsu8/Xr4EvvwSKF+cMySNH2nb/sbG23Z+K4siyjF69emH//v3YunWrYsoCAAwYMABHjhzBxYsXFZPBVvj6+mHLlnaIjlZWmQQArRa4fBn44ou0yxITE9G+fXuDyiQAEBFkWcapU6ewYsUKRa+PoCBWis+cAXr2BNzceAi1JOm9bl1fX2DUKODWLWDLFlWZVBJVoVRRySj2MPILAvDZZ/wEyZcPqFTJtvtXJyYyHc+fP8esWbMQFxdncPno0aOxZs0arFq1Cg0aNHCwdPq0bt0a/v7+WLBggaJy2IK5c4FjxzJuRbMVsgzMmgWkNv4OHToUBw4cMKhMpmTTpk12lM4yqlQBFi0CHj4ElizhiZ4qVYwrlu7uQJ06PDyGhPC2U6cCxYo5TGyV9FB4yl1FJfOzeLHtnYF0QTMdO3LOi7//5jBJW+2/TBll+0zFYqZNm0YA6N1336UXL17oLdNF9M6aNUsh6dLy9ddfk6enJ0VFRSktitVcv86pe5TwmTTWJImoWDGiV69YzoULF6bxl0yv+fj40OvXr5XtWBPExxOdP895Pv/4g2j3bs5pefWq8+Y2VSE1KEdFJcNs3mzbp0VAQHL5DYCoZk0Oqf3+e9sdo1YtpXtNxULq1q1LAEfsli9fnv75L+R3zZo1BIBGjRqlsIT63L9/nyRJonnz5iktitU0b+681VRFkRNBrF+/3uycmrr2yy+/KN21KlkQNShHRSWj3LnDvo62QJKA48cBf3/g7beBqCj+/rvvgGHDgLp1ef4tI7i4cIDPTz9lXF4VhxATEwNfX18k/jfvqtFokD9/fkyaNAl9+/ZFp06dsHz5cgi29uXNIB9++CGuXbuGCxcuOJ1sprh1CyhVitU3ZyVPHiBnzpK4f/+WweWSJEEUkz3bZFmGVqtFQEAArly54igxVbIJGqUFUFHJ9BQtCnh5AS9fZnxfY8cC1asDTZsmK5MA8NVXwAcfcER55coZC6pJSGCveJVMw8GDB5OUSYADMB49eoRevXohODgYixcvdkqFrX///mjUqBGOHDmCd999V2lxLGLhQk6yoHQgjjH+/RdYvvwiAgOvIi4uDjExMQbbq1ev9D4XUx0OVeyAaqFUUbEFjRoBf/2VMXNGYCBw6hSwfDnQt2/a5TVrAkePskf+8OHWHwfgUNHy5TO2DxWHMXjwYCxYsEBPqdSRI0cObN68GS1atFBAMuPIsoyAgABUrVoVa9euVVocs4mNBQoU0H+nc0YkCahVCzh0SGlJVFTUKG8VFdsQHMyje0a4cAFwdTWsTALJSdsyqkx6eABly2ZsHyoOZfv27QaVSQBISEhAq1atsGzZMgdLZRpRFNGvXz9s2rQJjx8/Vlocs9mxwzJlcvNmzvTl7Z3+OqtXA/HxnLQhVy5g8mROBfvqFRARAWzYYPk7nlbL0d537li2nYqKPVAVShUVW9Cli/PkFTGGRgN0755x5VfFJLLMykJUFFu8rDVe37lzB7dv3053OREl5aCcOnUqnG3SqUePHtBoNFiyZInSophNaKhlORGXLAFy5gQ6dza83MsLaNsW2L4dePIE2LYNGDKE6xW0aAGMGcOeLMePA0WKWC7vyZOWb6OiYmtUhVJFxRaULQvUr+/8ilpiom0r7qgAAF68AHbtAqZMYcWhYEG+FDw9AR8fVjZcXIAKFTjX3vz5rASY45+3Z88es/0jv/jiC6dLJu7r64tOnTph4cKFJnMkOgsnT7Krsbn88Qfwzz9Ar16Gl3fqxPkTlywBSpYE6tUDfvgBmDkTOHCArZfdurHi2a6dZbK6uHDlGRUVpVEVShUVWzF4sHN78EsSZwSuUEFpSbIMp09zpY/8+YFmzYDx44Hff+dky6nRaoErV4CVK4FBg9gltkgRTspsbDZ49+7depG6qZH+e4l59913sWbNGlSsWDGjp2Vz+vfvj3v37mHnzp1630dFRWH58uW4dctwlLISEFmuoMkysGIFUK0aJ2dITc+ewIMHrHjqFNXUU+ovXvBfS+PtEhJUC6WKk6BkziIVlSxFQgJRwYKcIE7pJHXptQ0blO6lLMHmzURVq3KXZjRPoSjyPjp14sTNKUlISCBPT880eQQlSSIA5OfnR6NHj6Zr164p0xEWUL16dWratCkREZ09e5b69OlDbm5uBIC+/vprhaVL5uZN637HkiU56fYPP+h/HxDA+506Nfm7LVuIIiKI3nuPyMODqGxZoj17iO7cIfLxsfzYuXIp22cqKkREqoVSRcVWaDSca0SWlZYkLRoNT8m3b6+0JJmaJ0+ADz/kaclz5/i7jLrOyjLvIyQEqFgR+PHHZEP3qVOnEBMTk7SuJEkQBAHvv/8+Nm3ahIcPH+Lbb79F6dKlMyaEA+jTpw927dqFSpUqoUqVKli2bBliY2MhCAK8jUWzOJh//7Vuu5s3Odq6Sxe+3XTopsGXLk3+rkMHDvzZvx+IiQH+/puDderVS7ZUWkJ0tHMOOyrZC1WhVFGxJS1bAl27Op8vpYsLsGwZ1whXsYoNG4By5YDffuPPtn6AJyZyFPCwYeyZcO0aMHfu3KTl/v7++PLLL3Hnzh3s3r0b7dq1g4slkSMKcfv2bYwZMwajRo0CgCQfT13UOhEhd+7cismXmoykeF2yBMibF2jVij9LEiuYhw4BN24kr7dgAb+YDBnCtQo6duTfft8+64JyMiq3iootUBVKFRVbM2sW4OfHWZGdhR9/5ATsKhZDBIwcCXz0EVuPHOEme/o0R/2+eVML5cqVw/bt23H//n1MmDABRazVOBzM3r170bRpU5QsWRIzZ87Ei/9Mb2QgCt3Pz8/R4qVLRt65Nm7ka6RnT/7cvDnns0wZ4N6kCfDpp5wdbNYsTvsTEsKpbHPnBiZMcLzcKiq2QK2Uo6Jia3x9OWyzWTOlJWGltkUL4H//U1qSTIksA/36cXoXwHFl+BIT+djbtg3Eli0D4YQ5y40SFxeHDz74ALH/mc1MRXc7k4XSzc36bWNjgbVrgT59WJHs1YsLaIWEJK9TuTL/PXVKf9uoKLZiGgrqMQdXV+u2U1GxFU5kQlFRyUI0asRKpZJmA0kCatfmuVrVfGExRDz9rFMmHY0sszW0XTv2tctMuLq6IiQkBG5ubklR6MZwJgtl3rwZ237JEvahHDmSLZTr1gFv3iQvf/CA/wYH62+XOzdQpgwnObcUb2/nmhBRyZ6opRdVVOzJ6tWcSFwXkOkoRBF4913OoJwrl+OOm4WYPx8YOFBpKfindHMDLl4ESpRQWhrLCAsLQ9OmTfH8+XOjVsrHjx8jX758DpQsfYhYQYuOtn4f585xgJUocnqolGl9PDyAS5c45+SUKcCZM4C/PyugFSsCDRsCBw+afyxBABo0APbutV5eFRVboL7TqKjYky5dgC1bOKuxxgEeJjpL5IcfcqZtVZm0ilu3Ml7h0lbIMgds9OiR+SJ5g4KCcPr0aZQqVcqopdLX19eBUhlHEIDq1TO2jyVLWJm8fDltjshXr9g6+euv7E6xcycwYwYnRq9TxzJlEuBhJaPyqqjYAtVCqaLiCO7fZ0/8PXvsdwxJYgVywQKOIFGnua1Cljl9y4kTzldNc+5c57CaWsqLFy/Qpk0bHD58GHIqrdjDw0MvNZIzMHYsV7Fxtt8/PUJC1IxgKsqjWihVVBxB4cJsMVyyhOe8bOnwpLN8tmzJCe0+/lhVJjPAggXAkSPOqUyMGMHW08yGj48P9uzZg84Gil37+PgoIJFxgoOd8/dPj5o1lZZARUVVKFVUHIcgcNjnrVvA5Mlc8BmwPmelKPK27drxPNmWLVwDUMVq4uKAr75SWor0SUwEpk9XWgrryJEjB1auXInx48frfe8sATkJCQk4fvw4iAjNmnGQjLMjSew/Wbiw0pKoqKgKpYqK48mXj+fU7t3jws9NmgCensnLXVz0LZiCwFZIndVRFIGyZYGJEzkkdP16zo6sWiUzzKZNwPPnSkuRPomJXAs8dR3ozIIgCJgwYQKWL1+eVJ/ckG9lTAxw9CjHlIWEcH7HHTs4P6e9EnhPmDABtWrVwqZNm5Ajx//bu/fwmO78D+Dvc2YGqfvjNjTVoK0Gi22oujWsWlYtVaUVdIXuWhbttrTd2j5auq3d7bqVtM2aaCpCy1p+olRZpGiCUbShdalrioYmDRJJ5vL749Mhkplk7ueQ9+t55olkzsl8JWnz9r18PsCECfrrT1CW3Q5Mnqz1KIgE91AS6YHTKb3brFZ5nDkjvzltNiAiQo6Etm8PdOoEdOggy+YUdF27yiEKPR9+URQpiH2rB4nNmzejX79+aNu2LT777CA++ki+9hkZMonv6TeTwQBER8syb7du0sYw0LNnhw8fRvv27WGz2dC4cWMcPXoUeXl1EBUV3uIMvjKb5X8V4TjvR1QZBkoiIgAHD0pW1ztFAVq1ktaMt/KktNMJpKaewkcfNcSGDTXhcMjku7d7F41GmaGrUUM600yY4F9RcKfTidjYWOzatQt2ux2qqmLixIl455138PjjwLp1Tths+vtCK4qUHfrLX7QeCZHgkjcREaRHt7dLnKtXAwUFUq/Qk5QUKffjKq/YoAEwbx5w4oRMPp8/LyVjfK2Y43RKR5WjR327T0927pSai6NG3Y0NG2pe7wzky0EYm02+FoWFQGKifL7evSVo+yIlJQWff/759TqZDocDixYtwq5duxAVNQ82WwEAfc27GAzAvfdK4X0ivWCgJCKC7M/zdr3GYpGdCG4OLQOQHQpDhgBpacAPP0jh6sxMoH9/YNYsaaQ0YYIEw2rV/Buv1erffVoqKJAQ1LOnFCQAgnOa2vU5duyQYDl3rnc913Nzc/Hss89CKTPVqygKevfujQULpqJ//00A9DVD6XAAS5ey3SLpCwMlEREk8Hm7d3LDBilEPXas++dHjJBa9haLvJ+QIL/8H3oISEoCPv9cDuVPmQJcuOD7WE2mWy9QZmbKkvT8+RLcvQl8vrLZZFb4+eel6+jx4xVf//LLLyM/Px9ld345HA4UFxdj6tSp+OSTIejTRz/7FBUFeOkl4MEHtR4J0c24h5KIqrwffvC94tLf/ga88orMiH399c3PZWRIKRfX47vvgNdfl4P5wdKzJ5CeHrzPF0rr10t1K7s9NEHSHYNBtiT873/u98ZmZGSgW7du5cJkaREREThy5Aiczkh06gRcuhS+8btjMAAxMVIlrEYN7cZB5A5nKImoyjt82Pd7kpJkRrPsLKXrBHJysjzfs6ccNvn+eyA1VXpEFxYCW7fKjKW/srL8vzec1qwBBg8GSkrCG8bsdimv1LOn9MsuzWaz4ZlnnrleusiT4uJiTJ48GXfdJd+v2rW1lfBT+gAAHEhJREFUKyVkNEq1sI0bGSZJnxgoiajKu3rV93uOH5cZwlGjbl4OdQXMpCR5e+ed8vbtt2Xf5dChsveyfn2ZPfvFL/wbc2Ghf/eF0+bNUtbH4dCm/I7dLvs2H3nkxp5NAFi4cCGysrKuH8TxfL8da9aswbp169CmjcwM1qsX/uVvgwFo2xbYts33Q1xE4cJASURVXlGRf/dZLECjRsCgQfK+wSABMz1dDtwAN2rUnz0rYXLTJtk/2b+/BK0XX/TvtYuL/bsvXM6fl/7SWoVJF7sdyM8HHntMvs9nz57FSy+95NW9rqLrL7/8MgApBZuRIWWbwlmy6ZFHJMw2ahS+1yTylU62GRMRacdk8u++VauAd96ROoirVwMDBkix6dJ55dIlebt5882Hfs6fBw4cAB54ILxjDgenE/jDH6TjjR6KxNvtUk5o5kzgk09+i2IPaVxVVdSrVw9msxnNmjVDs2bNYDab0a1bt+vX3HOPfN9mzgRmz5ZgGYqlfKNRDnItWCA/X7dyzVGqGhgoiajKi4jw775r14Dly4Hf/16C5NixMhu2cuWNaw4e9Hy/ovgfuPRcMmb5cmmbqCdOpwTA0aPH4c47N+Lhhx+G2WyG2WxGkyZNYDab0bBhQ7etIMuqXl0OZQ0ZIjPS334rs9PBCJaqaofDYcCvfgUsXsw+3XTr4ClvIqryTp0CoqL8uzcmRmpYzpkDTJoEfPABMH78jecVRT7/lStSNscVIJs2lWXx1FQJpP6+rt5cuiRLwvn5+mtbaDDIDGNWVvAO19hsEp4XLJA9jv4GS1WVR3T01/jmmyk4dSoVTZuagzNIojDgHkoiqvKaN6+4601FrFZZAn3uOSlS7qo96eJ0An/+s5zQXbtWlsWHDQM+/VT2Qb71lu+vaTLJSXI9sljkJLvewiQgQe/bb6WOaLAYjTJTuXWrHPyZPFm+1669s4oi3y+jUR4m083bFUwm2fYwc6bss01Pj0REhBVz5vwreIMkCgPOUBIRQbrXbNniXxCaPFlmqLKyPPeTHjQImD5dTnUXFckhi7/8xb+SRYAEN0+F1bVit8tM79mzWo/EM4MB6NNHAn0oXb0q/9CwWmUm+to1+QdE9epAzZpAmzYyy9y2bfn9sNOnT8f8+fNx6tQpNGjQILQDJQoSBkoiIkiR8n/+MzitAMNh/373Bbu1tH49MHCg1qPwzrFjsjSvRzk5OYiKisLUqVPx+uuvaz0cIq9wyZuICMCvf33rhMlGjWRmS28SErQr/O0LgwFITNR6FJ41atQI48ePx4IFC5Cfn6/1cIi8wkBJRAQgNha49179l2dRVWDiRP30lnZxOORQirsDKatXS4HxivappqTIknDjxsDo0XJS/Jtv5POdOOH+niVLZIuCp4enfaZ2O/DZZz7/FcNq6tSpKCgoQEJCwvWP2Ww2pKSkoF27dpg1a5aGoyMqj4GSiAgSJCdP1noU3vHnVHioHT0qodEdi0VKM8XFuX++Th052JKWJn3VR4+WGdjdu6UjkSezZkn7yrKPnBzZx7lnj+d7v/5a2kHqVbNmzRAfH485c+YgNzcXiYmJaNWqFUaPHo2srCwc9nfzLVGIMFASEf3s6af1Xd/RYJC+2K52jnpitXp+bsMGIDvb8yGiESOAO+64cUK+Xz/pSvP003Iq25PvvgMyM29+VK8uWwKWLKm4xmdJif77oU+ZMgUXL17E3XffjfHjx+PMmTMApIOPUW9T1FTlMVASEf2sbl3g2WdvlHzRG4dDTobrkdXquXuPwwEkJwOdOrk/BR8fD3z//Y1yPoEcFR03Tl7P1Uu9sjHrUX5+PmbPno0ePXrA6XTi8uXLAADXGVpFURgoSXd0+r9NIiJtzJgBtGihv8MlqgpMmwZ07qz1SNw7caLiQ01JSRL0ys5SRkfLXsfk5MDbNNapI/3Dt2wBTp6s+FqTqfJrwi03NxczZsxAZGQkpk+fjtzcXI/XMlCS3jBQEhGVEhEhB0T00IPaxWAAWrYE9FxBprCw4pnF48eB9HRpVVg6C7kCpjczipUpu3RemcLCwF8zmF588UXMnDkTly9fhqOCH0BFUbxqEUkUTgyURERlPPQQMHWqfk58OxwScmvU0HoknnlzwMVikf2NgwbJ+waDBMz0dKkLGahx44CLF4H//te76/V2KGf69Ono0KEDVC/2XHCGkvSGgZKIyI033pDalHrYT/n++/pttegSEVH5NatWAXl5smcSkDaUZrP3M4oV+cUvZDuAq/yQN/QW0KOiopCRkYFx48ZVei0DJemNDv5XSUSkP9WqSf3E7t21DZVvv63PMkFl1a5d+b7Ta9ekvmT//hIkx44F8vOBlSsDf31XBlu82LvrnU4Zs97UqFEDiYmJWLJkCapVq+ZxaZuBkvSGgZKIyIM77pCez/36hXf5W1Xl9RISgBdeCN/rBsJTD/OyLBbZQzltmsxQrlgR+F7GatVk6Twz0/tSQDab92PWwpgxY5CZmYnIyEi3oZKBkvSGgZKIqAIREcDatcD06RL0Qn0WwmAAGjaUvtgTJoT2tYIpJsZ9l5yyrFbgwAHgueckCLpb7o6OBoYOlYfZLMHe9X50dPnrH3sMaNDA+9nJ0mPWs44dO2L//v34zW9+c9PHHQ4HD+WQ7jBQEhFVwmSSriy7d4euPaNrWX3UKCnmXSZD6J4v4cxikb9vVpZ8TcsaPlz2W65aJbUrGze+8f7w4eWvHzcOuHJFZju9Vb8+EBnp/fVOpxM2DZq916tXD2vXrsXs2bOh/PyDZ7fbPQbKggIgI0O6Dq1cKV+ztDSZvfXUyYgoGBSnM5ASskREVUtRETB7NjBnjuz/U9XASgwZjbL8ev/9wL/+JcvAt6o775QC5XqnKEDfvrKdwcXpdCInJwcnT568/jh16hROnDiBY8eO4cyZM6hWrRouXrwIk6cK7iG2detWDB48GJcvX8bYsWNhsVhw6RLw8ccSGDMypAWmp59HVQXuu+9Gi8rhwyVYEwUDAyURkR8KC2VGbP58WcJ1TRh5s+xrMknJGqNRlnH/9CegRw/9lCny1wsvyNfDm6+B1iZMsCIr63koioIzZ84gOzsbRUVF1583GAxQVRUlpWoL3XPPPThy5Mj1mUItnD17Fg8//DCGDfsHfvjhCaSmys+SwVBxYfnSjEb5HlWrBowcKT9/DzwQ2nHT7Y+BkogoQFYr8L//yduMDODUKffXmUxyEKRLF1kiHjhQ9gjeLo4elRkwvatVCxgy5I9YuvR9r+8xGAyYMWMGXn311RCOrHJWKzB+vLx1zW4HwvU5HnwQSEwEOnQIzjip6mGgJCIKsrw8CZWFhTJ7VKOGhJh77vHc7/p20bcvsHWrfmcpjUZg0iRgzhwn4uPjkZyc7PW9x44dQ6tWrUI4Os+KiqQ26ptvykx2sL++rhn2V18FXnnl9v85peBjoCQioqBZu1ZOXevZN98ArVvLaen4+HgsXboUFf0qVBQFMTEx2LNnTxhHecOBA8BTT8lhrVD/xlYUoG1bqReq57JKpD885U1EREHz298CXbve3K9bLwwGORHeurW8r6oqkpKS8MQTT1S4L9LpdKJu3bo4ePBghcEzFLZvB7p1k+0E4XhppxM4fFi+hzt3hv716PbBGUoiIgqqo0elFWKpMy6aU1WgSROZnaxT5+bnSkpK8MQTTyAtLQ0ON0ekFUVBvXr1kJubi7Zt22LkyJEYMWIEoqKiQjrmLVvk1L/NFlglAX+oqix7b9oEPPxweF+bbk2coSQioqC6914praQnDgfwwQflwyQAmEwmfPzxx+jTpw/UMn02DQYD+vTpgwsXLiAtLQ3t27fHG2+8gRYtWqBHjx5ISEjAxYsXgz7ezEw5tKVFmATkNUtKpB6q1Rr+16dbD2coiYgo6BwOoFcvYNcu7Q/oKArwhz8A771X8XWFhYXo378/du7cCXupQX/44YcYPXr09fevXLmCtWvXIjU1FZ9++ikURUG/fv0QFxeHwYMHo2bNmgGNNy9P6pLm5GgTJkszGKQSweHD+ux9TvrBQElERCGRmwt07y5L4Bo0mQEgS7f9+wNr1nh3cvnKlSt45JFHsHfvXtjtdlSvXh05OTmo7SFN5eTkYOXKlVi2bBl27dqFmjVr4rHHHkNcXBz69u3rVxH0+Hhg6VLtg7iLqgK//33lgZyqNgZKIiIKmfPnZaby2LHwByRVBX71K+D//k96snvrp59+Qq9evbB//34MHz4cH330kVf3nThxAsuXL8eyZctw6NAhNGzYEE8++STi4uLQtWtXrwqib9ig325JmzcDffpoPQrSKwZKIiIKqUuXZJbQag3PSWWXoUOB1FTpCOOrH3/8ERMnTsS0adMQ40ujcsip8AMHDmDZsmVYvnw5srOzERUVhbi4OIwcORJt2rRxe9/Vq0CrVvpY6i5LVYGmTWW22ZdwTlUHAyUREYVccTEwa1boCnO7GI2ytD1vnizTat3O0uFwID09HcuWLcOqVauQl5eHjh07Ii4uDiNGjEBkZOT1axMTgT/+Mbyh21cffAD87ndaj4L0iIGSiIjCZt8+YNQoKd8DBC88qarM6vXuDSxZAtx9d3A+bzAVFRVhw4YNSE1Nxbp161BUVITY2FjExcVh6NAnEBtbH1lZ+g2Uqgp07MhT3+QeAyUREYVVcTGQkADMnw+cPBlYT2rXvb/8JfDCC0BcnPazkt7Iz8/H6tWrkZqaii1btkBVu8NmS9d6WF7Zswfo1EnrUZDeMFASEZEmHA4p3r1wIZCWJu+rasVL4kajPOd0yt7IkSOBiRNv7YBz7tw5DBqUi337WsPhMGg9nAoZjbLkvXix1iMhvWGgJCIiza1btwMWy5fo1m0yrFZg924pO1RUJAGzRg2gUSPgoYeAmBh5dOwIBFjyUTeaNpUT8WWtXi0Hmpo2BX76yf29KSnA8OFAZCTQr5+cEv/lL6XA/OnTQIsW5e/p0AH429+ko1GjRkBhofQKX7QIWLas4rE2bw6cOuX735FubwyURESkqezsbNx3330oKChAfn6+x5qPt6ucHKBxY/fPPfqozN5OnAi8+2755+vUAc6dAz79FHj8cWmVaDYD+/dL+DaZ3AfK2FjgqaeAHTuA7GwJ5iNHAiNGAH/9q4TNiuTmAvXq+f53pdsXAyUREWnmypUr6NatG7766isAwObNm9GnihU73LhRWhy6o6oyy3juHNC5c/nnx4+XguMDBwLr18tsruu3+rp1QLt27gOlJ198ATRrVvmhpi1bpMYnkQt7eRMRkSbsdjtGjBiBrKwsAICqqvjiiy80HlX4Wa3S4tAdhwNITpY9ou3alX8+Ph74/nspiA4EfkL84sXKD0ipKrB3b2CvQ7cfBkoiItLE1KlTsX79ejh+ruLtdDqxc+dOjUcVfidOVHwyPSlJguXYsTd/PDoa6NJFAqe/hdAVRcJsw4bAhAmyB/Pvf6/4HoOBeyipPAZKIiIKu3fffRfz5s1D6V1XTqcTu3btQlXbiVVYWPHM4vHjQHq61O80Gm983BUwk5L8f+2EBJmRzMkB5s4FpkyRAusVcTplzESlMVASEVFYbdy4EZMmTXL7XH5+Po4cORLmEWnLZqt8qdpikdPYgwbJ+waDBMz0dOmT7q8335Tl9AEDJJguXCj1PCvidAIlJf6/Jt2eGCiJiChsvvrqKwwdOrTCWchdu3aFcUTaq1698mLsq1YBeXmyZxKQAGg2S9AMxJkzsodzwwY5SZ6YCLz1liyBe+Iq40RUGgMlERGFxblz59CvXz8UFRV5DJRGo7HKHcypU0cOulTk2jVg+XKpSWk2y3J3fj6wcmVwx7J7t5QaatnS8zWKAlSxyk7kBQZKIiIKuatXr2LAgAHIycmB3VMbHAA2mw3p6bdGC8JgadfOu9aTFovsoZw2TWYoV6wI/l7G3r2lE9F333m+pqTE/YlzqtqMlV9CREQUmJSUFOzfvx9qZVNxAI4cOYKffvoJdevWDcPItBcT4125H6sVOHAAeO45mdF0t9wdHQ20aSN/NpuBO+4Ahg6V9w8dAg4flj+//77McO7eDVy4IEvcw4ZJsfN//EPKB1U2ZqLSWNiciIhCrrCwEB9//DG2bduGzZs34+zZsxVev2nTJvTt2zdMo9PWtWtArVqe+5eXNnkysGABkJXlfpZwxgzgtdfc3/vaa8Drr8ufx4yR/ZjR0dLx5soVCauLF1feetFkAq5elbdELgyUREQUdqdPn8bcuXMxb948NG3aFOfOnQMgeyhtNhtee+01zJgxQ+NRhk/79sDPzYJ074EHZLaUqDTuoSQiorBr3rw5FEVB8+bNkZ2djdOnT2Pp0qV4+umn0bJlyyqz3O3y5JOVH8zRA1WVZXGisjhDSUREmoiJiUG7du2QnJys9VA0d/48cNdd3h3O0ZLJJK0eKyorRFXTLfDvISIiut3k5eXhyy+/RGxsrNZD0QWzWQ7PGHV8VNZolNlJhklyh4GSiIjC7vPPP4fT6USvXr20Hopu/OlP+p6htNlkjETuMFASEVHYbd++HXfddRdatGih9VB0o0cP4NFH9TlLaTAAQ4YAXbpoPRLSKwZKIiIKu23btiE2NhZKZT0HqxBFkdaHERGVt2IMJ1dnnPfe03okpGcMlEREFFau/ZNc7i6vWTNg0SLvCp2Hi9MphdAbN9Z6JKRnDJRERBRWO3bsgMPhYKD0YNQoYPBgfZQRUlXpoDNsmNYjIb3T4U4NIiK6lV24ABw5AhQUAMXFQLVq0gKwdWuZ5dq+fTsiIyPRsmVLrYeqS4oi3Wr69AH27vWug04oGAxA165AcrK+luBJnxgoiYgoIHv3Ap98Im8zM4EffvB8bZMmQFHRQDRv3gX79insCe1BzZrAxo1A377Al1+GP1SqKtC5M7B+vezpJKoMC5sTEZHPCguBFSukr/T+/TKb5XQCDkfl9yqKA6qqwG5X8MADwJQpwPDhDC7uXLkCDBoEbNsW3n2VffsCa9bIzDKRNxgoiYjIazYb8PbbwFtvAfn5MpPlTYj0xHV/3brA9OnAn/+sz7I5WrLZgDlzgL/+VUJlqGpVGo3y/Zg9W0K+wRCa16HbEwMlERF5JSsLGD1aZiRD9ZujUydg6VLg/vtD8/lvZYcPy9ffapU9jcH6Hrg+V5cuwIcfAvfdF5zPS1WLDs6QERGRnjkcMmvVsSNw8GBol1737wfatwf++U99lc7Rg+hoICND6kG6Ql8gs7mue6OjgX//G9i5k2GS/McZSiIi8shmA8aMkVPH4TZmjAQdLoGX53QCO3ZIzcr//Ee+T6oqs42eDvCU3udqNEopoEmT5CQ3T3FToBgoiYjILZtNDsusWaPNbKGiSOhJTeV+vopcviwnwffuleXwzEzgxx+Ba9fka1i9OtCggSxpd+oExMTIbHPt2lqPnG4nDJRERFSOwwHEx8t+Ri1/SyiKjGPxYs6iEekZ91ASEVE5CQlyQEPrKQenE0hKkh7XRKRfnKEkIqKbHD8OtGsnS6Z6EREBHDoEREVpPRIicoczlEREdJ3DAfzud6GrdeivkhJZ+uYUCJE+MVASEdF1778v5WP0FihtNukWY7FoPRIicodL3kREBEDKzURFAWfPaj0S9xRFxnfsmJTIISL94H+SREQEANi4Ub9hEpDl7hMngC1btB4JEZXFQElERACAhQv1X+/RaJRxEpG+cMmbiIhw8iTQsuWtcehFUYDTp4HISK1HQkQunKEkIiJs2+ZbmFy9GigoAOrW9XxNSgpQXAw0bgzUqgXMny9L6teuAd9+C0yb5t9eSKcT2L7d9/uIKHQYKImICFYrYDJ5f73FIrUh4+LcP1+nDjBkCJCWBly6BHz2GTBqFPDmm8DAgcC6dcDs2cDcub6P1WSS8RKRfnDJm4iI8OCDwJ493l+vqrLsfO4c0Llz+efHjwfee0/CY61awIoVwOOPA//9741r3nsPeOYZoE0b4MgR38bbvTuwY4dv9xBR6HCGkoioirPZgIMHfbvH4QCSk4FOnaSrTlnx8cD33wMbNkj4czjkz6WlpckhoCFDfB/zvn3yOYlIHxgoiYiquAsXgKIi3+9LSpJQN3bszR+Pjga6dJHA6XAA1arJ25KSm69zvWb79r6/dmEhcPGi7/cRUWgwUBIRVXH+9uw+fhxIT5e9kUbjjY+7AmZSkrw9dEief+ihm+/v0UPeNmjg3+vrqdc4UVXHQElEVMWVnTn0hcUCNGoEDBok7xsMEjDT06WjDQAsWyYHcxITZa9m3brAU08BU6bI8/4uXRcX+z9uIgouBkoioiquenX/7121CsjLkz2TADBgAGA239xz+9IloH9/+XNmplz/zjvA88/Lx7Kz/XvtGjX8HzcRBRcDJRFRFVerlv/3XrsGLF8ugdFsluXu/Hxg5cqbr9u7F2jbVnpxt20LNGsGHD4sz6Wnh3/cRBRcDJRERFVcw4ZA/fr+32+xyB7JadNkhnLFCjk0486pU7KnsqQEeOEFmZ0sGz690agRUK+e/2MmouAyVn4JERHdzhRFakl+9pl/rRetVuDAAeC556Q+Zenlbpc33gC++krqVjZvLjOZXboAjz7q++EaVZV7iUg/OENJRETo3FkO1PjLYpGgl5UF7N5d/vn69YG//x3YtAmYM0f2VXbpIi0ffaWqUv+SiPSDnXKIiAhr1vhXYFwr69ZJFx4i0gcGSiIiwuXLQJMmnvc+6kmtWsD580DNmlqPhIhcuORNRESoXVtK/xh1vrPeaJT+3wyTRPrCGUoiIgIg+x/d9eXWm2++AVq31noURFQaZyiJiAiA1Ifs0SOwwzmhZDAAvXoxTBLpEWcoiYjoup07gZ49/SsfFGqKAmRkSPtGItIXzlASEdF13bvfqCepJ6oKvPgiwySRXnGGkoiIblJQIHspT58G7HatRyNL3S1bAgcPsn83kV7p7N+gRESktTvuAFJStB6FUBR5LF3KMEmkZwyURERUTrduEuIUReuRAKmpbLVIpHcMlERE5NaIEUBi4o1ZwnByvWZSEjBsWHhfm4h8xz2URERUoRUrgNGj5eR3OPZUGgwSJlesAIYODf3rEVHgOENJREQVeuopYM8e4P77Qz9TqShyIMhqZZgkupUwUBIRUaU6dgT27QNmzJAZxGC3aDQa5TFzpoTX9u2D+/mJKLS45E1ERD45eBB4801g1SpZBnc4/P9cqiqPYcOAV165NVo/ElF5DJREROSX8+cBiwVYtAg4d04+ZjIBJSWe7yn9fLNmwKRJwNixQJMmoR8vEYUOAyUREQXEbgeysmTfo9Uq7RG//RYoKpLwaDJJDcnWraX8T6dOQEwM0KaNfvuGE5FvGCiJiChknE591LIkotDioRwiIgoZhkmiqoGBkoiIiIgCwkBJRERERAFhoCQiIiKigDBQEhEREVFAGCiJiIiIKCAMlEREREQUEAZKIiIiIgoIAyURERERBYSBkoiIiIgCwkBJRERERAFhoCQiIiKigDBQEhEREVFAGCiJiIiIKCAMlEREREQUEAZKIiIiIgoIAyURERERBYSBkoiIiIgCwkBJRERERAFhoCQiIiKigDBQEhEREVFAGCiJiIiIKCAMlEREREQUEAZKIiIiIgoIAyURERERBYSBkoiIiIgCwkBJRERERAFhoCQiIiKigPw/fmgeTZ47+u4AAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "draw_graph(graph_true, highlight_nodes={\"X\", \"Y\"})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 255
    },
    "id": "xoFVii0MdmFC",
    "outputId": "019c2901-050b-4496-d202-7a676506b000"
   },
   "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>V1</th>\n",
       "      <th>V2</th>\n",
       "      <th>V3</th>\n",
       "      <th>V4</th>\n",
       "      <th>V5</th>\n",
       "      <th>V6</th>\n",
       "      <th>V7</th>\n",
       "      <th>V8</th>\n",
       "      <th>V9</th>\n",
       "      <th>V10</th>\n",
       "      <th>...</th>\n",
       "      <th>X</th>\n",
       "      <th>V15</th>\n",
       "      <th>V16</th>\n",
       "      <th>V17</th>\n",
       "      <th>V18</th>\n",
       "      <th>V19</th>\n",
       "      <th>V20</th>\n",
       "      <th>Y</th>\n",
       "      <th>V21</th>\n",
       "      <th>V22</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>-0.288362</td>\n",
       "      <td>-0.170262</td>\n",
       "      <td>0.036313</td>\n",
       "      <td>0.561619</td>\n",
       "      <td>-0.477592</td>\n",
       "      <td>0.026122</td>\n",
       "      <td>0.235572</td>\n",
       "      <td>0.193228</td>\n",
       "      <td>0.219485</td>\n",
       "      <td>0.224018</td>\n",
       "      <td>...</td>\n",
       "      <td>-0.334338</td>\n",
       "      <td>-0.274177</td>\n",
       "      <td>0.069914</td>\n",
       "      <td>-0.425716</td>\n",
       "      <td>-0.244209</td>\n",
       "      <td>0.489664</td>\n",
       "      <td>0.597804</td>\n",
       "      <td>-0.410179</td>\n",
       "      <td>0.024974</td>\n",
       "      <td>0.362765</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>0.172179</td>\n",
       "      <td>0.498972</td>\n",
       "      <td>-0.637497</td>\n",
       "      <td>-0.403630</td>\n",
       "      <td>-0.510466</td>\n",
       "      <td>0.803194</td>\n",
       "      <td>-0.629094</td>\n",
       "      <td>0.042610</td>\n",
       "      <td>-0.364143</td>\n",
       "      <td>0.421791</td>\n",
       "      <td>...</td>\n",
       "      <td>0.177786</td>\n",
       "      <td>0.222926</td>\n",
       "      <td>0.336022</td>\n",
       "      <td>1.071612</td>\n",
       "      <td>0.246272</td>\n",
       "      <td>-0.468573</td>\n",
       "      <td>-0.761310</td>\n",
       "      <td>-0.179248</td>\n",
       "      <td>1.007144</td>\n",
       "      <td>-1.160541</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>-0.148207</td>\n",
       "      <td>0.089956</td>\n",
       "      <td>-0.371467</td>\n",
       "      <td>-0.015495</td>\n",
       "      <td>-0.305202</td>\n",
       "      <td>-0.524397</td>\n",
       "      <td>-0.312337</td>\n",
       "      <td>0.700555</td>\n",
       "      <td>0.219698</td>\n",
       "      <td>0.031795</td>\n",
       "      <td>...</td>\n",
       "      <td>-0.008598</td>\n",
       "      <td>0.869035</td>\n",
       "      <td>0.767204</td>\n",
       "      <td>-0.080318</td>\n",
       "      <td>-0.410839</td>\n",
       "      <td>0.171747</td>\n",
       "      <td>-0.962758</td>\n",
       "      <td>0.183980</td>\n",
       "      <td>0.891160</td>\n",
       "      <td>-1.165926</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>0.392713</td>\n",
       "      <td>0.466647</td>\n",
       "      <td>-0.414478</td>\n",
       "      <td>0.146909</td>\n",
       "      <td>0.437273</td>\n",
       "      <td>0.467879</td>\n",
       "      <td>0.124036</td>\n",
       "      <td>-0.259591</td>\n",
       "      <td>0.100210</td>\n",
       "      <td>-0.001243</td>\n",
       "      <td>...</td>\n",
       "      <td>-0.220841</td>\n",
       "      <td>0.620429</td>\n",
       "      <td>0.497164</td>\n",
       "      <td>0.094684</td>\n",
       "      <td>-0.674582</td>\n",
       "      <td>0.418828</td>\n",
       "      <td>-0.232072</td>\n",
       "      <td>0.204818</td>\n",
       "      <td>-0.205904</td>\n",
       "      <td>0.029435</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>-0.382722</td>\n",
       "      <td>-0.243545</td>\n",
       "      <td>-0.381480</td>\n",
       "      <td>-0.179152</td>\n",
       "      <td>-0.413537</td>\n",
       "      <td>-0.195249</td>\n",
       "      <td>0.376594</td>\n",
       "      <td>0.571088</td>\n",
       "      <td>-0.329178</td>\n",
       "      <td>0.052103</td>\n",
       "      <td>...</td>\n",
       "      <td>0.121882</td>\n",
       "      <td>0.093363</td>\n",
       "      <td>0.942534</td>\n",
       "      <td>0.472701</td>\n",
       "      <td>-0.514804</td>\n",
       "      <td>-0.204782</td>\n",
       "      <td>0.434075</td>\n",
       "      <td>0.746962</td>\n",
       "      <td>0.583769</td>\n",
       "      <td>-1.445666</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>5 rows × 24 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "         V1        V2        V3        V4        V5        V6        V7  \\\n",
       "0 -0.288362 -0.170262  0.036313  0.561619 -0.477592  0.026122  0.235572   \n",
       "1  0.172179  0.498972 -0.637497 -0.403630 -0.510466  0.803194 -0.629094   \n",
       "2 -0.148207  0.089956 -0.371467 -0.015495 -0.305202 -0.524397 -0.312337   \n",
       "3  0.392713  0.466647 -0.414478  0.146909  0.437273  0.467879  0.124036   \n",
       "4 -0.382722 -0.243545 -0.381480 -0.179152 -0.413537 -0.195249  0.376594   \n",
       "\n",
       "         V8        V9       V10  ...         X       V15       V16       V17  \\\n",
       "0  0.193228  0.219485  0.224018  ... -0.334338 -0.274177  0.069914 -0.425716   \n",
       "1  0.042610 -0.364143  0.421791  ...  0.177786  0.222926  0.336022  1.071612   \n",
       "2  0.700555  0.219698  0.031795  ... -0.008598  0.869035  0.767204 -0.080318   \n",
       "3 -0.259591  0.100210 -0.001243  ... -0.220841  0.620429  0.497164  0.094684   \n",
       "4  0.571088 -0.329178  0.052103  ...  0.121882  0.093363  0.942534  0.472701   \n",
       "\n",
       "        V18       V19       V20         Y       V21       V22  \n",
       "0 -0.244209  0.489664  0.597804 -0.410179  0.024974  0.362765  \n",
       "1  0.246272 -0.468573 -0.761310 -0.179248  1.007144 -1.160541  \n",
       "2 -0.410839  0.171747 -0.962758  0.183980  0.891160 -1.165926  \n",
       "3 -0.674582  0.418828 -0.232072  0.204818 -0.205904  0.029435  \n",
       "4 -0.514804 -0.204782  0.434075  0.746962  0.583769 -1.445666  \n",
       "\n",
       "[5 rows x 24 columns]"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data = get_dataset(mixing_matrix=A, noise_vars=noise_vars,\n",
    "                   num_samples=1000, labels=labels)\n",
    "data.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "_REfaO6ID9N5"
   },
   "outputs": [],
   "source": [
    "sample_cov = data.cov().values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "-Q_8Qc_jdsm6",
    "outputId": "e909fad7-7577-44e0-a42f-79b2aef6a76b"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Total CI tests done: 1303\n",
      "Estimated ATE set: {-0.25886095726500385, 0.4173748925681714}\n"
     ]
    }
   ],
   "source": [
    "pc_alg = PCAlgorithm(use_ci_oracle=True, graph_true=graph_true)\n",
    "cpdag_pc, _, _ = pc_alg.run(data)\n",
    "print(\"Total CI tests done: %d\" % pc_alg.ci_test_calls[\"total\"])\n",
    "print(\"Estimated ATE set: %s\" % str(\n",
    "    set(get_ATE_using_cov(g, sample_cov) for g in generate_local_dags_from(cpdag_pc))))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "ocEqP4jDejM1",
    "outputId": "dfc5fbf3-7837-460e-c601-6d1787e6179b"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Total CI tests done: 749\n",
      "Estimated ATE set: {0.4173748925681714, -0.25886095726500385}\n"
     ]
    }
   ],
   "source": [
    "mb_by_mb_alg = MBbyMBAlgorithm(use_ci_oracle=True, graph_true=graph_true)\n",
    "result_mb_by_mb = mb_by_mb_alg.run(data)\n",
    "print(\"Total CI tests done: %d\" % mb_by_mb_alg.ci_test_calls[\"total\"])\n",
    "print(\"Estimated ATE set: %s\" % str(\n",
    "    set(get_ATE_using_nodes_and_cov(list(data.columns),\n",
    "                                    list(result_mb_by_mb[\"tmt_parents\"]) + par,\n",
    "                                    sample_cov)\n",
    "    for par in get_all_combinations(result_mb_by_mb[\"unoriented\"], \n",
    "                                    result_mb_by_mb[\"non_colliders\"]))))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "AhlAd4GUjgA5",
    "outputId": "0ab5914e-694d-45d2-a048-a4af846f5a75"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Total CI tests done: 883\n",
      "Estimated ATE set: {0.4173748925681714, -0.25886095726500385}\n"
     ]
    }
   ],
   "source": [
    "sd_alg = SequentialDiscoveryAlgorithm(use_ci_oracle=True, graph_true=graph_true)\n",
    "result_sd = sd_alg.run(data)\n",
    "print(\"Total CI tests done: %d\" % sd_alg.ci_test_calls[\"total\"])\n",
    "print(\"Estimated ATE set: %s\" % str(\n",
    "    set(get_ATE_using_nodes_and_cov(list(data.columns),\n",
    "                                    list(result_sd[\"tmt_parents\"]) + par,\n",
    "                                    sample_cov)\n",
    "    for par in get_all_combinations(result_sd[\"unoriented\"], \n",
    "                                    result_sd[\"non_colliders\"]))))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "TV5TXUenuTbq",
    "outputId": "8538dcf1-14ba-40dc-af5f-d53f9383edce"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Total CI tests done: 1618\n",
      "Estimated ATE set: {0.4173748925681714, -0.25886095726500385}\n"
     ]
    }
   ],
   "source": [
    "ldecc_alg = LDECCAlgorithm(use_ci_oracle=True, graph_true=graph_true, ldecc_do_checks=True)\n",
    "result_ldecc = ldecc_alg.run(data)\n",
    "print(\"Total CI tests done: %d\" % ldecc_alg.ci_test_calls[\"total\"])\n",
    "print(\"Estimated ATE set: %s\" % str(\n",
    "    set(get_ATE_using_nodes_and_cov(list(data.columns),\n",
    "                                    list(result_ldecc[\"tmt_parents\"]) + par,\n",
    "                                    sample_cov)\n",
    "    for par in get_all_combinations(result_ldecc[\"unoriented\"], \n",
    "                                    result_ldecc[\"non_colliders\"]))))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "bcQ_-A7eP5y9"
   },
   "source": [
    "The code below evaluates the various algorithms on randomly generated linear graphs.\n",
    "We use `ipyparallel` to execute the runs in parallel. The following command starts the `ipyparallel` engines:\n",
    "\n",
    "```\n",
    "ipcluster start -n <num_engines>\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "R74t7mjIfRDj"
   },
   "outputs": [],
   "source": [
    "import ipyparallel as ipp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "id": "dsMrRbQdhDnY"
   },
   "outputs": [],
   "source": [
    "# Verify that ipcluster is running and import the necessary Python packages.\n",
    "\n",
    "parallel_client = ipp.Client(debug=False)\n",
    "dview = parallel_client[:]\n",
    "# Execute an identity map in parallel.\n",
    "ar = dview.map(lambda x: x, (i for i in range(0, 20000, 2)))\n",
    "assert ar.get()[0] == 0\n",
    "\n",
    "# Import the required Python packages.\n",
    "with dview.sync_imports():\n",
    "  import numpy\n",
    "  import os\n",
    "  import pandas\n",
    "  import networkx\n",
    "  from itertools import combinations, permutations\n",
    "  import logging\n",
    "  from causallearn.utils.cit import CIT\n",
    "  import collections\n",
    "  import copy\n",
    "  import ipyparallel\n",
    "\n",
    "  from causal_discovery.utils import orient_colliders, get_dataset, get_connected_component_with_node, remove_directed_only_edges, apply_meek_rules, adj_matrix_to_nx_graph, dag_to_cpdag, get_neighbor_types, get_all_combinations, get_ATE_using_nodes_and_cov, get_ATE_using_cov, generate_local_dags_from, to_dict    \n",
    "  from causal_discovery.pc_alg import PCAlgorithm\n",
    "  from causal_discovery.mb_by_mb import MBbyMBAlgorithm\n",
    "  from causal_discovery.sd_alg import SequentialDiscoveryAlgorithm\n",
    "  from causal_discovery.ldecc import LDECCAlgorithm\n",
    "  \n",
    "  try:\n",
    "    from cPickle import dumps, loads, HIGHEST_PROTOCOL as PICKLE_PROTOCOL\n",
    "  except ImportError:\n",
    "    from pickle import dumps, loads, HIGHEST_PROTOCOL as PICKLE_PROTOCOL\n",
    "\n",
    "# Make sure ipyparallel is still able to execute functions.\n",
    "dview = parallel_client[:]\n",
    "ar = dview.map(lambda x: x, (i for i in range(0, 20000, 2)))\n",
    "assert ar.get()[0] == 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "mMjn9PkChGOp"
   },
   "outputs": [],
   "source": [
    "def execute_iteration(params, num_samples, seed, progress_log_filepath=None):\n",
    "  np = numpy\n",
    "  nx = networkx\n",
    "  \n",
    "  np.random.seed(seed)\n",
    "\n",
    "  # linearly scale down alpha.\n",
    "  alpha = (32000 / num_samples) * 0.01\n",
    "\n",
    "  result = {\n",
    "      \"graph_true_sample_cov_ate\": None,\n",
    "      \"graph_true\": None,\n",
    "      \"oracle\": {\n",
    "          \"ate\": None,\n",
    "          \"neighbors\": None,\n",
    "          \"parent_sets\": None,\n",
    "      },\n",
    "      \"oracle_sample_cov\": {\n",
    "          \"ate\": None,\n",
    "      },\n",
    "      \"pc-alg\": {\n",
    "          \"ci_tests\": None,\n",
    "          \"ate\": None,\n",
    "          \"neighbors\": None,\n",
    "          \"parent_sets\": None,\n",
    "      },\n",
    "      \"ldecc\": {\n",
    "          \"ci_tests\": None,\n",
    "          \"ate\": None,\n",
    "          \"neighbors\": None,\n",
    "          \"parent_sets\": None,\n",
    "      },\n",
    "      \"ldecc-checks\": {\n",
    "          \"ci_tests\": None,\n",
    "          \"ate\": None,\n",
    "          \"neighbors\": None,\n",
    "          \"parent_sets\": None,\n",
    "      },\n",
    "      \"mb-by-mb\": {\n",
    "          \"ci_tests\": None,\n",
    "          \"ate\": None,\n",
    "          \"neighbors\": None,\n",
    "          \"parent_sets\": None,\n",
    "      },\n",
    "      \"sd-alg\": {\n",
    "          \"ci_tests\": None,\n",
    "          \"ate\": None,\n",
    "          \"neighbors\": None,\n",
    "          \"parent_sets\": None,\n",
    "      },\n",
    "  }\n",
    "\n",
    "  A, noise_vars, labels, cov_matrix = params\n",
    "  graph_true = adj_matrix_to_nx_graph(A, labels)\n",
    "  data = get_dataset(mixing_matrix=A, noise_vars=noise_vars,\n",
    "                     num_samples=num_samples, labels=labels)\n",
    "  sample_cov = data.cov().values\n",
    "\n",
    "  result[\"graph_true\"] = graph_true\n",
    "  result[\"graph_true_sample_cov_ate\"] = get_ATE_using_cov(graph_true, sample_cov)\n",
    "\n",
    "  cpdag_oracle = dag_to_cpdag(graph_true)\n",
    "  result[\"oracle\"][\"ate\"] = [get_ATE_using_cov(g, cov_matrix) \n",
    "                             for g in generate_local_dags_from(cpdag_oracle)]\n",
    "  result[\"oracle\"][\"neighbors\"] = get_neighbor_types(cpdag_oracle, \"X\")                           \n",
    "  result[\"oracle\"][\"parent_sets\"] = [list(g.predecessors(\"X\")) for g in generate_local_dags_from(cpdag_oracle)]                     \n",
    "  \n",
    "  result[\"oracle_sample_cov\"][\"ate\"] = [get_ATE_using_cov(g, sample_cov)\n",
    "                                        for g in generate_local_dags_from(cpdag_oracle)]\n",
    "\n",
    "  pc_alg = PCAlgorithm(use_ci_oracle=False, alpha=alpha)\n",
    "  cpdag_pc, _, _ = pc_alg.run(data)\n",
    "  result[\"pc-alg\"][\"ci_tests\"] = to_dict(pc_alg.ci_test_calls)\n",
    "  result[\"pc-alg\"][\"neighbors\"] = get_neighbor_types(cpdag_pc, \"X\")\n",
    "  result[\"pc-alg\"][\"parent_sets\"] = [list(g.predecessors(\"X\")) for g in generate_local_dags_from(cpdag_pc)]\n",
    "  result[\"pc-alg\"][\"ate\"] = [get_ATE_using_cov(g, sample_cov)\n",
    "                             for g in generate_local_dags_from(cpdag_pc)]\n",
    "  \n",
    "  ldecc_alg = LDECCAlgorithm(use_ci_oracle=False, alpha=alpha, ldecc_do_checks=False)\n",
    "  result_ldecc = ldecc_alg.run(data)\n",
    "  result[\"ldecc\"][\"ci_tests\"] = to_dict(ldecc_alg.ci_test_calls)\n",
    "  result[\"ldecc\"][\"neighbors\"] = {\n",
    "      \"parents\": result_ldecc[\"tmt_parents\"],\n",
    "      \"children\": result_ldecc[\"tmt_children\"],\n",
    "      \"unoriented\": result_ldecc[\"unoriented\"], \n",
    "  }\n",
    "  tmt_par = list(result_ldecc[\"tmt_parents\"])\n",
    "  result[\"ldecc\"][\"parent_sets\"] = list(list(p) + tmt_par for p in get_all_combinations(result_ldecc[\"unoriented\"], result_ldecc[\"non_colliders\"]))\n",
    "  ate_ldecc = [get_ATE_using_nodes_and_cov(list(data.columns), tmt_par + par, sample_cov)\n",
    "               for par in get_all_combinations(result_ldecc[\"unoriented\"], \n",
    "                                               result_ldecc[\"non_colliders\"])]\n",
    "  result[\"ldecc\"][\"ate\"] = ate_ldecc\n",
    "\n",
    "  ldecc_alg = LDECCAlgorithm(use_ci_oracle=False, alpha=alpha, ldecc_do_checks=True)\n",
    "  result_ldecc = ldecc_alg.run(data)\n",
    "  result[\"ldecc-checks\"][\"ci_tests\"] = to_dict(ldecc_alg.ci_test_calls)\n",
    "  result[\"ldecc-checks\"][\"neighbors\"] = {\n",
    "      \"parents\": result_ldecc[\"tmt_parents\"],\n",
    "      \"children\": result_ldecc[\"tmt_children\"],\n",
    "      \"unoriented\": result_ldecc[\"unoriented\"], \n",
    "  }\n",
    "  tmt_par = list(result_ldecc[\"tmt_parents\"])\n",
    "  result[\"ldecc-checks\"][\"parent_sets\"] = list(list(p) + tmt_par for p in get_all_combinations(result_ldecc[\"unoriented\"], result_ldecc[\"non_colliders\"]))\n",
    "  ate_ldecc = [get_ATE_using_nodes_and_cov(list(data.columns), tmt_par + par, sample_cov)\n",
    "               for par in get_all_combinations(result_ldecc[\"unoriented\"], \n",
    "                                               result_ldecc[\"non_colliders\"])]\n",
    "  result[\"ldecc-checks\"][\"ate\"] = ate_ldecc\n",
    "\n",
    "  mb_by_mb_alg = MBbyMBAlgorithm(use_ci_oracle=False, alpha=alpha)\n",
    "  result_mb = mb_by_mb_alg.run(data)\n",
    "  result[\"mb-by-mb\"][\"ci_tests\"] = to_dict(mb_by_mb_alg.ci_test_calls)\n",
    "  result[\"mb-by-mb\"][\"neighbors\"] = {\n",
    "      \"parents\": result_mb[\"tmt_parents\"],\n",
    "      \"children\": result_mb[\"tmt_children\"],\n",
    "      \"unoriented\": result_mb[\"unoriented\"], \n",
    "  }\n",
    "  tmt_par = list(result_mb[\"tmt_parents\"])\n",
    "  result[\"mb-by-mb\"][\"parent_sets\"] = list(list(p) + tmt_par for p in get_all_combinations(result_mb[\"unoriented\"], result_mb[\"non_colliders\"]))\n",
    "  ate_mb = [get_ATE_using_nodes_and_cov(list(data.columns), tmt_par + par, sample_cov)\n",
    "            for par in get_all_combinations(result_mb[\"unoriented\"], \n",
    "                                            result_mb[\"non_colliders\"])]\n",
    "  result[\"mb-by-mb\"][\"ate\"] = ate_mb\n",
    "\n",
    "  sd_alg = SequentialDiscoveryAlgorithm(use_ci_oracle=False, alpha=alpha)\n",
    "  result_sd = sd_alg.run(data)\n",
    "  result[\"sd-alg\"][\"ci_tests\"] = to_dict(sd_alg.ci_test_calls)\n",
    "  result[\"sd-alg\"][\"neighbors\"] = {\n",
    "      \"parents\": result_sd[\"tmt_parents\"],\n",
    "      \"children\": result_sd[\"tmt_children\"],\n",
    "      \"unoriented\": result_sd[\"unoriented\"], \n",
    "  }\n",
    "  tmt_par = list(result_sd[\"tmt_parents\"])\n",
    "  result[\"sd-alg\"][\"parent_sets\"] = list(list(p) + tmt_par for p in get_all_combinations(result_sd[\"unoriented\"], result_sd[\"non_colliders\"]))\n",
    "  ate_sd = [get_ATE_using_nodes_and_cov(list(data.columns), tmt_par + par, sample_cov)\n",
    "            for par in get_all_combinations(result_sd[\"unoriented\"], \n",
    "                                            result_sd[\"non_colliders\"])]\n",
    "  result[\"sd-alg\"][\"ate\"] = ate_sd\n",
    "  \n",
    "  if progress_log_filepath is not None:\n",
    "    os.system(\"echo '%d' >> %s\" % (seed, progress_log_filepath))\n",
    "\n",
    "  np.random.seed(None)\n",
    "  return result"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "7ii8I2ZPlWYh"
   },
   "outputs": [],
   "source": [
    "def async_results_to_results(async_results, results):\n",
    "  for N, ar in zip(results[\"sample_sizes\"], async_results):\n",
    "    if N not in results:\n",
    "      results[N] = []\n",
    "    for r in ar.get():\n",
    "      results[N].append(r)\n",
    "  return results\n",
    "\n",
    "def execute_estimation_in_parallel(results, start_graph, end_graph, iters_per_graph=1):\n",
    "\n",
    "  num_threads = len(parallel_client.ids)\n",
    "  dview = parallel_client[:]\n",
    "\n",
    "  num_graphs = end_graph - start_graph\n",
    "  print(\"Executing %d iterations across %d engines\" % (\n",
    "      num_graphs*iters_per_graph, num_threads))\n",
    "\n",
    "  dview[\"get_dataset\"] = get_dataset\n",
    "  dview[\"get_noise_matrix\"] = get_noise_matrix\n",
    "  dview[\"execute_iteration\"] = execute_iteration\n",
    "\n",
    "  # The random seed below replicates the results in the paper.\n",
    "  RANDOM_SEED = 716697\n",
    "  np.random.seed(RANDOM_SEED)\n",
    "  new_param_list = []\n",
    "  \n",
    "  erdos_renyi_seed = RANDOM_SEED\n",
    "  for i_param in range(start_graph, end_graph):\n",
    "    N = 30\n",
    "    P = 0.08\n",
    "    MIN_SIZE = 20\n",
    "\n",
    "    A, noise_vars, labels, cov_matrix, erdos_renyi_seed = get_params(n=N, p=P, min_size=MIN_SIZE, seed=erdos_renyi_seed, get_next_seed=True)\n",
    "    params = (A, noise_vars, labels, cov_matrix)\n",
    "    graph_true = adj_matrix_to_nx_graph(A, labels)\n",
    "\n",
    "    for graph_iter in range(iters_per_graph):\n",
    "      new_param_list.append(params)\n",
    "      results[\"params\"].append(params)\n",
    "      results[\"ground_truth_ate\"].append(get_ATE_using_cov(graph_true, cov_matrix))\n",
    "      \n",
    "  sample_sizes = [8000, 16000, 24000, 32000]\n",
    "  results[\"sample_sizes\"] = sample_sizes\n",
    "  async_results = []\n",
    "  for N in sample_sizes:\n",
    "    print(\"%s: Starting iterations %d for N=%d\" % (\n",
    "        str(datetime.now()), len(new_param_list), N))\n",
    "    \n",
    "    def execute_estimation_iteration(i):\n",
    "      return execute_iteration(new_param_list[i], N, (RANDOM_SEED + N + i))\n",
    "    \n",
    "    async_result = dview.map(execute_estimation_iteration, range(len(new_param_list)))\n",
    "    async_results.append(async_result)\n",
    "\n",
    "  np.random.seed(None)\n",
    "  return results, async_results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "rvMPJwZfodBA"
   },
   "outputs": [],
   "source": [
    "results = {}\n",
    "results.update({\n",
    "    \"params\": [],\n",
    "    \"ground_truth_ate\": [],\n",
    "})\n",
    "\n",
    "results, async_results = execute_estimation_in_parallel(results, \n",
    "                                                        start_graph=0,\n",
    "                                                        end_graph=250,\n",
    "                                                        iters_per_graph=5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "Vfbe49lusOGd"
   },
   "outputs": [],
   "source": [
    "# Monitor progress of the tasks.\n",
    "for ar in async_results:\n",
    "  ar.wait_interactive(timeout=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "Fg8iX-g7sYBv"
   },
   "outputs": [],
   "source": [
    "results = async_results_to_results(async_results, results)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "FoTkEiC7QZ40"
   },
   "source": [
    "The code below plots the results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "lDvbGe-Y1cOP"
   },
   "outputs": [],
   "source": [
    "# For the color map:\n",
    "# https://gist.github.com/AndiH/c957b4d769e628f506bd\n",
    "\n",
    "# Tableau 20 Colors\n",
    "tableau20 = [(31, 119, 180), (174, 199, 232), (255, 127, 14), (255, 187, 120),  \n",
    "             (44, 160, 44), (152, 223, 138), (214, 39, 40), (255, 152, 150),  \n",
    "             (148, 103, 189), (197, 176, 213), (140, 86, 75), (196, 156, 148),  \n",
    "             (227, 119, 194), (247, 182, 210), (127, 127, 127), (199, 199, 199),  \n",
    "             (188, 189, 34), (219, 219, 141), (23, 190, 207), (158, 218, 229)]\n",
    "             \n",
    "# Tableau Color Blind 10\n",
    "tableau20blind = [(0, 107, 164), (255, 128, 14), (171, 171, 171), (89, 89, 89),\n",
    "             (95, 158, 209), (200, 82, 0), (137, 137, 137), (163, 200, 236),\n",
    "             (255, 188, 121), (207, 207, 207)]\n",
    "  \n",
    "# Rescale to values between 0 and 1 \n",
    "for i in range(len(tableau20)):  \n",
    "    r, g, b = tableau20[i]  \n",
    "    tableau20[i] = (r / 255., g / 255., b / 255.)\n",
    "for i in range(len(tableau20blind)):  \n",
    "    r, g, b = tableau20blind[i]  \n",
    "    tableau20blind[i] = (r / 255., g / 255., b / 255.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "U9mLTMRc1kkK"
   },
   "outputs": [],
   "source": [
    "SMALL_SIZE = 8\n",
    "MEDIUM_SIZE = 8\n",
    "BIGGER_SIZE = 8\n",
    "\n",
    "plt.rc('font', size=SMALL_SIZE)          # controls default text sizes\n",
    "plt.rc('axes', titlesize=SMALL_SIZE+6)     # fontsize of the axes title\n",
    "plt.rc('axes', labelsize=MEDIUM_SIZE+6)    # fontsize of the x and y labels\n",
    "plt.rc('xtick', labelsize=SMALL_SIZE + 4)    # fontsize of the tick labels\n",
    "plt.rc('ytick', labelsize=SMALL_SIZE+4)    # fontsize of the tick labels\n",
    "plt.rc('legend', fontsize=SMALL_SIZE+2)    # legend fontsize\n",
    "plt.rc('figure', titlesize=BIGGER_SIZE+20)  # fontsize of the figure title"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "2MAKzoHd1mJj"
   },
   "outputs": [],
   "source": [
    "plot_prop = {\n",
    "    \"pc-alg\": [\"dashed\", tableau20blind[0], \"o\"],\n",
    "    \"ldecc\": [\"dashdot\", tableau20blind[1], \"^\"],\n",
    "    \"ldecc-checks\": [\"solid\", tableau20blind[5], \"s\"],\n",
    "    \"mb-by-mb\": [\"dashdot\", tableau20blind[7], \"D\"],\n",
    "    \"sd-alg\": [\"dotted\", tableau20blind[3], \"v\"],\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "NmlANMMf2HlC"
   },
   "outputs": [],
   "source": [
    "def get_hausdorff_distance(arr_1, arr_2):\n",
    "\n",
    "  arr_1 = np.array(arr_1)\n",
    "  arr_2 = np.array(arr_2)\n",
    "\n",
    "  def get_distance_point_from_array(point, arr):\n",
    "    return np.min(np.abs(arr - point))\n",
    "\n",
    "  dist_1_to_2 = np.max([get_distance_point_from_array(p, arr_2) for p in arr_1])\n",
    "  dist_2_to_1 = np.max([get_distance_point_from_array(p, arr_1) for p in arr_2])\n",
    "\n",
    "  return max(dist_1_to_2, dist_2_to_1)\n",
    "\n",
    "def plot_mse(results):\n",
    "  plt.figure(figsize=(6, 4))\n",
    "  samples = results[\"sample_sizes\"]\n",
    "\n",
    "  dict_to_name = {\n",
    "      'pc-alg': \"PC\", \n",
    "      'ldecc': \"LDECC\", \n",
    "      'ldecc-checks': \"LDECC-checks\",\n",
    "      'sd-alg': \"SD\",\n",
    "      'mb-by-mb': \"MB-by-MB\",\n",
    "  }\n",
    "\n",
    "  errors = collections.defaultdict(list)\n",
    "\n",
    "  for N in samples:\n",
    "    for k in dict_to_name.keys():\n",
    "      errors[k].append([])\n",
    "\n",
    "    for i, res in enumerate(results[N]):\n",
    "      for k in dict_to_name.keys():\n",
    "        errors[k][-1].append(get_hausdorff_distance(\n",
    "            res[\"oracle\"][\"ate\"], res[k][\"ate\"]\n",
    "        )**2)\n",
    "\n",
    "      \n",
    "  samples = np.array(samples)\n",
    "  mses = {}\n",
    "  stds = {}\n",
    "  for k in dict_to_name.keys():\n",
    "    errors[k] = np.array(errors[k])\n",
    "    mses[k] = np.mean(errors[k], axis=-1)\n",
    "    stds[k] = np.std(errors[k], axis=-1)\n",
    "\n",
    "    plt.plot(samples, mses[k], label=dict_to_name[k],\n",
    "             color=plot_prop[k][1], linestyle=plot_prop[k][0],\n",
    "             marker=plot_prop[k][2])\n",
    "    ci = 1.96 * stds[k] / np.sqrt(errors[k].shape[-1])\n",
    "    plt.errorbar(samples, mses[k],\n",
    "                 yerr=ci, ls=\"none\", color=plot_prop[k][1])\n",
    "  \n",
    "  plt.title(\"MSE vs Sample size\")\n",
    "  plt.xlabel(\"Sample size\")\n",
    "  plt.ylabel(\"MSE (Hausdorff distance)\")\n",
    "  plt.legend()\n",
    "  plt.show()\n",
    "  \n",
    "plot_mse(results)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "7q0qKVkeAZd8"
   },
   "outputs": [],
   "source": [
    "def plot_accuracy(results):\n",
    "\n",
    "  dict_to_name = {\n",
    "      'pc-alg': \"PC\", \n",
    "      'ldecc': \"LDECC\", \n",
    "      'ldecc-checks': \"LDECC-checks\",\n",
    "      'sd-alg': \"SD\",\n",
    "      'mb-by-mb': \"MB-by-MB\",\n",
    "  }\n",
    "\n",
    "  def is_correct_adj_set(r, key):\n",
    "    return set(frozenset(p) for p in r[\"oracle\"][\"parent_sets\"]) == set(frozenset(p) for p in r[key][\"parent_sets\"])\n",
    "\n",
    "  N = len(results[results[\"sample_sizes\"][0]])\n",
    "  print(\"Total = %d\" % N)\n",
    "\n",
    "  same_count = collections.defaultdict(lambda: 0)\n",
    "  correct_count = collections.defaultdict(lambda: collections.defaultdict(lambda: 0))\n",
    "  \n",
    "  for s in results[\"sample_sizes\"]:\n",
    "    for i in range(N):\n",
    "      for k in dict_to_name.keys():\n",
    "        if is_correct_adj_set(r=results[s][i], key=k):\n",
    "          correct_count[k][s] += 1\n",
    "        \n",
    "  plt.figure(figsize=(6, 4))\n",
    "\n",
    "  accuracy = {}\n",
    "  samples = results[\"sample_sizes\"]\n",
    "  for k in dict_to_name.keys():\n",
    "    accuracy[k] = np.array([correct_count[k][s] * 100 / N for s in samples])\n",
    "\n",
    "    plt.plot(samples, accuracy[k], label=dict_to_name[k],\n",
    "             color=plot_prop[k][1], linestyle=plot_prop[k][0],\n",
    "             marker=plot_prop[k][2])\n",
    "  \n",
    "  print(accuracy)\n",
    "  plt.title(\"Accuracy vs Sample size\")\n",
    "  plt.xlabel(\"Sample size\")\n",
    "  plt.ylabel(\"Accuracy (%)\")\n",
    "  plt.legend()\n",
    "  plt.show()\n",
    "\n",
    "plot_accuracy(results)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "FOmTtLHaBBX2"
   },
   "outputs": [],
   "source": [
    "def plot_recall(results):\n",
    "\n",
    "  dict_to_name = {\n",
    "      'pc-alg': \"PC\", \n",
    "      'ldecc': \"LDECC\", \n",
    "      'ldecc-checks': \"LDECC-checks\",\n",
    "      'sd-alg': \"SD\",\n",
    "      'mb-by-mb': \"MB-by-MB\",\n",
    "  }\n",
    "\n",
    "  N = len(results[results[\"sample_sizes\"][0]])\n",
    "  \n",
    "  correct_count = collections.defaultdict(lambda: collections.defaultdict(lambda: 0))\n",
    "  \n",
    "  for s in results[\"sample_sizes\"]:\n",
    "    for i in range(N):\n",
    "      for k in dict_to_name.keys():\n",
    "\n",
    "        if frozenset(results[s][i][\"graph_true\"].predecessors(\"X\")) in set(frozenset(p) for p in results[s][i][k][\"parent_sets\"]):\n",
    "          correct_count[k][s] += 1\n",
    "  \n",
    "  plt.figure(figsize=(6, 4))\n",
    "  \n",
    "  coverage = {}\n",
    "  samples = results[\"sample_sizes\"]\n",
    "  for k in dict_to_name.keys():\n",
    "    coverage[k] = np.array([correct_count[k][s] * 100 / N for s in samples])\n",
    "\n",
    "    plt.plot(samples, coverage[k], label=dict_to_name[k],\n",
    "             color=plot_prop[k][1], linestyle=plot_prop[k][0],\n",
    "             marker=plot_prop[k][2])\n",
    "  \n",
    "  plt.title(\"Recall vs Sample size\")\n",
    "  plt.xlabel(\"Sample size\")\n",
    "  plt.ylabel(\"Recall (%)\")\n",
    "  plt.legend()\n",
    "  plt.show()\n",
    "\n",
    "plot_recall(results)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "aDhNjpg8g_sG"
   },
   "outputs": [],
   "source": [
    "def plot_number_of_tests(results):\n",
    "\n",
    "  plt.figure(figsize=(6, 4))\n",
    "\n",
    "  dict_to_name = {\n",
    "      'pc-alg': \"PC\", \n",
    "      'ldecc': \"LDECC\", \n",
    "      'ldecc-checks': \"LDECC-checks\",\n",
    "      'sd-alg': \"SD\",\n",
    "      'mb-by-mb': \"MB-by-MB\",\n",
    "  }\n",
    "\n",
    "  samples = results[\"sample_sizes\"]\n",
    "  samples = np.array(samples)\n",
    "  \n",
    "  tests_count = collections.defaultdict(lambda: collections.defaultdict(lambda: 0))\n",
    "  for k in dict_to_name.keys():\n",
    "    for s in samples:\n",
    "      tests_count[k][s] = np.mean([t[k][\"ci_tests\"][\"total\"] for t in results[s]])\n",
    "\n",
    "  for k in dict_to_name.keys():\n",
    "    plt.plot(samples, [tests_count[k][s] for s in samples], label=dict_to_name[k],\n",
    "             color=plot_prop[k][1], linestyle=plot_prop[k][0],\n",
    "             marker=plot_prop[k][2])\n",
    "\n",
    "  \n",
    "  plt.title(\"Number of CI tests\")\n",
    "  plt.xlabel(\"Sample size\")\n",
    "  plt.ylabel(\"Number of CI tests\")\n",
    "  plt.legend()\n",
    "  plt.show()\n",
    "\n",
    "plot_number_of_tests(results)"
   ]
  }
 ],
 "metadata": {
  "colab": {
   "provenance": []
  },
  "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.8.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}