diff --git a/doc/conf.py b/doc/conf.py
index 64ee47c8..c8b3b805 100644
--- a/doc/conf.py
+++ b/doc/conf.py
@@ -264,7 +264,7 @@
}
# prevent jupyter notebooks from being run even if empty cell
-nbsphinx_execute = "never"
+# nbsphinx_execute = "never"
nbsphinx_allow_errors = True
# Custom sidebar templates, maps document names to template names.
diff --git a/doc/tutorials/markovian/example-pc-algo.ipynb b/doc/tutorials/markovian/example-pc-algo.ipynb
index f49e4198..d5fa254e 100644
--- a/doc/tutorials/markovian/example-pc-algo.ipynb
+++ b/doc/tutorials/markovian/example-pc-algo.ipynb
@@ -7,7 +7,7 @@
"source": [
"# PC algorithm for causal discovery from observational data without latent confounders\n",
"\n",
- "In this tutorial, we will demonstrate how to use the PC algorithm to learn a causal graph structure.\n",
+ "In this tutorial, we will demonstrate how to use the PC algorithm to learn a causal graph structure and highlight some of the common challenges in applying causal discovery algorithms to data.\n",
"\n",
"The PC algorithm works on observational data when there are no unobserved latent confounders."
]
@@ -25,6 +25,7 @@
"\n",
"from pywhy_graphs import CPDAG\n",
"from pywhy_graphs.viz import draw\n",
+ "\n",
"from dodiscover import PC, make_context\n",
"from dodiscover.ci import GSquareCITest, Oracle"
]
@@ -78,7 +79,7 @@
"
0
\n",
"
1
\n",
"
1
\n",
- "
0
\n",
+ "
1
\n",
"
1
\n",
"
1
\n",
"
1
\n",
@@ -89,33 +90,33 @@
"
1
\n",
"
1
\n",
"
1
\n",
- "
1
\n",
- "
1
\n",
"
0
\n",
- "
1
\n",
- "
1
\n",
+ "
0
\n",
+ "
0
\n",
+ "
0
\n",
+ "
0
\n",
"
0
\n",
" \n",
"
\n",
"
2
\n",
"
1
\n",
"
1
\n",
- "
1
\n",
+ "
0
\n",
"
1
\n",
"
0
\n",
"
1
\n",
"
1
\n",
- "
0
\n",
+ "
1
\n",
"
\n",
"
\n",
"
3
\n",
"
1
\n",
- "
0
\n",
- "
0
\n",
"
1
\n",
"
0
\n",
+ "
1
\n",
"
0
\n",
- "
0
\n",
+ "
1
\n",
+ "
1
\n",
"
0
\n",
"
\n",
"
\n",
@@ -124,10 +125,10 @@
"
1
\n",
"
1
\n",
"
1
\n",
- "
0
\n",
"
1
\n",
"
1
\n",
- "
0
\n",
+ "
1
\n",
+ "
1
\n",
"
\n",
" \n",
"\n",
@@ -135,11 +136,11 @@
],
"text/plain": [
" A T S L B E X D\n",
- "0 1 1 0 1 1 1 1 1\n",
- "1 1 1 1 1 0 1 1 0\n",
- "2 1 1 1 1 0 1 1 0\n",
- "3 1 0 0 1 0 0 0 0\n",
- "4 1 1 1 1 0 1 1 0"
+ "0 1 1 1 1 1 1 1 1\n",
+ "1 1 1 0 0 0 0 0 0\n",
+ "2 1 1 0 1 0 1 1 1\n",
+ "3 1 1 0 1 0 1 1 0\n",
+ "4 1 1 1 1 1 1 1 1"
]
},
"execution_count": 2,
@@ -165,27 +166,11 @@
]
},
{
- "attachments": {},
"cell_type": "markdown",
"id": "84323abd",
"metadata": {},
"source": [
- "We'll use the PC algorithm to infer the ASIA network. The ASIA network is a case study of an expert system for diagnosing lung disease from Lauritzen and Spiegelhalter (1988). Given respiratory symptions and other evidence, the goal is to distinguish between tuberculosis, lung cancer or bronchitis in a given patient. The ground truth causal DAG is as follows:\n",
- "\n",
- "\n",
- "\n",
- "The variables in the DAG have the following interpretation:\n",
- "\n",
- "* T: Whether or not the patient has **tuberculosis**.\n",
- "* L: Whether or not the patient has **lung cancer**.\n",
- "* B: Whether or not the patient has **bronchitis**.\n",
- "* A: Whether or not the patient has recently visited **Asia**.\n",
- "* S: Whether or not the patient is a **smoker**.\n",
- "* E: An indicator of whether the patient has either lung cancer or tuberculosis (or both).\n",
- "* X: Whether or not a chest X-ray shows evidence of tuberculosis or lung cancer.\n",
- "* D: Whether or not the patient has **dyspnoea** (difficulty breathing).\n",
- "\n",
- "Note we have three kinds of variables, diseases (B, L, T, and E which indicates one or more diseases), symptoms (X and D), and behaviors (S and A). The goal of the model is to use symptoms and behaviors to diagnose (i.e. infer) diseases. Further, note that diseases are causes of symptoms, and are effects of behaviors."
+ "We'll use the PC algorithm to infer the ASIA network. The ASIA network is a case study of an expert system for diagnosing lung disease from Lauritzen and Spiegelhalter (1988). Given respiratory symptions and other evidence, the goal is to distinguish between tuberculosis, lung cancer or bronchitis in a given patient."
]
},
{
@@ -217,6 +202,55 @@
"ground_truth = nx.DiGraph(ground_truth_edges)"
]
},
+ {
+ "cell_type": "markdown",
+ "id": "de591d7c",
+ "metadata": {},
+ "source": [
+ "The ground truth DAG can be visualized and is seen as follows:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 4,
+ "id": "d7815b62",
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "image/png": "iVBORw0KGgoAAAANSUhEUgAAApQAAAHzCAYAAACe1o1DAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAA9hAAAPYQGoP6dpAABhYUlEQVR4nO3dd1gUZ6MF8LMFUBALWBKjaIztC6LGiL1hFDUqWCNRUBG7SBTFgl0jdk3sxhJjL1GjJhpLxK6IlQULEFRskaZIkbK7c/8wejWigLvL7Oye3/Pc537uLDMHAsPhnXfekQmCIICIiIiI6APJxQ5ARERERNLGQklEREREOmGhJCIiIiKdsFASERERkU5YKImIiIhIJyyURERERKQTFkoiIiIi0gkLJRERERHphIWSiIiIiHTCQklEREREOmGhJCIiIiKdsFASERERkU5YKImIiIhIJyyURERERKQTFkoiIiIi0gkLJRERERHphIWSiIiIiHTCQklEREREOmGhJCIiIiKdsFASERERkU5YKImIiIhIJyyURERERKQTFkoiIiIi0gkLJRERERHphIWSiIiIiHTCQklEREREOmGhJCIiIiKdsFASERERkU5YKImIiIhIJyyURERERKQTFkoiIiIi0gkLJRERERHphIWSiIiIiHTCQklEREREOmGhJCIiIiKdsFASERERkU5YKImIiIhIJ0qxAxARfai0TDXuJKYhS62FpVKOivY2sLHiaY2IqKDxzEtEkhL1OAWbQ2IRfCsOsUnpEF7bJgPgYGcNl2ql0au+A6qUsRUrJhGRWZEJgiDk/jYiInHdS0pH4B4VTkUnQCGXQaN996nr5famlUsiqLMTyttZF2BSIiLzw0JJREZvW2gspuyLgForvLdI/pdCLoNSLsM0N0d4ODsYMCERkXljoSQio7Y0OArzD0fqvJ/RrlXh61JFD4mIiOi/eJc3ERmtbaGxeimTADD/cCS2h8bqZV9ERPQmFkoiMkr3ktIxZV/Ee9+TcvkP3J3dAY9+8c/TPifvi8C9pHR9xCMiotewUBKRUQrco4I6l/mSaRHHoShWBlmPIpH95GGu+1RrBQTuUekrIhER/YuFkoiMTtTjFJyKTnjvDTjZT/9B5oMbsGvpA7l1MaRFHM91vxqtgFPRCYiOS9FjWiIiYqEkIqOzOSQWCrnsve9JizgOeaEiKFzZGdbVGuepUAIv7vzedJ5zKYmI9ImFkoiMTvCtuFyXB0q7fhzWVRtBprCAzefNoH7yEJmPcr+BR6MVEBwZp6+oREQEFkoiMjKpmWrE5nLjTOY/0VAn3of1580AAFblHKGwLZnnUcrYxHSkZap1jUpERP9ioSQio3I3MQ25LY6bFhEMuU1xFHJwAgDIZDLY/K8p0m6chKDV5HoMAcCdxDTdwxIREQAWSiIyMllq7Xu3C1oN0m+cQiGHmlAnP0b2k4fIfvIQlmWrQZv2FBl3r+nlOERElHdKsQMQEb3OUvn+v3Mz7oZBk5qE9BsnkX7j5Fvb0yKOo/CndXQ+DhER5R0LJREZlYr2NpAB77zsnRZxHHLr4rBzHfzWtvTIc0iPPAdtdibkFlbvPIbs3+MQEZF+sFASkVGxsVLCwc4ad3O4MUebnYn0yLOwrt4ENtWbvLVdUcQe6ddP4Hl0CGz+1+ydx3Cwt4aNFU9/RET6wms+RGR0XKqVznEdyufRIRCynsO6cv0cP87qk2q5LnKukMvgUrW0vqISERFYKInICPWq75DjOpRpEcchU1qi0Ke1c/w4mUyOwp8543nMZWieP8vxPRqtAM8GDvqMS0Rk9mSCIOS2QgcRUYHzWhuCszGJuS5wnh8KuQyNKtljo0/OI5xERPRhOEJJREYpqLMTlLk8fjG/lHIZgjo76XWfRETEQklERqq8nTWmuTnqdZ/T3RxR3s5ar/skIiIWSiIyYh7ODhjtWlUv+wpwrYYezpw7SURkCJxDSURGb1toLKbsi4BaK+RrTqVCLoNSLsN0N0eWSSIiA2KhJCJJuJeUjsA9KpyKToBCLntvsXy5vWnlkgjq7MTL3EREBsZCSUSSEvU4BZtDYhEcGYfYxPQ3nqgjw4tFy12qloZnAwdULm0rVkwiIrPCQklEkpWWqcZntRvAo6cnBvh4o6K9DZ+AQ0QkAp55iUiybKyUyI67jY8tM+FYtpjYcYiIzBbv8iYiSVOr1VAq+bcxEZGYWCiJSNJYKImIxMdCSUSSxkJJRCQ+FkoikjSNRsNCSUQkMhZKIpIsQRBYKImIjAALJRFJlkajAQAWSiIikbFQEpFkqdVqACyURERiY6EkIslioSQiMg4slEQkWSyURETGgYWSiCSLhZKIyDjwLExEksVCSURiSctU405iGrLUWlgq5ahobwMbK/M9F5nvZ05EksdCSUQFKepxCjaHxCL4Vhxik9IhvLZNBsDBzhou1UqjV30HVCljK1ZMUfAsTESS9bJQKhQKkZMQkSm7l5SOwD0qnIpOgEIug0YrvPUeAcDdpHRsDLmL9efuoGnlkgjq7ITydtYFH1gEnENJRJLFEUoiMrRtobFotegEzsYkAkCOZfJ1L7efjUlEq0UnsC001uAZjQHPwkQkWSyURGRIS4OjMP9w5Ad9rEYrQKMVMG63CgmpmfB1qaLndMaFI5REJFkslERkKNtCYz+4TP7X/MOR2G7iI5U8CxORZLFQEpEh3EtKx5R9ETluSw07isQDP/z/CwoLKArbwqJUBRT+zBlFnFpBbvX2vMnJ+yLQ6LOSJjunkiOURCRZLJREZAiBe1RQ5zJXsljTXrDvMAr2bYbC9ssOAIAnR1fj4VpfZMXdfuv9aq2AwD0qg+Q1BjwLE5FksVASkb5FPU7BqeiEXN9XuFJdWH38//MiizX8Bs/vXEP8r9MR9+sMlB2wAnILq1fbNVoBp6ITEB2XgsqlTW9JIY5QEpFksVASkb5tDomFQi77oI8tXLEWijXuAc2zOKRFBL+1XSGXYdN505xLyUJJRJKl0WgAsFASkf4E34rLdWmg97FxbAkAyLh95a1tGq2A4Mi4D963MWOhJCLJ4gglEelTaqYasUnpOu1DWbQkZFY2UD/9J8ftsYnpSMtU63QMY8RCSUSSxUJJRPp0NzENHz42+f/kloWgzcq5mAoA7iSm6eEoxoWFkogki4WSiPQpS63Vy360WRmQW757eSB9HceYsFASkWSxUBKRPlkqda9F6mcJEDLToCzxsUGPY2xM7zMiIrPxslAqFAqRkxCRKahob4MPu7/7/6VFHAMAFPq0To7bZf8ex9SwUBKRZHGEkoj0ycZKCQcdnmTz/M41JJ/ZDmWxMiji2CLH9zjYW8PGyvTOWab3GRGR2WChJCJ9c6lWGhtD7ua6dNDzmIvITrwPaDXQpD9Fxt1ryLh9FYpipVGq2yTIlJZvfYxCLoNL1dKGii4qnoWJSLJYKIlI33rVd8D6c3dyfV/yqc0v/odCCUUhW1iUqogSrQa881newIt1KD0bOOgxrfHgWZiIJItzKIlI36qUsUXTyiVxNiYxx1HKIjVboUjNVvner0IuQ6NK9ib52EWAhZKIJEytVkMul0Mu53RwIvowKSkpiI6ORlRUFKKionD06FGkCFZQNv1Opyfm/JdSLkNQZye97c/YsFASkWSp1Wpe7iaiD3LgwAH06dMHCQkJr16TyWQQBAElS5bEknmOGLdbpbfjTXdzRHkdbvgxdvyznogki4WSiD5U0aJF3yiTACAIL0Ykt2/fDg9nB4x2raqXYwW4VkMPZ9OcO/kSCyURSRYLJRF9qCZNmsDPz++N1xQKBdq0aYOWLVsCAHxdqmB2FydYKeVQyPO3QqVCLoOVUo45XZwwzKWy3nIbKxZKIpIsFkoi+lBXrlzBkSNHIJPJXs3D1mq1mD9//hvv6+RUGtVjfsVnRTQAkGuxfLm9USV7HB3Z3ORHJl9ioSQiyWKhJKL80mg0mDVrFurXrw9LS0vs2LHjVaHs27cvatSo8eq9ly5dQrVq1bBvyzrYXlyPIyOawat+BVSwt37riToyABXsreFVvwKOjmyGjT71TXrO5H/xTExEksVCSUT58ffff6N37944f/48xo4di6lTp8LS0hJ3797F999/j+nTpwMAsrOzMXPmTMyYMQNarRYA4OzsjCplbDHVzRFT4Yi0TDXuJKYhS62FpVKOivY2JvkEnLwy38+ciCRPo9GwUBJRrgRBwJo1azBy5EiUKVMGJ0+eROPGjV9tHzVqFIYPHw5LS0uEh4ejZ8+eCA8Pf3WTjkwmQ9GiRd/Yp42VEo5lixXo52HMeMmbiCRLrVZzUXMieq/Hjx/Dzc0NAwcORM+ePXH16tU3yuRLlpYvHpXo5uYGlUr1qkwCLwqpjY1NgWWWIhZKIpIsXvImovf57bffUKNGDVy4cAH79u3DTz/9BFvb9z+pZsOGDXByensB8iJFihgqpklgoSQiyWKhJKKcPHv2DP369UPnzp3RpEkThIeHo2PHjnn62CZNmuDYsWOwtrZG4cKFX73OQvl+LJREJFkslET0XydPnkStWrXw66+/Yt26ddi9ezdKlSqVr33MnDkTcrkcERERmD59OipWrIhq1aoZKLFpYKEkIslioSSilzIzMzF27Fi0aNEC5cuXx7Vr1+Dt7Q2ZLH8LkkdHR2PZsmUYP348Pv30U0yaNAm3b99GxYoVDRPcRPBMTESSxUJJRAAQFhYGLy8v3LhxA7Nnz8aoUaM++Ia9sWPH4qOPPsLIkSP1nNK0cYSSiCSLhZLIvGk0GsybNw/Ozs4QBAGhoaEYM2bMB5fJU6dOYffu3QgKCnpj/iTljoWSiCSLhZLIfN25cwctW7bE2LFj8d133yE0NBS1atX64P1ptVr4+/ujbt266Nmzpx6TmgeeiYlIslgoicyPIAj45Zdf4OfnBzs7OwQHB6N58+Y673fr1q24ePEiTpw48epRjJR3/IoRkWSxUBKZl/j4eHTt2hXe3t7o2rUrwsLC9FImnz9/jvHjx6Nz585o1qyZHpKaH56JiUiyWCiJzMfvv/8OHx8faLVa7N69G507d9bbvn/44Qc8evQIc+bM0ds+zQ1HKIlIslgoiUxfamoqBg4ciI4dO8LZ2RkqlUqvZTIuLg6zZs3CsGHDUKVKFb3t19zwTExEkqVWq2FlZSV2DCIykLNnz8LLywuPHz/GqlWrMGDAgHyvK5mbKVOmQKFQYPLkyXrdr7nhCCURSZZarYaFhYXYMYhIz7KysjBhwgQ0bdoUZcqUwbVr1zBw4EC9l8mIiAj89NNPmDx5Muzs7PS6b3PDEUoikiy1Wv3B680RkXGKiIiAl5cXVCoVZsyYgTFjxhhsaktAQAA+/fRTDBs2zCD7NycslEQkWZxDSWQ6tFotFi9ejHHjxuGzzz5DSEgI6tSpY7DjHTlyBAcPHsSvv/4KS0tLgx3HXPCSNxFJlkajYaEkMgGxsbFo3bo1Ro4ciSFDhuDixYsGLZMajQajRo1CkyZN0KVLF4Mdx5zwTExEksURSiJpEwQBmzdvhq+vL2xtbXH06FF89dVXBj/u+vXroVKpEBISovd5meaKI5REJFkslETSlZiYiB49esDLywsdO3aESqUqkDKZmpqKiRMnomfPnqhXr57Bj2cueCYmIslioSSSpj///BP9+vVDRkYGtm/fjm+++abAjj137lw8efIEQUFBBXZMc8ARSiKSLBZKImlJS0vD0KFD0a5dO9SsWRPh4eEFWibv37+P+fPnY+TIkahQoUKBHdcc8ExMRJLFQkkkHSEhIfDy8sL9+/exbNkyDBkypMDnL06cOBFFihTB+PHjC/S45oAjlEQkWSyURMYvOzsbU6ZMQePGjVG8eHFcuXIFQ4cOLfAyefnyZWzYsAHTpk1D0aJFC/TY5oBnYiKSLBZKIuN28+ZNeHl54cqVK5g0aRICAwNFebqVIAgYNWoUqlevjgEDBhT48c0Bz8REJFkslETGSavVYvny5QgICECFChVw7tw5ODs7i5Zn//79OH78OP744w+eMwyEl7yJSLJYKImMz4MHD9C2bVsMHz4c/fv3x+XLl0Utk9nZ2QgICECrVq3Qrl070XKYOp6JiUiyWCiJjMv27dsxZMgQFC5cGIcOHYKrq6vYkbBq1SpERUVhx44dXMTcgDhCSUSSpVaroVAoxI5BZPaePHmCnj17wsPDA66urlCpVEZRJp8+fYqpU6fC29sbtWrVEjuOSeOf9kQkWRyhJBLf0aNH0bdvX6SlpWHLli349ttvxY70ysyZM5GRkYEZM2aIHcXkcYSSiCRLo9GwUBKJ5Pnz5/juu+/QunVrVK9eHSqVyqjK5O3bt7F48WKMGTMGZcuWFTuOyeOZmIgkiyOUROK4dOkSPD09cefOHfz444/w9fWFXG5cY1Tjxo1DyZIlMWrUKLGjmAXj+q9PRJRHgiBwhJKogKnVasyYMQMNGjSAtbU1Ll++DD8/P6Mrk+fOncOOHTswc+ZM2NjYiB3HLPBMTESSpNFoAICFkqiAREVFwcvLC6GhoQgMDMSkSZNgaWkpdqy3CIIAf39/1K5dG7179xY7jtngmZiIJEmtVgNgoSQyNEEQsGrVKowaNQply5bF6dOn0bBhQ7FjvdOOHTtw/vx5/PXXX0Y3cmrK+JUmIklioSQyvEePHqF9+/YYMmTIq0coGnOZzMjIwLhx49CxY0e0bNlS7DhmhWdiIpIkFkoiw9q1axcGDRoEpVKJ33//He3btxc7Uq6WLFmCe/fu4eDBg2JHMTscoSQiSWKhJDKM5ORk9O7dG926dUPz5s0RHh4uiTKZkJCAmTNnYvDgwahevbrYccwOz8REJEkslET6d/z4cfTp0wdPnz7FL7/8Ai8vL8k8rnDatGkQBAFTpkwRO4pZ4gglEUkSCyWR/mRkZGDUqFFo2bIlKlWqhLCwMPTu3VsyZfLmzZtYsWIFJkyYgFKlSokdxyzxTExEksRCSaQfV69ehaenJ6KiojB//nyMGDFCcndHjxkzBuXLl4efn5/YUcyWtL5jiIj+9bJQKhQKkZMQSZNGo8Hs2bNRr149KJVKXLp0Cf7+/pIrk8HBwdi/fz9mz56NQoUKiR3HbEnru4aI6F8coST6cDExMWjevDkCAwMxatQohISEoEaNGmLHyjetVotRo0ahQYMG+Oabb8SOY9Z4JiYiSeKTcojyTxAErF27FiNHjkSpUqVw8uRJNGnSROxYH2zjxo24cuUKzpw5I5n5nqaKI5REJEkcoSTKn8ePH8Pd3R0DBgxAjx49cO3aNUmXybS0NAQGBqJ79+5o1KiR2HHMHs/ERCRJLJREeffbb79h4MCBAIC9e/fCzc1N5ES6W7BgARISEjB79myxoxA4QklEEsVCSZS7Z8+ewcfHB507d0bDhg0RHh5uEmXy0aNHmDt3Lvz8/FCpUiWx4xA4QklEEsVCSfR+p06dQu/evZGQkIC1a9fC29vbZOYZTpo0CYUKFcKECRPEjkL/4gglEUkSCyVRzjIzMzF27Fg0b94c5cqVw7Vr19CvXz+TKZPXrl3DunXrMGXKFBQvXlzsOPQvnomJSJJYKIneplKp4OnpiRs3bmDWrFkYPXq0Sa3VKggCRo8ejSpVqmDw4MFix6HX8ExMRJLEQkn0/zQaDRYtWoQJEyagatWquHDhAmrXri12LL07ePAgjh49ir1798LCwkLsOPQaXvImIklioSR64e7du/jqq68wZswY+Pn5ITQ01CTLpFqtxujRo9GiRQt07NhR7Dj0HzwTE5EksVCSuRMEARs2bMDw4cNRokQJBAcHo3nz5mLHMpg1a9bg5s2b2LRpk8nMBzUlHKEkIklioSRzlpCQgG7duqFv377o0qULwsLCTLpMPnv2DJMnT4aXlxfq1KkjdhzKAc/ERCRJLwulKd1wQJQXf/zxB3x8fKBWq7Fr1y506dJF7EgGN2vWLKSmpmLmzJliR6F34AglEUkSRyjJ3KSmpmLQoEHo0KEDvvzyS4SHh5tFmbx79y4WLVqE0aNHo1y5cmLHoXfgmZiIJEmj0QBgoSTzcPbsWfTu3RuPHj3CqlWrMGDAALOZRxgYGIgSJUpgzJgxYkeh9+AIJRFJEkcoyRxkZWVhwoQJaNq0KUqVKoVr165h4MCBZlMmQ0NDsWXLFsyYMQNFihQROw69B8/ERCRJarUaMpkMcjn/LibTdP36dXh6ekKlUmH69OkYO3asWf0BJQgC/P394eTkBG9vb7HjUC7M5zuTiEyKWq02q1+uZD60Wi0WL16McePGoVKlSjh//jy+/PJLsWMVuN27d+P06dM4dOgQb76TAP5pT0SSxEJJpujevXto3bo1Ro4cicGDB+PSpUtmWSazsrIwduxYtGvXDq6urmLHoTzg2ZiIJCMpKQlubm5ISUlBYmIiMjMz0bhxY1hYWGDAgAHo1auX2BGJPoggCNiyZQuGDRsGW1tbHDlyBK1atRI7lmiWLVuG27dvY+/evWJHoTxioSQiyVAoFLhy5QrS09NfvXb27FkAQLNmzcSKRaSTpKQkDBkyBDt27EDPnj2xdOlSlChRQuxYoklKSsKMGTMwYMAAODo6ih2H8oiXvIlIMooVK4bhw4e/NZ/KxsYGI0aMECcUkQ4OHTqEGjVq4MiRI9i2bRs2b95s1mUSAGbMmIHs7GxMmzZN7CiUDyyURCQp/v7+b8ydlMvlGDt2LOzs7ERMRZQ/6enp8PX1Rdu2beHk5ASVSoUePXqIHUt0UVFRWLp0KcaPH48yZcqIHYfyQSYIgiB2CCKi/PDz88PSpUshCAKKFy+O2NhY2Nraih2LKE8uXLgALy8v3Lt3D/PmzcPQoUPNZl3J3HTp0gUXL17ErVu3ULhwYbHjUD5whJKIJOf1J2ZMmDCBZZIkITs7G1OnTkWjRo1QrFgxXLlyBcOGDWOZ/NfJkyexZ88ezJo1i2VSgjhCSUSS1KhRI1y8eBHJycn85UNG79atW/D09MSVK1cwadIkBAYGwsLCQuxYRkOr1aJ+/foAgJCQED6wQIJ4lzcRSU5aphrLt+xD4tNkxDzJQkW5BWyseDoj4yMIApYtW4YxY8agfPnyOHv2LOrVqyd2LKOzdetWXLx4ESdPnmSZlCiOUBKRJEQ9TsHmkFgE34pDbFI6Xj9xyQA42FnDpVpp9KrvgCpleAmcxPfgwQP069cPhw8fxrBhwzB37lxYW1uLHcvoPH/+HNWqVUPdunWxe/dusePQB2KhJCKjdi8pHYF7VDgVnQCFXAaN9t2nrJfbm1YuiaDOTihvx1/eJI7t27djyJAhKFSoEH7++We0adNG7EhGKygoCFOnTkVERASqVKkidhz6QCyURGS0toXGYsq+CKi1wnuL5H8p5DIo5TJMc3OEh7ODARMSvenJkyfw9fXFli1b0L17d6xYsQL29vZixzJajx8/RuXKldG/f38sWrRI7DikAxZKIjJKS4OjMP9wpM77Ge1aFb4uHPUgwzt69Cj69u2L1NRULFu2DD179uQd3LkYPHgwduzYgejoaK4lK3Gc+UpERmdbaKxeyiQAzD8cie2hsXrZF1FOnj9/jhEjRqB169aoVq0aVCoVevXqxTKZi/DwcKxevRqTJk1imTQBHKEkIqNyLykdrRadQKZa+8brd2d3yNPHl/k2CIUq1HzjNSulHEdHNuecStK7S5cuwcvLCzExMZgzZw6GDx/Ou5TzqF27doiKisL169dhaWkpdhzSEdfZICKjErhHBXUO8yXtO4x6499p4ceQcefKW69blCz/1seqtQIC96iw0ae+fsOS2VKr1ZgzZw6mTp0KJycnXL58GZ9//rnYsSTj8OHD+PPPP7Fr1y6WSRPBEUoiMhpRj1PQ+oeTeXpv0uEVSLn8ByqM+z3P+z86shkql+aSQqSb6OhoeHl54cKFCxg/fjwmT57MUpQPGo0GX3zxBYoXL44TJ05waoCJ4Lg8ERmNzSGxUMgN88tFIZdh03nOpaQPJwgCVq1ahVq1aiE+Ph6nT5/G999/zzKZTz///DNUKhUWLFjAMmlCWCiJyGgE34rL1/JA+aHRCgiOjDPIvsn0PXr0CB06dMDgwYPh6emJq1evomHDhmLHkpyUlBRMnDgRPXv2hLOzs9hxSI84h5KIjEJqphqxSekGPUZsYjrSMtV8TCPly65duzBo0CAolUr8/vvvaN++vdiRJGvu3Ll4+vQpgoKCxI5CesYRSiIyCncT02DoCd0CgDuJaQY+CpmK5ORk9OnTB926dUOzZs2gUqlYJnVw//59LFiwAP7+/qhQoYLYcUjP+Gc6ERmFrP8sEyT145C0HT9+HH369MGTJ0+wfv169O7dm/P9dDRhwgTY2tpi3LhxYkchA+AIJREZBUtlwZyOCuo4JE0ZGRkYPXo0WrZsiYoVKyIsLAx9+vRhmdTRpUuXsGHDBkybNg1FixYVOw4ZAEcoicgoVLS3gQww6GVv2b/HIcrJ1atX4eXlhcjISMydOxcjR46EQqEQO5bkCYKAUaNG4fPPP0f//v3FjkMGwj/Vicgo2Fgp4WDgJ9k42Fvzhhx6i0ajwZw5c1CvXj3I5XJcvHgRo0ePZpnUk3379uHEiROYN28elEr+/JkqFkoiMhou1UobdB1Kl6qlDbJvkq6YmBi0aNEC48ePh7+/Py5cuAAnJyexY5mM7OxsjBkzBq1bt0a7du3EjkMGxEJJREajV30Hg65D6dnAwSD7JukRBAFr165FrVq1cP/+fZw4cQKzZ8+GlZWV2NFMysqVKxEVFYX58+dzHqqJ46MXicioeK0NwdmYRL0WS4VchkaV7PksbwIAxMXFYcCAAdi3bx/69euHRYsW8UYRA3jy5AmqVKmCTp06Yc2aNWLHIQPjCCURGZWgzk5Q6vmyt1IuQ1BnXsYkYO/evahRowbOnTuH3377DWvXrmWZNJCZM2ciIyMDM2bMEDsKFQAWSiIyKuXtrDHNzVGv+5zu5ojyBr7hh4xbSkoKfHx80KlTJzRs2BDh4eFwd3cXO5bJiomJwZIlSzBmzBh8/PHHYsehAsBL3kRklJYGR2H+4Uid9xPgWg3DXCrrIRFJ1enTp9G7d2/Ex8fjhx9+QL9+/Tifz8C++eYbnD17Frdu3YKNDZfqMgccoSQio+TrUgWzuzjBSinP953fCrkMVko55nRxYpk0Y5mZmRg3bhyaNWuGsmXL4tq1a/Dx8WGZNLCzZ89i586dmDlzJsukGeEIJREZtXtJ6Qjco8Kp6AQo5LL33qwjlwFaAWhauSSCOjvxMrcZCw8Ph6enJ65fv47p06cjICCA60oWAEEQ0LBhQ2RlZeHixYuQyzluZS64wigRGbXydtbY6FMfUY9TsDkkFsGRcYhNTH/jiToyAOrkf1C9BLB0xDeoXNpWrLgkMq1Wi0WLFiEwMBBVqlTBhQsXULt2bbFjmY3t27cjJCQEx44dY5k0MxyhJCLJSctU405iGrLUWlgq5ahob4NBPn2hUqlw7do1seORSO7evYs+ffrg5MmTGDlyJGbOnIlChQqJHctsZGRkoHr16qhZsyb27dsndhwqYByhJCLJsbFSwrFssTdec3d3x+bNm3Hnzh1UrFhRnGAkCkEQsHHjRgwfPhzFihXDX3/9BRcXF7FjmZ3FixfjwYMHOHTokNhRSAQcjyYik9C2bVtYWlpyZMTMJCQkoFu3bujTpw86deoElUrFMimC+Ph4zJw5E4MHD0a1atXEjkMi4CVvIjIZ7dq1Q1ZWFv766y+xo1ABOHDgAPr164fs7GysWrUK3bp1EzuS2Ro2bBg2b96M6OholCxZUuw4JAKOUBKRyXB3d8eJEyfw5MkTsaOQAaWmpmLw4MFo37496tSpg/DwcJZJEd24cQOrVq3ChAkTWCbNGEcoichkPHjwAOXKlcOmTZvQq1cvseOQAZw7dw5eXl549OgRFi5ciIEDB3JdSZF17NgR4eHhuHHjBm+CMmMcoSQik/HJJ5+gbt262Lt3r9hRSM+ysrIwceJENGnSBKVKlcLVq1cxaNAglkmRHTt2DL///jvmzJnDMmnmOEJJRCbl+++/x9y5cxEfHw8rKyux45AeXL9+HV5eXggLC8OUKVMwbtw4KJVcpERsGo0GdevWReHChXHmzBmWezPHEUoiMinu7u5ISUnB8ePHxY5COtJqtfjxxx9Rp04dpKen4/z585g4cSLLpJHYsGEDrl69igULFrBMEkcoici0CIKAzz77DG3btsXy5cvFjkMf6N69e/D29sZff/0FPz8/zJ49G4ULFxY7Fv0rLS0NVapUQdOmTbF9+3ax45AR4AglEZkUmUwGNzc37Nu3D/x7WXoEQcCWLVvg5OSEmzdv4siRI/jxxx9ZJo3M/PnzkZiYiNmzZ4sdhYwECyURmRx3d3c8ePAAly5dEjsK5UNSUhI8PDzQq1cvfP3111CpVGjVqpXYseg/Hj58iLlz5+K7777Dp59+KnYcMhIslERkcpo2bYoSJUrwbm8JOXToEJycnHD48GFs3boVW7ZsQYkSJcSORTmYNGkSChcujMDAQLGjkBFhoSQik6NUKtG+fXs+hlEC0tPT4evri7Zt28LR0RHh4eHw8PAQOxa9w9WrV/Hzzz9j6tSpKF68uNhxyIjwphwiMkm//vorunfvjpiYGF6WM1KhoaHw9PTEvXv3MHfuXAwdOhRyOcc5jJUgCGjdujXu378PlUoFCwsLsSOREeFPLhGZpDZt2sDS0pKjlEYoOzsbU6dORcOGDVGsWDFcuXIFvr6+LJNG7sCBA/jrr78wb948lkl6C0coichktWvXDpmZmTh27JjYUehft27dgpeXFy5fvoyJEydiwoQJLCcSoFarUbNmTXz00Uf466+/uO4kvYV/DhKRyXJ3d8fJkyeRlJQkdhSzJwgCli1bhi+++AJPnz7F2bNnMXXqVJZJiVi9ejVu3rzJRczpnVgoichkubm5QaPR4MCBA2JHMWsPHz5Eu3bt4OvrC29vb1y5cgX16tUTOxblUXJyMiZPnozevXvjiy++EDsOGSkWSiIyWWXLloWzszPnUYpox44dqFGjBsLCwnDw4EEsW7YMNjY2YseifJg1axbS0tIwc+ZMsaOQEWOhJCKT5ubmhoMHDyIzM1PsKGbl6dOn8PT0RI8ePdCqVSuoVCq0bdtW7FiUT3fu3MEPP/yAgIAAfPLJJ2LHISPGm3KIyKSpVCrUrFkTBw8eZKEpIH/99Rf69u2LlJQULFu2DD179uS8OyOVlqnGncQ0ZKm1sFTKUdHeBjZWylfbe/bsieDgYERFRaFIkSIiJiVjp8z9LURE0lWjRg18+umn2Lt3LwulgT1//hzjx4/Hjz/+iJYtW2L9+vUoX7682LHoP6Iep2BzSCyCb8UhNikdr48qyQA42FnDpVppOBV+iq1bt2LNmjUsk5QrjlASkckbOXIkduzYgXv37nGtQwO5fPkyPD09ERMTg9mzZ8PPz49fayNzLykdgXtUOBWdAIVcBo323b/+X25XxEfhr1k+qFjKtgCTkhTxp52ITJ67uzsePnyIy5cvix3F5KjVasycORP169dHoUKFcOnSJYwYMYJl0shsC41Fq0UncDYmEQDeWybf2F66CtosPo1tobGGjkgSx594IjJ5TZo0QYkSJbB3716xo5iU6OhoNGvWDJMnT8aYMWNw/vx5ODo6ih2L/mNpcBTG7VYhU63NtUj+l0YAMtVajNutwtLgKAMlJFPAQklEJk+pVKJ9+/YslHoiCAJWrVqFWrVqIS4uDqdOncLMmTNhaWkpdjT6j22hsZh/OFIv+5p/OBLbOVJJ78A5lERkFn799Vd0794dMTEx+PTTT8WOI1n//PMPfHx8cODAAQwcOBALFizgDRtG6l5SOlotOoFMtfatbalhR5F44Ic3XpNbF4NFSQcUq98VhT+rm+M+rZRyHB3ZHOXtrA0RmSSMI5REZBbatGkDS0tLjlLqYPfu3ahRowYuXbqE/fv3Y9WqVSyTRixwjwrqXC5xF2vaC/YdRsG+gz+K1u8KbXoy4nZORXr0hRzfr9YKCNyjMkRckjgWSiIyC7a2tvjqq69YKD9AcnIy+vbti65du6JZs2ZQqVTo0KGD2LHoPaIep+BUdEKucyYLV6qLIjVcUKRGSxSr3wVlPOcCciXSrp/I8f0arYBT0QmIjksxRGySMBZKIjIbbm5uOHXqFJKSksSOIhknTpxAzZo1sXv3bqxfvx67du1CqVKlxI5FudgcEguFPP+LycutbCCzsIRMrnjnexRyGTad51xKehMLJRGZDTc3N2g0Ghw4cEDsKEYvIyMDAQEBcHFxQcWKFREWFoY+ffrwiTcSEXwrLk93dGsz06BJT4YmPRlZ8XeRdGgZhKwM2Di6vPNjNFoBwZFx+oxLJoBPyiEis1G2bFk4Oztj79698PT0FDuO0bp27Ro8PT0RGRmJuXPnYuTIkVAo3j1iRcYlNVON2KT0PL03btvEN19QWMD+6+9Q+NMv3vtxsYnpSMtUv/GYRjJv/E4gIrPi7u6O2bNnIzMzE1ZWVmLHMSoajQbz58/HpEmT8L///Q+hoaGoWbOm2LEon+4mpiGvy7fYuQ6B0u4TAIAm7QnSIo4j8eBiyC0Lw7pao3d+nADgTmIaHMsW0z0wmQRe8iYis+Lu7o7U1FQcO3ZM7ChG5fbt22jRogXGjx+PkSNH4sKFCyyTEpWVwzJB72L5cVUUrlgbhSvWRhFHF5TuPgUW9uWRdGQlBE223o5Dpo+FkojMiqOjIypVqsS7vf8lCALWrVuHmjVr4v79+zh+/DjmzJnD0VsJs1R++K92mUyOQhVqQpOahOykhwY7DpkefjcQkVmRyWRwc3PDvn37oNWa9whLXFwcOnfuDB8fH3Tv3h3Xrl1Ds2bNxI5FOqpobwOdbp3SagAAQnbGO98i+/c4RC+xUBKR2XF3d8ejR49w6dIlsaOIZt++fahRowbOnj2L3377DevWrUPRokXFjkV6YGOlhMMHPslG0Kjx/PYVQKGEhX35d77Pwd6aN+TQG1goicjsNGnSBHZ2dmZ52TslJQX9+/eHu7s7GjRoAJVKBXd3d7FjkZ65VCudp3Uon8dcRGp4MFLDg/Hswh78s3E01E8eoqhzJ8itci6lCrkMLlVL6zsySRyf5U1EZql37964cuUKVCrzeYzc6dOn0bt3b8TFxeGHH36Aj48P15U0UVGPU9D6h5Pv3J7Ts7xlSkso7crB9ou2KFK73Xu/N46ObIbKpW31FZdMAMericgsubu7Y+PGjYiJiUGlSpXEjmNQWVlZmDJlCubMmYOGDRviyJEj+Oyzz8SORQZUpYwtmlYuibMxiTkucF6kZisUqdkq3/tVyGVoVMmeZZLewkveRGSW2rRpAysrK5O/7B0eHo569ephwYIFCAoKwsmTJ1kmzURQZycoP+Dxi++jlMsQ1NlJr/sk08BCSURmqUiRIvjqq69MtlBqtVosXLgQX375JbKzsxESEoJx48bxiTcmLCkpCQcPHsSUKVNQo0YNVP3EHv7N331jzYeY7uaI8h94ww+ZNl7yJiKz5ebmhqFDhyIxMRH29vZix9Gbu3fvom/fvjh+/Dj8/f0xc+ZMFCpUSOxYZADPnj3DqFGjcOzYMcTExAB4sTSWIAgoVKgQ+jatgmyFFeYfjtT5WAGu1dDD2UHn/ZBp4gglEZmtjh07QqvV4sCBA2JH0QtBELBhwwbUrFkTf//9N44dO4YFCxawTJqwlJQUbNiw4VWZBF58HwDAtGnTYGVlBV+XKpjdxQlWSnme7vx+nUIug5VSjjldnDDMpbJes5Np4V3eRGTW6tevj/Lly+PXX38VO4pOEhISMHjwYOzatQteXl5YvHgxihcvLnYsKgCbNm2Cl5fXG69ZW1vj4cOHKFbs/5+1fS8pHYF7VDgVnQCFXJbjzTovvdzetHJJBHV24mVuyhUveRORWXN3d0dQUBAyMjIkO5J38OBB9OvXD1lZWdi5cye6desmdiQqQCVKlICFhQWys188e1uhUGDw4MFvlEkAKG9njY0+9RH1OAWbQ2IRHBmH2MR0vF4rZXixaLlL1dLwbODAu7kpzzhCSURmLSIiAjVq1MAff/yBr7/+Wuw4+ZKamorRo0dj1apVaNeuHdauXYuPP/5Y7FhUQARBwKxZszBx4kS0b98eKSkpOHHiBBQKBW7fvo3y5XO/ISctU407iWnIUmthqZSjor0Nn4BDH4TfNURk1j7//HN89tln2Lt3r6QK5fnz5+Hl5YWHDx9ixYoVGDRoEBcpNyNpaWno168fduzYgcmTJ2PKlClITk6Gs7MzWrRokacyCbx4TKNj2WK5v5EoFxyhJCKzN2rUKGzZsgUPHjyAXG7c9ypmZ2dj+vTpCAoKgrOzMzZu3IgqVaqIHYsK0J07d9CpUydER0djw4YN6NKly6tt2dnZkMvlXB6KCpxxnzmJiAqAm5sb/vnnH4SGhood5b1u3LiBBg0aYPbs2Zg6dSpOnz7NMmlmjh8/jrp16+LZs2c4d+7cG2USACwsLFgmSRQslERk9ho3bgw7Ozvs27dP7Cg50mq1WLx4MerUqYP09HScO3cOkyZNglLJWUvmQhAELFmyBK1atULt2rURGhoKJyc+sYaMBwslEZk9pVKJDh06GOVTc+7fvw9XV1d89913GDhwIC5fvoy6deuKHYsKUGZmJvr37w8/Pz/4+fnhzz//NKmF+Mk0sFASEeHF8kERERH4+++/xY7yytatW+Hk5ISbN2/i8OHD+PHHH1G4cGGxY1EBevToEVq0aIHNmzfjl19+wcKFCzkyTUaJhZKICICrqyusrKyMYpQyKSkJHh4e6NmzJ9q1aweVSoXWrVuLHYsKWEhICOrWrYvY2FicPHkSvXv3FjsS0TuxUBIRAShSpAhatWqFHTt2YPny5XB1dcWCBQsKPMfhw4fh5OSEQ4cOYevWrdiyZQtKlChR4DlIXOvXr0ezZs1QoUIFXLx4EfXq1RM7EtF7cdyciMxedHQ0tm7dimvXruH+/fu4cOECBEHARx99VGAZ0tPTMXbsWCxduhStW7fGunXrUK5cuQI7PhkHtVqN0aNH48cff4SPjw+WLVsGKysrsWMR5YqFkojMXvv27REZGfnq34IgQKFQwNHRsUCOHxoaCi8vL9y9exeLFy/GsGHDjH49TNK/xMREfPPNNzh58iSWLl2KoUOHcrF6kgyesYjI7C1atAhKpfKNX94ajQaff/65QY+bnZ2NadOmoWHDhihSpAiuXLmC4cOHs0yaobCwMDg7OyMsLAxHjhzBsGHDWCZJUnjWIiKz9/XXX2Pbtm1vvW7IEcrIyEg0adIEM2bMwIQJE3Du3DlUr17dYMcj4/Xrr7+iYcOGKFasGC5evIgWLVqIHYko31goiYgAdO3aFRs2bHj1bwsLC1SsWFHvxxEEAcuXL0ft2rXx5MkTnDlzBtOmTYOFhYXej0XGTavVYtKkSejevTs6duyIM2fOoEKFCmLHIvogLJRERP/y9PTEihUrAAC2trZ6v/T88OFDtGvXDsOGDUPfvn1x5coV1K9fX6/HIGl49uwZOnXqhJkzZ2LWrFnYunUrrK2txY5F9MF4Uw4R0WsGDx6M6OjoV3fWpmWqcScxDVlqLSyVclS0t4GNVf5PnTt37sTgwYNhaWmJAwcOoF27dvqOThIRGRmJTp064eHDh/j999/x9ddfix2JSGcyQRAEsUMQERmTqMcp2BwSi+BbcYhNSsfrJ0kZAAc7a7hUK41e9R1QpYzte/f19OlT+Pr6YvPmzejWrRtWrlzJx+aZsT///BMeHh746KOPsHfvXlSrVk3sSER6wUJJRPSve0npCNyjwqnoBCjkMmi07z49vtzetHJJBHV2Qnm7ty9XHjt2DH369MGzZ8+wbNky9OrVi3fumilBEDBv3jyMGzcO7du3x6ZNm1CsWDGxYxHpDedQEhEB2BYai1aLTuBsTCIAvLdMvr79bEwiWi06gW2hsa+2PX/+HCNHjsRXX32FKlWqQKVSwdPTk2XSTKWnp6NXr14YO3YsAgMDsXfvXpZJMjkcoSQis7c0OArzD0fm/sZcjHatisbFU+Hp6Ym///4bs2bNwnfffcd1Jc1YbGwsOnXqhFu3bmH9+vXo3r272JGIDIKFkojM2rbQWIzbrdLb/p78uRSV8A82bdpUYE/aIeN08uRJdOvWDTY2Nvjtt99Qq1YtsSMRGQwLJRGZrXtJ6Wi16AQy1dq3tqWGHUXigR/e+bEfec2H1Sf/XYhcgAICDo9ohs/K8JKmuRIEAStXroSfnx+aNm2KHTt2oGTJkmLHIjIoLhtERGYrcI8K6lzmShZr2gvKYh+99bqyxMc5vFsGyOWY+vtNbPTh+pLmKCsrC76+vli9ejX8/Pwwf/58LlpPZoGFkojMUtTjFJyKTsj1fYUr1YXVx1XyvF+NVsCp6AREx6Wgcun3LylEpuWff/5B165dcfHiRaxbtw7e3t5iRyIqMJwpTkRmaXNILBRyw9x1rZDLsOl8bO5vJJNx8eJF1K1bF7dv38aJEydYJsnssFASkVkKvhWX69JAAKDNTIMmPfnN/3v+7L0fo9EKCI6M01dUMnIbN25EkyZNUK5cOVy8eBENGjQQOxJRgeMlbyIyO6mZasQmpefpvXHbJr79osICFQL2vPfjYhPTkZap/qDHNJI0qNVqjB07FgsXLoS3tzeWL1+OQoUKiR2LSBQ80xGR2bmbmIa8Lm9h5zoESrtP3nhNJsv94o4A4E5iGhzL8m5vU5SUlAQPDw8cO3YMixcvhq+vLxeuJ7PGQklEZicrh2WC3sXy46r5uinnQ49D0hEeHg53d3ckJyfj8OHDaNmypdiRiETHOZREZHYslQVz6iuo41DB2bNnDxo0aIAiRYogNDSUZZLoXzzbEZHZqWhvA0NfnJT9exwyDVqtFlOnTkWXLl3w9ddf4+zZs/j000/FjkVkNFgoicjs2Fgp4WBnbdBjONhb84YcE5GSkoIuXbpg+vTp+P7777F9+3bY2PCPBaLX8WxHRGbJpVppbAy5m+vSQc9jLiI78f5br1uV+x8sir/9BB3gxTqULlVL6yUniSs6Ohru7u64d+8e9u3bhw4dOogdicgosVASkVnqVd8B68/dyfV9yac25/i6/dcj3lkoNVoBng0cdIlHRuDQoUPw8PBA6dKlceHCBVSv/t9ntxPRSyyURGSW/r5yBsLD65CXrQ4hh9k/RWq2QpGarfK9X4VchkaV7PnYRQkTBAELFizA2LFj0aZNG2zZsgXFixcXOxaRUeMcSiIyK0+fPkW/fv3Qrl07VEoMgYVCodf9K+UyBHV20us+qeA8f/4cXl5eCAgIwJgxY7B//36WSaI8YKEkIrPx+++/w9HREbt27cLq1atxbP9OTHevoddjTHdzRHkD3/BDhnHv3j00bdoUu3fvxtatWzFr1iwo9PwHB5GpYqEkIpOXlJQELy8vdOzYEbVq1UJ4eDj69+8PmUwGD2cHjHatqpfjBLhWQw9nzp2UotOnT6Nu3bqIj4/HmTNn4OHhIXYkIklhoSQik7Znzx58/vnn+P3337F+/Xr88ccfKF++/Bvv8XWpgtldnGCllEMhz98KlQq5DFZKOeZ0ccIwl8r6jE4FZNWqVWjZsiWqV6+Oixcv4osvvhA7EpHksFASkUmKj4+Hh4cHunTpgnr16iEiIgJ9+vR55/OWPZwdcHRkcziWfHGvYm7F8uX2RpXscXRkc45MSlBWVhaGDBmCwYMHY+DAgTh69ChKlSoldiwiSeJd3kRkcnbu3Ilhw4ZBo9Fg8+bN+Pbbb99ZJF8SBAHL5k7H/jlzMCzwe5Ss3wnBkXGITUzH6ytVyvBi0XKXqqXh2cCBd3NL1OPHj9GtWzeEhIRg9erV6N+/v9iRiCSNhZKITMbjx48xbNgw7Nq1C126dMGyZcvw0Uc5rxX5uqysLAwaNAjr168HAJSzVWCcmyOmwhFpmWrcSUxDlloLS6UcFe1t+AQcibt06RI6deqE7OxsHD9+HI0aNRI7EpHk8axIRJInCAK2bt0KPz8/yGQybN++Hd27d891VBJ4ccNOp06dcObMGQCAXC5HRkbGq+02Vko4li1msOxUsLZs2QIfHx84OTlh9+7dKFeunNiRiEwC51ASkaQ9evQInTp1Qq9evdCqVStcv34d33zzTZ7KZHR0NOrWrYuzZ89Cq9UCeFEonz59auDUVNA0Gg0CAgLQq1cvfPPNNzh58iTLJJEecYSSiCRJEARs2LABI0aMgKWl5avL3Pkxa9Ys3L59+639Jicn6zMqiezJkyf49ttvceTIESxatAjfffddnv7gIKK84wglEUnO/fv30aFDB/Tt2xft27fH9evX810mAWDBggWYN2/eqzt75XI5NBoNRyhNyPXr11GvXj1cuHABhw4dwogRI1gmiQyAhZKIJEMQBKxduxaOjo64cuUK9u3bh02bNsHe3v6D9le8eHGMHj0azZs3R9myZdG0aVMAQGZmpj5jk0j27t2L+vXro1ChQggNDUWrVvl/NjsR5Q0LJRFJQmxsLNq2bYv+/fujS5cuiIiIQMeOHXXe7/3797Fnzx6MGzcOx48fx82bN7FixQo9JCaxaLVaTJ8+HZ06dYKrqyvOnTuHzz77TOxYRCaNcyiJyKhptVr89NNPCAgIQPHixXHgwAG0a9dOb/tfuXIlChcujD59+gAAqlWrprd9U8FLTU1Fnz59sHv3bkyfPh0TJkyAXM6xEyJDY6EkIqN1+/Zt+Pj4IDg4GAMGDMC8efNQrJj+lvDJyMjATz/9BG9vbxQtWlRv+yVxxMTEwN3dHXfu3MFvv/0Gd3d3sSMRmQ3+2UZERker1WLp0qVwcnLC33//jcOHD+Onn37Sa5kEgB07diA+Ph6+vr563S8VvKNHj8LZ2RkZGRk4f/48yyRRAWOhJCKjEh0dDRcXFwwfPhy9e/dGeHg4WrdurffjCIKAxYsXo02bNqhatare908FQxAELFq0CG3atEHdunVx4cIFODo6ih2LyOywUBKRUdBoNFi0aBFq1qyJe/fu4dixY1i+fDlsbQ3zrOzz58/j0qVLGD58uEH2T4aXkZGBvn37wt/fH6NGjcKBAwdQokQJsWMRmSWZIAiC2CGIyLzdunUL3t7eOHfuHPz8/BAUFAQbGxuDHrNnz564cOECIiMjedOGBD148ACdO3eGSqXCmjVr0KtXL7EjEZk1nkWJSDRqtRpz585FrVq1EB8fj5MnT+LHH380eJl89OgRdu7cCV9fX5ZJCTp79izq1q2LR48e4fTp0yyTREaAZ1IiEkVERAQaNWqEcePGYdiwYbh27dqrhcUNbdWqVbCysoK3t3eBHI/0Z82aNWjRogUqV66Mixcv4ssvvxQ7EhGBhZKICpharUZQUBDq1KmDlJQUnDlzBgsWLIC1tXWBHD8rKwsrV65E79699X7XOBlOdnY2fH19MWDAAPj4+OCvv/5CmTJlxI5FRP/iOpREVGDCwsLg7e2Nq1evIiAgAFOnTkWhQoUKNMPOnTvx+PFjLhUkIfHx8ejevTvOnDmDlStXYtCgQWJHIqL/4E05RGRwWVlZmDVrFmbOnImqVavi559/hrOzsyhZGjRoAFtbWxw5ckSU41P+XL16Fe7u7sjIyMCuXbvQpEkTsSMRUQ54yZuIDOrKlSuoV68eZsyYgbFjx+LSpUuilckLFy4gJCSESwVJxPbt29GoUSOULFkSFy9eZJkkMmIslERkEJmZmZg0aRKcnZ0hCAIuXLiAGTNmwMrKSrRMS5YsQcWKFdG+fXvRMlDuNBoNxo8fDw8PD3Tp0gWnT59G+fLlxY5FRO/BOZREpHehoaHw9vbGrVu3MGnSJIwfPx6WlpaiZnr8+DG2b9+OoKAgKBQKUbOYurRMNe4kpiFLrYWlUo6K9jawscrbr5unT5+iZ8+eOHToEObNm4dRo0ZBJpMZODER6YqFkoj0JiMjA1OnTsW8efNQu3ZtXLp0CTVr1hQ7FgDgp59+goWFBXx8fMSOYpKiHqdgc0gsgm/FITYpHa9PzpcBcLCzhku10uhV3wFVyuT89KObN2/C3d0dcXFxOHDgANq0aVMg2YlId7wph4j04ty5c+jXrx9iYmIwZcoUBAQEwMLCQuxYAF4sOVOxYkV06NABq1atEjuOSbmXlI7APSqcik6AQi6DRvvuXykvtzetXBJBnZ1Q3u7/l4r6/fff0atXL5QrVw579+5F5cqVCyI+EekJ51ASkU6eP3+O0aNHo3HjxrC1tcXly5cRGBhoNGUSAHbv3o2HDx9yqSA92xYai1aLTuBsTCIAvLdMvr79bEwiWi06gW2hsRAEAUFBQXBzc4OLiwvOnz/PMkkkQRyhJKIPdvr0afTr1w+xsbGYPn06/P39oVQa30yaxo0bw9LSEsHBwWJHMRlLg6Mw/3Ckzvspl3gZZ1ZPxuTJkzFlyhQ+CpNIoozvzE9ERi8tLQ2BgYFYsmQJGjRogH379qF69epix8rR5cuXcfbsWezevVvsKCZjW2isXsokANy3r4OAlXswbVAnveyPiMTBQklE+XL8+HH4+Pjg0aNHWLBgAfz8/Iz6ruklS5bAwcEBHTt2FDuKSbiXlI4p+yLeuT0r7g6Sz2xB5qMoaNKeQlHYFhYlHVC4cn0UrZvzf4O9960wPCn9jTmVRCQtvLZARHmSkpKCoUOHwsXFBZ988gmuXbuGkSNHGnWZjI+Px9atWzF06FCjvBQvRYF7VFC/Y65kxv0bePTLCGTF3UaRWm1g5zoYRWq1AWQypFzc9859qrUCAveoDBWZiAoAz7BElKujR4+if//+iI+Px+LFizFs2DBJzHVbvXo1ZDIZ+vfvL3YUkxD1OAWnohPeuf3Zue2QW9ng4z6LIC9U5I1tmrSn7/w4jVbAqegERMeloHLpnJcUIiLjZvy/EYhINMnJyRg4cCBat26NSpUqQaVSYfjw4ZIok2q1GitWrECvXr1gb28vdhyTsDkkFgr5uxcZz37yDyxKOrxVJgFAYVP8vftWyGXYdD5W14hEJBLj/61ARKL4888/UaNGDWzduhUrVqzA0aNHUalSJbFj5dlvv/2G+/fv87ndehR8K+69SwMpi5VC1j/RyIq/k+99a7QCgiPjdEhHRGLiskFE9IYnT57A398f69evR+vWrbF69WpUqFBB7Fj51rx5cwiCgJMnT4odxSSkZqrhNPUQ3vcL4/ntK4jbMQUAYFW2KqzKOaJQxVoo5FATMkXuM6xkAMKntsnzYxqJyHjwp5aIXtm/fz8GDRqEtLQ0rF69Gj4+PpJ8jvK1a9dw8uRJ7NixQ+woJuNuYtp7yyQAFP70C3zUez6Sz+1Exu3LyHxwE89CdkFuXQz27fxgXaX+ez9eAHAnMQ2OZYvpLTcRFQxe8iYiJCYmwtPTE25ubqhduzYiIiLQv39/SZZJ4MVSQZ988gk6deokdhSTkaXW5ul9Vh9XRekuE1B+xDZ81GchijbsDiHrOeL3zEJWQu5zJPN6HCIyLiyURGZuz549cHR0xB9//IH169fjjz/+QLly5cSO9cESExOxefNmDB061Kge/yh1lsr8/bqQKSxg9XFVlGjeB3auQwGtGuk3T+v9OERkHPiTS2Sm4uPj4eHhgS5duqB+/fqIiIhAnz59JDsq+dLatWshCAIGDBggdhSTUtHeBh/6nWH58Ytnc2tSk977Ptm/xyEi6WGhJDJDO3fuhKOjI44cOYLNmzfjt99+Q9myZcWOpTO1Wo1ly5bBw8MDpUqVEjuOSbGxUsIhlyfZZNwNQ073eT7/+yIAwMLu/SPfDvbWvCGHSKL4k0tkRh4/foxhw4Zh165d6NKlC5YvX44yZcqIHUtv9u/fj9jYWC4VpGcv75bPvHMZ8qLVoH3HWGXSkZUQsjNRuGpDWNiXAzRqZDy4gfQbp6AoVgZFarZ65zEUchlcqpY21KdARAbGZYOIzIAgCNi6deurRcmXLVuG7t27S/7y9n+1bNkSmZmZOHPmjNhRTEJWVhZ27tyJhQsX4vLly6herwWetxz9zvc/j7mEtJunkfngBjQpiRA02VAWLYXCleqiWKMeuS5ufnRkMz4ph0iiOEJJZOIePXqEwYMHY9++fejRoweWLFlikpeDw8PDERwcjK1bt4odRfKSkpLw008/YenSpXjw4AFcXV3x559/wtXVFb3XXcDZmMQcFzgvXOlLFK70Zb6Pp5DL0KiSPcskkYRxhJLIRAmCgA0bNmDEiBGwsrLCihUr0LlzZ7FjGcygQYOwf/9+3L17l3d3f6CoqCj8+OOP+Pnnn6HRaODp6YkRI0agRo0ar95zLykdrRadQKYel/exUspxdGRzlM9ljiYRGS/elENkgu7fv48OHTqgb9++6NChAyIiIky6TD558gSbNm3C4MGDWSbzSRAEnDhxAu7u7qhWrRp27NiBgIAA3L17F2vWrHmjTAJAeTtrTHNz1GuG6W6OLJNEEsdL3kQmRBAErFu3Dv7+/ihSpAj27duHjh07ih3L4NatW4fs7GwMGjRI7CiS8d/5kZ9//jlWr16NXr16oVChQu/9WA9nBySkZmL+4UidcwS4VkMPZwed90NE4uIlbyITcffuXQwcOBCHDx9G3759sXDhQpQoUULsWAan0WhQpUoVNG7cGBs3bhQ7jtHLaX6kv78/XF1d832T1rbQWEzZFwG1VshxTuW7KOQyKOUyTHdzZJkkMhEcoSSSOK1Wi59++gkBAQEoXrw4Dhw4gHbt2okdq8D88ccfuH37NrZt2yZ2FKOWl/mR+eXh7IDGn5VE4B4VTkUnQCGXvbdYvtzeqJI9gjo78TI3kQnhCCWRhMXExKB///4IDg7GgAEDMG/ePBQrVkzsWAWqdevWePbsGUJCQsSOYnRerh+5cOFC7N+/HyVLlsTQoUMxZMgQva8/GvU4BZtDYhEcGYfYxHS8/otFhheLlrtULQ3PBg68m5vIBLFQEkmQVqvFsmXLMG7cOJQsWRJr1qxB69atxY5V4G7cuIHPP/8cGzduhKenp9hxjEZO8yP9/f3zND9SH9Iy1biTmIYstRaWSjkq2tvwCThEJo4/4UQSExUVBR8fH5w6dQpDhgzBnDlzYGtrniM+S5cuRZkyZdC9e3exoxiF960fWZCL2NtYKeFY1rxGyonMHQslkURoNBosXrwYEyZMwEcffYRjx47BxcVF7FiiSU5Oxi+//IJRo0bByspK7DiiMsT8SCKi/GChJJKAmzdvol+/fjh//jyGDx+OoKAg2NjYiB1LVD///DMyMzPNdqmgnOZHBgQEGGR+JBFRblgoiYyYWq3GwoULMXnyZDg4OODkyZNo0qSJ2LFE93IOabdu3VC2bFmx4xQoXdaPJCIyFBZKIiMVEREBb29vXLx4Ef7+/pg+fTqsrbnMCgD8+eefiI6OxoYNG8SOUmCMZX4kEVFOeJc3kZHJzs7G3LlzMX36dFSqVAnr1q1Dw4YNxY5lVNq1a4f4+HiEhoaafJni/EgikgKOUBIZkbCwMHh7e+Pq1asICAjA1KlTeRnzP27duoU///wT69evN9kyyfmRRCQ1LJRERiArKwuzZs3C999/j2rVquH8+fNwdnYWO5ZRWrZsGUqWLIkePXqIHUXvOD+SiKSKhZJIZJcvX4a3tzciIiIwfvx4TJw40eyXwXmXlJQUrF+/Hn5+fiZVsDg/koikjoWSSCSZmZmYMWMGZs+ejRo1aiA0NBRffPGF2LGM2i+//IL09HQMHjxY7Ch6wfmRRGQqeFMOkQhCQ0Ph7e2NW7duYeLEiRg/fjwsLS3FjmXUtFot/ve//6FWrVrYsWOH2HE+WEE+X5uIqKBwhJKoAGVkZGDq1KmYN28eateujUuXLqFmzZpix5KEI0eOIDIyEmvWrBE7ygfh/EgiMmUcoSQqIOfOnUO/fv0QExODKVOmICAgABYWFmLHkowOHTrgwYMHuHz5sqTmFeY0P9Lf35/zI4nIpHCEksjA0tPTMWnSJCxatAjOzs64fPkyHB0dxY4lKdHR0Thw4ADWrFkjmRLG+ZFEZE5YKIkM6NSpU+jXrx/u3buHOXPmYOTIkVAq+WOXX8uWLYOdnR2+/fZbsaO8F9ePJCJzJRc7AJEpSktLw3fffYfmzZujVKlSrxYqZ5nMv9TUVKxbtw79+/dH4cKFxY6To6ysLGzevBl169ZFixYtEB0djdWrVyM2NhZTp05lmSQik8ffbkR6dvz4cfj4+ODRo0dYsGAB/Pz8oFAoxI4lWRs3bkRqaiqGDh0qdpS3cP1IIqIXeFMOkZ6kpKRg7NixWLFiBZo2bYq1a9eiSpUqYseSNEEQ4OjoiP/973/YtWuX2HFe4fxIIqI3cYSSSA+OHj2K/v37Iz4+HkuWLMHQoUMhl3NGia7++usv3LhxA8uXLxc7CudHEhG9B0coiXSQnJyM0aNHY82aNXBxccGaNWtQqVIlsWOZDHd3d8TExCAsLEy0S8g5rR/p7+/P9SOJiF7DEUqiD3Tw4EEMHDgQT58+xYoVKzBw4ECOSurR7du3sX//fqxcuVKUMsn5kUREecdCSZRPT548gb+/P9avX4/WrVtj9erVqFChgtixTM7y5ctRvHhx9OrVq0CPy/mRRET5x0JJlA/79+/HoEGDkJaWhjVr1qBfv34crTKAl1/f/v37w8bGxuDH4/xIIiLd8PocUR4kJibC09MTbm5u+OKLLxAREQEfHx+WSQPZvHkzkpOTDb5UENePJCLSD96UQ5SL3bt3Y+jQocjMzMSPP/4ILy8vFkkDEgQBNWvWRKVKlbB3716DHIPP1yYi0i9e8iZ6h/j4ePj6+mLHjh1wc3PDihUrULZsWbFjmbwTJ04gPDwcP/zwg973zfmRRESGwRFKov8QBAE7d+7EsGHDoNVqsWTJEnz77bccudKztEw17iSmIUuthaVSjor2NrCxUqJr1664efMmwsPD9fI1z2l+5LBhwzBkyBCULl1aD58JERFxhJLoNY8fP8bQoUOxe/dudO3aFcuWLeM8Oj2KepyCzSGxCL4Vh9ikdLz+16wMwCfFrHAjpRRGDGqb5zIpCEKO781p/cjVq1dz/UgiIgPgCCURXpSSLVu2wM/PD3K5HMuXL0f37t3FjmUy7iWlI3CPCqeiE6CQy6DRvvu0I4MAATI0rVwSQZ2dUN7OOsf3CYKAvn37IjIyEmfPnn1VKv87P7JNmzbw9/dH69atOcpMRGQgLJRk9h4+fIjBgwdj//796NGjB5YsWYJSpUqJHctkbAuNxZR9EVBrhfcWyf9SyGVQymWY5uYID2eHt7Z///33mDRpEgDgyJEjqFChAudHEhGJhIWSzJYgCNiwYQNGjBgBKysrrFixAp07dxY7lklZGhyF+Ycjdd7PaNeq8HWp8urfe/bsQZcuXQAAcrkc9vb2SEhI4PxIIiKRsFCSWbp//z4GDhyIgwcPwtPTEz/88APs7e3FjmVStoXGYtxuld72N6eLE3o4OyAsLAz169dHRkbGG9unT5+OgIAAzo8kIhIBFzYnsyIIAtasWQNHR0dcu3YN+/fvx8aNG1km9exeUjqm7It45/bUsKO4O7sDMh9F5Xmfk/dF4LwqCk2bNn2rTCoUCty7d49lkohIJCyUZDbu3r2LNm3aYMCAAejatSsiIiLQoUMHsWOZpMA9KqjzMV8yL9RaAUPXn8azZ8/e2qbRaLBp0ybwggsRkTi4bBCZPK1Wi1WrVmHMmDEoXrw4Dh48iLZt24ody2RFPU7BqegEve9XoxWQZFkaB89eQZlCAlJTU5Gamoq0tDSkpqbCzs6Od3ETEYmEhZJMWkxMDPr374/g4GAMGDAA8+bNQ7FixcSOZdI2h8TmujTQh1LIZTgfb4Gpbo563zcREX04XvImk/TyCTdOTk6IiYnBkSNH8NNPP7FMFoDgW3EGKZPAi1HK4Mg4g+ybiIg+HAslmZyoqCi0aNECfn5+6Nu3L1QqFVq1aiV2LLOQmqlGbFK6QY8Rm5iOtEy1QY9BRET5w0JJJkOj0WDhwoWoVasWHjx4gODgYCxbtgy2trZiRzMbdxPTYOjbYgQAdxLTDHwUIiLKDxZKMgk3b95E06ZNMXr0aAwcOBBhYWFo0aKF2LHMTpZaa1LHISKivGGhJElTq9WYO3cuateujYSEBJw8eRI//PADbGxsxI5mliyVBXNKKajjEBFR3vCsTJIVERGBRo0aYfz48fD19cXVq1fRpEkTsWOZtYr2NjD0wj2yf49DRETGg4WSJCc7OxszZ85EnTp1kJKSgjNnzmD+/PmwtrYWO5rZs7FSwsHOsP8dHOytYWPFFc+IiIwJz8okKWFhYejbty+uXbuGMWPGYMqUKXzcnpFxqVYaG0Pu5mnpoNSwI3gec+mt14vWdYPc6u1iqpDL4FK1tF5yEhGR/rBQkiRkZWUhKCgIM2fORPXq1RESEoK6deuKHYty0Ku+A9afu5On96ZeOZDj60WcWuVYKDVaAZ4NHHSJR0REBiAT+PBbMnKXL1+Gt7c3IiIiMH78eEycOBFWVlZix6L38FobgrMxiXpd4Fwhl6FRJXts9Kmvt30SEZF+cA4lGa3MzExMnDgR9erVg0wmQ2hoKGbMmMEyKQFBnZ2glOv39hylXIagzk563ScREekHCyUZpdDQUHz55ZeYO3cuJk+ejAsXLuCLL74QOxblUXk7a0zT8/O2p7s5oryBb/ghIqIPw0JJRiUjIwNjx45FgwYNYGVlhUuXLmHy5MmwtLQUOxr9K6+zZDycHTDatapejhngWg09nDl3kojIWHEOJRmNc+fOwdvbG7dv38bUqVMREBAApZL3jUndttBYTNkXAbVWyNecSoVcBqVchulujiyTRERGjiOUJLr09HSMGjUKjRs3RrFixXDlyhWMHz+eZdJEeDg74OjI5mhUyR7Ai6L4Pi+3N6pkj6Mjm7NMEhFJAEcoSVSnTp1Cv379cO/ePcyYMQMjR45kkTRhUY9TsDkkFsGRcYhNTMfrJx8ZXixa7lK1NDwbOKByaVuxYhIRUT6xUJIo0tLSMH78eCxduhQNGzbEunXrUK1aNbFjmT1BECCTGfrhiS+kZapxJzENWWotLJVyVLS34RNwiIgkioWSClxwcDB8fHzwzz//ICgoCMOHD4dCoRA7FhEREX0gzqGkApOSkoKhQ4eiZcuWKF++PMLCwjBixAiWSSIiIonj9SUqEEeOHMGAAQOQkJCAJUuWYOjQoZDL+fcMERGRKeBvdDKo5ORkDBgwAK6urvjss8+gUqng6+vLMklERGRCOEJJBnPw4EEMHDgQycnJWLlyJQYOHFhgN3wQERFRweEwEendkydP0LdvX3z99df4/PPPER4ejkGDBrFMEhERmSiOUJJe7d+/H4MGDUJ6ejrWrl0Lb29vFkkiIiITxxFK0ovExER4enrCzc0NderUQUREBPr168cySUREZAY4Qkk62717N4YOHYrMzEz88ssv8PLyYpEkIiIyIxyhpA8WHx+PHj16oGvXrmjQoAGuX7+O3r17s0wSERGZGY5QUr4JgoAdO3bA19cXgiBgy5Yt8PDwYJEkIiIyUxyhpHx5/PgxunXrBg8PDzRv3hwRERH49ttvWSaJiIjMGEcoKU9ejkT6+flBoVBgx44d6N69u9ixiIiIyAhwhJJy9fDhQ7i7u8PT0xOurq6IiIhgmSQiIqJXOEJJ7yQIAjZs2IARI0bAysoKe/bsQadOncSORUREREaGI5SUo/v376N9+/bo27cvOnbsiOvXr7NMEhERUY5YKOkNgiBgzZo1cHR0xLVr17B//35s2LABdnZ2YkcjIiIiI8VCSa/cvXsXbdq0wYABA9C1a1dERESgQ4cOYsciIiIiI8dCSdBqtVixYgVq1KiBGzdu4ODBg1i3bh2KFy8udjQiIiKSABZKMxcTE4NWrVph6NCh6NmzJyIiItC2bVuxYxEREZGEsFCaKa1WiyVLlsDJyQkxMTE4evQoVq1ahaJFi4odjYiIiCSGhdIMRUVFoUWLFvDz84O3tzdUKhW++uorsWMRERGRRLFQmhGNRoOFCxeiZs2aePDgAYKDg7F06VLY2tqKHY2IiIgkjIXSTNy8eRNNmzbF6NGjMWjQIISFhaFFixZixyIiIiITwEJp4tRqNebMmYPatWsjISEBp06dwg8//AAbGxuxoxEREZGJYKE0YeHh4WjUqBECAwMxfPhwXLt2DY0bNxY7FhEREZkYFkoTlJ2djZkzZ6JOnTpITU3F2bNnMW/ePBQuXFjsaERERGSClGIHIP0KCwtD3759ERYWhoCAAEyZMgWFChUSOxYRERGZMI5QmoisrCxMnToVX375JbKzs3H+/HnMmjWLZZKIiIgMjiOUJuDy5cvw9vbG9evXMX78eEyYMAFWVlZixyIiIiIzwRFKCcvMzMTEiRNRr149yGQyhIaGYvr06SyTREREVKBkgiAIYoeg/AsNDYW3tzciIyMxadIkjBs3DhYWFmLHIiIiIjPEEUqJycjIwNixY9GgQQMUKlQIly5dwqRJk1gmiYiISDScQykh586dg7e3N27fvo3vv/8eAQEBUCr5n5CIiIjExRFKCUhPT8eoUaPQuHFjFC9eHFeuXMH48eNZJomIiMgosJEYuVOnTqFfv364f/8+5s6di5EjR0KhUIgdi4iIiOgVjlAaqbS0NPj5+aF58+YoU6YMrl27htGjR7NMEhERkdHhCKURCg4Oho+PD/755x8sWrQIvr6+LJJERERktDhCaURSUlIwZMgQtGzZEuXLl0dYWBi+++47lkkiIiIyahyhNBJHjhxB//79kZiYiKVLl2LIkCGQy9n3iYiIyPixsYgsOTkZAwYMgKurKypXrgyVSoVhw4axTBIREZFkcIRSRAcPHsTAgQORnJyMVatWYcCAAZDJZGLHIiIiIsoXDoOJ4MmTJ+jbty++/vprODo6Ijw8HAMHDmSZJCIiIkniCGUB279/PwYNGoT09HSsXbsW3t7eLJJEREQkaRyhLCCJiYnw9PSEm5sb6tSpg4iICPTr149lkoiIiCSPI5QFYPfu3Rg6dCiysrKwYcMGeHp6skgSERGRyeAIpQHFx8ejR48e6Nq1Kxo2bIiIiAh4eXmxTBIREZFJ4QilAQiCgB07dsDX1xeCIGDr1q3o0aMHiyQRERGZJI5Q6iAlJQXPnz9/47XHjx+jW7du8PDwgIuLC65fvw4PDw+WSSIiIjJZZj9CmZapxp3ENGSptbBUylHR3gY2Vrl/WdRqNerXrw9LS0uEhoZCqVRiy5Yt8PPzg0KhwM6dO9GtW7cC+AyIiIiIxGWWhTLqcQo2h8Qi+FYcYpPSIby2TQbAwc4aLtVKo1d9B1QpY5vjPlauXIkbN25AJpNh4sSJuHHjBvbv349vv/0WixcvRsmSJQvkcyEiIiISm0wQBCH3t5mGe0npCNyjwqnoBCjkMmi07/7UX25vWrkkgjo7obyd9attiYmJqFSpEp49e/bqNXt7e6xduxbu7u4G/RyIiIiIjI3ZzKHcFhqLVotO4GxMIgC8t0y+vv1sTCJaLTqBbaGxr7ZNmTIFqampr/4tk8lQtmxZtG/f3gDJiYiIiIybWYxQLg2OwvzDkTrvZ7RrVbQolYmaNWsipy/bggUL4O/vr/NxiIiIiKTE5OdQbguN1UuZBID5hyOx6OLWt8qkXC5HmTJloFSa/JeTiIiI6C0mPUJ5LykdrRadQKZa+8brCfvnI+3mGZT1WQoLu0/e2JZ8bieenvgFpbpNhnXlem/tUwEt6sUfwtfN66NChQqoUKECPvnkE5ZJIiIiMlsm3YIC96igzmGuZImW/fH874tI/HMZPuoZ9Or17Kf/IPnMNlhXa5RjmQQAyBVQNvSCl1d9Q8UmIiIikhSTvSkn6nEKTkUn5HjzjcKmOIq36IvM2DCkqv569XrS4RWAQoESrQa+c78arYBT0QmIjksxSG4iIiIiqTHZQrk5JBYK+bufTlOkVhtYlfscT46theb5M6RdP4GMmEso3tQLStv3ryGpkMuw6Xzse99DREREZC5MtlAG34p779JAMpkMdm2GQZuZjqRDy/HkrzWw/KgKbL/MfekfjVZAcGScPuMSERERSZZJFsrUTDVik9JzfZ9lqQooWr8z0m+ehiY9GXZth0Emy9uXJDYxHWmZal2jEhEREUmeSRbKu4lpyOut64rCRV/8f1t7WJaqkOdjCADuJKblPxwRERGRiTHJQpn1n2WC3kX9LB5PT2+BRakK0DyLR/L5XQY5DhEREZEpM8lCaanM26eVdHglAKB092mwrt4Ez87tQPbTf/R+HCIiIiJTZpKNqKK9Dd59f/cL6bfO4nl0CIo39YSyaEmU+GoAoFC+WDooD2T/HoeIiIjI3JlkobSxUsLBzvqd27WZ6Ug6+hMsy3wG2y87AACUtvYo3tQTGTGXkHbzdK7HcLC3ho2VSa8LT0RERJQnJlkoAcClWul3rkP59ORGaFKTXtzVLVe8et22TntYlvkMT47+BG3mu+8SV8hlcKlaWu+ZiYiIiKTIZAtlr/oOOa5DmflPNFIu/wHbL76G1cdV39gmkytg13YYNGlP8fTkxnfuW6MV4NnAQe+ZiYiIiKRIJghCXlfYkRyvtSE4G5P43gXO80shl6FRJXts9OGzvImIiIgAEx6hBICgzk5Qvufxix9CKZchqLOTXvdJREREJGUmXSjL21ljmpujXvc53c0R5d9zww8RERGRuTHpQgkAHs4OGO1aNfc35kGAazX0cObcSSIiIqLXmfQcytdtC43FlH0RUGuFfM2pVMhlUMplmO7myDJJRERElAOzKZQAcC8pHYF7VDgVnQCFXPbeYvlye9PKJRHU2YmXuYmIiIjewawK5UtRj1OwOSQWwZFxiE1Mx+tfABleLFruUrU0PBs4oHJpW7FiEhEREUmCWRbK16VlqnEnMQ1Zai0slXJUtLfhE3CIiIiI8sHsCyURERER6cbk7/ImIiIiIsNioSQiIiIinbBQEhEREZFOWCiJiIiISCcslERERESkExZKIiIiItIJCyURERER6YSFkoiIiIh0wkJJRERERDphoSQiIiIinbBQEhEREZFOWCiJiIiISCcslERERESkExZKIiIiItIJCyURERER6YSFkoiIiIh0wkJJRERERDphoSQiIiIinbBQEhEREZFOWCiJiIiISCcslERERESkExZKIiIiItIJCyURERER6YSFkoiIiIh0wkJJRERERDphoSQiIiIinbBQEhEREZFOWCiJiIiISCcslERERESkExZKIiIiItIJCyURERER6YSFkoiIiIh0wkJJRERERDphoSQiIiIinbBQEhEREZFOWCiJiIiISCcslERERESkk/8DlIH8+cL3FkcAAAAASUVORK5CYII=\n",
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "pos = nx.spring_layout(ground_truth, seed=1234)\n",
+ "nx.draw(ground_truth, with_labels=True, pos=pos)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "ec0242ea",
+ "metadata": {},
+ "source": [
+ "The variables in the DAG have the following interpretation:\n",
+ "\n",
+ "* T: Whether or not the patient has **tuberculosis**.\n",
+ "* L: Whether or not the patient has **lung cancer**.\n",
+ "* B: Whether or not the patient has **bronchitis**.\n",
+ "* A: Whether or not the patient has recently visited **Asia**.\n",
+ "* S: Whether or not the patient is a **smoker**.\n",
+ "* E: An indicator of whether the patient has either lung cancer or tuberculosis (or both).\n",
+ "* X: Whether or not a chest X-ray shows evidence of tuberculosis or lung cancer.\n",
+ "* D: Whether or not the patient has **dyspnoea** (difficulty breathing).\n",
+ "\n",
+ "Note we have three kinds of variables, diseases (B, L, T, and E which indicates one or more diseases), symptoms (X and D), and behaviors (S and A). The goal of the model is to use symptoms and behaviors to diagnose (i.e. infer) diseases. Further, note that diseases are causes of symptoms, and are effects of behaviors."
+ ]
+ },
{
"cell_type": "markdown",
"id": "b769b7f1",
@@ -229,7 +263,7 @@
},
{
"cell_type": "code",
- "execution_count": 4,
+ "execution_count": 5,
"id": "30a29f1d",
"metadata": {},
"outputs": [],
@@ -274,7 +308,7 @@
},
{
"cell_type": "code",
- "execution_count": 5,
+ "execution_count": 6,
"id": "902d126f",
"metadata": {},
"outputs": [],
@@ -298,7 +332,7 @@
},
{
"cell_type": "code",
- "execution_count": 6,
+ "execution_count": 4,
"id": "89a4f0af",
"metadata": {},
"outputs": [],
@@ -318,7 +352,7 @@
},
{
"cell_type": "code",
- "execution_count": 7,
+ "execution_count": 14,
"id": "e5cba859",
"metadata": {},
"outputs": [],
@@ -337,25 +371,156 @@
},
{
"cell_type": "code",
- "execution_count": 8,
+ "execution_count": 56,
"id": "918c6f5f",
"metadata": {},
"outputs": [
{
"data": {
- "image/svg+xml": "\n\n\n\n\n",
+ "image/svg+xml": [
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n"
+ ],
"text/plain": [
- ""
+ ""
]
},
- "execution_count": 8,
+ "execution_count": 56,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"graph = pc.graph_\n",
- "draw(graph)"
+ "\n",
+ "draw(graph, direction='TB')"
]
},
{
@@ -363,29 +528,182 @@
"id": "162d98be",
"metadata": {},
"source": [
- "Compare this against the ground truth CPDAG."
+ "Compare this against the ground truth CPDAG, which should match exactly."
]
},
{
"cell_type": "code",
- "execution_count": 9,
+ "execution_count": 57,
"id": "67d47367",
"metadata": {},
"outputs": [
{
"data": {
- "image/svg+xml": "\n\n\n\n\n",
+ "image/svg+xml": [
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n"
+ ],
"text/plain": [
- ""
+ ""
]
},
- "execution_count": 9,
+ "execution_count": 57,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
- "draw(ground_truth_cpdag)"
+ "draw(ground_truth_cpdag, direction='TB')"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 58,
+ "id": "0569e7c4-a571-484c-a07b-e6935c55fd39",
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "The oracle learned CPDAG via the PC algorithm matches the ground truth in directed edges: True \n",
+ "and matches the undirected edges: True\n"
+ ]
+ }
+ ],
+ "source": [
+ "match_directed = nx.is_isomorphic(ground_truth_cpdag.sub_directed_graph(), graph.sub_directed_graph())\n",
+ "match_undirected = nx.is_isomorphic(ground_truth_cpdag.sub_undirected_graph(), graph.sub_undirected_graph())\n",
+ "\n",
+ "print(f'The oracle learned CPDAG via the PC algorithm matches the ground truth in directed edges: {match_directed} \\n'\n",
+ " f'and matches the undirected edges: {match_undirected}')"
]
},
{
@@ -393,18 +711,29 @@
"id": "db1ff244",
"metadata": {},
"source": [
- "Now, we will show the output given a real CI test, which performs CI hypothesis testing to determine CI in the data. Due to finite data and the presence of noise, there is always a possibility that the CI test makes a mistake."
+ "Now, we will show the output given a real CI test, which performs CI hypothesis testing to determine CI in the data. Due to finite data and the presence of noise, there is always a possibility that the CI test makes a mistake. In order to maximize the chances that the graph is correct, you want to ensure that the CI test you are using matches the assumptions you have on your data. \n",
+ "\n",
+ "For example, the G^2 binary test is a well-suited test for binary data, which we validated is the type of data we have for the ASIA dataset."
]
},
{
"cell_type": "code",
- "execution_count": 10,
+ "execution_count": 9,
"id": "5b80aeba",
"metadata": {},
"outputs": [],
"source": [
- "ci_estimator = GSquareCITest(data_type=\"discrete\")\n",
- "pc = PC(ci_estimator=ci_estimator)\n",
+ "ci_estimator = GSquareCITest(data_type=\"binary\")\n",
+ "pc = PC(ci_estimator=ci_estimator, alpha=0.05)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 10,
+ "id": "64a5d4a9-1e07-4282-813c-a080322eaff0",
+ "metadata": {},
+ "outputs": [],
+ "source": [
"pc.fit(data, context)"
]
},
@@ -416,9 +745,104 @@
"outputs": [
{
"data": {
- "image/svg+xml": "\n\n\n\n\n",
+ "image/svg+xml": [
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n"
+ ],
"text/plain": [
- ""
+ ""
]
},
"execution_count": 11,
@@ -428,7 +852,7 @@
],
"source": [
"graph = pc.graph_\n",
- "draw(graph)"
+ "draw(graph, direction='TB')"
]
},
{
@@ -436,7 +860,7 @@
"id": "ff0b314b",
"metadata": {},
"source": [
- "The resulting graph captures some of the graph but not all of it. The problem here is a violation of the [faithfulness assumption](https://plato.stanford.edu/entries/causal-models/#MiniFaitCond); in the Asia data, it is very hard to detect the edge between E and X.\n",
+ "The resulting graph captures some of the graph but not all of it. The problem here is a violation of the [faithfulness assumption](https://plato.stanford.edu/entries/causal-models/#MiniFaitCond); e.g. in the Asia data, it is very hard to detect the edge between E and X. This highlights a common problem with causal discovery, where the inability to detect certain edges may lead to incorrect orientations.\n",
"\n",
"Beyond faithfulness violations, in general causal discovery algorithms struggle when there is no user-provided causal knowledge to constrain the problem. A core philosophy of dodiscovery is that causal domain knowledge should be provided to constrain the problem, and providing it should be easy. For example, for this data, we know that smoking (S) causes both lung cancer (L) and bronchitis (B) and not the other way around.\n",
"\n",
@@ -448,9 +872,9 @@
],
"metadata": {
"kernelspec": {
- "display_name": "venv",
+ "display_name": "pywhy-discover",
"language": "python",
- "name": "venv"
+ "name": "pywhy-discover"
},
"language_info": {
"codemirror_mode": {
@@ -462,11 +886,11 @@
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
- "version": "3.10.9"
+ "version": "3.9.13"
},
"vscode": {
"interpreter": {
- "hash": "83bc06ec5be13f4230f46a3f77f7cefbb44c2fa59a857bbc121b0c3cdb0063f8"
+ "hash": "f69a7104467f431c4bacbebec40c4cb5787ef707a55bea5c5fb34f2af39396ab"
}
}
},
diff --git a/doc/whats_new/_contributors.rst b/doc/whats_new/_contributors.rst
index 14a98c41..4d8bbc02 100644
--- a/doc/whats_new/_contributors.rst
+++ b/doc/whats_new/_contributors.rst
@@ -22,3 +22,4 @@
.. _Adam Li: https://adam2392.github.io
.. _Chris Trevino: https://py-why.github.io
+.. _Robert Osazuwa Ness: https://py-why.github.io
\ No newline at end of file
diff --git a/doc/whats_new/v0.1.rst b/doc/whats_new/v0.1.rst
index 2bb8cd56..269c8bfe 100644
--- a/doc/whats_new/v0.1.rst
+++ b/doc/whats_new/v0.1.rst
@@ -37,7 +37,7 @@ Changelog
- |Feature| Implement FCI algorithm, :class:`dodiscover.constraint.FCI` for learning causal structure from observational data with latent confounders under the ``dodiscover.constraint`` submodule, by `Adam Li`_ (:pr:`52`)
- |Feature| Implement Structural Hamming Distance metric to compare directed graphs, :func:`dodiscover.metrics.structure_hamming_dist`, by `Adam Li`_ (:pr:`55`)
- |Fix| Update dependency on networkx, which removes a PR branch dependency with pywhy-graphs having the MixedEdgeGraph class that was causing a dependency conflict, by `Adam Li`_ (:pr:`74`)
-- |Chore| Add tutorial for PC algorithm with Asia data, by `Robert Osazuwa Ness`_ (:pr:`67`)
+- |Enhancement| Add tutorial for PC algorithm with Asia data, by `Robert Osazuwa Ness`_ (:pr:`67`)
Code and Documentation Contributors
-----------------------------------
@@ -47,3 +47,4 @@ the project since version inception, including:
* `Adam Li`_
* `Chris Trevino`_
+* `Robert Osazuwa Ness`_
diff --git a/dodiscover/constraint/pcalg.py b/dodiscover/constraint/pcalg.py
index e4b79603..ec9e6020 100644
--- a/dodiscover/constraint/pcalg.py
+++ b/dodiscover/constraint/pcalg.py
@@ -1,5 +1,5 @@
import logging
-from itertools import combinations, permutations
+from itertools import combinations
from typing import Optional
import networkx as nx
@@ -134,37 +134,36 @@ def orient_edges(self, graph: EquivalenceClass) -> None:
A skeleton graph. If ``None``, then will initialize PC using a
complete graph. By default None.
"""
- node_ids = graph.nodes
-
# For all the combination of nodes i and j, apply the following
# rules.
idx = 0
finished = False
while idx < self.max_iter and not finished: # type: ignore
change_flag = False
- for (i, j) in permutations(node_ids, 2):
- if i == j:
- continue
- # Rule 1: Orient i-j into i->j whenever there is an arrow k->i
- # such that k and j are nonadjacent.
- r1_add = self._apply_meek_rule1(graph, i, j)
-
- # Rule 2: Orient i-j into i->j whenever there is a chain
- # i->k->j.
- r2_add = self._apply_meek_rule2(graph, i, j)
-
- # Rule 3: Orient i-j into i->j whenever there are two chains
- # i-k->j and i-l->j such that k and l are nonadjacent.
- r3_add = self._apply_meek_rule3(graph, i, j)
-
- # Rule 4: Orient i-j into i->j whenever there are two chains
- # i-k->l and k->l->j such that k and j are nonadjacent.
- #
- # However, this rule is not necessary when the PC-algorithm
- # is used to estimate a DAG.
-
- if any([r1_add, r2_add, r3_add]) and not change_flag:
- change_flag = True
+ for i in graph.nodes:
+ for j in graph.neighbors(i):
+ if i == j:
+ continue
+ # Rule 1: Orient i-j into i->j whenever there is an arrow k->i
+ # such that k and j are nonadjacent.
+ r1_add = self._apply_meek_rule1(graph, i, j)
+
+ # Rule 2: Orient i-j into i->j whenever there is a chain
+ # i->k->j.
+ r2_add = self._apply_meek_rule2(graph, i, j)
+
+ # Rule 3: Orient i-j into i->j whenever there are two chains
+ # i-k->j and i-l->j such that k and l are nonadjacent.
+ r3_add = self._apply_meek_rule3(graph, i, j)
+
+ # Rule 4: Orient i-j into i->j whenever there are two chains
+ # i-k->l and k->l->j such that k and j are nonadjacent.
+ #
+ # However, this rule is not necessary when the PC-algorithm
+ # is used to estimate a DAG.
+
+ if any([r1_add, r2_add, r3_add]) and not change_flag:
+ change_flag = True
if not change_flag:
finished = True
logger.info(f"Finished applying R1-3, with {idx} iterations")
diff --git a/pyproject.toml b/pyproject.toml
index 20d97b43..1f39eb94 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -48,7 +48,6 @@ importlib-resources = { version = "*", python = "<3.9" }
networkx = "^2.8.8"
pywhy-graphs = { git = "https://github.com/py-why/pywhy-graphs.git", branch = 'main', optional = true }
pygraphviz = { version = "*", optional = true }
-bnlearn = { git = "https://github.com/erdogant/bnlearn.git", branch = 'master', optional = true }
[tool.poetry.group.style]
optional = true
@@ -89,6 +88,7 @@ sphinx_rtd_theme = { version = "^1.0.0" }
graphviz = { version = "^0.20.1" }
ipython = { version = "^7.4.0" }
nbsphinx = { version = "^0.8" }
+bnlearn = { git = "https://github.com/erdogant/bnlearn.git", branch = 'master', optional = true }
dowhy = { version = "^0.8" }
typing-extensions = { version = "*" } # needed in dowhy's package
joblib = { version = "^1.1.0" } # needed in dowhy's package
@@ -97,7 +97,6 @@ tqdm = { version = "^4.64.0" } # needed in dowhy's package
[tool.poetry.extras]
graph_func = ['pywhy-graphs']
viz = ['pygraphviz']
-data = ['bnlearn']
[tool.portray]
output_dir = ['site']