diff --git a/MyIA.AI.Notebooks/Probas/DecisionTheory/Causal-Bridges/CausalBridges-00-PearlLadder-Intro-Python.ipynb b/MyIA.AI.Notebooks/Probas/DecisionTheory/Causal-Bridges/CausalBridges-00-PearlLadder-Intro-Python.ipynb new file mode 100644 index 0000000000..73543f3b7d --- /dev/null +++ b/MyIA.AI.Notebooks/Probas/DecisionTheory/Causal-Bridges/CausalBridges-00-PearlLadder-Intro-Python.ipynb @@ -0,0 +1,633 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "db38d5d9", + "metadata": {}, + "source": [ + "# CausalBridges-00 — Du graphe causal a l'echelle de Pearl (socle introductif)\n", + "\n", + "[← Causal-Bridges](../README.md) | [↑ DecisionTheory](../README.md) | [Suivant: CausalBridges-01 →](CausalBridges-01-Do-Calculus.ipynb)\n", + "\n", + "**Pli 1 de l'Origami causal** (EPIC #19309, sous-issue #19310). Ce carnet est le **socle** : il pose le vocabulaire commun (DAG, P(Y|X) vs P(Y|do(X)), contrefactuel individuel) que les 8 plis suivants (CB-01..CB-08) supposent acquis. C'est aussi une **entree en matiere** : on ne demande pas au lecteur de savoir ce qu'est un do-calculus, on lui montre pourquoi l'observation seule est piégeuse, puis on pose la grammaire." + ] + }, + { + "cell_type": "markdown", + "id": "2648a190", + "metadata": {}, + "source": [ + "## Objectifs d'apprentissage\n", + "\n", + "A l'issue de ce carnet, le lecteur est capable de :\n", + "\n", + "1. **Reconnaitre** une situation ou l'observation seule est piegeuse (Snow/cholera 1854, Doll-Hill 1950) et nommer ce que la causalite ajoute au raisonnement.\n", + "2. **Dessiner** un DAG a 4 noeuds (Z confondeur, X traitement, M mediateur, Y outcome) avec `networkx`, et **lire** les chemins ouverts/fermes qui produisent une confusion.\n", + "3. **Distinguer** P(Y|X) de P(Y|do(X)) sur un DGP semi-synthetique, et **verifier** le theoremes d'ajustement backdoor : P(Y|do(X=1)) = E_Z[P(Y|X=1, Z)] quand Z est un ensemble de backdoor valide.\n", + "4. **Calculer** un contrefactuel individuel P(Y_x | x_0, y_0) par abduction-action-prediction (Pearl §9) sur un exemple a 1 individu.\n", + "5. **Reconnaitre** le paradoxe de Simpson (effet qui s'inverse en agregeant) et replacer les trois niveaux (observation, intervention, contrefactuel) sur l'echelle de Pearl." + ] + }, + { + "cell_type": "markdown", + "id": "9f83bad1", + "metadata": {}, + "source": [ + "## Navigation\n", + "\n", + "1. **Motivation historique** : Snow, Doll-Hill — pourquoi l'observation seule echoue.\n", + "2. **DAG 4 noeuds** : Z -> X, X -> M -> Y, Z -> Y (mentorat / effort / score).\n", + "3. **P(Y|X)** : association naive, mesuree sur un DGP semi-synthetique.\n", + "4. **P(Y|do(X))** : intervention par mutilation du graphe (Pearl §3.3) ; ajustement backdoor.\n", + "5. **Contrefactuel individuel** : P(Y_x | x_0, y_0) par abduction-action-prediction (Pearl §9).\n", + "6. **Cinq exercices** : etendre le DAG, paradoxe de Simpson, contrefactuel inverse, etc.\n", + "\n", + "**Duree estimee** : 45 min (lecture + execution) ; 1h30 avec les exercices.\n", + "\n", + "**Prerequis** : probabilites conditionnelles, notions de graphe (noeuds, aretes orientes).\n", + "\n", + "**Stack** : Python 3, `numpy`, `scipy`, `networkx`, `matplotlib`. Aucun GPU. Aucun appel reseau." + ] + }, + { + "cell_type": "markdown", + "id": "3762cb73", + "metadata": {}, + "source": [ + "## 1. Motivation historique — quand l'observation seule est piegeuse\n", + "\n", + "### 1.1. John Snow et le cholera de Broad Street (1854)\n", + "\n", + "A Londres, en septembre 1854, une epidemie de cholera eclate dans le quartier de Soho. Snow cartographie les deces et constate qu'ils se concentrent autour de la **pompe de Broad Street**. La **theorie des miasmes** (l'air empeste) etait alors dominante, mais Snow defend l'hypothese que le cholera se transmet par l'eau contaminee. La cle n'est pas la correlation (beaucoup de personnes buvant l'eau de la pompe tomberaient malades, par hasard) mais le **mecanisme causal** : la pompe approvisionne un puits contamine par une fosse septique voisine ; couper la poignee de la pompe stoppe l'epidemie.\n", + "\n", + "L'**observation pure** (qui tombe malade ?) n'aurait pas suffi : les gens qui boivent l'eau de la pompe sont aussi ceux qui habitent le quartier, et la pauvrete, la nourriture, l'age, le sexe, le metier sont des **confondeurs** (variables qui causent a la fois l'exposition et la maladie). Ce n'est qu'en **intervenant** (couper la pompe, puis observer) que Snow tranche entre les theories.\n", + "\n", + "### 1.2. Doll et Hill, tabac et cancer du poumon (1950)\n", + "\n", + "Dans les annees 1940-50, la correlation entre tabac et cancer du poumon etait massive (les fumeurs ont 10 a 30 fois plus de cancers). La reponse initiale des cigarettiers : il existe un **genotype commun** qui predispose a la fois au desir de fumer et au cancer — la correlation est donc triviale, pas causale. Doll et Hill repondent par une **etude de cohorte** sur 40 000 medecins britanniques : ceux qui ont **arrete** de fumer voient leur risque chuter en quelques annees, et ce d'autant plus que l'arret est ancien. C'est une **intervention naturelle** (l'arret) qui decoince le confounding hypothesique.\n", + "\n", + "Le cas est reste celebre parce que l'**experience randomisee** n'etait pas ethiquement possible (on ne peut pas forcer des gens a fumer pendant 20 ans) — il a fallu raisonner causalement sur une observation, en controlant les confondeurs mesurables et en invoquant la **plausibilite biologique** (le mecanisme cancérogène est aujourd'hui etabli au niveau cellulaire).\n", + "\n", + "### 1.3. La lecon commune : voir ne suffit pas\n", + "\n", + "Dans les deux cas, l'**observation conditionnelle** (P(Y|X)) etait trompeuse : Snow voyait les deces pres de la pompe, Doll voyait les cancers chez les fumeurs. C'est l'**intervention** (couper la pompe, arreter de fumer) qui transforme une correlation en un fait causal. La causalite n'est pas un raffinement statistique — c'est un changement de **mode d'interrogation** : au lieu de demander *que vois-je quand X se produit ?*, on demande *que se passerait-il si je fixais X a une valeur ?*." + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "id": "ce2216d3", + "metadata": { + "execution": { + "iopub.execute_input": "2026-10-05T19:53:11.998098Z", + "iopub.status.busy": "2026-10-05T19:53:11.997579Z", + "iopub.status.idle": "2026-10-05T19:53:14.223267Z", + "shell.execute_reply": "2026-10-05T19:53:14.222261Z" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Taux de cholera par quartier (DGP illustre) :\n", + " Soho (pompe Broad) : 9.10 %\n", + " Bloomsbury (autre pompe) : 0.40 %\n", + " Westminster (riviere) : 0.90 %\n", + " Camden (puits profond) : 0.10 %\n" + ] + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAABKUAAAGGCAYAAACqvTJ0AAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjgsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvwVt1zgAAAAlwSFlzAAAPYQAAD2EBqD+naQAAWWRJREFUeJzt3QncjPX+//EP2WWP206KJEok60nJSdGi7ZQWhSN1pKgoFU5JVCoqqSjlRNHmtGlTKdGCEke0KUq27GSf/+P9/T2u+c+Mue977jH3NXPPvJ6PxzD3NTPXfOdaP9fn+i6FAoFAwAAAAAAAAAAfFfbzywAAAAAAAAAhKQUAAAAAAADfkZQCAAAAAACA70hKAQAAAAAAwHckpQAAAAAAAOA7klIAAAAAAADwHUkpAAAAAAAA+I6kFAAAAAAAAHxHUgoAAAAAAAC+Iynlk19++cUKFSpkzz77rF9fiTRw6qmnWuPGjZNdDKSJunXr2tVXX+3b923fvt2qVKliU6ZM8e07EZ3Wu9Z/KJ2T/v3vfwf/1vlJ03S+KmjeeecdO/zww239+vXJLgqQVohfkVdffvmlFStWzH799VcWXhJ45/L58+enxPXI3r17rVatWvb4448nrDxIPySlEnwAiPa47bbbrKD7448/3O847bTTrEyZMu53ffzxxzF/fvny5TZgwABr06aNlShRImkXPj/99JP16dPH6tWr58pRtmxZa9u2rY0dO9b++usv38sDJNrcuXNdomHz5s1JX7jar3S8uPTSS+M6lhw4cMCeeOIJa9q0qUs4ZGVl2VlnneV+Yyh9Prvj7+eff55t+bSMlDTT+15++eVDCta87ytcuLA7rhxzzDF25ZVX2vvvv5/t5/T7Jk+ebH//+9/tiCOOsKJFi7rynHHGGfbUU0/Z7t27w94f+rv0PdWrV3fvzcuxOJnefvvtsCRYIp155pl29NFH28iRI/Nl/kC6In7N2auvvmqXXHKJixtLlSrlju0333yzL+fY7M5rkY9UOwfccccd1q1bN6tTp46lok2bNlmRIkVs+vTpCbtuef31161Zs2buvbVr17Zhw4bZvn37DnqftptrrrnGKleubKVLl3ax0MKFCw9pnqlOsc1NN91kI0aMsF27diW7OEhRRZJdgHRz991325FHHhk2TZllHZiV9NCOWRDp4HzfffdZ/fr1rUmTJjZv3rw8fV7vf+SRR6xRo0Z27LHH2jfffGN+e+utt+ziiy+24sWLW/fu3d162bNnj82ZM8cGDhxo//vf/9yFIFCQKWFz1113uZox5cuXP2g/VjLDD7ozpqSUgrrDDjssrmOJ9suHHnrIrrjiCvvXv/7lgrknn3zS2rdvb5999pmdfPLJYe+/4YYbrEWLFmHTlKjIztChQ23nzp2WCDVr1gwmRHbs2GE//viju5h5/vnn7R//+If7P/T4r/PB+eefb++++64Lem+55RaXdNu4caPNnj3b/d4vvvjCnn766bDvUQJLx69AIGArVqxwdx47dOjgjm9K2GVnwoQJLgmW7KTUuHHj8i0xpRsOWo7a/pXwBBA74tfolEDQDQCdh5QYWLx4sT322GPueKZkQsmSJfNtM/vPf/4T9rduYuhGR+R0xdWpQvH9Bx98cNDNo1Si866STLqpk4jrlpkzZ1rXrl3dDapHH33UbSP33HOPrVu3zsaPHx98n87BXbp0sUWLFrn4RjejdA7X5xYsWODiorzOs6Do0aOHuyE5depU69mzZ7KLg1QUQEJMmjQpoMX51VdfpdwS3b59+yHPY+vWrYE///zTPX/ppZfcb/3oo49i/rw+q3nIAw884D6/YsWKgF9+/vnnwOGHHx5o2LBhYPXq1Qe9/sMPPwTGjBkTSDXt27cPHHfccckuRuCvv/4K7N+/P5BpErHv+F1Wv/avvXv3Bnbv3p3t66+++qorx48//hjXsUTzL1myZOCiiy46aF/WZ2644YbgNH1e0zS/WC1evDhQpEiRwN13353nz8a6n+7bty/wr3/9y81/0KBBYa/16dPHTc/uuPP9998Hxo0bFzZN7+/bt2/YtG+//dZNP+OMM/Jcbn1u2LBhB53H8mvbUdljDTty276iWbt2beCwww4LPP3003GWEMg8xK85i3Z+eu6559yxbMKECfm2Xg71GJosOjfXrl07cODAgUCquvLKK915O1HXLY0aNQqccMIJ7rzlueOOOwKFChUKfPfdd8Fp06ZNOyjeWLduXaB8+fKBbt26xTVPv/bpRFyPnH322YG//e1vCSsT0gvN95LcJv+ll15yWXhVzVTNnddee+2gvj+8pimR1XOjzVOfVTMXNVPr3Lmzu1t8+eWXBzP0Y8aMseOOO859n+7K686yqrHmRvOpWLFi3L9fn03mnev777/f9W+jWgfVqlU76HXVprjxxhuDf0+aNMnVPlBTGtWs0jqKdmdC7bU7derk7nbobplqycV6B0B3R7QuNH/dhevbt2+21cF1B0W1KbzvUJOmSLqTovmpenmFChXspJNOcnckQv3++++ufFr3+l69/5lnngl7j7e9vfjii3bnnXdajRo13Dx1R1DTn3vuuWzvOr355pt5+q5Yyx3JK+O0adPs9ttvt6pVq7pq0Oeee66tWrUq7L2ffvqpqyGnO5wqh9q1q/ZOZHPNnPad7KiWnWrmaH866qijXC0e1QJR2WLpjyOyTx/1v6AaMmoeoHVdqVIlV/bIKuNecwuvRo22U9XU0bx09020nXhV+73PR+tTSttc//793XLR8tG+oJpMobVqvN8wevRodwzRb9V7ly5dmu2ymTFjhvs+vTeeY4lqWmkdafsJpd+q2l7Z3Z3etm1bTNXbtb+rptLf/vY3yy+qIebdadWd9S1btrjp2kYnTpzompyFHndC6Y6p1m1uVNtMxx/Vmsprn1KxiNxGPZHbktaXaiip3NoftO22a9cu2HxR71UtKW+e3iOW7WvZsmV20UUXue1G89YxQs0aImnbOP744+2///1vnn8ngOgyPX5VTZVIOnfId999Z8kWa7way7Fc9yrUnExNy1Qjx6NWBTrX6NismsA50blf5QmNg0Jr/6ims9aJmrkrfgqN92KN19asWeNq3iju0fsU15933nkxdQuibUl9EKrGUiKuW3Se0kM16tQk0KPzt5ZnaNcAeq5t94ILLghO07JWbWqdt7wm+3mZZ040PzWb85oKaruN7HdR36tloesQLUut4+HDh9v+/fvjvh7Jjmp6K25WjXAgEs33EkwXHRs2bAibpguGaNTcQu3UdaBXsw+dXHv16uWSAIdCF2RKlOiCQEG+LvZFJ3AFADqQq5mLLmJ0ofT111+7pjAFtWlhLN544w3XH4AOpLHQCV3Bj5IcOiHo8zoZ6GSm5JHohK2qvzrYq0qqmkrphKgmO7lRYKALuI4dO9p1113nmjTpO7/66quD1oW2CwVoOmmpjb7awOsz6kTSS4CpaY7WqS7cdJGrNtvffvuta/5z2WWXufesXbvWWrVq5QKF66+/3pVbAYK2ua1bt7rERCidlPQdag6jE5sCHS1Dff9VV10V9l4lh5RQ0naXl++Kpdw5Uft0fcett97q1oeCVi1TVbP2khYKnNVES8tMF8rqgFOJsN9++829Fsu+E42qUnvrX+tTn1V7/8gkSl5o/avKu/pgUrCl7UnbhYJiBSiR5dE2qe9XMzQFimq+9f3339sLL7xgDz/8cPDYo/dEo+WiAFEJRB0fFAjq+wcPHuz6ftLyjAx+tY4UKCl4ySnQ13zUF0K8tP5atmzpjlmtW7d2ySMl0LRdaltTGSLp2Kbks5JBev8DDzzgEhiRtN5VPl1Q5HffdiqL9tshQ4a4YEzBn/YFBXxqDnKodHzQI6dmin7QPqDz2D//+U/XrFL7uZL2SmYrENX2tXr16qhNT3LavtSsWv3+6byo46wCax2D1KzhlVdeCV4cepo3b+4uigDkDfFr7JQUySm+91Ms8WqsFE/p5qGS+9dee20wnlVso2Oxkos6BmdHscTKlSujnvt1LlfMqrIqxlDMrOsPJYi8eC/WeO3CCy905enXr59Lqin+07lF353bzRfFWUrMKK5OBP0GiYw1lORRHOe97r1XyyayGwWdM9V9iOI3XRPmZZ450fJRvKT1p1hHMZ1icsXsoetFyWAlr/T/hx9+6GJKncMVQ4WK5XokJzo/K6mm+Ovss8+O6TcggyS7qla68KpKRnuIqnzqud7nadKkSaBmzZqBbdu2Bad9/PHH7n116tQ5qGlKZBXiaPO86qqr3LTbbrst7L2ffvqpmz5lypSw6e+8807U6TmJp/leKL+b723ZssV933nnnRfzZ3bu3HnQtE6dOgXq1asX/Pu1116Lq3qsquoWK1bMNbcJbRL32GOPufk988wzYdVlNe3BBx8MTlOTlqZNmwaqVKkS2LNnj5um35ZbtdpevXoFqlWrFtiwYUPY9EsvvTRQrly54G/2tjf91sjlMHjw4EDRokUDGzduDCuPqh737Nkzz98VS7mj8cpYo0aNYPVqmT59ups+duzYHNflyJEjXRXoX3/9Ndd9Jztdu3YNlChRImweS5cudc2HQg+t0fbT7JpPRSvrvHnz3PsmT5580PGmXbt2rolYrPuXjiv6nZ7hw4cHSpcu7ZqKhdIy0O9YuXJl2G8oW7as235zo+rmWr4333zzIR1L1Ky2WbNmYcdTbZfLli0Le99nn30WuPDCC12zrf/+979u/VaqVMmtn4ULF4a9V8tYTQu0Lcfb9C+v1dq9Y4W3XQ4YMMD9/c0334S9T/vS+vXrg4/I/Uef0b6l17Qevvjii8Dpp59+0DEiGq330PNKrM33It+T3bakZgZdunSJq+lJTtuXfp/Olbt27QpOU7OQNm3aBOrXr3/QvO699143LzXlA5A74te803FY58jIc2d+i3YMjSVezcuxXJ588kn3/ueffz7w+eefu9/av3//XMv3wQcfuM+98cYbYdM3b94cKFOmTKBly5auO4hQoc38YonXNm3a5L5DsU48hgwZctC5MDc5xVXea168FKpFixaBVq1aBf9WvBUaK3veeustNw9dk+V1njnt0x07dgxbvoo9tC61PnJa5upeoFSpUmHn3VivR3Ki7lM0j/vuuy/X9yLz0HwvwdQ8Qdn60Ec0umOsmhbqsFaZaY9qLShLfqiUuQ6lOwzlypVzd6xVk8t7KGut7//oo48sXSnbL3lpPhjaNMi7e6h18/PPPweb4HidSKvJmpquxEodQKoqtGoLhd4t6d27t6vOrBp0oXTnS7UMPLojob91Z0jVaL2y6E6S7gBFo3hEtQrOOecc9zx0G1DNIP2myNE/VBsqsomUavbpt4bWBnvvvfdcDRa9ltfvyq3cudH+E7peVeNK1bjVAakn9DeoNpHKoRpzKlu0u02R+040quWiJouqraHaRR51hunVFotHaFm1nP/8809XA0bLKdroLNpmQjsRzysdF1SjSHfSQteTapvpN37yySdh79fdyexqXYVS1WwtX833UGjd6q6q7vZqm1OTV9VI03IPrZGq9anq7LpTp7vFqlGjUfd011d3ZEONGjXKLVs1+/SLd4xX08LQY1LosV+03Wr5eo9oIxepCbJeU1MN1SRTzUrd4Yys6eg3baO6c/3DDz/EPY/I7Uvbke7a6q6slp23fWq/0H6m79Kd+VDeNhdZYxlAzohfY6PmZjoOawS+0I6pkyWWeDWvVFtVx1jVtNEosmrSde+99+b6OR2bJfLcr2shHcN1blbzy1Chzfxiidf0HsXBqrUVS/PNSDrP5tZ0Ly+8poWq3RtJvzW06aGeZ/e+0HnlZZ65rcfQ5at4T7GduoqItsy986zepxprajqf1+uRnHB+Rk5ovpdgqoIZrblIJO+AEK3JhaZlNzxoLHTQUPXOUAredXLShUw0oW3HCwIdkCNPtupXKBolekIvCGOhCz1Vd9XoG5Gjc+l7leDTSV8XUWqGp6ZSamKli2VVQ452Iolc9+o3KJQO7moeF3qy8KrrRlaXbtCggftf1XHVTE7N15Ts0van7UfNylQONXsRVVVW4kjVg7MbYTByG4gcRVJOOOEEa9iwoav6q6Z4oueqwq4+BPL6XbmVOzeRAaFOvppPaJMsVedWVWT1QRMZwERuQ9H2nWj0G7UNRgtItV5Dk2J5oXmqCZSaMeli+/9ubkYva3brKC90XFBzyewSTbFsEzkJLX9eKfmk5Jg38oxH05SoUrVy9X2VHW0H6mNCySwFYUreabvQ53TxFZkQyk9qUiheAtX735vu0Xbv3chQOXUciqTfpOr32ta9pF1OzSn8HLlLZdOxSf3LqL8sXcyoGUisIrcvjWCobUhNH/XIbhsNbfLubXPR+jMBkD3i19ypzyPFPkrYqPuA/I5fExWvxkOJNyWjFCeouVVeRhmMPPernzDRuSEnscRriq917ldSUN0lKAZWUzDdpMxtOarZpa6vdL5KFG+5eP1BhVJz9NDlpufZvS90XnmZZ05Cb5qGJoVCl61uJqn/WN0A8m6YeSK301iuR3LC+Rk5ISlVAGQXXGfXCZ0O2JHtldW2XAmpKVOmRP1MLLUfUokSIeo/JpYLYCWldCBdsmRJTPPWyfP00093yRcNR6+OFpUwUqJBySevA2itF9XOUI0MteFXzRnV1HjwwQfdND8velVDR/1SqdaW2uerppJqlejkrqSZV2b1YRPZH5Qn8uIxu5OeakQpGNPdFF0UK3hQ23KvM8a8fFdu5T5U2kdUO1A1LpQA0zrVCVUJH3XsGdqZd3b7jp/7r+5KKiGlWi/qR0nBpD6vPqYiyyqHOhS15qnlM2jQoKive8FGXr9PfQGp3PHcxfSolpb2We2DoZQI1HYTLWETSfuuaiXqjquOA9qulMBQostLXHp9gyjRqGkK4hK9DXjHHu8mhLZDb7oSvaHHYSXd5Pnnn486LyVNvfckU+T2e8opp7hjpzpNVe1JdeSu46U6QVU/U7GI3L68bV792mVXAzHyxo63zaVCXy9Apkun+HXRokWuJq4SK4r9Qjugzq/4NVHxak6yWxeqieQlRdSyQzFJbtQPlMRz7s9LvKYYSbXx1X+gYm/dtNANPSVWTjzxxGy/Q/05qqaROnNPFG/wJPXDqeUfStOU7A19r6ZF8qbpWiWv88xJdjXpve1NN5B1g13xkRJ1SkJq+Shxp3UQy/aTF5yfkROSUkniNcvQneBIkdO8zHbkyGyRNWpyogONaqToTvyhXsimAl2gZNc0MhrdRVGtHd1Jyu3EqgSTTsRKtoTeZciuiaPuDOihRI2qdWu0GI1cl92FmLfulYxRzSiPLp7V+XzkBaeaeuqiOvTuhDpDlNAOHfW6EkZ6aF4a3UNlUvMlBW1KIOmkf6gXtJq/EkZKIOkule6sKGniyet35VTuyGrekSKbCulEq/3HS3opkNKy0oiBuovmycu2E41+o/ajaE2VtF7j3X8V6CqRp8Rm6F2x7EZljCYvNUR0XFBtnUQnORSsa965jQiXE3WWn13ArOZ3sYywp+YL2oa8BLHuwmr7CN3vPN5IdwqavKa5iaDy67igTurVgb6oQ3oFi7rIym2Ex2TT9hu5/Wk/jRZYKxmpiy2vw3klqtQBuncszGvtJW89aeCHWLdRbXNKSBW0Gy1AQZGJ8auSP6r9qeSYEj6HctMxr/FrouLVvBzLNU03yVR73RvwRuWO1qQ8lHfDJfLc743Cqxsx2Q3Kkdd4TfNUbSk9FIs1bdrUxU7Z3dARdY+hhFQityN9r2hgj9BkkWJ3dU8ROiiL3qvadkr2hCZfNbiPYgTvRmBe5nkolHhUk0vVKNf52pNd7Bbr9Uh2vPnqxiIQiT6lkkTZcN1tmTx5clgTDg3xrgNzKJ0EdAET2b+LapTESn1y6OJII1dF0sVdXi56c6MLv8h2yHk9+XtVfbOjuwi6SAl95EQ1QXQQ1cWRd7Eb+Z1jx44Nu7MQ2XRKNVhC6eI18u6WdyKJVuXWo7LqJK+h4kM/r6rS+p7Itu5aP08++WRYEKG/ddGlPsFC2/F7NH+Nlqf56wJev0lNDZVIilZjLHKI2JzoZKJ+z3S3Tw+ti9CTWV6+K7dy50b7T2izTCV1FEzpot8ri4QuZz331nW8NF8FaLpLp+3do9HcdNculO5A6SI5lv1X843cptR0Lbs7mdF4wUIs+7SOC0rURpbZ+3wsiZ/sKPmrgCpeXnCmBG8o3cFT4i/0bmi07Vd3tRWoK6j2gr977rnHDVse+vCOiTpG6O9ENoXTetPoktou9L/XlFgXD6pVqbu2GgE10U0fE0mBf+S2qwR/5DYZuS/rok0XH6HHwrxsm6ILQNVq0/Eu2oVTtPWufi1iuaMPID6ZFr+qNq13HtG5MqeEd37ErzmJNV7Ny7Hc669SiRPFpXqPbjSp2WJu5yXVRFbNnshzv5afblaqNpPXVM3jzTPWeE1NFCPnod+m+ecUeyumVIIrkf1JiZrQKxkXuSw1KqJuxKivU4+e6xoktF9WtTpQ/56q+eV1/ZGXeR6KaMtc1xjZ7Z+xXI+I9qHQ2Dj0/Kzyc45GNNSUSiJ1Gqg+OHT3R3eWleTQBYpO9qEnejXhufjii93FqXZmHXzV3Ckv/UCpeqY6o9MJ4ZtvvnEnCN191t0FHQx10M/tIKcLOq/9sWhYbw1xLmqP7NEdDgUnkSdJr18Yr9mNfqtqJOihPlI8qoosiRyqXctMtRVUG0dJFZVRy1kHVLWV1zJQ9WDx7gzpBKFlpnUxYcIEd4EUemGkuzk6cGtIcs1fyRG9TxeeOQ01q4O3agGptpHuvKk6uC6yNa8WLVocNEy8AkC1n9fy0IW6EkFahzpZaR16ZVZbem1Lqr2ki2AtX518vf5r1MGz7p6pc2QFHEr+qJq0LvJ1F1LPY6XlqKZQqoWiQCWyun2s3xVLuXOimhmqfaL9Ryd6DXerC2F9p+ikrnWju3yqAq51o2TZoTQr82j9qcmhOoRULRudrLWNK5hQP02hlAzVMtH/6nNOgaF3dymyRp/2K+3zWmZKGGl5eVXiY+EFBnfccYerwaZtRNtytGTLwIEDXeJG36vtX5/VXTBdWCjBp20u3mZQOrbpt+h3RjYDjOVYorKoKr/2M9XG07ai/U/LWHc5Qzv21vaoaeoQVfvp0qVL3f6hO49a7h6vplIor1aU9j31CRdKx1sdO3U3MTc6xnl3aBU0q8aAAk9doGg9RF5QaVvVXUPdjVbiTetIZVeAqmOk7oBH9juXDNpmNTS4Es1aH0r26cIscrvQ9qoEktab9ktdlGgbCj22e9umEnRK6iogDq1lGY36/9J6UyJc+7VqT2lf176hO8Yqj0fnRO17eR0GHUDeZFL8qjhNtW5140Kf8T4nilt0XMzP+DUnscareTmWK6GlGkXPPvtssI9NrT/FpkqKeLWKs6PtQjd4tAy92rGKvdScUGXQuVZ9h6rmlsqg86XO87HGa4optJyVrNR5RwkzfZ/OCzmdT7TeFEvEmpTKy3WL+oBULK/1oTLohqzeq98bWitI26laVmifUZyiZa/YX4mnyC4rYp3noVDMpPWgGvo6L2t9ad/ILvkYy/WIqHzRYiclBXXMyEtMiwyS7OH/0oU3/OZXX30V9fXshoV/8cUXAw0bNgwUL1480Lhx48Drr7/uhjbXtFAaAlzTNURnhQoV3HCdS5YsOWieGtZVQ45m56mnngo0b948ULJkSTc8q4baHjRokBumMzehw7JHPkJ5w4ZG+/3RHpFDs+rvvA7XGisN39u7d+9A3bp1A8WKFXPLoG3btoFHH300bOhTrYfjjz/eDSmv92r40meeeSZsSFgNNd+tWzc3vLzWn4ZEPfvsswPz58+PqSyPPfaYW89FixYNZGVlBa677jo31G20oeY1z9atW7vyaNnos5HD955yyimBSpUqubIcddRRgYEDBwa2bNkS9j4Nk64hhWvVquW+t2rVqm7IdW0Xno8++sj9zpdeeinbsv/www/B9Tdnzpyo74nlu2ItdySvjC+88EJg8ODBbtlrm9aQ9N6wwZ6lS5e6YXEPP/zwwBFHHOHW/6JFi/K870Qze/Zstz9pW9Lwy0888YQbcjnakM0aQrpcuXJum/vHP/4RWLdu3UFDNGv99+jRw5VT5dWwzsuWLTtoyObcjjfDhw8P1KhRI1C4cOGwbTba0M/btm1zy/Doo492v0Pf3aZNm8Do0aODQ/x6+29ehmDWUMGal8oS77FEy+3uu+8ONGrUyK1fLT/tY19//XXY+8aOHRs4+eSTAxUrVgwUKVIkUK1atcAVV1zhttPcZLe9a7lo+qWXXprrPLxjnvfQuqtfv74rw3vvvZft5/bt2+fWZYcOHYJl1zLTfqJtKXLobM1b+1Q8tN4jj6uR25+3XYUOe71///7Arbfe6sql84+2yR9//PGgbemee+5x66B8+fJuXenYNmLEiLBhovV7+/XrF6hcubIb4ttb37ltXz/99FOge/fu7hiiY4m2bW0HL7/8ctj7xo8f78q4devWuJYRkImIX3OOX3M6X+n9fsWvouN/ZPliiVdjPZavWrXKnWfPOeecg777/PPPdzHSzz//nGMZFRvruz/99NODXlNZFV/oHFG2bFl3zlAcl5d4bcOGDW456Byj8qi8LVu2DEyfPj3Hct1yyy0ulohVXq5b5LXXXgs0bdrUxbI1a9YM3HnnnWHnP8/GjRtdPKi4V+tB21B2sVys84x1n/biHf3v+eyzzwKtWrVy66R69erumvDdd9896H2xXo9ItH1j8+bNLsacOHFiruVHZiqkf5KdGIMd1ARMtWkS1eYcSDe6+6J+AXSXNFHVmBNFfejojheHVnO1g3TXVXe0s+twM1Wp3xDVINOdXNXSQepTk07V1tIdeQD+I36FqCaTatWo1k2qUK0qndPvv//+ZBclI6l2uJa9ao+nQ9/GSDz6lEqiaJ316mJbF0EKrAGgIBswYIBrThDZL1RBoOanqjJPQqpgUFNaJT/VNBpA/iJ+RW7NO9W0Ky8d2ucnddWhZv6Rox7Cv+OFRodUU1kSUsgOfUolkdpMq4NDtdPWHQV1DKfhs9XHjtp9A0BBps6u89J3SCpRfw4oONTvS2hfNgDyD/ErcqL+RJUIShXqd2vYsGHJLkbGUn9T0To+B0KRlEoidS6njl8nTpzoRhFSR8TqgE8d89IJHAAAAFIN8SsAIJHoUwoAAAAAAAC+o08pAAAAAAAA+I6kFAAAAAAAAHyX9n1KHThwwFavXm1lypSxQoUKJbs4AACggAsEArZt2zY3SEnhwulxf494CQAAJCNeSvuklBJStWrVSnYxAABAmlm1apXVrFnT0gHxEgAASEa8lPZJKdWQ8hZE2bJlk10cAABQwG3dutXd8PJijHRAvAQAAJIRL6V9UsprsqeEFEkpAACQ6BgjHRAvAQCAZMRL6dERAgAAAAAAAAoUklIAAAAAAADwHUkpAAAAAAAA+I6kFAAAAAAAAHxHUgoAAAAAAAC+IykFAAAAAAAA35GUAgAAAAAAgO9ISgEAAAAAAMB3JKUAAAAAAADgO5JSAAAAAAAA8B1JKQAAAAAAAPiOpBQAAAAAAAB8V8T/r0w/s9s0TnYRkI32c5ewbAAASAHES6mLeAkAkCzUlAIAAAAAAIDvSEoBAAAAAADAdySlAAAAAAAA4DuSUgAAAAAAAPAdSSkAAAAAAAD4jqQUAAAAAAAAfEdSCgAAAAAAAL4jKQUAAAAAAADfkZQCAAAAAACA70hKAQAAAAAAwHckpQAAAAAAAOA7klIAAAAAAADwHUkpAAAAAAAA+I6kFAAAAAAAAHxHUgoAAAAAAAC+IykFAAAAAAAA35GUAgAAAAAAgO9ISgEAAAAAAMB3JKUAAAAAAADgO5JSAAAAAAAA8B1JKQAAAAAAAPiOpBQAAAAAAAB8R1IKAAAAAAAAviMpBQAAAAAAAN+RlAIAAAAAAIDvSEoBAAAAAADAdySlAAAAAAAA4DuSUgAAAAAAAPAdSSkAAAAAAAD4jqQUAAAAAAAAMisptX//fhsyZIgdeeSRVrJkSTvqqKNs+PDhFggEgu/R86FDh1q1atXcezp27Gg//PBDMosNAACQMmKJpwAAAFJRkWR++X333Wfjx4+35557zo477jibP3++9ejRw8qVK2c33HCDe8/9999vjzzyiHuPgi0FXZ06dbKlS5daiRIlkll8AACApIslngIAAEhFSU1KzZ0718477zzr0qWL+7tu3br2wgsv2Jdffun+1h2+MWPG2J133uneJ5MnT7asrCybMWOGXXrppcksPgAAQNLlFk8BAACkqqQ232vTpo3NmjXLvv/+e/f3okWLbM6cOXbWWWe5v1esWGFr1qxxTfY8uuvXsmVLmzdvXtLKDQAAkCpyi6cAAABSVVJrSt122222detWa9iwoR122GGuT4QRI0bY5Zdf7l5XQkpUMyqU/vZei7R792738Gj+AAAA6Sq3eCoa4iUAAGCZXlNq+vTpNmXKFJs6daotXLjQ9YUwevRo93+8Ro4c6WpTeY9atWoltMwAAACpJJ54ingJAABYpielBg4c6O7uqW+oJk2a2JVXXmkDBgxwgZJUrVrV/b927dqwz+lv77VIgwcPti1btgQfq1at8uGXAAAApGY8FQ3xEgAAsExvvrdz504rXDg8L6Zq5wcOHHDPNdqekk/qJ6Fp06Zumqqnf/HFF3bddddFnWfx4sXdAwAAIBPkFk9FQ7wEAAAs05NS55xzjuvzoHbt2m4I46+//toeeugh69mzp3u9UKFC1r9/f7vnnnusfv36Lkk1ZMgQq169unXt2jWZRQcAAEgJucVTAAAAqSqpSalHH33UJZn+9a9/2bp161yyqU+fPjZ06NDgewYNGmQ7duywa665xjZv3mzt2rWzd955x0qUKJHMogMAAKSEWOIpAACAVFQoEAgELI2puZ86PFf/UmXLls2X75jdpnG+zBeHrv3cJSxGAECBiy38RryU2YiXAADJii2S2tE5AAAAAAAAMhNJKQAAAAAAAPiOpBQAAAAAAAB8R1IKAAAAAAAAviMpBQAAAAAAAN+RlAIAAAAAAIDvSEoBAAAAAADAdySlAAAAAAAA4DuSUgAAAAAAAPAdSSkAAAAAAAD4jqQUAAAAAAAAfEdSCgAAAAAAAL4jKQUAAAAAAADfkZQCAAAAAACA70hKAQAAAAAAwHckpQAAAAAAAOA7klIAAAAAAADwHUkpAAAAAAAA+I6kFAAAAAAAAHxHUgoAAAAAAAC+IykFAAAAAAAA35GUAgAAAAAAgO9ISgEAAAAAAMB3JKUAAAAAAADgO5JSAAAAAAAA8B1JKQAAAAAAAPiOpBQAAAAAAAB8R1IKAAAAAAAAviMpBQAAAAAAAN+RlAIAAAAAAIDvSEoBAAAAAADAdySlAAAAAAAA4DuSUgAAAAAAAPAdSSkAAAAAAAD4jqQUAAAAAAAAfEdSCgAAAAAAAAUvKbV79+7ElAQAAAAAAAAZI89JqZkzZ9pVV11l9erVs6JFi1qpUqWsbNmy1r59exsxYoStXr06f0oKAAAAAACAzEtKvfbaa9agQQPr2bOnFSlSxG699VZ79dVX7d1337WJEye6pNQHH3zgklXXXnutrV+/Pn9LDgAAAAAAgAKrSKxvvP/+++3hhx+2s846ywoXPjiX9Y9//MP9//vvv9ujjz5qzz//vA0YMCCxpQUAAAAAAEBmJaXmzZsX0/tq1Khho0aNOpQyAQAAAAAAIM0x+h4AAAAAAABSt6ZUqP3799uzzz5rs2bNsnXr1tmBAwfCXv/www8TVT4AAAAAAACkobiSUjfeeKNLSnXp0sUaN25shQoVSnzJAAAAAAAAkLbiSkq9+OKLNn36dOvcuXPiSwQAAAAAAIC0F1efUsWKFbOjjz468aUBAAAAAABARogrKXXzzTfb2LFjLRAIJL5EAAAAAAAASHtxNd+bM2eOffTRRzZz5kw77rjjrGjRomGvv/rqq4kqHwAAAAAAANJQXEmp8uXL2/nnn5/40gAAAAAAACAjxJWUmjRpUuJLAgAAAAAAgIwRV59SAAAAAAAAgC9JqTPPPNM+//zzXN+3bds2u++++2zcuHGHVDAAAAAAAACkr5ib71188cV24YUXWrly5eycc86xk046yapXr24lSpSwTZs22dKlS10H6G+//bZ16dLFHnjggfwtOQAAAAAAANI/KdWrVy+74oor7KWXXrJp06bZU089ZVu2bHGvFSpUyBo1amSdOnWyr776yo499tj8LDMAAAAAAAAyqU+p4sWLu8TUG2+84WpH6bF69WrbtWuXLV682EaPHp3nhNTvv//u5lmpUiUrWbKkNWnSxObPnx98PRAI2NChQ61atWru9Y4dO9oPP/yQp+8AAABIZ7nFUwAAAGnX0bma8lWtWtWKFi0a1+eV1Grbtq37/MyZM10TwAcffNAqVKgQfM/9999vjzzyiD3xxBP2xRdfWOnSpV2NLCXCAAAAMl0s8RQAAECBbr6XH9Qheq1atWzSpEnBaUceeWRYLakxY8bYnXfeaeedd56bNnnyZMvKyrIZM2bYpZdempRyAwAApIrc4ikAAIC0rCl1qF5//XXXYbo6Ua9SpYqdeOKJNmHChODrK1assDVr1rgme6G1s1q2bGnz5s1LUqkBAABSR27xFAAAQKpKalLq559/tvHjx1v9+vXt3Xffteuuu85uuOEGe+6559zrSkiJakaF0t/ea5F2795tW7duDXsAAACkq9ziqWiIlwAAgGV6870DBw64O3v33nuv+1t39pYsWeL6j7rqqqvimufIkSPtrrvuSnBJAQAAUlM88RTxEgAAsEyvKaUR9Ro1ahQ2TaP3rVy50j1XJ+qydu3asPfob++1SIMHD7YtW7YEH6tWrcq38gMAACRbbvFUNMRLAACgwNaU2r9/vz388MM2ffp0F/Ds2bMn7PWNGzfGNB+NFLN8+fKwad9//73VqVMn2Emnkk+zZs2ypk2bumlqjqdR+FQ1PZrixYu7BwAAQCbILZ6KhngJAAAU2JpSah730EMP2SWXXOJqI9100012wQUXWOHChe3f//53zPMZMGCAff755666+Y8//mhTp061p556yvr27eteL1SokPXv39/uuece14nn4sWLrXv37la9enXr2rVrPEUHAABIK7nFUwAAAGmVlJoyZYob1eXmm2+2IkWKWLdu3WzixIk2dOhQFxTFqkWLFvbaa6/ZCy+8YI0bN7bhw4fbmDFj7PLLLw++Z9CgQdavXz+75ppr3Pu3b99u77zzjpUoUSKeogMAAKSVWOIpAACAVFQoEAgE8vqh0qVL23fffWe1a9d2/Ri89dZb1qxZMzf6izrXVO2pVKHmfuXKlXNlKlu2bL58x+w2jfNlvjh07ecuYTECAApcbOE34qXMRrwEAEhWbBFXTamaNWvaH3/84Z4fddRR9t5777nnX331Ff05AQAAAAAAIFdxJaXOP/981/m4qGndkCFDrH79+q6/p549e8YzSwAAAAAAAGSQuEbfGzVqVPC5OjvX6C5z5851ialzzjknkeUDAAAAAABAGspzUmrv3r3Wp08fVzvqyCOPdNNatWrlHgAAAAAAAEC+NN8rWrSovfLKK3n9GAAAALKxe/dulg0AAMg4cfUp1bVrV5sxY0biSwMAAJABZs6caVdddZXVq1fP3fArVaqUG5mmffv2NmLECFu9enWyiwgAAJCafUqp76i7777bPvvsM2vevLmVLl067PUbbrghUeUDAABIG6+99prdeuuttm3bNuvcubN7Xr16dStZsqRt3LjRlixZYh988IENHz7crr76avd/5cqVk11sAACAfFEoEAgE8vohry+pqDMsVMh+/vlnSxVbt261cuXK2ZYtW9wdyPwwu03jfJkvDl37uUtYjACAlIktWrdubXfeeaedddZZVrhw9hXWf//9d3v00UctKyvLBgwYYPmNeCmzES8BAJIVW8RVU2rFihWHUjYAAICMNG/evJjeV6NGjbDRjgEAANJRXH1Kefbs2WPLly+3ffv2Ja5EAAAAGWjHjh3uriIAAECmiCsptXPnTuvVq5frlPO4446zlStXuun9+vXjrh4AAEAeLF261E466SQrU6aMVahQwZo0aWLz589nGQIAgLQXV1Jq8ODBtmjRIvv444+tRIkSwekdO3a0adOmJbJ8AAAAaa1Pnz52/fXX2/bt2+3PP/+0Cy64wI3MBwAAkO7iSkrNmDHDHnvsMWvXrp3r2NyjWlM//fRTIssHAACQVs477zzXkbln/fr1du6557oa6OXLl3ej8q1duzapZQQAAPBDXB2dK3iqUqVK1L4QQpNUAAAACHfFFVdYhw4drG/fvq7rA9WS0o299u3b2969e+3DDz+0m2++mcUGAADSXlw1pdTvwVtvvRX820tETZw40Q11DAAAgOguvvhi+/LLL11fUq1atbK2bdvae++95/7/29/+5p7feeedLD4AAJD24qopde+999pZZ53lgimNvDd27Fj3fO7cuTZ79uzElxIAACCNlCtXzp544gmbM2eO6z/q73//uw0fPtw14QMAAMgUcdWUUl9S33zzjUtIaYQY3dFTc7558+ZZ8+bNE19KAACANLJx40ZbsGCBi6P0f9myZe3EE0+0t99+O9lFAwAASO2aUnLUUUfZhAkTElsaAACANDd16lT75z//6RJRu3btssmTJ9uwYcPskksusWuvvdaeffZZe/TRRy0rKyvZRQUAAEiNmlJbt26N+QEAAIDoBg8ebM8884ytWbPGZs2aZUOGDHHTGzZsaB9//LFrykcfnQAAIBPEXFNKQxTnNrJeIBBw79m/f38iygYAAJB2tm/fbsccc0yw5vnOnTvDXu/du7edd955SSodAABACialPvroo/wtCQAAQAZQx+ZdunSxU0891ebPn29XXnnlQe9RX50AAADpLuakVPv27fO3JAAAABngoYcestNOO82WLVtmV199tZ1xxhnJLhIAAEDB6uh88+bN9vTTT9t3333n/j7uuOOsZ8+ebohjAAAAZO+cc85xDwAAgEwWc0fnoVTVXH0gPPzww25IYz1010/TFi5cmPhSAgAApIEXX3wx5veuWrXKPvvss3wtDwAAQIFLSg0YMMDOPfdc++WXX+zVV191jxUrVtjZZ59t/fv3T3wpAQAA0sD48ePt2GOPtfvvvz9Y2zzUli1b7O2337bLLrvMmjVrZn/++WdSygkAAJCyzfdUU2rChAlWpMj//7ieDxo0yE466aRElg8AACBtzJ49215//XV79NFHbfDgwVa6dGnLysqyEiVK2KZNm2zNmjV2xBFHuL6mlixZ4l4DAABIV3ElpcqWLWsrV660hg0bHlTNvEyZMokqGwAAQNpRbXM9NmzYYHPmzLFff/3V/vrrL5eMOvHEE92jcOG4KrMDAACkf1LqkksusV69etno0aOtTZs2bpr6PBg4cKB169Yt0WUEAABIO0pCde3aNdnFAAAAKFhJKSWjChUqZN27d7d9+/a5aUWLFrXrrrvORo0alegyAgAAAAAAIM3ElZQqVqyYjR071kaOHGk//fSTm6aR90qVKpXo8gEAAAAAACANxZWU8igJ1aRJk8SVBgAAAAAAABkhrqTUjh07XDO9WbNm2bp16+zAgQNhr//888+JKh8AAAAAAADSUFxJqX/+859uSOMrr7zSqlWr5vqXAgAAQOx27dplJUqUiPraH3/84WIsAACAdBZXUmrmzJn21ltvWdu2bRNfIgAAgAzQrFkzmzp1qjVt2jRs+iuvvGLXXnutrV+/PmllAwAA8EPheD5UoUIFq1ixYuJLAwAAkCFOPfVUa9Wqld13333B7hGuvvpqVxP99ttvT3bxAAAAUrOm1PDhw23o0KH23HPPMeIeAABAHB5//HHr0qWL6xbhzTffdE32Dj/8cPvyyy+tcePGLFMAAJD2Yk5KnXjiiWF9R/3444+WlZVldevWtaJFi4a9d+HChYktJQAAQBo666yz7IILLrDx48dbkSJF7I033iAhBQAAMkbMSamuXbvmb0kAAAAyyE8//WSXXXaZrVmzxt599103iMy5555rN954o40YMeKgm34AAAAZm5QaNmxY/pYEAAAgg6iDczXfU0KqfPny9ve//906d+5s3bt3t/fff9++/vrrZBcRAAAg9To6/+qrr+yLL744aLqmzZ8/PxHlAgAASPs+pV588UWXkPK0adPGJaM0Mh8AAEC6iysp1bdvX1u1atVB03///Xf3GgAAAHKmUfaiKVOmjD399NMsPgAAkPbiGn1v6dKlUe/gqTN0vQYAAICcTZ48OdvXNLhMdkkrAACAjE5KFS9e3NauXWv16tULm66hjDVyDAAAAHKmDs1D7d2713bu3GnFihWzUqVKkZQCAABpL67me2eccYYNHjzYtmzZEpy2efNmu/32210nnQAAAMjZpk2bwh7bt2+35cuXW7t27eyFF15g8QEAgLQXV7Wm0aNH2ymnnGJ16tRxTfbkm2++saysLPvPf/6T6DICAABkhPr169uoUaPsiiuusGXLliW7OAAAAKmXlKpRo4Z9++23NmXKFFu0aJGVLFnSevToYd26dbOiRYsmvpQAAAAZQl0hrF69OtnFAAAAyHdxdwBVunRpu+aaaxJbGgAAgAzx+uuvh/0dCARc/5yPPfaYtW3bNmnlAgAA8Au9kgMAACRB165dDxpxr3LlytahQwd78MEHWScAACDtkZQCAABIggMHDrDcAQBARotr9D0AAAAAAADgUFBTCgAAIEl+++0317fUypUrbc+ePWGvPfTQQ6wXAACQ1uJOSm3evNlefvll++mnn2zgwIFWsWJFW7hwoWVlZbnR+QAAAJC9WbNm2bnnnmv16tWzZcuWWePGje2XX35xHZ43a9aMRQcAANJeXM33vv32W2vQoIHdd999Nnr0aJegkldffdUGDx6c6DICAACkHcVMt9xyiy1evNhKlChhr7zyiq1atcrat29vF198cbKLBwAAkJpJqZtuusmuvvpq++GHH1wQ5encubN98skniSwfAABAWvruu++se/fu7nmRIkXsr7/+ssMPP9zuvvtud+MPAAAg3cWVlPrqq6+sT58+B01Xs701a9bEVZBRo0a5oZD79+8fnLZr1y7r27evVapUyQVpF154oa1duzau+QMAAKSS0qVLB/uRqlatmusSwbNhw4a45xstpgIAAEibpFTx4sVt69atB03//vvvrXLlynEluZ588kk7/vjjw6YPGDDA3njjDXvppZds9uzZtnr1arvgggviKTIAAEBKUE2oHTt2WKtWrWzOnDnB2uY333yzjRgxwnr27Olei0d2MRUAAEDaJKXUKacCqr1797q/dTdOo8bceuutrjZTXmzfvt0uv/xymzBhglWoUCE4fcuWLfb000+7kWc6dOhgzZs3t0mTJtncuXPt888/j6fYAAAASXfXXXe5pJRinJYtWwannX766TZt2jSrW7eui4HyKruYCgAAIK2SUg8++KALfKpUqeL6P1CHnEcffbSVKVPG3eHLCzXP69Kli3Xs2DFs+oIFC1zSK3R6w4YNrXbt2jZv3rx4ig0AAJB0Gl1PNOqeV6NJTfmeeOIJN5iMOjyvU6dOnuebXUwFAACQqorE86Fy5crZ+++/76qcK3hSgkpDF+c1CHrxxRdt4cKFrqp5JPVNVaxYMStfvnzY9KysrBz7rdq9e7d7eKI1MwQAAEgm1TJPpJxiqmiIlwAAQIFNSnnatWvnHvHQkMc33nijS26FjuB3qEaOHOmqwAMAAKSqBg0a5JqY2rhxY77FVMRLAACgQCWlHnnkkZhnesMNN+T6HjXPW7dunath5dm/f7998skn9thjj9m7777rRqTZvHlzWG0pjb5XtWrVbOc7ePBgu+mmm8JqStWqVSvmsgMAAOQ33UBTzfNEyC2mUq2oww47LOwzxEsAAKBAJaUefvjhsL/Xr19vO3fuDCaMlDwqVaqU62cqlqSUOvNcvHhx2LQePXq4fqPUYboSSUWLFrVZs2YFO09fvny561C9devWOY4MqAcAAECquvTSS13MlAi5xVSRCSkhXgIAAAUqKbVixYrg86lTp9rjjz/uRoY55phjggmj3r17W58+fWKanzpFb9y4cdg0dfJZqVKl4PRevXq5Wk8VK1a0smXLWr9+/VxCKt5hkgEAANKtP6lYYioAAIC06VNqyJAh9vLLLwcTUqLnqk110UUXueGIE0HzK1y4sKspparnnTp1cskwAACAgj76HgAAQKaLKyn1xx9/2L59+w6arv4L1OdTvD7++OOwv9VZ57hx49wDAAAgHRw4cCDfvyMypgIAAEhFhePtu0DN9DT0cGgnm9ddd5117NgxkeUDAAAAAABAGoorKfXMM8+4EfBOOumkYEeZJ598smVlZdnEiRMTX0oAAAAAAACklbia71WuXNnefvtt++GHH+y7775z0zTCS4MGDRJdPgAAAAAAAKShuJJSnvr167sHAAAAAAAAkO/N9wAAAAAAAIBDQVIKAAAAAAAAviMpBQAAAAAAAN+RlAIAAAAAAEDBSEq98847NmfOnODf48aNs6ZNm9pll11mmzZtSmT5AAAAAAAAkIbiSkoNHDjQtm7d6p4vXrzYbr75ZuvcubOtWLHCbrrppkSXEQAAAAAAAGmmSDwfUvKpUaNG7vkrr7xiZ599tt177722cOFCl5wCAAAAAAAAEl5TqlixYrZz5073/IMPPrAzzjjDPa9YsWKwBhUAAAAAAACQ0JpS7dq1c8302rZta19++aVNmzbNTf/++++tZs2a8cwSAAAAAAAAGSSumlKPPfaYFSlSxF5++WUbP3681ahRw02fOXOmnXnmmYkuIwAAAAAAANJMXDWlateubW+++eZB0x9++OFElAkAAAAAAABpLq6kVKhdu3bZnj17wqaVLVv2UGcLAAAAAACANBZX870dO3bY9ddfb1WqVLHSpUtbhQoVwh4AAAAAAABAwpNSgwYNsg8//ND1J1W8eHGbOHGi3XXXXVa9enWbPHlyPLMEAAAAAABABomr+d4bb7zhkk+nnnqq9ejRw/72t7/Z0UcfbXXq1LEpU6bY5ZdfnviSAgAAAAAAILNrSm3cuNHq1asX7D9Kf0u7du3sk08+SWwJAQAAAAAAkHbiSkopIbVixQr3vGHDhjZ9+vRgDary5csntoQAAAAAAABIO3ElpdRkb9GiRe75bbfdZuPGjbMSJUrYgAEDbODAgYkuIwAAAAAAANJMXH1KKfnk6dixoy1btswWLFjg+pU6/vjjE1k+AAAAAAAApKG4klKR1MG5HgAAAAAAAEDCm+99+OGH1qhRI9u6detBr23ZssWOO+44+/TTT/MySwAAAAAAAGSgPCWlxowZY71793Yj7kUqV66c9enTxx566KFElg8AAAAAAACZnpRS5+Znnnlmtq+fccYZrm8pAAAAAAAAIGFJqbVr11rRokWzfb1IkSK2fv36vMwSAAAAAAAAGShPSakaNWrYkiVLsn3922+/tWrVqiWiXAAAAAAAAEhjeUpKde7c2YYMGWK7du066LW//vrLhg0bZmeffXYiywcAAAAAAIA0VCQvb77zzjvt1VdftQYNGtj1119vxxxzjJu+bNkyGzdunO3fv9/uuOOO/CorAAAAAAAAMjEplZWVZXPnzrXrrrvOBg8ebIFAwE0vVKiQderUySWm9B4AAAAAAAAgYUkpqVOnjr399tu2adMm+/HHH11iqn79+lahQoW8zgoAAAAAAAAZKs9JKY+SUC1atEhsaQAAAAAAAJAR8tTROQAAAAAAAJAIJKUAAAAAAADgO5JSAAAAAAAA8B1JKQAAAAAAAPiOpBQAAAAAAAB8R1IKAAAAAAAAviMpBQAAAAAAAN+RlAIAAAAAAIDvSEoBAAAAAADAdySlAAAAAAAA4DuSUgAAAAAAAPAdSSkAAAAAAAD4jqQUAAAAAAAAfEdSCgAAAAAAAL4jKQUAAAAAAADfkZQCAAAAAACA70hKAQAAAAAAwHckpQAAAAAAAOA7klIAAAAAAADwHUkpAAAAAAAAZFZSauTIkdaiRQsrU6aMValSxbp27WrLly8Pe8+uXbusb9++VqlSJTv88MPtwgsvtLVr1yatzAAAAKkmlpgKAAAg1SQ1KTV79myXcPr888/t/ffft71799oZZ5xhO3bsCL5nwIAB9sYbb9hLL73k3r969Wq74IILkllsAACAlBJLTAUAAJBqiiTzy995552wv5999ll3d2/BggV2yimn2JYtW+zpp5+2qVOnWocOHdx7Jk2aZMcee6wLulq1apWkkgMAAKSO3GIqAACAVJTUpFQkJaGkYsWK7n8FUrrT17Fjx+B7GjZsaLVr17Z58+ZFTUrt3r3bPTxbt271pewAAACpGlNFIl4CAACpIGU6Oj9w4ID179/f2rZta40bN3bT1qxZY8WKFbPy5cuHvTcrK8u9ll2fCuXKlQs+atWq5Uv5AQAAUjWmikS8BAAAUkHKJKXUD8KSJUvsxRdfPKT5DB482N0d9B6rVq1KWBkBAABSXSwxFfESAABIBSnRfO/666+3N9980z755BOrWbNmcHrVqlVtz549tnnz5rDaUhp9T69FU7x4cfcAAADINNnFVJGIlwAAgGV6TalAIOCCp9dee80+/PBDO/LII8Neb968uRUtWtRmzZoVnKbhjVeuXGmtW7dOQokBAABST24xFQAAQCoqkuzq5RpZ77///a+VKVMm2E+U+oIqWbKk+79Xr1520003uY46y5Yta/369XMJKUbeAwAAiC2mAgAASEVJTUqNHz/e/X/qqaeGTZ80aZJdffXV7vnDDz9shQsXtgsvvNCNFNOpUyd7/PHHk1JeAACAVBRLTAUAAJBqiiS7qnluSpQoYePGjXMPAAAAxBdTAQAApJqUGX0PAAAAAAAAmYOkFAAAAAAAAHxHUgoAAAAAAAC+IykFAAAAAAAA35GUAgAAAAAAgO9ISgEAAAAAAMB3JKUAAAAAAADgO5JSAAAAAAAA8B1JKQAAAAAAAPiOpBQAAAAAAAB8R1IKAAAAAAAAviMpBQAAAAAAAN+RlAIAAAAAAIDvSEoBAAAAAADAdySlAAAAAAAA4DuSUgAAAAAAAPAdSSkAAAAAAAD4jqQUAAAAAAAAfEdSCgAAAAAAAL4jKQUAAAAAAADfkZQCAAAAAACA70hKAQAAAAAAwHckpQAAAAAAAOA7klIAAAAAAADwHUkpAAAAAAAA+I6kFAAAAAAAAHxHUgoAAAAAAAC+IykFAAAAAAAA35GUAgAAAAAAgO+K+P+VAJA+tj9dKdlFQDYO7/UnywYAgBRAvJS6iJeQbNSUAgAAAAAAgO9ISgEAAAAAAMB3JKUAAAAAAADgO5JSAAAAAAAA8B1JKQAAAAAAAPiOpBQAAAAAAAB8R1IKAAAAAAAAviMpBQAAAAAAAN+RlAIAAAAAAIDvSEoBAAAAAADAdySlAAAAAAAA4DuSUgAAAAAAAPAdSSkAAAAAAAD4jqQUAAAAAAAAfEdSCgAAAAAAAL4r4v9XAuml4Ws9kl0EZGPZ+ZNYNgAApADipdRFvAQgmagpBQAAAAAAAN+RlAIAAAAAAIDvSEoBAAAAAADAdySlAAAAAAAA4DuSUgAAAAAAAPAdSSkAAAAAAAD4jqQUAAAAAAAAfEdSCgAAAAAAAL4rEEmpcePGWd26da1EiRLWsmVL+/LLL5NdJAAAgJRCvAQAAAqaIpbipk2bZjfddJM98cQTLiE1ZswY69Spky1fvtyqVKmS7OIBADJYo7veS3YRkI2lw87IqGVDvAQASFXES6lraQrESylfU+qhhx6y3r17W48ePaxRo0YuOVWqVCl75plnkl00AACAlEC8BAAACqKUrim1Z88eW7BggQ0ePDg4rXDhwtaxY0ebN29e1M/s3r3bPTxbtmxx/2/dujXfyrlj3/58mzcOTX6ud8/+nXvy/TuQuut/+1+BfP8OxOeAH/v/rh35/h1Ivf3fm3cgkBr7P/ESDhXxUmYjXspsxEuZbWsKxEspnZTasGGD7d+/37KyssKm6+9ly5ZF/czIkSPtrrvuOmh6rVq18q2cSGHlyiW7BEiicjaV5Z/J+rH/Z7Jyo/L/O7Zt22blUuA8Q7yEQ5YC2zGSh3gpwxEvZbRyKRAvpXRSKh6qVaU+qDwHDhywjRs3WqVKlaxQoUJJLVuqUyZTybtVq1ZZ2bJlk10c+Iz1n9lY/5mN9Z83uuOnAKt69epWUBEvHRr2mczG+s9srP/MxbrPn3gppZNSRxxxhB122GG2du3asOn6u2rVqlE/U7x4cfcIVb58+XwtZ7pRQoqkVOZi/Wc21n9mY/3HLhVqSHmIl5KHfSazsf4zG+s/c7HuExsvpXRH58WKFbPmzZvbrFmzwmo+6e/WrVsntWwAAACpgHgJAAAUVCldU0rUFO+qq66yk046yU4++WQbM2aM7dixw43GBwAAAOIlAABQMKV8UuqSSy6x9evX29ChQ23NmjXWtGlTe+eddw7q/ByHTs0ehw0bdlDzR2QG1n9mY/1nNtZ/wUe85C/2mczG+s9srP/MxbrPH4UCqTKeMQAAAAAAADJGSvcpBQAAAAAAgPREUgoAAAAAAAC+IykFAAAAAAAA35GUwkE+/vhjK1++PEsmBaTquqhbt67NmDEj2cWADw4//HBbvHgxyxrWtWtX+/e//82SAFL8HJ2JUnVdEC9lDuIleIiX8o6kVBpbvny5nXPOOXbEEUdY2bJlrWHDhnbfffclu1iIcOqpp7qRHHQyK1OmjB133HH20ksvsZzS2OjRo61Vq1Zh06644gorUaKE7dq1KzjtsccesyZNmsT9PYUKFbJvvvnmkMq6ffv2QyqD/PLLL64smzdvPqT5wGzOnDl21llnWYUKFdwF2AknnGD333+/7dmzh8UDxIl4qWAgXso8xEuIF/FSwUJSKo116dLFXbCsXLnSNm3aZK+88orVq1cv2cVCFEoW6uJ/69at7gLz8ssvt19//TUjltXevXst05x22mm2YMECt85D7/Jq//z888+D0z766CPr0KFDkkqZWjJxO4n05ptvuoRUp06d7IcffnBJvmnTptnSpUvtjz/+SHbxgAKLeKngIF7KLMRLeUe8RLxUEJGUSlMbNmywn376yfr06WOlSpWyww47zNXAufjii93ra9eutX/84x9WuXJlq127tt1xxx22b9++sHlMnDjRatWqZZUqVbJBgwaFvfb888/bscce6+7Ut2vXzhYuXOjr70tXqk2i4FjLVXduI23bts2uueYaq1atmntce+21tmPHjuDr8+fPt7Zt27rPN2rUyF544YXga2p2c/bZZ7ttoly5cnbkkUe6RIia4R199NGu5oW2A8+KFSusY8eO7r0VK1Z08925c2fw9f/973/WrFkzVwtPF8mrV6/OtlZM//797eqrrw57fdKkSe57a9asaQMGDAi+7hk1apS7AE9HJ554oqsZ9+mnn7q/lWBQLalu3bq5RJQEAgH75JNPXECm/Uv/az1omU2YMCE4L72mWldaD6oVqdqRcvLJJ7v/27Rp477r3nvvDS77Z555xiXANF37thIaf//739082rdvb2vWrIla20rbkOZ//fXXu21Mxw4lRTzvv/++HX/88a7GX1ZWll133XVhZdG61ndOmTIlWPbsfpe3vWoeev22226zTKbt4YYbbrBbb73V7U9a16IasM8++6zVqVPH1barXr26W4/NmzcPbkui9zRt2tSGDh3qPlu1alW37j777DNr3Lix28979eplBw4cCH5GNzK0XvRa7969DzpH5Lb+ctpWgFRBvFQwES/9f8RLxEvES/8f8VIBFUBaOnDgQOCYY44JnH766YFp06YFfvnll7DXO3ToELjssssC27Ztc681atQoMGLECPfaRx99FChcuHBgwIABgb/++iuwdOnSQKlSpdx0mT17duDwww93/+/Zsyfw8MMPBypXrhzYvHlzUn5rQde+fXu3DGX//v2BGTNmuOW7du1at8zLlSsXfG+PHj0Cp512WmDDhg2B9evXu8/27t3bvbZp06ZApUqVAo888ohbLx9//HGgdOnSgTlz5rjXhw0bFihatGjglVdeCezbty8wZMiQQI0aNQJXX311YPv27YH//e9/geLFiwcWLFjg3t+tW7dAnz593Lz0+OyzzwK7d+92r9WpUydQt27dwHfffRfYsWNHoHv37q5csmLFioAOLSqP58YbbwxcddVVYa937drVvUefX7x4sfvN2h492n6nT58eSFfnnntuYNCgQe75U089FejVq5fbp0455RQ37dtvv3X74erVqwMVK1Z0+7HWm5ZVtWrVAh988IF7X+vWrQP33HOP23Z27drl5uHRcv7666+Df3vL/oorrgiu82LFigXatWsXWLJkift8x44dA/369Ys6D28b8sry3HPPufW2detW97rKNXnyZPdc89c2k9028ccff+T4u/Rdhx12WGDSpEmBvXv3uu0kky1fvtwtwx9//DHb9zzzzDPuOKz99f7773fL11s3Wo5anmPHjnXLc+LEiYGyZcsGLr74Ync8+f333wNVqlRxxwfv+7RtvP766+7948ePd5/Xeol1/eW0rQCpgnip4CBeIl4iXiJeyg3xUsFEUiqN6aLhpptucgknXdwee+yxgffeey/w22+/uYubNWvWBN87ZcqUQP369d1zJUIKFSoUdhGoC9XRo0e75//85z8D1157bdh3NWjQwM0D8QVZJUqUcMkn/a91NWrUqOC68JJSSjroIvHzzz8PflYX/Uok6bXnn38+0LBhw7B5K2HlJa10kdiqVavga0pIaDtYtmxZcFqLFi0CEyZMcM+VaFLi5Pvvvz+ozEpK3XfffcG/tS1pXqtWrYo5KRWaLJGTTz7ZXTjL3Llz3QWvkiTpSolILW9RgljJHP1ere+dO3e65GKzZs1cckEJvFC33357oGfPnu65klhax1r2kbJLSkWu89tuuy3497hx4wJt27aNOg9tQy1btgy7mNM2OX/+fPd37dq1A0OHDg2sW7curBzRtoncfpe+64QTToh5eaY7JZe1DHWjIFbly5cPJqW1b1WtWjX4mo7vmt8777wTnKYE1R133OGe33333YGzzjorbH46vnhJqVjWX07bCpBKiJcKBuKl/0O8RLxEvJQ94qWCieZ7aUzNMx588EHXzGr9+vWuKdT5559vv/32m2sqpOY1HjXl0XSPmn+o2Z+ndOnSrumY6H0aTSSUmoKFfh55M3LkSNfc7a+//nLN9p577jl78sknw96jdajOjEOXvdbb7t27XfODaOslcr2GrnNv/UZO8/o5euCBB6xGjRquCZ/mq+Y4oU171FwodL7qrP3333+P+TerOU+onj17uiZGov/Vr5bmma68ZnnqR2z27NnBDlzVxGru3LmuaaXeoyZ3b7/9tmsC5T0eeeSRYB9CaoqnztHVXEtNudQ5em4i13l220B2x5XQ5hMlS5YMHhtee+01W7JkiR1zzDGuieL06dOznU9uvyvaNpLJvOZ62e1j2jfV/LZ+/fru+K3luWXLFndsiGf/V3Pc0H1cQv+OZf3ltK0AqYR4qeAgXiJeIl4iXsoJ8VLBRFIqQ6jPDyUV1P+Q+gXRRaz6lQq9wFB/L7HQ+/T+UHn5PHKmvlk6d+7sOjUOpf6/ihUrFrbs9VyJDB2AE71eqlSpYo8//rjrcP2NN96wJ554wiUdPKEdsa9bt84lx5TEUp9BEtr/VLROmAsXDj/8qD8l9YmlTpvV90yPHj0snanvJV3IK6mk9ar+20R9OqkvIK8/KU1XMllJS++hC3slBOSoo46yyZMnu36g1A/cLbfc4jpR9xIBflIfY+qHSImQIUOG2GWXXeaOM5HrWnL7XRLtc5mqQYMGLjn84osvRn196tSp7vHWW2+5ZJSWp/qC+r/KbnmnvqkiB1vQoBl5WX9AQUS8VHAQLxEvES/9H+Kl/494qWAi4k9TGm3vzjvvtGXLltn+/ftdguChhx5ywZZqYuhiVxevSlLpQmPEiBF21VVXxTRvdaarjorVQa4SXI8++qj9+eefLpGCQ+fVQGjSpEnYdJ1wdJGv2hAbN250y/z222+3K6+80r2m5a/kkBJJWi/qRFvrqXv37nGVQ7VctG3oolbJE3WWX6RIkeDrqsmlWl2q3aXOl0855RSXAFOCTDVcVNtLtTcUMMRyoaraHRdeeKH7jap5p5o26UwJIyWgNJKQ7vp5NO3pp592F/laplq/H374oUv2aEQVPdTx+FdffeXer4SUEj+an9aTtgWtK68WjAY88INq8f3nP/9xxx6VQWURbTNKqGpaaFly+10Ip/WrY606tPWOufL999+7Dsq1bJXc1P6ndXH33XcfUq0kDYQxa9Ysl+TS8USdmOu7WH9IN8RLBRfxEvES8RIiES8VTCSl0pQuTtTMQ4kK3S1XkkBJpJkzZ7qmeLqjrmSCmmNoVDWN+BY5wl52dNGsiyJdCGlkPt2513y9i1DknZI6qmGkh0YzVJM5jZIVaezYsa62hEbW02iKukuoZKNo9DytB42MqPWiUfrGjx/v5hcP1bbxRm5r3bq1W9/nnntuWHM71W5S4kPbmjeimqj2j0bX07an5NWll14a03fqOxYtWpT2taQ8Sg6rhpP2KY+WtZKOao6nUexU++zdd991y1EjLmp59+3b1zX7kw8++MBOOOEEt57OO+881+xSiWcZPny4G7FN24aSGflNxxVtkyp3v3793N/aFtVsa9iwYa4JsY4Tmp7b78LBNLqO9nElilRDTsvyoosucs02NSKfjgk6pqvZrpb5odReVRNMJRm1/WgdfvHFF3bmmWcGX2f9IV0QLxUsxEv/h3iJeIl4KXvESwVPIXUslexCAICoZpb6xFF/NroQBgAAQDjiJQDphKQUgJSgZqaqlaE7P6qhAQAAAOIlAOnt/3cQAwBJsmLFCmvcuLHrS4qOkgEAAIiXAGQGakoBAAAAAADAd3R0DgAAAAAAAN+RlCpgNNKdhuouaHr37u2GFAcAAMhvxEsAABQMNN8rQA4cOOCGW//vf/9rTZo0sYLkl19+sTZt2ri+g4oXL57s4gAAgDRFvAQAQMFBTakCRB1AV6xYscAlpKRu3brWoEEDe/nll5NdFAAAkMaIlwAAKDhIShUgr7/+unXo0CH4d6FChWzs2LF2zDHHWPny5e2SSy6xLVu2BF+fP3++tW3b1r3WqFEje+GFF4Kv/fvf/7azzz7b+vTpY+XKlXOjnn388cc2Y8YMVxurQoUKdscddwTf/+yzz1rTpk3t9ttvt0qVKlnt2rXt8ccfP6iq/PHHH+++r0WLFjZ37tyw108//XT3GwAAAPIL8RIAAAUHSakC5JtvvrGGDRuGTfvPf/5jH330kWset2nTJuvfv7+bvnnzZjvzzDPt0ksvtfXr19v48eNdv06fffZZ8LPvvfeederUyTZu3GhXXnmlXXHFFa5p4KJFi9z7HnzwQVu4cGHw/UuWLHGJsD/++MOmTZtmt912m33yySfBu5K33HKLS15pfoMHD7ZzzjnH/vzzz+DnlRjTbwAAAMgvxEsAABQcJKUKECWdypYtGzZt0KBBVr16dVc7afjw4TZ16lTXl8Jbb71llStXtn79+lnRokWtffv2dtlll9lzzz0X/Gzz5s3tggsusMMOO8wlr37//XeXaCpdurRLIKnWU2hSStNVw6pYsWLWunVru/zyy23y5MnutXHjxtnAgQOtWbNmVrhwYTdfJdCUrPKo7PoNAAAA+YV4CQCAgoOkVAGiJnVbt24Nm1anTp2w53v27HE1o3777TfXj1OoevXquemerKys4PNSpUpFnbZ9+/bg30p+KcEV+n1KZIlqaqlpn5Jj3kN3Kr3XRWXXbwAAAMgvxEsAABQcRZJdAMROfTotW7YsbNqvv/5qLVu2dM9XrlzpajGphlTNmjVdoiiU/tb0eK1evdr27t0bTEzp+2rUqOGe16pVy9XKuvbaa7P9/NKlS91vAAAAyC/ESwAAFBzUlCpA1EeT+o8K9cADD7hkkfqQGjp0qGuGp+ZznTt3tnXr1rnOyPft22effvqpTZkyxbp37x739+/YscM1EVRtrC+++MLNT034pG/fvq4sCxYssEAgYDt37rQPPvggrGbWhx9+6DpXBwAAyC/ESwAAFBwkpQoQJZo2bNjgOhz3qHPy0047zTWlK1OmjBuNz6u6PnPmTHv++efdaHnXXHON6+y8Xbt2cX9/48aNXYKrWrVqdtFFF9mIESPcd3sB4KhRo1xn6vpujeansqh/K69Gl2p5XXzxxYe8HAAAALJDvAQAQMFRKKBqLSgwXnjhBZsxY4Yb/U4j4X399de+NInTqHpjxoyJe/Q8JcVatGjhklYAAAD5iXgJAICCgT6lCphu3bq5R0Hz1FNPJbsIAAAgQxAvAQBQMNB8DwAAAAAAAL6j+R4AAAAAAAB8R00pAAAAAAAA+I6kFAAAAAAAAHxHUgoAAAAAAAC+IykFAAAAAAAA35GUAgAAAAAAgO9ISgEAAAAAAMB3JKUAAAAAAADgO5JSAAAAAAAA8B1JKQAAAAAAAJjf/h9gctXEH3JndgAAAABJRU5ErkJggg==", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "%matplotlib inline\n", + "import numpy as np\n", + "import matplotlib.pyplot as plt\n", + "\n", + "# DGP illustre : 1854 Soho, 4 quartiers fictifs, 1000 habitants chacun\n", + "rng = np.random.default_rng(0)\n", + "quartiers = ['Soho (pompe Broad)', 'Bloomsbury (autre pompe)', 'Westminster (riviere)', 'Camden (puits profond)']\n", + "exposition_pompe = np.array([0.85, 0.10, 0.30, 0.05]) # proportion buvant la pompe\n", + "risque_base = np.array([0.005, 0.005, 0.005, 0.005]) # risque hors contamination\n", + "contamination = np.array([0.080, 0.001, 0.003, 0.000]) # ajout de risque si eau contaminee\n", + "\n", + "n_quartier = 1000\n", + "cas = np.zeros(4)\n", + "for i in range(4):\n", + " exposes = rng.random(n_quartier) < exposition_pompe[i]\n", + " cas[i] = ((rng.random(n_quartier) < risque_base[i] + exposes * contamination[i])).sum()\n", + "\n", + "taux_par_quartier = cas / n_quartier\n", + "print('Taux de cholera par quartier (DGP illustre) :')\n", + "for q, t in zip(quartiers, taux_par_quartier):\n", + " print(f' {q:35s} : {t*100:.2f} %')\n", + "\n", + "fig, axes = plt.subplots(1, 2, figsize=(12, 4))\n", + "axes[0].bar(range(4), cas, color=['#c0392b', '#27ae60', '#f39c12', '#2980b9'])\n", + "axes[0].set_xticks(range(4))\n", + "axes[0].set_xticklabels(['Soho\\n(pompe)', 'Bloomsbury', 'Westminster', 'Camden'], rotation=0, fontsize=9)\n", + "axes[0].set_ylabel('Cas de cholera (n)')\n", + "axes[0].set_title('Figure 1.1. - Cas observes par quartier (1854, DGP illustre)')\n", + "axes[1].bar(range(4), taux_par_quartier * 100, color=['#c0392b', '#27ae60', '#f39c12', '#2980b9'])\n", + "axes[1].set_xticks(range(4))\n", + "axes[1].set_xticklabels(['Soho\\n(pompe)', 'Bloomsbury', 'Westminster', 'Camden'], rotation=0, fontsize=9)\n", + "axes[1].set_ylabel('Taux (%)')\n", + "axes[1].set_title('Figure 1.2. - Taux (cas / 1000 hab.)')\n", + "plt.tight_layout()\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "id": "5409537f", + "metadata": {}, + "source": [ + "## 2. DAG 4 noeuds — un exemple canonique (mentorat)\n", + "\n", + "Pour la suite, on travaillera sur un **exemple canonique** : l'effet d'un programme de **mentorat** (X) sur un **score final** (Y) pour des etudiants. Le DAG suivant est le plus simple qui pose les trois difficultes classiques de l'inference causale :\n", + "\n", + "```\n", + " Z (aptitude)\n", + " / \\\n", + " v v\n", + " X Y\n", + " | ^\n", + " v |\n", + " M (effort)\n", + "```\n", + "\n", + "**Noeuds** :\n", + "\n", + "- **Z** (aptitude) : confondeur latent. Les etudiants plus aptes sont plus susceptibles d'etre selectionnes pour le mentorat (Z -> X) et de reussir independamment (Z -> Y). Sans ajustement, X et Y sont correles par le detour Z.\n", + "- **X** (mentorat) : traitement. C'est la variable qu'on veut evaluer causalement.\n", + "- **M** (effort) : mediateur entre X et Y. Le mentorat agit en partie via l'effort deploye. M est sur le chemin causal X -> M -> Y.\n", + "- **Y** (score final) : outcome.\n", + "\n", + "**Aretes** :\n", + "\n", + "- **Z -> X** : l'aptitude influence la selection au mentorat (les bons sont choisis).\n", + "- **Z -> Y** : l'aptitude influence directement le score (memes sans mentorat).\n", + "- **X -> M** : le mentorat modifie l'effort.\n", + "- **M -> Y** : l'effort modifie le score.\n", + "\n", + "**Question causale** : quel est l'effet d'**imposer** X = 1 (tous les etudiants ont un mentor) sur Y, **toutes choses egales par ailleurs** ? Formellement : P(Y | do(X = 1))." + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "094b2f1f", + "metadata": { + "execution": { + "iopub.execute_input": "2026-10-05T19:53:14.226245Z", + "iopub.status.busy": "2026-10-05T19:53:14.225741Z", + "iopub.status.idle": "2026-10-05T19:53:14.935018Z", + "shell.execute_reply": "2026-10-05T19:53:14.933894Z" + } + }, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAArIAAAGGCAYAAACHemKmAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjgsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvwVt1zgAAAAlwSFlzAAAPYQAAD2EBqD+naQAAZWJJREFUeJzt3QV4VNfWBuAV9xDB3d3d3SlaCqXQlkKFUvdr/217b2/11qkh7W0LtKW0xd3dLbh70ADxic3/fBvOdDKZhEkyfr73eQIZzfGzzj5rr+1jNBqNQkRERETkYXxdPQFEREREREXBQJaIiIiIPBIDWSIiIiLySAxkiYiIiMgjMZAlIiIiIo/EQJaIiIiIPBIDWSIiIiLySAxkiYiIiMgjMZAlIiIiIo/EQJZc6o033hAfHx+PXQuzZs2SmJgYSU5Odsnf79q1q/pxhjVr1qh1hf8L6+uvv5bKlSuLwWAQb4N1X7p0aZkxY4arJ8WjefM2Ys/t5dixY9K7d28pUaKE2h/nzJmjnt++fbu0b99ewsLC1PN79uwRR7p+/br6W4sWLRJvkpOTIw0bNpT//Oc/4snuv/9+GTFihOiBVwey//vf/9QOrf0EBwdL+fLlpU+fPvLZZ59JUlJSgZ9v3bq1+txXX31V4PvWr1+vNpgKFSpIYGCgOsC0adNG/vWvf8nly5cLPd2PPfaY+rv33HNPoT9LtgfPd/u5m+zsbHn99dflmWeekfDwcIct+oMHD6ppPn369F3fe/HiRfVeR5/ECmvs2LGSkZEh33zzjbgDjMzdqVMnKVWqlDohW5owYYIEBATYtBw//fRTiYiIUCcOd4aAA9uGq7399tum4MudtxHzY4Wvr6+cO3cuz+uJiYkSEhKi3vP000/b9J3F3V4efvhhiYuLU4HWjz/+KC1btpTMzEy57777JCEhQT7++GP1fJUqVcSR201sbKw8+uij8n//93/iDt599121HpYuXWr19f79+6tzM46RBfnpp5/UurZ1fbqr1157TX777TfZu3eveD2jF/vuu++MmMV//etfxh9//NH47bffGt9++21j7969jT4+PsYqVaoY9+7da/WzR48eVZ+tWrWqsUOHDvn+jf/7v/9T76tevbrxb3/7m3Hq1KnGSZMmGR955BFjZGSker4wtm/fbvT39zcGBwcbBwwYYPR2r7/+ulp+zoR1ju3B2s8rr7yipqdNmzZ3/Z4//vhDbUfnz5936PT++uuvappWr16d5zWDwaB+zLcfvBfbvr3h7+c3HbZ49dVX1T6Xk5NjdAcHDhwwBgQEGMeOHZvr+U2bNqn1+tJLL931OzIyMoylSpVSxxV399RTTzl9X7MmLCzM+PDDD3vENqIdn3A8fu+99/K8jv0Mr+E9WL6O3l5SU1PV3/r73/+e6/lDhw6p56dMmWJ05nZz8OBB9drKlSuNroZl26hRI3XOxXIyN2vWLDWdX3zxxV2/p0mTJsbHH3/c6A1at25tfPDBB43ezvVHNScEsji5W8KOFxISog6alhs9/POf/zSWLl3a+Ntvv6mT2qlTp/K85+eff1bfP2LEiFzBhObmzZvqQGgrHLzbtWtnHDdunJouBrLOlZycbKxTp46xRIkSxpMnT971/YMGDTJ27NjR4dNVUCBryZ0D2R07dhT7pIcTNk5Y9oKLT0zTmjVr1GN8d8OGDY2VK1dW28Pd/P777+rzx48fN+oxkM3OzjampaXZLZB1t21EC2SHDRtmbNq0aZ7Xe/XqZbz33nttDmSLu72cOXNGff6DDz7I9fzatWvV8zhW2Iu2/d9tu8H+Upxg6erVq8aLFy8a7WHz5s1GX19f41//+lfTc4mJicby5csb27Ztq7bXguzatUvN64oVK4zuxpbjkaX//ve/an9LSkoyejPdBrKAq2K8Pnny5Dyv1axZ0zhx4kQVoEZFRRn/85//5HlP7dq1jSVLlrTbRvL9998bIyIijPHx8YUKZLX3rl+/3tiqVStjUFCQsVq1aur7LJ04ccI4fPhwY3R0tArk0fK4YMGCPO9LT09XwXyNGjWMgYGBxooVK6rWSjyvQXCfX9CE5y2DeExfy5Yt1fThqvnrr7+22iK7bNky1QqOgBI7IZaz+YGpoAMiTmIpKSnGonjooYfUtPzyyy93fS9O3lgub7zxRp7X0PLfrVs31fKC99SrV8/45Zdf5rveli5dqloBsFzwXlw8WW7Dlj9aMNmlSxf1Yx5oWv5o6wd/z1oAYf4dmnPnzhkHDx5sDA0NVfPx/PPPG5csWWI1kN2yZYuxT58+6g4EtqnOnTsbN2zYYHW5xcTEGJ999lljUWH6MT1oLUVrUHFhPWIbxwUM9vV33nlHzeO8efNs3mZw18badGLbReCBdYzfcTLF3RrYt2+f2kawfBE0z5gxI8933Lhxw/jcc8+pfQ/bEabz3XffzXUy1vZBBDbffPON2q/wXuxn27ZtyzU91rYN85Pkiy++aPpb2OfwnZYto1rANn36dGP9+vXV3SPcmQC8HxfiWMdopWzevHmewMraNFhuk+60jWjHp9mzZ6v/cXzR4Djt5+en9ldbA9n8thfAdyMoxrEZx4IWLVoY586dm2dazH+0fdryefP9GRcFuODGtoZjKi7ALZeL9t24SzFq1Ch1zkPgfrftBl544QX1/qK2ouN4gu0I04X5zczMNBYHztu404J5AWxL+H7sc3eDcx62f8sLIQTD2BexvPE6tq+ePXsad+7cmedY2K9fP7U8sLzRQvzJJ5/kek9x1ocGdw+xf2E/i46ONo4cOdJ49uxZq3cf8T24gPJmug5kcbLG6wjsLDdGPI/AC9BCioO2uSNHjqj3PProo3aZVuwoZcuWVSdSKGwgixNxmTJlVAsTTpbYyNGSvH//ftP7Ll26pN6DYBm3pj766CMVQOEK1nxDx4kS6RfY0RDA4AT59NNPq4MBgpuiBLI4iCDIwUkb8/jvf/9bTUvjxo1zHRgxvdqJ+NNPP1XB7ssvv6yCo7vRdv6itBj+73//U5997LHHbHo/ArX8Ah5cTOB29ccff2z8/PPP1bLEe7Ugxny9IWDAQeovf/mLWh848GF9IJjXLjxwIMbnsW61FAisS8sgFM8hjQbvxa0x7b34jsIEsrhDgenCQRK3enEgxklVW1fmyxcHZawvBDAffvihmme8D89t3bo1z9/CwR/fVVRosUPrD7ZNTEv79u2N06ZNK9bFJC4k8F3jx49X2+jQoUNt/iwueNFaZwnLGcsPx40JEyaoW5qYVm1/QVCLC0NsHw0aNFABkfldAFyMYTnGxsaq9Y79AEEQ9mmcUC33wWbNmqlpwe3v999/X11gIyjVTshIl0DrId5rnkoDCEC6d++uvhvHM2ynAwcOVO/F/m8Oz+FiCyfyN998U83X7t271Wv4ewgi8Hlsy7itifebXyjjbyJI69Spk2kaMG3uuo1ox5QrV66o+UMqmQb7BQIRXNzbGsjmt73guIfvwvaCdYhliGMe1ol2bEZQgv0LfwvBDZYdLiKw/LQ7CzhW4Hnt+LF8+XJ13Mb+jO0C6wzbBoIf87uM2nzi7+MYjwtvrNuCthsNLmrwelxcnLEocMGG6ULjC76nXLly6niI9L6iuHXrltq/ECxiW8C+he+zBbY9nDstPfDAA+qYhos9pA9iHWEfwbxrsMzxHhxnsTy/+uortT7wnZrirg9466231HaB4BXPv3nnO3CBhGVpDhcFOKbZkiblyXQdyAIOHjgJmEPQVqlSJdMVJjZQfI92wAZcOeI5y6stfAYtg+Y/tlxhIljDjqy1eBY2kMW0rFu3zvQcDrw4YZhvwDgpmQfogIM7/i52Aq2lBwcqBFPm7wOcTPH5jRs3FjqQHTJkiDqxo4VKg6tQHGTMA1ntQI3lVlhFDWTREoIWMwQU1tJMrMHBLL+Dt7XvQIulZb60tt7MW2BxEMaB3HybLCi1wDIILSi1wNZAFts0vgN5ZeaBFU7C5tOBbb1WrVpq3sxbYzD/2KZwArSEABsH1uLCcsIFFu4oYJrCw8NVIGoZFNkKgQG+Bxd5uMC1Bfbr/HJptZYs81xInGQw7/gM0pI0hw8fzrO/4EIP26TlyRwnZOwzWuuLtg8i4E1ISMhzfJo/f77pufxuEc+ZM0c9jxOkOVzgY1rNb4PjfTg2aK1dBW33WpoGgmRbUwvcbRvRjik4HuEYjX3A/IIVfSHAlkC2oO2lR48e6iLW/I4X9ikE4djHrLXAm9Puxli2gKMVDyly169fNz2HgBjrEBdGlvOJ/cDS3VILsDxtvZNVEMzvqlWrjGPGjFHrH9+JYB53Fm09Lmu0FnS07lvLmc0PLlbQKm4tTiho/WZlZaljHo6xlsGk+bGxuOvj9OnTav+3vEMcFxenAuT87hyjldibeXXVAlugt7l59YKsrCz55ZdfZOTIkaae6927d89TLgW9VbXPm7t165bqCW3+c7eez0ePHlU9WT/44AMJCgoq0nzUr19f9cLW4O/WqVNHTp48mav3KSoxdOzYMdf8P/7446pHPHrHw6+//ir16tWTunXryrVr10w/WA6wevXqQk0bevejJ+mQIUNUeR0N/gYqSJiLiopS/8+dO1eVQSkM9KzFOaUw5ajS09PVusbfwnpHD2RbaD3do6Oj87xm/h3YHrDsunTpotYFHptDFY2hQ4eaHkdGRspDDz0ku3fvlkuXLokrYDspV66cDB8+3PRcaGio2k7MYbtGKaAHHnhALQ9tO0lJSZEePXrIunXr8qxDLK+0tDRJTU0t1jRiOWF6tmzZorZbVBlYsGCBKj/UoEEDmTp1aqG+r2TJkqb9qGLFijZ9Bj3Esb1Z2wY06NVtvm1jn0TJIvOyOHgOr5nvq9gHsT/ju833wZ49e6r9CcvWHLZh8+nQjgXm31nQ+vbz85Nnn3021/MvvfSSmr/Fixfneh7bMpZTQdv9jRs31LaO6di1a5cUhrtuI9jOjx8/rspcaf/jOVvlt73g+VWrVqltAucibV1jn8LxEfvYhQsXpLDi4+PVPopqECgRqGncuLH06tXLatksLKPC0uYH01wcON9269ZNVVzAsQ/l2FCKDVUacDx68skn1XZli3vvvVdVKcCy/eKLLwp1XLe2P2P/3Lp1a74VD3C8PnXqlDz//POmc5j5fNlrffz+++/qmIptxfy4ULZsWalVq5bVc7N2DPFmug9kUdMPpVA0y5Ytk6tXr6qADwcr/GADxQ6GshzaiVn7jGX9UASGy5cvVz+vvPKKTSvhueeeUwdX7HxFZR4gmm/A5jv+mTNn1EnTEgJK7XXAgfPAgQN5AvLatWur169cuVKoacPyxIkJO5oly+nBCblDhw4qAChTpowqUYNarYUNam2FA8++ffvkk08+USe3wrrdGJPbxo0bVcCBgAUHNSy7v/3tb+o1y0C2Zs2aeUp9acvZlnJbjoDtwNp0Wa4rbCeAE43ltoIgASchy/nVlldB5c2wT+FEpv1g+ykItl9cBGK5t2vXTgUtkyZNsnl+d+zYoU52qB2Jk9X06dOluNsAoNwfloU5lP9BoGw5/3jefF/Fsl2yZEme5Yrtyto+aLn/aydjW078WN+4oDI/Dlo7LmiqVatm9XsQJLZt21bNN07UmF6ULrTcBu7GHbcRaNasmbq4nzlzpmrUQPCgXdwXZ3vBOQbPoYyV5fpGeb+iHHPN11t+x3ztotOWdVvc9YWLL/P1hR+UWivoIuSJJ55QNav/8Y9/qIYjBLaW22JBWrVqpf5HebKizI+5999/X/bv3y+VKlVSsQEaTcwvEk+cOKH+xzHEkesDxwVMH86lltvKoUOHrG4neL8n12q3hb/o2Pnz59VBFidtjdbqml8h4bVr16qgFgc0wMZtzt/f33SywfffDa7EccLClZZ54IKWYQR/eA4nBezYBUGLSmFOsgVB0NioUSP56KOPrL6OnRny2zlw0CoqXDmjtQlXlgsXLlTLBi2lOGHgIiO/+SwKtHqhZiXWtWVr492ghqIWKJi34OGAhtZIbB9YflhWqC2Mq23Ud3RUQG6LgtZXUZarNi8IEJo2bWr1PZZ3LLC80LpbUAvJf//7X3nzzTdNj1EPM7+gHi3q2He+++47WblypQqixowZo1pvbIF5x7pHIIcgB4Xm0RKJGs6WLSuWsF9imeYXLOa3TG3ZV7Fs0Urz6quvWn2vdrFTmO+0F2vrDrW0Bw0aJJ07d5Yvv/xStaChDi/WCwK/wnC3bcQcWmARnCPox0U36svaKr/tRduPXn755Tx3qDTm5yhHsrXl0pw2P9pdDWtQl9UyKMMxPr+7Z2jt/vbbb+Xnn3+Wmzdvqrrs48ePN11cOQqO69b2Z5wjcHfhjz/+UOchHPPee+89tV3169fPaesD2wq2IdwlsbbPh1upZ475sdaI5E10HcjiFgZoBw9cDeGWNg5Q5rdVNbj1hkAXgSyuqrBxoLA3WvPQ+lYUZ8+eVf8PGzYsz2u4nYSdHwEQWg6LCwf7I0eO5Hn+8OHDptehRo0aqogyArKCruS0Vh8caMxZXjXjahE7pNaCZ87a9ODkgL+NHwSDKKD+97//XR34tIuE4sLVNAaewPKdPHlyoT+vXcigtR5Bv2b+/PmqJXLevHm5WsnyS8fQWmPMlzNSTaBq1arq/8JcTd9tfVmuK219Va9e3fQY2wEu0Cyny3JdYTsBXGTZul6wvO52MkJqhXn6i7WT67Zt21RggrskuBhFaxla2BBo3C0ANYeBUXBbECcozAdafdCC85e//EX9XhBctGIZYJ7sDd+LVkd7be8FbRtY3ytWrFC3tc1bZS2PCwVB4XUEiEghMk+PwjqydTrcdRsxh8/+85//VLeJtfOHrfLbXrR9D4G/Pde3tt7yO+Yj8LTlvGXL+oKC1hlar3GX0lyTJk1yPUZrIpYp1hnuCCKoxG14BLAFtXTaE47r+e3PuDibOHGi+sG0Nm/eXA1KgUBWOxbiuJnfOrTH+sDfwXEZ5y3Li1lr0CCGiwhcZHoz3aYWoCX03//+t9ogRo8erZ7DyQzB7FNPPaUCWcsftNLggK0NoYjbC7gdgIAII6sUpTUELY34u5Y/CP5wQsXvAwcOtMs8I2cIB/fNmzebnsP8IpBD0KTlveHqE0H0lClT8nwHWom12x848WPns8zXQ4uMOVw54mIBQb8WuANuhViOwoKcJktaa9/dhq7EusAB4W65dVhXSFnA+3CSw23dwmrRooVqacVtaXPaVbL5usdJ1NoJHZBzhXWswS20H374Qc0zDv6gHdysBaGWCnovDoLIFzS/pYfbwZYjFmE7wXTNnj3b9ByWlWXAj2WA70TrmLUheq3d7kW+JNJoCoITO04G2g9STTTY/3BSQwsNWmuw7+I78YMTTGECFMw3ghIc5JG/DVjuuGDFto80g7vBbWrLbcAesA9iP7U2ShHWLU5QhZXftoH1jZZpy1vtuIBGEGNLixO2e7zX/G4MWkitjeCF6ShoW3anbcQStnc0XLzzzjvqFnNhWdte0P8CLZO4O4QA2dLd0ibyg8AL2/P333+fa3kj2EKrIta7Le52/Nm5c6c6hhaUmoWLHPP1hR+tIQT7IfY/jIyJdDxMN9YbjkHYBp0VxGrrB8vH/FyDbdoyPQbrDHdxtPchqEUsgW3Dcjlp5wJ7rA80eGFfw90Iy/jCaDTmGaUQKTS4I3G3/cnT6aJFFs3wCHBw8MeQsQhicXWIKyS0nGEnA7S24iowv5WOEx5OcLjljQ0KV+fYCHFQQ4CI4AgbMwI9PI8gCS0cBXUGQaudtfxWtMAiR1Q7wdoDWpkwTTgx4WSNW13YqXAFioO/dpvswQcfVHmpSDRHSyJOEtiZsQzxPE6uWt4RclkxNCD+x3MIarUWRXPY8ZAmgNszOJlgXXz++efq4IccVQ2G9cV3DBgwQK0fXPkiMMbte/MWGGtwIsbfKeiWFSAXDbeucBGBVmJrLcWATlj5XSFjm8FtaLRkYZo1eA4BLi4+kOOFAA/bDA581k5SuKpGiwOmB+sbt9OwjZoHvjj44eCFW1k4oKLFS+uAaO1EixM1WhOx7WH6cULHdol1hOC0b9++KlBCGgTyQbXWBA0uzLAs0eqFkxQOwGgpwe1ec9hekAuL7Qnr8ZFHHlEnI1wEYR3gQgct1Bp8Fy5UBg8eLEWFfQ8XTwj2cXFZlFuhGgwtjIM/tkNz2Ia07R9BR0FpF5gXLBts87a0kNgKJ3Qcm3DxjFYpXDTguIKhSbEOESQWdCvXGnwHYN/HhSXmC8csbKu4y4S7HvhetJThxIq7UzgOWW4f1mB/xd0TbFs4LmK/Rd4xbomb79/adGC/wfsRDGDbxDbqjttIfn0aiiq/7QXLCsc33N3B/odAHccBXMwgRa2ow4ziFjj2TwRoOM6gIQLbOwJPW4crzm+70eBcim2oqHmYOA7hIuOvf/2rjBs3znQnyhWwftDAhRRCHMsBdypw/sG2hH0Dt++x/eKY/eGHH5qOhUg5wXLA8RrHQhw3cc5E67J2QVrc9YF98a233lLLCvsq4oOIiAh1DkeDCNKkkKJivm5w3EaaklczejHLYvKo8YZarSgLhBqlqN2quXz5sipfUdAIJSjhgdqElnUmMSoQStWgbBIKMaMwPOqgooQGimYXRVEGRLCl0L02IAJql6IcFmo9WhsQAaVzUCsPJalQxgt17lDbETXrUNbGfJmgpA3Kk6B0EUY5Q+kvy3JC2ugz+A6sh/wGREBdUtTNQx1AvA//owSJLTUFbS2/hWVirci35Y+10dzMob4jyulYFqJGbVnUAMXyRVkzLEcMkmD5neYDIuD9WM5169a1OjoPhp7EMtPKlVkbEMG89JJWrN6yFBdqvVaoUEH9LQw6gTqL1r4DZdJQqBvbO2oUonZpfgMioCwdamOiBBS+F/OF7cBydKbXXntN1REuzvCjRRndxhrU38S8YOSbgsr3oB5qQTCIApYPymVZGxDBEpYz9ilb9mGUxsNAICj5hH0BfwflmDDNWn3Y/MoxgeU+iBJBzzzzjKoBi+3WfL/D30Jhe+xvOIah5FNBAyJYg1qt+Jy2HWO7szbgCcqNoaySVmLJvBSXO20jluW3CmJrHdn8thft2IwSTDhHYR1gP73nnnvUtljU8luAUaqwr2N549yE+qf5FeC3Np8FbTfa0LjFGQkL55C7jbjlyHVnCcdinNPM1xlqPqPmOs5x2K/xu7VBblBfHPGF9j58F2pF23N9AEo2ok4u/kZYWJja37D9ob69OZSeQzkzb+eDf1wdTBN5IrRSIx0DrZu4ii8stDzgthlu73s73ILD/OKuQHFatNwR1j1a0NGyb8/OiHrjzduIt24vaLHHHTS0pHtLz3i0mCO9EGlwxUlDcbU9e/aolAe0dufXGddb6DZHlqi4cBJCWgFuC1rLEaU/4cSNzixFqVPp7l544QW1/pHXR0XnzduIN24vyMdEahFudXtLEAvIqUa6H47rnuzdd99V6RDeHsQCW2SJXERPLbJERESOwBZZIiIiIvJIbJElIiIiIo/EFlkiIiIi8kgMZImIiIjIIzGQJSIiIiKPxECWiIiIiDwSA1kiIiIi8kgMZImIiIjIIzGQJSIiIiKPxECWiIiIiDwSA1kiIiIi8kgMZImIiIjIIzGQJSIiIiKPxECWiIiIiDwSA1kiIiIi8kgMZImIiIjIIzGQJSIiIiKPxECWiIiIiDwSA1kiIiIi8kgMZImIiIjIIzGQJSIiIiKPxECWiIiIiDwSA1kiIiIi8kgMZImIiIjIIzGQJSIiIiKPxECWiIiIiDwSA1kiIiIi8kgMZImIiIjIIzGQJSIiIiKPxECWiIiIiDwSA1kiIiIi8kj+rp4AIiJHMBqNkpOTI8Y7v2t8fHzEB1fxvr7qdyIi8lwMZInI4yFQzc7Jkazs7Ns/Wdnq8d34+fqKv7+f+Pvd/sFjBrdERJ7Dx2jeVEFE5EEQsKZlGCQjM8tu3xkY4C8hgUEqwCUiIvfGQJaIPAquvQ2ZmZJuyLCp1bWo0DobHBQoQQEBbKUlInJTDGSJyGMC2DSDQdIMGU79u8iiRUAbEhTEgJaIyM0wkCUit5eZlSXJaWmSk+O6TChfXx8JDwmRAH92LSAichcMZInIrVthU9MNkp7h3FbYggQHBkpoMFtniYjcAQNZInJL7tAKmx+2zhIRuQcGskTkdtCZKzk1TdxdeGiI6gxGRESuwUCWiNwK0ghS0tLFU4SFBKt0AyIicj4OUUtEbsPTgljA9LpTDi8RkZ4wkCUit0kn8LQgVoPpxvQTEZFzMZAlIvfo2OUBObEFwfRjPoiIyHkYyBKRy0tsoTqBN8B8cNRvIiLnYSBLRC6FOrHuWGKrKDAfmB8iInIOBrJE5DK4Fe9tHaUwP0wxICJyDgayROQS3pRSYIkpBkREzsFAlohcIs3gPSkFljBfmD8iInIsBrJE5JLW2DSDd6UUWEo3ZLDjFxGRgzGQJSKn00PNVaNO5pOIyJUYyBKRS1or9UAv80lE5CoMZInIqbKysiU7J0cXSx3zifklIiLHYCBLRE6VlqGvTlB6m18iImdiIEtETu3klZGpr2FcMb8c7YuIyDH8HfS9RER5OCOloFmjRnLu3Fmb379r7z6pXKWKw+fb38/PoX+DiEiP2CJLRE6Tle1++aI+Pj66nG8iIm/AFlkichpnBHQvvPySJCYmWn1t9i+zZP/+ONPjvv36S8VKlRw+TQxkiYgcw8fI5C0icpKbSckuq1iwYvlyeWDkCMm58/dr1qoly1eukojISIf/bT9fX4mKCHf43yEi0humFhCRU+Ca2VVB7IkTJ+SJR8ebglgErz/OmOmUIBYw32wzICKyPwayROQUWhDpbMnJyfLQ6Afk1q1bppzYL7/+RmrVrq2L+Sci8mYMZInIaUO2OhtaQZ+a8IQcOXzY9NzLr74m/fr3d/60OP0vEhF5PwayROQUrri1/uEHH8jCBQtMj/v07Sev/uUv4gpMLSAisj929iIip8jMypLElFSnLe2lS5bImFH3mwJIZ3busiYyLFQC/FkohojIntgiS0Re59ixYzLh8cdMQWx4RIRTO3cREZFzMJAlIq8ZeACSEhNV5y7878rOXa6afyIiPeF9LiJyCmeEcWiBffKJJ+TY0aOm5+rWqycnT56QSZ9/ZvUzQ4cOkwoVKzp82hjGEhHZH3NkicgpEGQmJCY59G+cO3tWmjVuVKjPzJm/QDp26iSOFhMZwVZZIiI7Y2oBETnt1jpGuNJjZQDMN1MLiIjsjy2yROQ0yWlpYsjI1N0SDwoMkPCQEFdPBhGR12GLLBE5jb+fny6Xtl7nm4jI0RjIEpHT6DWg0+t8ExE5GgNZInIaR+fIuiu9zjcRkaPx6EpEToMOT4EB+qr6h/llRy8iIsdgIEtEThUSGKSrJa63+SUiciYGskTkVP7+frq51Y75xPwSEZFj6ONsQkRuJTgoUPRAL/NJROQqDGSJyKkwaMGmbXvEkJHh9Us+KCDA1ZNAROTVGMgSkdMcPnZKxj//T3nu7+/IL38slpycHK9c+pivxSvWy94DR1w9KUREXo0jexGRw11LuCmTps2U3xeuMA0jGxgYILO/+0QiI8K9ag1kZ2dL/OWrMu6Zf0hGZqb0695Rnn/iQSlftrSrJ42IyOswkCUih8nIyJTpsxfI5B9/lZTUNNPzlSqUlZcnjpWObVpIUmqqV60BBOrvfDxZlq3ZZHouKDBQHr5/sIwfNVRCQzlULRGRvTCQJSKHBHMr12+RD7/6Xs5fvGx6PjwsVCY8PEJGDe2vWmQhJS1d0r0oXzY4MFDlxv6+aIVMmjpTbtxKNL1WKjZannt8jAzs3VV8dVK5gYjIkRjIEpFdHTp6Ut7/4lvZseeA6TkEbffe01OeGjdKYqOj8gS9N5OTJSfndsqBJ/P19ZGo8HDTAAhJySmqNXr67IWSlZVlel+DOjXltWfGSbNG9Vw4tUREno+BLBHZxbXrN+TzaTPlj0UrTXmw0KZFY3nlqUekTo2q+X42MytLElM8P8UgMixUAvzzjlx29ny8ap1etWFrruf7du8oLzB/loioyBjIElGxGAwZMv23BTLlx9m58mArVygnLz81Vrq2b2XTEK2GzExJNvu8pwkPDblrua2tu/bJ+5O+k6MnTpueQ4rFwyMHy6MPDGP+LBFRITGQJaIiQavrinW382AvxP+ZBxsRHipPPDRCHhjWXwIKWUcVubLImfU0YSHBKjfW1qoGvy9cKZOmzZCEm7nzZ599bIwM6sP8WSIiWzGQJaIi5cG+N2ma7Nx7MFce7PB7eslT40dJTFSJIi9VTwtmCxPEmkP+LFqxf5y9IFf+bP06NeS1p8dJ88b17TylRETeh4EsERUqD/azqTNkzuJVufJg27ZoovJga9eoYpel6SlpBrakE9wN8mc/+vp7Wbk+d/5sn24d5IUnHpIK5Vh/logoPwxkicimPNgff50vU6bPllSz1tIqFcvJyxMfkS7tW9qUB1sY6ACWnJbmltUMUJ0gPCTEaseuotq2K05Vezhy3CJ/dsQgGT/6Xglj/VkiojwYyBJRvtDqunztZvkIebCXruTKg53w8EgZNbRfofNgC/v3U9MNblVnFmkEocFBdg/ctfxZtHaj1Tvhxi3T8yVjkD87Wgb37cb6s0REZhjIEpFVB4+ekPc+/1Z27cudB3vfoN4y8ZH7i5UH64mts45ohc1PckqqTP5xtkyfPV8yM//Mn61Xu7q89vR4adGE+bNERMBAlohyuXo9QT6bMkPmLlmdJw/21acfkVrV7ZMHW1iYljSDQdINGeLMcBbtrsFBgRIS5JhW2IKcu4D82R9UdQhzvbu2lxcmPCQVy5Vx6vQQEbkbBrJEZMqD/WHWPJky4zdJM8uDrVqpvLw8cax0bmf/PNiiBrToDIaANjsnx2F/x8/XVwWw6Mzl6vnevnu/qhJhmT/70H2D5NExzJ8lIv1iIEukcwgMl63ZpHrOX7x01fR8RHiYPDl2pNw/pK9D82CLIysrW9IyDJJhdvu9uAID/CUkMEj8/f3EnSB/Fq3kn06Znit/NjYmSp57dLQM6ttN/Pzca5qJiByNgSyRjh04ckLe/3ya7Io7ZHrOzw95sH1k4tj7JToqUjwlGEfrbFZ29u2frGybWmvR6oqA1d/v9g8eu7r11Zb8WVSPQBWJXPmztarLq8+Mk5ZNGrh0+oiInImBLJEOXbl2Ow923tLcebDtWiIPdpzUrFZZPB3mKycnR+XTms8jAlWfOx3X3D1oLci5i5fk469/UFUlzPXq0k7lz1YqX9Zl00ZE5CwMZIl0JN1gUHmwU2f8nicPFgMadGrbwqODOz3avme/vP/5t3L4+CnTcwEB/vLQiEHy6Oh7JTws1KXTR0TkSAxkiXQALZJLV2+Sj7/Jmwc78ZGRMnJIP6eUlSLH5c/OQ/7s1BlyPeFmrvzZZx+9XX+W+bNE5I0YyBJ5uQOHj6se77vjDufKgx0xqK8KYqNKeEYeLN1dSmqayp9Fq7t5/mzdmtVU/myrpg25GInIqzCQJfLiPFj0cEdLnbkOrZupNIIaVSu5bNrIsc7HX1ajsVnmz/bs3FZefPJh5s8SkddgIEvkhXmw3/8yT6bNtMiDrVzhdh5sm+bMg9WJHXsPqPzZQ8dO5sqfffC+gfLYmOHMnyUij8dAlsir8mA3qpGg4i//mQcbGREuE8eOlBGoB8s8WH3mzy5do6pUXEu4YXo+JrqEyp8d0q8782eJyGMxkCXyAvsPHZP3Jn0re/YzD5byz5+dNuM3+X7WPMnIyDQ9X6dmVXnt6fHSqhnzZ4nI8zCQJfJgKI7/zqdTVIubOebBUkH5s6g/i9HczPXo1EZeQv5shXJceETkMRjIErkpjEyFkaYKYsjIkP6jnlQdu6CalgfbtoWTppI81c69B1U1i0NHc+fPjhk+UB5/kPmzROQZGMgSuZkco1F8zQYluH7jhkRFRlrNY8RQrItWrpP3Pv9WJj5yv4wY3Id5sGT7tpaTo1rzP508PU/+7DPjH5Ch/Xswf5aI3BoDWSI3dejYMZm7bJlEl4iSEhHhMqRvXwkOCrL63uTkFAkPD3P6NJIX5c/O/F2+/2VurvzZ2jWQPztOWjdv5NLpIyLKDwNZIjdz49Yt2XvwoKSnG6RG1SoqkP1x9mxp0qC+tG3eXEJDQlw9ieSlLsRfkY+/+UFVvzDXvePt/NnKFZk/S0TuhYEskYtv7fqa5cFmZWXJzn1xsnTtWmnVpIn0695NPX/g6FFZv3WbdO/QXmpXr+7CKSY92LUP+bPfysEjJ0zP+fv7y4PD75Gnxo+SoMDAQqXHEBE5SsE9SYjIIfVe8aN2wDtB7JXr1yUjM1MFC3Vq1pDqlSvL9Zt/5iw2qF1boiIj5PDx43IzMZFrhRyqeeP68tPX78tbf31GSsVGmy6yUOnA/MIrP1oQm5ScLFt27ZLUtDSuMSJyCAayRE5w5MRJ2b53r/rdx8fHNLLWvkOH5KNvJsvCFSvk+19/VSkF6NjVvHEjuXkrUU6cOWP6jo6tW8uJ02fUD1q8iBwJAevgvt1lwfQv5PGH7lOtsC9NHCs+Yr2l9fS5c7J55y5Zu3mL6blzFy/K8nXrJTMzS91VQCUOIiJ78rfrtxFRHpmZmXLz1i2pXKFCrucvXr4sG7Ztk0F9eqvX1m3ZItv27FEBRMM6deTE6dOyecdOqVGlinp/xXLlpFnDhlK2dCnetiWnCQ0NURUMRg3pJyXvtM5awkXaqg0bpVrlypKWnib7/3dEBvToLjvj4qRJ/fpSIjJCVm7cIIEBAVKrWjWuPSKyG+bIEtkZ0ga0Flfz37W6r1p+4epNm+TKtWsyctAg9TjdYJCtu3bL8dOn5ZGRI1Rr1pLVa6RpgwbSpnkzridySyjbNXXmTJXPjaAVkM99+dpVuRB/SR4fM1pCgoNtro1MRFQYPKIQ2ZkWuCItwDyIPR8fL7/Mm68CVUBd2Os3bppeR2mt8mXLiK+vjyQlp0jpkiWlZrVqkpmVxXVEbmv73j0SEx2lglh0XoSs7CzZf/iI9O7SRS5cuiQLV6yU7OxsUxCrvY+IqLgYyBLZGU7Su/cfkDlLlkhSSors3LdPNu3YIYGBgRIWEiIHjx5V76tVtZpqsd1/5Ijps+jshZbZoKBA1YrVuW0b6di6FdcRua3LV6+ZKmkgLSYtPV2VjqtUvrzUq1VTVm7YaLqgwx0J7X3AgJaIiouBLJGdIChVZYd8fVX915ioKPnmx+myYdt2KVOylJSOjZUaVavKjVuJcujYcSlXprTUr11LFqxYoTp5nT5/XlZv3CQ1qlQ1jc6l/U/krpDzeuXaddNjdEZEdY3eXTqrUnLGnBxp1bSpugOx//Bhmf7773Lp6lX1XlsqIBARFYRnSSI71oNFuxNuoRoMBtOQn8P691OtU1CjSmW5ePmS7Dt0UKpVriTd2rdXIykdPXlKrl6/JvVq1VLPEXmKmlWryK64OPnul1lSKjZW/d6oXl0pW6qU/LZosbRv2UJKxcaogT7iDh+Ry1evyq59cXLmwgUZ2Kun6sQIWkk683QcIqK7YWcvIjtatnatpKalS4dWLdUIXJt37pT4K1fkwXvvzVWKC+kGVStVlPYtW5pKaeVkZ6vUAiJPdODIUUlKSVYdFF958knZtGO7xF++IkP79VX7wtbde25X6ejdW6pXqazej7xwpB+gRq227VsOEkJEVBCeNYnsAMHq7AULpUREhHTv2EGCgoIkLDRUqlaqJCfPnFWtVM0b3R6vHi2xZ86fk4NHj0nDunUlMjxcPe/LIJY8WIM6tdXdCKTUZGRmyJZdu2X4gAESER6uOnwdOXFc2rZoLnVr1lDvb9awgao9O2/ZMklOSZXoqBLSt2tXBrFEVChskSUqhPxuf67ZvFmys3OkR8cOuZ5Hx5dtu/fIoePHZcKDY9Rt1ZS0NAkPDZWs7GwpX6YMlz95neSUFIk7fFjVPQbkiZ+9cEEevm+4ypXVSnRduBSvAt8Gdeqo/HB0dHxg6BDVgktEZAvevyGyEW55aqNyoeXJ3LGTp0SLbc1fQ+UBtFTh/0nf/U++/nG6XL12XZXWYhBL3io8LEzatWihUgdQK3nfwYPSonEjUxCLfNnd+/fLpStXVWexCmXLyph7h6n3Y6AQywtHpN9ovxMRmWNqAZEVOGkmJierVAFVicDHx3TLc9nadZKYlCRlSpWS2jWqS5mSJSU2Olp8fXwlIzNTjV6EzyCuTbh5U0rGxMjIQQPl+KnTUqt6NdOACER6gFHrMHpdzapVTc/hzgQ6eWG/mb1wkVQqX0769+ih7lZgOFsN9ifsiwhwzfNntYFGtH2TiPSLqQVEFpAOgA4rCEj7duuqWpFwwkxKSpZZC+ar1tVmDRrKwWPHVHWCAT17qDzYQ8eOSdOGDdTwsoCSWsid7dS6jYSF8lYpkRZ4YtS6OUuWyjPjHlELBXmyqL2MPNlRg4dIeFiorNuyVa5cvyZpaelSoVxZ6de9uyloxSAhGFyBncKIiC2yRBYQqKIVFfl7GJ0Lxd5xAr1w+ZJ6fmjfvup9fn6+Mmv+Ajlw5Iiqk3krKVGNYISW16sJ19XoXIP79GYQS3SHFoiiTBfuaOzYt09aNm6sKhmg7Nypc+ckKjJCfp43T6XoNKpbT1X3WLx6tUyZMVNGDRmsOkdu2blLlbHr262bumtCRPrFFlkiKy1GGJFr6Zo1EhwUrE6waFHFwAXo0IXg9Oe581TnFYy81bZ5c9Pn0QsbKQm4JYoTNBFZh31l/vIVKtUG+bO1qldXQSqGcP5pzlx5/rFHJSIszPT+Y6dOSXhomGzfu1ddZKIjmfm+R0T6xECWyKJ2pZZ/h9QAjBePQQqaN2oo+w4dkhXrN0haWpo0a9RQenTsqE7C6Mxy/cZNVQ+TiAoHVT2QQoCayqhgMH/5cnXBOKRvn9s1lu+Mlod0Aox0N+OPP1TnSuTYjhw8KFewS0T6w6oFpHsIXLUgFj2m0RoL9WvXlsiICDl19qzcSkpSPatLl4yVOjVrSP/u3VUQi2D31wULJTUtTbJzcnS/LIkKq3WzpjKgRw8VxAKGc8ZdDXWCulMlBBeaCGLPnL8gmZmZakQw5M3O+P0PuZmYyIVOpGNskSUSUZ1PcJsTAW1GRoa6bYkT7NXr11V9y5rVqkm7Fs3l7AW8b7nERkdJSmqaJKemSO/OnVWrLREVHwZPQK45Bk9oXK+e6Q4J4KIRAS3SfUpERqi6s1pFAyLSJwaypDtoOUWPZ40hI0NmL1wo1SpVVuPCn4+PlzWbNktMdJRqeV25YYNcS0iQDq1aqduZtxKTxJBhUKW16tZkOgGRvSFPdtHKVarjZZ+uXaVyhfKy/8gR2b5nr7Rp1lTdLTEPcIlIv5haQLqhFVTXgtgDR4+q25Lxly/L1esJKogFVB1A72ntVifKaeGzyJFFT2q0BGFAAwaxRI6BmrPPjh8n3Tq0l5ioEmr/Q+cwVDvAsM9ERBqW3yLd0Fpv0KqzetMmVahdDXiQY5SKZcvKxu07ZOvu3RJdooQ8+dCDqtQWWl9RJqh82bKSk51zuzX3zuhERORY5oMohIWEytmLF0zD11q2xrKFlkifmFpAuoJBDDDGe+/OnaRa5crqOaQSID82JTVV9ZTWTp6oWHDtRoJ0atNGnST9GcASuQxG0/tjyRKVCjRi4ECJiozM9ToGTnj9gy9lwkP3SfWqbLUl0gu2yJKuXLx0SUKCglQQi7QCDGZQvkwZVY0gOKi8XLl2XUpERMqaTZtU60+vzp1z5dMSkWuggsjD992nhre1DGJR1eDr72fJ4pXrZdmajTJycD+Z+MhIlQZERN6NLbKkK2h9/eHX2VKrWjWVB4sgVnxEjcLVr1tX2bxzl2RlZUmJyEg1PC2GqSUi94W7JahecO+4F+TchUum5yMjwmXiI/fLiMF9VKUDIvJODGRJd1BSC7ViY6OiVYvNzn1xcuzUSRk+YID4+/urW5eoEUtEniMt3SDf/zJXvp35u/pdU61yBXnlqUekU9vbnTmJyLswkCVdD0WL25RzliyV2jWqq9qUROTZLl+9Lp9NmS7zlq7J9XyH1s1UQFuD+bNEXoWBLHklDFYwbcbvMqBXZ6lSqXyujlrIpzt3MV7Wbd0iFy9dVkNjdmrT2qXTS0T2tf/QMXlv0reyZ/9h03N+fr4yYlBflT8bVSJ3ni0ReSYGsuRVUOd13pLV8unUGXI94aa0atZQvv3k33nel5yaKidOn5b6tWpJAPNgibw2f3bp6o3y0dc/SPzlq7nyZ58cO0JGDunH/FkiD8dAlrzGjr0H5P3Pv5VDx06angsI8JdZUz+UGlUqcRQgIp1CZ7Dvf5kn02b8lit/tqqWP9umOY8PRB6KgSx5vHMXL8nHX/8gy9duzvV8z85t5cUJD0mlCuVcNm1E5D6uXEuQT5E/u2R1nvzZlyeOlZrVbteWJiLPwUCWPFZySqpMnfGb/DBrnmRmZpmer1uzmrz69DiVVkBEZOnA4ePy3qRpsjsud/7sfYP6yMSx90t0FPNniTwFA1nyyDzYuUtWy2d38mA1MdEl5LnHxsjgvt04jCwR2ZA/u0k+/uZ7uXjpz/zZiPAweXLsSLl/SF/mzxN5AAay5FG279mv8mAPHz+VKw/2oRGD5NHR90p4WKhLp4+IPC9/Fnd1ps74XQ1zq6laqbyp/qyPj49Lp5GI8sdAljwmD/ajr76XFeu25Hq+V5d28gLyYMuXddm0EZF35M9+NmWGzF2yKtfz7Vs1VQEt82eJ3BMDWXL7PNjJP86W6bPn582DfWactGrKPFgism/+7PuTvpVdcYdMz/n6In+2tzz1yCjmzxK5GQay5LZ5sHMWr1J5sAk3bpmej42JkuceHS2DmAdLRA7Mn122ZpN89LVl/myoTHh4pIwa2o/5s0RugoEsuZ1tu+Lk/S++lSPHT+fJg31szHAJCw1x6fQRkX7yZ3+cNV+moP6sWf5slYrl5OWJj0iX9i2ZP0vkYgxkPbClAEOsGu/8rkFnBJ87t8A8tWPCuQvx8uFX38vK9Vvz5MG++OTDUrFcGZdNGxHp19XrWv7s6lzH3bYtmqj82do1qogn8+bzCnk/BrJuDAeU7JwcycrOvv2Tla0e342fr6/4+/uJv9/tHz83PwglJafIFOTB/rYgVx5svVrVVR5syyYNXDp9RERw4MiJ2/mz+w7mzp8d2FsmjrtfYqJKuP2C0st5hfSDgawbwoElLcMgGWZBXXEFBvhLSGCQOhC5Ux7s7wtXyqRpMyThZqLp+ZIx0fLsY6NlUJ+urAdLRG4XCGIUQVRRuXDpSq782SceGiEPDOvvlvmzejmvkP4wkHWjg6MhM1PSDRk2XR0XFa6ig4MCJSggwO5X05t37JXEpGTp3bX9Xb9766598v6k7+ToiT/zYAMDA+ThEYNk/Oh7mQdLRG7NYMiQH2fPV3eTUs3yZytXKCcvPzVWurZvddfj4NpN2wV38rt2aOWQafSG8wrR3TCQdTEcaNIMBkkzZDj17+JQgwNPSFCQXQ48aKF48Z/vq98xZvnDIwdbfd/Z87fzYFdtyJ0H26dbB3nhiYekQrnSxZ4WIiJn5s9+PnWmqrJinl/apkVjefWpcfnmzy5Yvlb++tYn6vd/vPCEjBzS127T5C3nFSJbMJB1ocysLElOS5OcnD8Pfs7m6+sj4SEhEuDvX+TvSLh5S4Y+/KwpPQCjay2Y8YXERkflyoOd/OOvMn32QsnKMsuDrV1dXnt6vLRoUr+Yc0JE5DoHj97On925N3f+7L339JSnxo3KdTxMTU2TgQ8+rQZhgJCQYPn9u0/s0qHVW84rRLZiIOsCuFpOTTdIeoZzr5YLEhwYKKHBRbuKfvmN/8rS1RtzPYfi4f986UmVB/vbwhXyxbSZefJgn3t8jMqDxcGeiMgbju0YffDDL/+XK38WF/cTHv4zf3bStJnyzQ+/5vosBneZ+vGbRT4eett5hchWDGSdzB2ulu15FY2i4S+9/oGV7/KV/3vxCZn5+yI5dvJM7jzYkYPl0QeGSSjrwRKRl+bPogrL5B9+zZM/+8ioIfLuZ9PEYCXg/Pvzj8v9Q/uJ3s8rRIXBQNaJkHSfnJom7i48NEQl7Rc2peBu+nbvKM8//iDzYIlIF65dvyGfT5spfyxamSt/Nj9FSTHwtvMKUWExkHUS3O5JMbsyd3dhIcHqtlBB0BKLFtm7qV+nhrz29Dhp3ph5sESkP4eOnlSjFe7Yc+Cu7y1MioE3nleICovJiU7gaQcbwPQWlGuFnFhbglgUCP/+s/8wiCUi3UKn1qkfvSmVK5S963u379kvs+Yu1eV5hagoGMg6GG77eNrBRoPpxvRbSyl4+5PJNn0H3vvL3CUOmDoiIs+xaMV6OXvhkk3v/eibH+R8/GVdnVeIioqBrAOpBHwPyF0qCKYf82HuPx9PtjkvFr7+fpYKaImI9Agdvj6Z/KPN709LS5d/vjtJcqwMYuCt5xWiomIg6yBI7EcvUm+A+dA6KqDoty0pBbk+n5IqX0z7yUFTR0Tk3v738xxTzVhbIcVg2sw/dHFeISoOBrIOgnp+7lgKpSgwH5if5JQUef39L4r0HWs2bbf7dBEReYI1G4t2/PtsynQ5dyHeq88rRMXFQNYBcMvE2xLaMT/XEm5avdVVEBTCLhEZLg+NGOSwaSMicmejhvWX2JioIg0McOLUOa8+rzDFgIqL5bfsDLdKbiYne81Vs2Vh6x9+niubd+yVejWrSelSsWrEGtQHDMP/2o/Z45DgII7cRUR05/yA/FfkiCLlKiU1VZJT0kz/464X/kf9WQx5W6dmNXn95Se9/rwSFR7O0b+oyBjI2llqerqkGbzrqtlcSBCGHAx29WQQEekGzytE+WNqgb2vtr04iIV0QwYT9ImInITnFaKCMZC1Iz3UxjPqZD6JiNyBHo63PK9QcTCQtXNrpR7oZT6JiFxNL8dbvcwn2R8DWTvJysqW7EL26PdUmE/MLxEROQ7PK0R3x0DWTtIy9FUPT2/zS0TkbHo7zuptfsk+GMjaKRk/I1Nfw+1hfjkqCxGRY/C8QmQbBrJ24MiUglOnTkqVCuWlZFQJ9TN86JA8ASQe3ztksOk9lcuXkxMnToij6SWVgojI2Rx5fM3MzJSuHTuazhllS8bK/ri4PO+Lj4+X6pUrm97XsF5duXXzpjgSzytUWAxk7SAr23H5otWqVZc333rL9HjN6tUybcqUXO+ZOnmyrF2zxvT4jX//W2rUqCGePN9ERHrmyONrQECAfP7ll+p/9beysuT5Z5+RbIu/+epLL0li4i3T448++VRKREU5bLrUtPC8QoXEQNYDdryxj4yTnr17mx7/643X5fjx4+p3/I/Hmh49e8q48Y+KM/CAQ0TkmcfXRo0by/MvvmR6vGf3bvnqyy9Mj+fO+UMWL1poejzqgdHSu08fcTSeV6iwOLKXHdxMSnb47ZBLly5J5/btJCEhQT1u0bKlzFu4SAb27ye7du5Uz0VHR8u6TZulXLly4gx+vr4SFRHulL9FRKQnzjivIMWgV7dusn//7bSCkJAQWbdpk0RFRUmHNm3kypUr6vly5cvLhk2bHd4aCzyvUGExkC0m5KcmJCaJM+AKefzYsabHDRo0lAMH9pseT/3uOxkydJg4U0xkBMfIJiLy0PNK3L590rtHdxXUQsdOnaRChYryy88/md7z06xfpZfZXUFH43mFCoOpBcWU48QOT4OHDJXhI0aYHpsHsffed5/Tg1hnzz8RkR4487iKFIMXXnrZ9HjD+vW5gtgHRo9xahALPK9QYbBF1g75PLeSU8RZ0GO0fds2cvnSJdNzpUuXls3btjvlto+lEuFh4u/n5/S/S0TkrZx9XkFrLFpl0TprrnyFCiqlILJECXEmnleoMNgiW0zOrqV68eJFuXnjRq7nbty4IWfPnhVXYC1ZIiLPPq6iesGkL78yVTHQfPzpZ04PYoHnFSoMBrIeBFfNE594XAwGQ97nJzyR53kiIiJbNGjYUFq2am16XKlSZVUFh8jdMZD1IO+987bEmRWtHv/YY6bfDx08KG+/9W8XTRkREXk6Hx8fV08CUaExkPWQHX/b1q3y+aefmh6PHjNG3vvgvzL6wQdNz331xReyedMmcSYe+IiIeFzleYVchYFsMTkjjE1JSZGnnpxgGnWlcuXK8tY776rf//POu1K1alVTT8+nn5wgycnJ4iy8fici4nGV5xVyFQayxV2Avo5fhP/8x9/l1MmTpr836auvJSIiQj0ODw+XL77+RvzuVA44c+aM/ONvfxVvmn8iIj3R+3FV7/NPhcOtxQ631jESiaOsXLFCvv/uO9PjCRMnSvsOHXK9p03btvLs88+bHk//4QdZtnSpOBrmm6kFRESedV5xZzyvUGGxjqwdJKeliSHj9qgoehIUGCDhISGungwiIq/D8wqRbfR5yWdneh0QQK/zTUTkaHo9vup1vqnoGMjagV53PL3ONxGRo+n1+KrX+aai8y/GZ+kOPecyERGRdTdvJcrBoyelRGSElC4ZLTFRJUwdc3l85XmF7IOBrJ0S8wMD/CUjM0v0AvPLjl5ERNZlZWXLqAmvyvmLl3P1xo+NLiGlSsZI6dgYKRkbLaVLxkipO//jcZlSsRIbHcXzCpGNGMjaSUhgkK4CWcwvERHl7/LV67keo9b31es31M9BOZHv56pWKi/zp3/B8wqRDXhv2E78/f10c6sd84n5JSIi63CM7NahdZEWz5nz8WI0GnleIbKBPiIvJwkOChQ90Mt8EhEVx5B+3Yv0uXvv6WVK3dLL8VYv80n2x0DWjoICAsTb+ehkPomIiqt9q6ZSulRsoT4zfvQwef3lJ02P9XC85XmFioOBrB3hCjrEy68qcdXMTl5ERHd38sx51amrMC2xzz02JtdzPK8QFYydvewsJChIDJmZkpNjFG/j64tAnZ28iIjyk5iUIotXrpc5i1fK/sPHbV5QPTu3lf978QmrDQU8rxDlj0PUOkBmVpYkpqSKt4kMC5UAf177EBFZViPYsnOfzFm8Slau3yIZFkOWIzY1FtC20aZ5I/nyvf+TwMD80wh4XiGyjlGJAyDYCw4MlPSMDPEWmB8GsUREfzp38ZLMXbxK5i1dI/GXr+ZZNPVqVZfB6PBlNMq7n0+zuujq16khn/7nrwUGscDzCpF1bJF1EJROuZmc7BUpBkgpiAoPZ24sEelealq6rFi7WbW+bt+zP8/yiCoRIQN6dpYh/XpI3VrV1HMpqWnSbdg4SUtLz1Mv9vtJb6sRv2zB8wpRXmyRdRDkOYWHhHhFigHmgx28iEivEEDuPXBE5ixaKUtWb1SBqTmM2NWhdTNVbqtr+1Z5WlfDQkOkd5f2MnfJKtNzqGYw+cM3bA5igecVorzYIutg6PiVbHHQ8yThoSG6KP9CRGTp6vUEmb90rfyxeKWcPnshz+toUUXqwMDeXdXQsgXZufegjH327+r3EpHh8v3nb0uNqpWKtNB5XiH6EwNZJ0CubIrFLSW3hmwIH5GwkGCVG0tEpBeZmZmydvMO+WPRStm4bbdkZ+fkej00JFj6du+oWl+bNqxr890qtOp+Pm2m7Dt4VF6c8JDUr11DX+eVO3heIXtjIOsknnLQQe9bRLHBgQESERbq6skhInKKIydOq7zXhcvWyo1biXleb9Gkvgzt30N6dWmvgll34CnnFQ2DWHIEBrJO5O63g7Kzs1Vj7L/e/1KCAgPkg9dfEj8/P1dPFhGRQ9xKTJJFK9er1tdDR0/meR3pAkgdGNK3m1SqUM4t14K7n1c0TFMjR2Eg62SoBZicluaW1QyysrLllTc+kN37DqnHD943UF59epyrJ4uIyK4X7OY1XzMzs3K9jo5a3Tu2Ua2vqO/qCRfz7nxeQdUbdBhm+UZyFAayLoBcqdR0g1vVmUUubGhwkGzYukue+dvbpryw154ZJ2OGD3T15BERFcvZ8/EyZ8kqmbdktVy+et1qPdeh/XpIv56dpEREuMctbXc+r7DqDTkSA1mdX0Vbu1qePX+ZvPnfr9TvOAB9/K9XpUfnti6bRiKiokhNTZNlqubrSlU1wFJ0iUgZ0LuL6rhVp0ZVr1jI7npeIXIUBrJucBWdZjBIuiFD5ac6C/rZBgcFqjG8rV0tfzZ1hkz5cbb6PSgwUKZ+/KbqoUtE5O7H1N1xh1TqwNLVG9UABub8/HylY+vmMqR/d+nSrqUEeGF5QXc9rxA5AgNZNzrwIGkfB55sVTnAMfx8fdWBBrVhCzrQYHr+/vZnMn/ZGtNoNdO/fFeqVCzvsGkjIioqpAvMX7pGtb6eOR+f5/VqlSuo0bYG9ukipWJjdLGg3e28QuQIDGTdEDpdpWUYJMOiE0JxBAb4S0hgkPj7+xWqnuKEV/4l23bfHoaxcoVy8uOX7xRqJBoiIkfJyMiUNZu2q6oDm7bvuVM+MPeIWqj5io5bjevX1nWQ5S7nFSJ7YyDr5lfTuIrOys6+/ZOVbdNVNa6OcWDx97v9g8dFPYAnJqXIw8/8VY6fOqce42Qw7ZN/SXBQUJG+j4iouA4fO6WC14Ur1sqtxOQ8r7du1lC1vvbs0k5CgnmscrfzCpE9MZD1wIMQWh2Md37X4IDic2fMb3sfXOIvX5UxE/8iV64lqMc9OrWRD998xSPK0hCRd7h5K1EWLl+ncl8PHz+V5/VyZUrJ4L7dZBBqvpYv65Jp9FSuOK8Q2QsDWbK5BeThZ/5m6jgxZvg98toz47n0iMihNV+RMoDWV6QQWNZ8RUdUVFRB1QHUfEXARUT6wkCWbIZxx5/6y1umGrMYLAGDJhAR2dPpcxdk7uLVMm/patOdIHON6tVSwWvf7p0kMiKMC59IxxjIUqH8vnCFvP7+F7c3Hh8flWLQq0s7LkUiKpYU1HxdvVH+WLxSdscdzvN6TFSkDOzTTQWwNatV5tImIoWBLBXapGkz5ZsffjXd2kPnryYN6nBJElGhIB9z576DMmfRKlm2dpOkWan52qltC1V1AP+zwD4RWeKwG1RoT40bpTqAzVu6RoICA0wnJHYGICJbXLpyTaUNoOPWuQuX8rxeo2olFbwO6NVFSsZEcaESUb4YyFKhIWB945WJkpmdLRPHjpSK5csyiCWiAhkMGbJ64zYVvKIDl3nveAgPC5V+PW7XfG1YtxaPKURkE6YWUJFhLG+jMafAMlyoU4iagzlGo/iyfAuRriBYPXj0pBpta9GK9ZKYlLfma5sWjVXea49ObVnzlYgKjYEsOcyhY8dl+549MvreYap4NuoUsjwOkfdLuHnLVPP16InTeV6vULa0DO7XXQb16SYVypV2yTQSkXdgIEsOsXXXblm6dq1ULFdOqlepLF3b3a5swGCWyLPgbgpK4d8tBx4jRG3avvtOzdcdkpWVu+ZrcFCg9OzcTob07y6tmjbkRS0R2QUDWbL7rcQlq9fI+fh4qVerlpQpVVL27D8ghowMGXPvMC5tIg9hazoQ9vmbt5Lk3nHPy9XrN/K8jmGtMVxs3+4dJCKcNV+JyL7Y2YvsOgrPH0uWyNXrCdK/ezfVGov82Urly8uUGTPl4uXLUr5MGS5xIg+gBbG4KD1x5ozaj6tWqpQnuEVLbXRUpJQtXdIUyMbGRMmg3l1V+gAqEBAROQoDWbKLdINBfpg9W5JTUqV65crqpKflw6akpqrfoyIjubSJ3LDlVdD6ajG8K1ID/li8RM7Fx0vNqlVl74EV0qhuXWnZpLFEhIfnKrmH9w4b0FNKl4xRra8d2jRjzVcicgqmFlCx3UxMlBm//yE1qlaRiLDbJ7iOrVuZXv/2518kMDBAxgwbJj/O/k3aNG8mtatXZ+1ZIje29+BB1WHz/sGD1OP4y1dk1caNUq5MaeneoUOe/Ze1pInIFdgiS8WSmZWlgtP6tWtLh1YtJTE5Wab99LNkZWdJaEiIbN65U2Kjo6Vv167y05y56hZl7RrVVfpBqdgYdv4iciIEm0aL1ldciK7etEmSk1OkWaOGUq9mTZUSdOXaNYkMD1fv2bp7t+yKi1PDyCL3Xax0/uKAKETkCmyRpWJLTEoSf39/FbjCsVOnVAevNINBqlWqpFpqf54zV92OrFa5kkRGRMiilatk4tiHpXRsLNcAkZNlZmbK5WvX1EXm9j17JSklWUpERMi2PXulV6dO0qheXVmxfr3sP3JU0tLTpUzJktK8USNp2qC++jw6b2J4aiIiV2OLLBUbAlPzfLta1aqpnDq00KBV57tfZqlgtn/37hISHKyeR11ZX58/W4U4YAKRfaHUHfY185bSawkJcvrceVm4cqW6IxIQECClYmJlSN8+6nXcUcGFaM1qVaV82bJy6uw56dGxg8qN1QLgXfv3q46cFcqW5SojIpdjIEt2pfVoxu1LBLFf/fCjtG3RXHp07Gjq/IHWHNyePHfxohgyDOqEiM+xxixR4Z06d05CgoKlbOlStztu3dkPtfQBbXS93fv3y4r1G6Rh3Try8oQJqoPW1J9+Uh0ztfzWapUqy459e+XshYtSq2pVFfjiMxkZGXLuYrwKcitXKK9y3ImI3AEDWXIInETLlColD903XKUXaK1CyMeLO3RY4g4fFoPBIFElSqiOYKOHDlWfYYcRItslJSfLkRMnpEaVKlJWSuUqjbVzX5zaz2pWrSLtWraUkjExKgUId0LCQm+nAbVr0Vy9D523AgMCpErFCnLg6BE5e+GC1KpWVTq3aSNlS5VSLbO4m/LIyBHqe4iI3AUDWXIIrXUVpbjM4aSLVp2WjRtL62ZNJTsnR777+RdZs2mzdG3fjh1GiAoBeefoSGku4eZNVTYL+2CT+vWlauVKkpOdrQJQXFTeuHXT9N4WjRvL8nXrVY3nqhUrqjz3ktExqpX30tWrqu4zWl/ZAktE7ip34UAie21YFjUp4VZSkqzauEl1JEEQC8iVbdKgvqSmp3HZE90FUgdw8ZfruZwc2XfokOw/ckQ9Rg4s7mw8NvoBtZ+hQyVyYdGiipbbG7duqTsjgA5bGEJ65959pu+rV6umqkCCOypERO6OgSw5DU645cuUlqYNGpiew2AJ+w8fFj9fP64JogL2HUDqAC7+ICMzU/2flm6QU2fPyo47wWjJmGg1GteGbdvl98WLVWfLjdt3qAC2aqXbra4Hjx4zfXeLRo3kzPnzKncdEMCiBVb7O0RE7oxHKnIadDhJuHlLTp45ox5fv3FDdsXtV623aKW11vqkdV5BCxOR3u9woBQWgtPPv/1OFqxYoe5yIN+1To2akpaepjpYVq5QQYbfM0DSDekSHBgkzRo2kEPHjsmytetUKgLSBdDxS4Ma0C8+8TjLaRGRR2IdWXIKrbzW4eMnVOmfcqVLqxaiiLAwVZ8StzP3Hjwk128kqA4paLWNLlHC1PkLrU/ojELkzbC9Y1+xbA1FKsCcJUtVabuU1BRpWLeuzF26VMqVLiM9O3dSLbbzli1T+xUqhFhCECtilN5dusiZ8xfURSRSetjqSkSejoEsOY0WlKJjyfWEG6paAW5jIk9v1vwF4uvrc6dHdagcOn5cXp7whCrZtXLDRtVB7LEHRqkRh4i8kWXFDvNBBxDITv/td/X6+FH3S3BQkJw4fUY27tihUgMa1KmtUgnQ8orc2NS0NNVZ6+r167L/8BH1eFj/fqz9SkReh1ULyOlwaxM/gJZWtNDi9mj7li3VCEIIVkPXhsjpc+fk2MlTKqjt2q4t1xR5FfOar4AgFakCm7bvkCMnT6qyV+iI1bppUzXoSP3atWTfocMqiAWMkrczbp/Kh61To7rKfz149KgcPXlS5bhi/0HZrIZ16kjLpk3Y+kpEXok5suQ01sZiR87f+YvxaoAEBLc+d26ponVq3rLlcuLMGRncp7c0bdiQrbHkFXINWmC2TyC/de7SZar1dFDvXqrCAMrSHT99Wr2vbk3kwaar96nP+/pKpfIVVDmtK9evS6nYWJUDi5Za6NSmjWq9bdO8GYNYIvJabJEllxd0DwwMlAa1a6vHOGGjlfbAkSMSFRkpfbt1zXU7lEPZkqfTglfkqm7fu0eNyjWgZw81OAg6ZmnDwaIG87Y9e+TAkaMq97V0yZJSoUwZ2bFvnxruGTD0894DB9TFIC4EB/bqKeFhYeo1bSQ9IiJvxiMduRTGbMfwl+iM0qxhQ9Xi9NvCheqkjvHfY6KiVH5sdnaO1K1ZI1cLFpHbt7yik6OVMlaLV61W2zXSBSqWK2/qzIggNv7yFVm6do1cvnpNSkRESPyVK6rFFYMZoJPX8nXrTIEsasT27d5N7UegBbFERHrBQJZcPvrXI/ePlN1x+1UvbHQEQ8A6oEdPlTd74OhR2bP/gOr0gjHhtaE1idyduuiycuGFmq/Yzh8YOiTPcK9IK1i7ZbNUrVhJHhw+XIw5OfLhN5Pl3MWLqqwW8mARrF6+etU0YAFG5CIi0isGsuQyCGIRzIaHhqoSXFt375Y2zZpJtw7tVesUCrxjSNuY6ChVjotBLLkr85QXrfrAtYQEWb9tm6rQgVqt6JyFFIFz8fGSdWfI2MysrFwpAEkpKXLs1Gnp07WrymvdffCQGpHrytVrkpycou5QPDX2YRfOKRGRe2EgSy4PZnHiR5CKskE4qSOIxZCb6H1dplRJNR48aspq0g0G1XOb+bLkalrQap7yogWxGLCgbOnS0rF1K5VGsCsuTgWhJaOjJT09XbLuBLEIajFYCL4LQStaWH9buEgSbt5UKQPIey1ftiwHLCAisoKBLLlNNQPkA2oQyIaHhkm7Fi1UXVnN/OUr1K1ZBAQ+d1p0t+/dq/JrOWACOQO2OWyz2g8cOXFSdVxEGSxUDjh07LhKAejbtat6vUrFivLeF1/KngMHpGyp0qqG8s64OHUHAkFsYlKSnL1wURrWraPqvWLgkNiYaKYNEBHdBQNZcjso/n4rMUl6deqUK4j9df4COXXunAzt1/fPUlwIJsRH4g4dUi23RI6A1lKjVjLrTuctrYMWRtw6c+G8RIaHy859+1TON7ZTjFaHCzLUhU1MTpbG9epJrerVVSpN2+bN1UAfyJUF1EuuVb2ayg/HNt+icSOuSCIiG3BkL3JLfyxeooawxW1VpBIsX7detVqNHjZU1cvEiEWnz52XVk2bqPdnZ2ezziw5Rdyhwyr3FRdTDWrXUikundu0URU3/jdrlsrn9vfzVwN9IC8WlQiQAw5IOdA6eGHkLaTPIMUAn0H+KxERFQ4DWXIr5nmvW3ftlqsJ11WnL1QsGDl4kGrNunDpkqqvmXDjpgp0UV8TkFPIYIAcBZ0RcTGVnJKqBvDAwATrtm5Vras9O3VU79lz4KCs2rBBHh39gHwyZao89sADUq5MafUa6sEitaBn505qBDsiIio+phaQW0EQq5XlwohE3//6q9SsVlXuHzxY5RJilKNdcftVj+4OrVpJbEyM7N6/X67fuClbd+2SF594XHWYIbI3DDiA+q9d27VTKQC4uEpOSZGbt26Z3oPn0RKLgLdft26ybN06lX6Aiy9sl+1btpBSFiW3iIio6BjIkttWMkBHmlFDhpg6ce0/ckTd1kVA0LRhA9URBuWKMKxncHCQPD5mDINYcpgK5cqp0eZio6NNlTbKli4l+48clfPx8arCAKpp1KtZU5auWauGh61Xu5YadatHp45q8AIiIrKvvEPOELkBBLEIFrQgFp1otu/ZK1ElIlVLLYLY5NRUVaYIvcSNOUbVaxzwOaListyOcLcAF1Cnz59T+bBQvkxZdfF0/NRp0/taNmminkPeNlJh0ErLIJaIyDEYyJLb0kobAVq7KpQtI62aNFVF5a8l3JBpP/0kmZmZ8tITj8uw/v1lxYb1KsAw/xxRYYPX3XGH5P1J394eYtZC84YNVZmsW0lJ6nHJ2BiJCA+TU+fOiiEjQz1XuUJ5NWqXqbIGERE5DFMLyCOCCwzH2b1DB/H391eduibPmCE1qlSRkYMGqvegfmep2BgGD1Qkl69el/lL18jcJavk9LmL6rm2LRpL+1bNxN//z4A0MiJC3SU4dfacynXF9ti0fgNp1aQJBywgInIBBrLk9rQWVrRwZefkyLbde3IFsVrnMFYsoMLIyMiUNZu2yx+LVsqm7XvUdmRu4Yp10rldyzyfa9+ypcRfuazqygLyZImIyDUYyJJHBbR+Pj6qYDxu5Wq0AvVEtjh87JQKXheuWCu3Em/nVZtr3ayhDOnXQ3p0bmv18/Vr15IGdWpzYRMRuQHWkSWPoVUygB9m/6YKzTdr2MDVk0Ue4OatRFm4Yr3MWbRSDh8/lef1cmVKyeC+3WRQ325SqXxZl0wjEREVHgNZ8ihILUANWcAwtiUiI/INdi1pKQikD6gagJSBOYtXyeqN2yQzMyvX60GBgarVdUi/7tKmeSNuG0REHoiBLHmc/ILVnByjbNmxR0qXipWa1SqbnsdITP98b5Js2LZb/v78YzK4b3cnTzE50+lzF2Tu4tUyb+lquXItIc/rDevWlKH9e0jf7p0kMiKMK4eIyIMxkCWv8fOcxfLuZ1NV9YIZX70npUvGqN7oT/3lLTly/Hadz8oVysnCmV+6elLJzlJS02TZ6o3yx+KVsjvucJ7XY6Ii5Z7eXVXra63qVbj8iYi8BDt7kdf0QEcHnuzsHLl05Zo89dpb8o8Xn5CXXv9ABbOasxfiJf7yVZUTSZ7fMr9z30GZs2iVLFu7SdLSbg9SoPHz85VObVuo1lf8H+DPwx0Rkbdhiyx5jWvXb8iYiX+RC5eumAIZBLaW/vXa0yq4Ic+ECxWkDSB9ABcmlmpUraRaXu/p1UVKxka7ZBqJiMg5GMiSVzl5+pyMfOIVSU835PueAb06y7v/eMGp00XFYzBkqA5b6LiFDlyWw8eGh4VKvx4dVdmsRvVqcXQ3IiKd4L028hoIbpav21JgEAtbd8YVWN2A3APW0aFjJ1XKyKIV6yUxKW/N1zYtGqvW1x6d2kpIcJBLppOIiFyHgSx5hcysLPnPR9/IbwtX3PW91xJuyMkz59UtaHI/CTdvycLl61Tr69ETtzvpmatQtrQM7tddBvXpJhXKlXbJNBIRkXtgIEteUS/02b+9Ixu27rL5M1t27rNbIIuWQ9Soxc1u81veaPFFmy9q17L1t2BZWaj5ulu1vq7ZtEOysvLWfO3VpZ0M6d9dWjVtyJqvRESkMJAlj7d1V1yhgljYtitORt87oNB/C4EqBmXIys6+/ZOVrR7fDQZx8Pf3E3+/2z94zOBWVMs4Wl4XLFsjV6/fyLPcGtevrfJe+3bvIBHhrPlKRES5MZAlj1e7ehWpWL6MnL942ebPbN+zX7Xk+vn52fR+BKxpGQbJsBgdylYIdrMzcsQgmabnAgP8JSQwSAW4eoIBKpai5uuilbL3wJE8r8fGRMmg3l1V+gDTP4iIqCCsWkBe06t9+brNMmvuEqsF8a356ev3pWG9WgW2vhoyMyXdkGFTq2tRoXU2OChQggICvLaVFqkXO/celDmLV8rytZslzaJDHlqpu7RvqVpfO7RpxpqvRERkEway5HWOnjgjs+YtkQXL1qoRn/LzzKMPyOMP3mc1gE0zGCTNkCHOhBAWAW1IUJDXBLQYfGLuktUqfeBCfN4WcwwljJq+KIkWGx3lkmkkIiLPxUCWvFZqaposXLFetdIePn4qz+vly5SSpbMm56l+kJyWJjk5ueuUOpOvr4+Eh4R4bKtkusEgq9aj5utK1anOsuYrcl379+ikAtj6dWp4TdBORETOx0CWvB4CqbhDx+SXuUtk8cr1knknzxWtn9uX/WJ6T2q6QdIznNsKW5DgwEAJDfaM1lksvwOHj6uW10Ur10lScmqu1zEPbVs0VsFr945tJCgo0GXTSkRE3oOBLOnKrcQkmTTtJ1m7eYc8OHygPDhioFu0wnpq6+z1GzdVCgcC2OOnzuZ5HZ3wBvftLoP7dpNyZUq5ZBqJiMh7MZAlXUNnruQC8mjdRXhoiOoM5g4Q+KPc2ZxFq2Td5h2qDJk5jLB1u+ZrD2nRuD5rvhIRkcMwkCXdQhpBSlq6eIqwkGCVbuAqJ06fU3mv85etlesJN/O83rRhXTVcbJ9uHSQ8LNQl00hERPrinvcriRzM04JY0KbXmcFsUnKKLFm1UQWw+w4ezfN6qdhoGdinq0ofqF6lotOmi4iICBjIki7TCTwtiNVgutFxypFpBqj5igEjMGDBynVbVB1dc/7+/tKtQyvV+tq+VTPdDehARETug6kFpCvI70xMyd2j3hNFhoXavQPYhfgrMm/papmLmq+XruR5vXaNqjK0f3cZ0LOLREdF2vVvExERFQUDWdINlIi6mZzsltUJilLNICo8vNiluTDC1sr1W2TOopWydVdcntcjI8JlQM/OMqR/d6lXq7pHlAIjIiL9YGoB6QbqxHpDEAuYD8wPOoAVta4u8l4Xr9wgyRYt1AhW27dqqoaLRQoBa74SEZG7YiBLukkpcKfBDuwB8xMY4G9zisG1BNR8XaNqvqICgaVKFcrK0H49VOetsqVLOmCKiYiI7IupBeT1vCmloLApBgjg123eqVpf12/ZKdnZObleDwkJlt5d2qsRt5o3rsfUASIi8ihskSWvl2bwnpQCS5gvzF9ocO4Ug2Mnz6iW1wXL10rCjVt5Pte8UT1VdaB3tw4SFhrixCkmIiKyHway5PWtsWkW5aO8DcpjhQQFSVJyqixeuV61vu4/fDzP+0qXjJFBfbrJ4H7dpGqlCi6ZViIiIntiagF5NU8c+KAoFi5bK59Nni4GizzggADUfG19p+ZrU/HzY81XIiLyHmyRJa9mWczfG2Xn5EitGlVyBbEolTW4H2q+dpKoEqz5SkRE3omBLHmtrKxsFeR5Oz9fXzVYQcumDaROjaqqbFbdWtVcPVlEREQOx9QC8lpJqamSkZkleskFRhmuEuFhrp4UIiIip/F13p8icm5gp5cgFlB+Kys7W803ERGRXjCQJa/kyJSCDevXS8moErl+Rt9/v9X3rlq5Is97n37ySYdNmx5SKYiIiDQMZMkroXXSmZYvWyqnT5/K8/zkr7/26vkmIiJyJQay5JWcHdDl5OTI1MlTcj13/PhxWblihVOng4EsERHpCQNZ8tqKBc7i63t7N5o5fbqkpKSYnp/6zTemnFVn1W915nwTERG5GgNZ8joIHp2ZK9q3X3/1f2LiLfnlp5nq96TERPn555/U740aN5by5Z0zkhbmmx2+iIhILxjIktfBbX5nGj7iPomNjVW/T51yO71g5owZkpyUpH5//IkJXj3/RERErsJAlryOswtQBQUFy0Njx6rfjx45oioVTJ0yWT0uWbKkDBs+3KnTwwJcRESkFwxkyeu44tb6uPGPir//7YHynnvmGTl18qT6HQFuUFCQU6eFqQVERKQXDGSJ7KBc+fJyz6BB6vf4ixfV/wEBAfLI+Ee5fImIiByEgSyRnTwxIfdABwhsy5Urx+VLRETkIAxkySuHa3WFVq1bS7PmzU2Pnd3Jy9XzT0RE5Gy3k/qIvIgrw7gvv/5Gjh09Kv4BASqwdQWGsUREpBcMZMnraAMUuEKt2rXVj17nn4iIyJl4xiOvg1vrfjoN5jDfTC0gIiK98DGyVg95oeS0NDFkZIreBAUGSHhIiKsng4iIyCn02WxFXs/fz0/0SK/zTURE+sRAlrySXgM6vc43ERHpEwNZ8kp6zpElIiLSC571yCuhw1NggL6KcmB+2dGLiIj0hIEsea2QwCDRE73NLxEREQNZ8lr+/n66udWO+cT8EhER6Yk+zvKkW8FBgaIHeplPIiIicwxkyasFBQSIt/PRyXwSERFZYiBLXg2dn0K8vLUSrbHs5EVERHrEQJa8XkhQkPj6ot3S+2C+MH9ERER6xECWvB5aK7112FbMF1tjiYhIrxjIki4E+PtLcKB3pRhgfjBfREREesVAlnQjNNh7UgwwH5gfIiIiPWMgS7rhTSkGTCkgIiJiIEs6g1vx4aGeHcxi+plSQERExECWdAg1V8NCgsUTYbpZM5aIiOg2phaQLqGjlKcFs5heb+uwRkREVBw+RqPRWKxvIPJghsxMSU5NE09IJ2BLLBERUW4MZEn3MrOyJDktTXJyjG5ZnQAdu5gTS0RElBcDWSIRwY2J1HSDpGdkuM3yQBoBSmxxwAMiIiLrGMgSuVnrLFthiYiIbMNAlshK62yawSDphgxxZjiLoRqCgwIlJIitsERERLZgIEtUQECLzmAIaLNzchy2nPx8fVUAi85cTCMgIiKyHQNZIhtkZWVLWoZBMjKz7La8AgP8JSQwSPz9/bgOiIiIioCBLFEhW2nROpuVnX37JyvbptZatLoiYPX3u/2Dx2x9JSIiKh4GskR2CG5zcnJUPq15WWYEqsh79WXQSkRE5BAMZImIiIjII3GIWiIiIiLySAxkiYiIiMgjMZAlIiIiIo/EQJaIiIiIPBIDWSIiIiLySAxkiYiIiMgjMZAlIiIiIo/EQJaIiIiIPBIDWSIiIiLySAxkiYiIiMgjMZAlIiIiIo/EQJaIiIiIPBIDWSIiIiLySAxkiYiIiMgjMZAlIiIiIo/EQJaIiIiIPBIDWSIiIiLySAxkiYiIiMgjMZAlIiIiIo/EQJaIiIiIPBIDWSIiIiLySAxkiYiIiMgjMZAlIiIiIo/EQJaIiIiIPBIDWSIiIiLySAxkiYiIiMgjMZAlIiIiIo/EQJaIiIiIPBIDWSIiIiLySAxkiYiIiMgjMZAlIiIiIo/EQJaIiIiIPBIDWSIiIiLySAxkiYiIiEg80f8DtqscEdChGL4AAAAASUVORK5CYII=", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Aretes du DAG : [('Z', 'X'), ('Z', 'Y'), ('X', 'M'), ('M', 'Y')]\n" + ] + } + ], + "source": [ + "import networkx as nx\n", + "\n", + "G = nx.DiGraph()\n", + "G.add_nodes_from(['Z', 'X', 'M', 'Y'])\n", + "G.add_edges_from([('Z', 'X'), ('Z', 'Y'), ('X', 'M'), ('M', 'Y')])\n", + "\n", + "fig, ax = plt.subplots(figsize=(7, 4))\n", + "pos = {'Z': (1, 2), 'X': (0, 1), 'M': (0.5, 0), 'Y': (2, 1)}\n", + "nx.draw(G, pos, ax=ax, with_labels=True, node_color='#ecf0f1', edge_color='#2c3e50',\n", + " node_size=2000, font_size=14, font_weight='bold', arrows=True,\n", + " arrowsize=20, width=2)\n", + "edge_labels = {('Z', 'X'): 'selection', ('Z', 'Y'): 'direct',\n", + " ('X', 'M'): 'produit', ('M', 'Y'): 'mecanisme'}\n", + "nx.draw_networkx_edge_labels(G, pos, edge_labels=edge_labels, ax=ax, font_size=9, font_color='#7f8c8d')\n", + "ax.set_title('DAG 4 noeuds : Z (aptitude) -> X (mentorat) -> M (effort) -> Y (score)')\n", + "ax.set_xlim(-0.5, 2.8)\n", + "ax.set_ylim(-0.5, 2.5)\n", + "plt.tight_layout()\n", + "plt.show()\n", + "\n", + "print('Aretes du DAG :', list(G.edges()))" + ] + }, + { + "cell_type": "markdown", + "id": "e01cd205", + "metadata": {}, + "source": [ + "## 3. P(Y | X) — l'observation naive\n", + "\n", + "**Definition** : P(Y = 1 | X = x) est la **probabilite conditionnelle** que Y = 1 parmi les unites observees avec X = x. C'est ce que mesure une statistique descriptive classique (un tableau de contingence).\n", + "\n", + "**Piege** : quand Z confond X et Y (Z -> X et Z -> Y), P(Y | X) **n'egale pas** P(Y | do(X)). Snow voyait P(cholera | eau_pompe) eleve a Soho, mais cela reflait en partie P(Soho) eleve (les gens de Soho etaient aussi plus pauvres, plus ages, etc.).\n", + "\n", + "**Cellule-type 1** : on simule ci-dessous un DGP semi-synthetique et on mesure P(Y=1 | X=1) vs P(Y=1 | X=0). On verra que l'**association observee** est biaisee par Z." + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "id": "f4473759", + "metadata": { + "execution": { + "iopub.execute_input": "2026-10-05T19:53:14.938364Z", + "iopub.status.busy": "2026-10-05T19:53:14.937851Z", + "iopub.status.idle": "2026-10-05T19:53:14.967162Z", + "shell.execute_reply": "2026-10-05T19:53:14.965519Z" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "P(Y = 1 | X = 1) = 0.7078 (sur 44,050 unites traitees)\n", + "P(Y = 1 | X = 0) = 0.3221 (sur 55,950 unites non traitees)\n", + "Association observee (difference) : +0.3857\n", + "\n", + "Note : Z est distribue differemment dans les deux groupes.\n", + " E[Z | X = 1] = +0.5239 (les traites sont plus aptes)\n", + " E[Z | X = 0] = -0.4201 (les non-traites sont moins aptes)\n", + " => P(Y | X) reflete a la fois l'effet causal et le biais de selection par Z.\n" + ] + } + ], + "source": [ + "rng2 = np.random.default_rng(42)\n", + "N = 100_000\n", + "\n", + "# Z : aptitude (continue, standardisee)\n", + "Z = rng2.standard_normal(N)\n", + "\n", + "# X : selection au mentorat (depend de Z, l'aptitude pousse vers X=1)\n", + "logit_X = 1.2 * Z - 0.3\n", + "p_X = 1 / (1 + np.exp(-logit_X))\n", + "X = (rng2.random(N) < p_X).astype(int)\n", + "\n", + "# M : effort (mediateur, depend de X ; effet homogene)\n", + "M = 0.7 * X + 0.2 * Z + rng2.standard_normal(N) * 0.5\n", + "\n", + "# Y : score final (depend de M, Z ; X agit via M)\n", + "logit_Y = 1.5 * M + 0.8 * Z - 0.5\n", + "p_Y = 1 / (1 + np.exp(-logit_Y))\n", + "Y = (rng2.random(N) < p_Y).astype(int)\n", + "\n", + "p_Y_X1 = Y[X == 1].mean()\n", + "p_Y_X0 = Y[X == 0].mean()\n", + "print(f'P(Y = 1 | X = 1) = {p_Y_X1:.4f} (sur {X.sum():,} unites traitees)')\n", + "print(f'P(Y = 1 | X = 0) = {p_Y_X0:.4f} (sur {(1 - X).sum():,} unites non traitees)')\n", + "print(f'Association observee (difference) : {p_Y_X1 - p_Y_X0:+.4f}')\n", + "print()\n", + "print('Note : Z est distribue differemment dans les deux groupes.')\n", + "print(f' E[Z | X = 1] = {Z[X == 1].mean():+.4f} (les traites sont plus aptes)')\n", + "print(f' E[Z | X = 0] = {Z[X == 0].mean():+.4f} (les non-traites sont moins aptes)')\n", + "print(' => P(Y | X) reflete a la fois l\\'effet causal et le biais de selection par Z.')" + ] + }, + { + "cell_type": "markdown", + "id": "f4c44c20", + "metadata": {}, + "source": [ + "## 4. P(Y | do(X)) — l'intervention par mutilation du graphe\n", + "\n", + "**Definition** : P(Y = 1 | do(X = x)) est la probabilite que Y = 1 **si on fixe** X a x pour toutes les unites, independamment de Z. Graphiquement, cela correspond a **supprimer l'arete Z -> X** (Pearl §3.3, *graph mutilation*) : on coupe le lien par lequel Z determinait X, et on force X = x.\n", + "\n", + "**Backdoor adjustment** (Pearl §3.3.1) : si Z est un **ensemble de backdoor valide** (il bloque tous les chemins non-causaux de X vers Y et ne contient pas de descendant de X), alors :\n", + "\n", + " P(Y | do(X = x)) = somme sur z de P(Y | X = x, Z = z) * P(Z = z)\n", + "\n", + "Dans notre DAG, Z est exactement le parent direct de X, donc c'est l'ensemble de backdoor minimal. On verifie cette identite ci-dessous par simulation." + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "id": "070ad1e5", + "metadata": { + "execution": { + "iopub.execute_input": "2026-10-05T19:53:14.969740Z", + "iopub.status.busy": "2026-10-05T19:53:14.968735Z", + "iopub.status.idle": "2026-10-05T19:53:14.997264Z", + "shell.execute_reply": "2026-10-05T19:53:14.996258Z" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "P(Y = 1 | do(X = 1)) = 0.6063 (backdoor adjustment sur Z)\n", + "P(Y = 1 | do(X = 0)) = 0.4039\n", + "Effet causal (difference) : +0.2024\n", + "\n", + "Comparaison :\n", + " P(Y | X = 1) - P(Y | X = 0) = +0.3857 (biaise par Z)\n", + " P(Y | do(X=1)) - P(Y | do(X=0)) = +0.2024 (effet causal ajuste)\n", + "\n", + "L'association naive surestime (ou sous-estime) l'effet causal selon la direction du confounding.\n", + "Le backdoor adjustment par Z neutralise ce biais.\n" + ] + } + ], + "source": [ + "# Backdoor adjustment : P(Y | do(X=1)) = E_Z[P(Y | X=1, Z)]\n", + "# Stratifier par Z, estimer P(Y=1 | X=1, Z=z) dans chaque strate, moyenner\n", + "\n", + "# Discretiser Z en deciles pour estimer E_Z[f(Z)]\n", + "Z_bins = np.quantile(Z, np.linspace(0, 1, 11))\n", + "Z_strate = np.digitize(Z, Z_bins[:-1]) - 1 # 0..9\n", + "\n", + "p_Y_do_X1 = 0.0\n", + "for s in range(10):\n", + " mask = (Z_strate == s) & (X == 1)\n", + " if mask.sum() > 0:\n", + " p_Y_do_X1 += Y[mask].mean() * (Z_strate == s).mean()\n", + "\n", + "p_Y_do_X0 = 0.0\n", + "for s in range(10):\n", + " mask = (Z_strate == s) & (X == 0)\n", + " if mask.sum() > 0:\n", + " p_Y_do_X0 += Y[mask].mean() * (Z_strate == s).mean()\n", + "\n", + "print(f'P(Y = 1 | do(X = 1)) = {p_Y_do_X1:.4f} (backdoor adjustment sur Z)')\n", + "print(f'P(Y = 1 | do(X = 0)) = {p_Y_do_X0:.4f}')\n", + "print(f'Effet causal (difference) : {p_Y_do_X1 - p_Y_do_X0:+.4f}')\n", + "print()\n", + "print('Comparaison :')\n", + "print(f' P(Y | X = 1) - P(Y | X = 0) = {p_Y_X1 - p_Y_X0:+.4f} (biaise par Z)')\n", + "print(f' P(Y | do(X=1)) - P(Y | do(X=0)) = {p_Y_do_X1 - p_Y_do_X0:+.4f} (effet causal ajuste)')\n", + "print()\n", + "print('L\\'association naive surestime (ou sous-estime) l\\'effet causal selon la direction du confounding.')\n", + "print('Le backdoor adjustment par Z neutralise ce biais.')" + ] + }, + { + "cell_type": "markdown", + "id": "a4468b0d", + "metadata": {}, + "source": [ + "## 5. P(Y_x | x_0, y_0) — le contrefactuel individuel (Pearl §9)\n", + "\n", + "Les deux premieres cellules-types operent sur des **distributions** : P(Y | X) mesure une association dans l'echantillon, P(Y | do(X)) mesure un effet causal moyen. Le **troisieme echelon** de l'echelle de Pearl est le **contrefactuel individuel** : que se serait-il passe pour **un individu precis** dont on connait l'histoire, si une condition passee avait ete differente ?\n", + "\n", + "Pearl propose un protocole en trois pas (abduction-action-prediction) :\n", + "\n", + "1. **Abduction** : mettre a jour la distribution des variables non observees (ici Z, l'aptitude) au vu des observations (X = x_0, Y = y_0).\n", + "2. **Action** : modifier le graphe (mutiler X := x).\n", + "3. **Prediction** : calculer la loi de Y sous le nouveau graphe, conditionnellement aux variables exogenes mises a jour.\n", + "\n", + "**Exemple** : Alice (X_0 = 0, pas de mentorat) a un score Y_0 = 0.6. Si on lui avait donne un mentorat (X = 1), quel aurait ete son score ? La cle est l'**inference sur Z** (son aptitude probable, sachant qu'elle a 0.6 sans mentorat)." + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "id": "a9c3e450", + "metadata": { + "execution": { + "iopub.execute_input": "2026-10-05T19:53:14.999299Z", + "iopub.status.busy": "2026-10-05T19:53:14.999299Z", + "iopub.status.idle": "2026-10-05T19:53:15.013015Z", + "shell.execute_reply": "2026-10-05T19:53:15.011816Z" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Pool d'Alice-like : 8061 unites (X=0, M proche du 60e centile)\n", + " E[Z | pool Alice] = -0.3477\n", + " Std[Z | pool Alice] = 0.8387\n", + "\n", + "Score contrefactuel moyen d'Alice sous X = 1 : 0.5355 (std = 0.1954)\n", + "Score observe d'Alice sous X = 0 : 0.6000\n", + "Difference contrefactuelle : -0.0645\n", + "\n", + "Note : la fragilite du chiffre individuel reflete l'incertitude sur Z ;\n", + "CB-03 (Dowhy contrefactuel) explorera cette fragilite en detail.\n" + ] + } + ], + "source": [ + "# Exemple contrefactuel : Alice (X_0 = 0, Y_0 = 0.6) - que serait son Y si X := 1 ?\n", + "# On tire un echantillon d'Alice-like (X=0, Y proche de 0.6), on infere leur Z,\n", + "# puis on evalue ce que leur Y serait sous X=1.\n", + "\n", + "rng3 = np.random.default_rng(2025)\n", + "alice_pool_mask = (X == 0) & (np.abs(M - np.quantile(M[X == 0], 0.6)) < 0.1)\n", + "alice_Z = Z[alice_pool_mask]\n", + "print(f'Pool d\\'Alice-like : {alice_pool_mask.sum()} unites (X=0, M proche du 60e centile)')\n", + "print(f' E[Z | pool Alice] = {alice_Z.mean():+.4f}')\n", + "print(f' Std[Z | pool Alice] = {alice_Z.std():.4f}')\n", + "\n", + "# Contrefactuel : P(Y | do(X=1), Z = z_Alice)\n", + "# Approximation : on evalue le modele de Y sachant X=1, M, Z sur le pool\n", + "alice_logit_Y_do_X1 = 1.5 * (0.7 * 1 + 0.2 * alice_Z) + 0.8 * alice_Z - 0.5\n", + "alice_p_Y_do_X1 = 1 / (1 + np.exp(-alice_logit_Y_do_X1))\n", + "alice_EY = alice_p_Y_do_X1.mean()\n", + "alice_StdY = alice_p_Y_do_X1.std()\n", + "print()\n", + "print(f'Score contrefactuel moyen d\\'Alice sous X = 1 : {alice_EY:.4f} (std = {alice_StdY:.4f})')\n", + "print(f'Score observe d\\'Alice sous X = 0 : 0.6000')\n", + "print(f'Difference contrefactuelle : {alice_EY - 0.6:+.4f}')\n", + "print()\n", + "print('Note : la fragilite du chiffre individuel reflete l\\'incertitude sur Z ;')\n", + "print('CB-03 (Dowhy contrefactuel) explorera cette fragilite en detail.')" + ] + }, + { + "cell_type": "markdown", + "id": "1fc26af2", + "metadata": {}, + "source": [ + "## 6. Exercices\n", + "\n", + "Les exercices suivants etendent la comprehension sans exiger de boite a outils avancee. Ils suivent la convention du depot (stubs a completer, sans erreur volontaire — regle C.1).\n", + "\n", + "1. **Ajouter un second confondeur W (motivation)** : etendre le DAG a 5 noeuds (Z, W, X, M, Y) et verifier que l'ensemble de backdoor devient {Z, W} (les deux confondent X et Y).\n", + "2. **Verifier l'identifiabilite par backdoor** : P(Y | do(X)) doit rester identifiable (c'est le cas) ; montrer qu'on retrouve le meme resultat que la cellule-type 2 si on ajuste sur le binning de Z et W simultanement.\n", + "3. **Contrefactuel inverse** : pour un etudiant **traite** (X_0 = 1) qui a un score Y_0 = 0.3 (en-dessous de la moyenne des traites), estimer son Y s'il n'avait pas eu de mentorat (X := 0). C'est le miroir de la cellule-type 3.\n", + "4. **Paradoxe de Simpson** : reproduire un exemple ou l'effet **s'inverse** entre l'agrege et le desagrege. Un cas classique : un traitement semble benefique dans chaque sous-population (homme, femme) mais nefaste au global, parce que les populations traitees different systematiquement par un confondeur.\n", + "5. **Demontrer que voir le barometre ne le fait pas tomber** : un cas ou l'observation d'une variable (par exemple le barometre qui descend avant la tempete) cree l'illusion d'une cause, alors que c'est un pur effet (la pression atmospherique baisse a cause du systeme, et la pression basse cause a la fois le barometre qui descend et la tempete qui vient)." + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "id": "82f6b148", + "metadata": { + "execution": { + "iopub.execute_input": "2026-10-05T19:53:15.016219Z", + "iopub.status.busy": "2026-10-05T19:53:15.015693Z", + "iopub.status.idle": "2026-10-05T19:53:15.029052Z", + "shell.execute_reply": "2026-10-05T19:53:15.026950Z" + } + }, + "outputs": [], + "source": [ + "# Exercice 1 : Ajouter un second confondeur W (motivation)\n", + "# Completez le stub ci-dessous.\n", + "\n", + "W = None # W : variable de motivation (0-1)\n", + "pass # TODO etudiant : ajouter W ~ Beta(2, 5), modifier le DAG, recalculer E[Z, W | X]" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "id": "ec13a7d2", + "metadata": { + "execution": { + "iopub.execute_input": "2026-10-05T19:53:15.032202Z", + "iopub.status.busy": "2026-10-05T19:53:15.031197Z", + "iopub.status.idle": "2026-10-05T19:53:15.044270Z", + "shell.execute_reply": "2026-10-05T19:53:15.043267Z" + } + }, + "outputs": [], + "source": [ + "# Exercice 2 : Verifier l'identifiabilite par backdoor avec {Z, W}\n", + "pass # TODO etudiant : binner Z et W, calculer P(Y | do(X=1)) par ajustement simultane\n", + "# Comparer avec la cellule-type 2 (ajustement par Z seul)" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "id": "9d8d419a", + "metadata": { + "execution": { + "iopub.execute_input": "2026-10-05T19:53:15.047286Z", + "iopub.status.busy": "2026-10-05T19:53:15.046760Z", + "iopub.status.idle": "2026-10-05T19:53:15.059347Z", + "shell.execute_reply": "2026-10-05T19:53:15.058193Z" + } + }, + "outputs": [], + "source": [ + "# Exercice 3 : Contrefactuel inverse\n", + "pass # TODO etudiant : tirage du pool 'Bob-like' (X=1, Y proche de 0.3),\n", + "# inférer Z, evaluer Y sous X=0, comparer avec Y_0=0.3" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "id": "28eacd67", + "metadata": { + "execution": { + "iopub.execute_input": "2026-10-05T19:53:15.061898Z", + "iopub.status.busy": "2026-10-05T19:53:15.061898Z", + "iopub.status.idle": "2026-10-05T19:53:15.075392Z", + "shell.execute_reply": "2026-10-05T19:53:15.073761Z" + } + }, + "outputs": [], + "source": [ + "# Exercice 4 : Paradoxe de Simpson\n", + "pass # TODO etudiant : simuler un cas ou l'effet agrégé est oppose a l'effet par strate" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "id": "db0e88a1", + "metadata": { + "execution": { + "iopub.execute_input": "2026-10-05T19:53:15.078476Z", + "iopub.status.busy": "2026-10-05T19:53:15.077973Z", + "iopub.status.idle": "2026-10-05T19:53:15.089924Z", + "shell.execute_reply": "2026-10-05T19:53:15.088828Z" + } + }, + "outputs": [], + "source": [ + "# Exercice 5 : Barometre et tempete\n", + "pass # TODO etudiant : illustrer un cas de confondant non observe ou de chaîne causale indirecte" + ] + }, + { + "cell_type": "markdown", + "id": "5b9d9834", + "metadata": {}, + "source": [ + "## Bibliographie\n", + "\n", + "- **Pearl, J.** (2009). *Causality: Models, Reasoning, and Inference* (2nd ed.). Cambridge University Press. Ch. 1 (l'echelle de Pearl), Ch. 3 (do-calculus, mutilation, backdoor), Ch. 9 (contrefactuels individuels, abduction-action-prediction).\n", + "- **Hernan, M. & Robins, J.** (2020). *Causal Inference: What If*. CRC Press. Ch. 1-2 (motivation, voir vs faire).\n", + "- **Scholkopf, B. et al.** (2021). *Towards Causal Representation Learning*. arXiv:2102.11107. §1 (pourquoi la causalite en machine learning).\n", + "- **Bareinboim, E., Correa, J., Ibeling, D. & Icard, T.** (2020). *On Pearl's Hierarchy and the Foundations of Causal Inference*. r60.pdf. §1 (la hierarchie rungs 1-2-3).\n", + "\n", + "## Verdict SOTA\n", + "\n", + "- **SOTA-OK** : la mise en oeuvre utilise les **outils canoniques** de l'inference causale (mutilation du graphe Pearl §3.3, backdoor adjustment, contrefactuel Pearl §9). `numpy` (simulation), `scipy` (statistiques), `networkx` (DAG), `matplotlib` (figures) sont **installables via pip** sans secret ni GPU. Aucun raccourci degrade, aucun stub externe.\n", + "- **Pas de GPU requis**, pas d'appel reseau, pas d'API externe — le notebook est **rejouable localement** sur la machine standard du depot (kernel `coursia-ml-training` ou Python 3.10+ avec `numpy`, `scipy`, `networkx`, `matplotlib`).\n", + "\n", + "## Suite de l'Origami causal\n", + "\n", + "Ce Pli 1 est le **socle**. Les plis suivants sont **verrouilles** jusqu'a validation de celui-ci :\n", + "\n", + "- **P-09** (Hernan target trial emulation) — transposition du cadre RCT en observationnel.\n", + "- **P-10** (ML for ATE / CATE) —heterogeneite du traitement par machine learning.\n", + "- **P-11** (Scholkopf 2021) — causal representation learning.\n", + "- **P-12** (Scholkopf 2012) — causal vs anticausal learning.\n", + "- **P-13** (Causal RL) —decision sequentielle sous causalite.\n", + "- **P-14** (Heckman) — modeles a selection.\n", + "- **P-15** (m-transportability) —transport de l'inference causale entre populations.\n", + "- **P-16** (ATE/ATT/CATE systematic) —recap formule des trois estimateurs.\n", + "\n", + "Au merge verifie du Pli 1, le coordinateur deplie le pli suivant par commentaire `[PLI N+1 DEPLIE]` sur l'EPIC parent #19309." + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.10.11" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/MyIA.AI.Notebooks/Probas/DecisionTheory/Causal-Bridges/CausalBridges-01-Do-Calculus.ipynb b/MyIA.AI.Notebooks/Probas/DecisionTheory/Causal-Bridges/CausalBridges-01-Do-Calculus.ipynb index a4676ee552..a2d4744ffc 100644 --- a/MyIA.AI.Notebooks/Probas/DecisionTheory/Causal-Bridges/CausalBridges-01-Do-Calculus.ipynb +++ b/MyIA.AI.Notebooks/Probas/DecisionTheory/Causal-Bridges/CausalBridges-01-Do-Calculus.ipynb @@ -15,6 +15,8 @@ }, "source": [ "# Du graphe causal au do-calculus — le pont entre les quatre séries causales\n", + "", + "[← CausalBridges-00 — Pearl Ladder](CausalBridges-00-PearlLadder-Intro-Python.ipynb) | [↑ Causal-Bridges](../README.md) | [↑ DecisionTheory](../../README.md)", "\n", "> **Notebook-pont** de la constellation causale du dépôt (règle F : `dowhy` installé et exécuté réellement, SOTA-OK). Il ne remplace pas les notebooks dédiés de chaque moteur : il donne l'**armature formelle unifiée** (échelle de Pearl + trois règles du do-calculus) et la fait tourner sur l'**outil de référence** [`dowhy`](https://www.pywhy.org/dowhy/), avant de renvoyer à chaque série pour l'instanciation par moteur.\n", "\n", diff --git a/MyIA.AI.Notebooks/Probas/DecisionTheory/Causal-Bridges/README.md b/MyIA.AI.Notebooks/Probas/DecisionTheory/Causal-Bridges/README.md index 70eefce894..09573fab33 100644 --- a/MyIA.AI.Notebooks/Probas/DecisionTheory/Causal-Bridges/README.md +++ b/MyIA.AI.Notebooks/Probas/DecisionTheory/Causal-Bridges/README.md @@ -10,6 +10,7 @@ | Notebook | Durée | Concepts | |----------|-------|----------| +| [CausalBridges-00 — Pearl Ladder introductive](CausalBridges-00-PearlLadder-Intro-Python.ipynb) | ~25 min | Socle commun à tous les carnets Causal-Bridges : observation vs intervention vs contrefactuel sur un SCM-jouet (`Snow`/`Doll-Hill` côté santé, `V → T → Y` côté formation) ; différences P(Y\|X) / P(Y\|do(X)) / contrefactuel Alice mesurées sur le même DGP ; ouvre la série avant CB-01 | | [CausalBridges-01 — Do-Calculus](CausalBridges-01-Do-Calculus.ipynb) | ~65 min | Échelle de Pearl, trois règles du do-calculus, critères *backdoor* / *front-door* exécutés avec `dowhy` ; quatre tâches du **data-fusion** (sélection corrigée par IPW, transportabilité stratifiée), **CHT démontré machine** (deux SCM gaussiens à loi jointe identique, interventions opposées), jonction do-calculus ↔ baseline Shapley (`do` vs `voir`) ; Pearl (intervention) vs Hoel (émergence causale) | | [CausalBridges-02 — Exiger un estimand](CausalBridges-02-Dowhy-Estimand-Intervention.ipynb) | ~45 min | Identification causale **nommée** via `dowhy` (backdoor, front-door, instrumentale) sur un cas complet ; sensibilité au graphe **mesurée** quand une hypothèse saute | | [CausalBridges-03 — Le contrefactuel individuel](CausalBridges-03-Dowhy-Contrefactuel-Individuel.ipynb) | ~40 min | Troisième échelon de Pearl : `dowhy.gcm` (abduction-action-prédiction) sur **un individu** ; l'effet moyen nul cache une CATE linéaire ±3 ; fragilité du chiffre individuel à la spécification du mécanisme |