Python/Learning & Inference in Discrete (Multi-valued RVs) Bayesian Networks.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "de1cd439",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<style>.container { width:100% !important; }</style>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from IPython.display import display, HTML\n",
    "display(HTML(\"<style>.container { width:100% !important; }</style>\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "c3bc0918",
   "metadata": {},
   "outputs": [],
   "source": [
    "import math\n",
    "import numpy as np\n",
    "import numpy.matlib\n",
    "import pandas as pd\n",
    "import networkx as nx        # for graphical model visualization\n",
    "from scipy import linalg\n",
    "from numpy.linalg import inv\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from ismember import ismember"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7e2117de",
   "metadata": {},
   "source": [
    "## 1. [Modeling] Construction of a Bayesian Network Graphical Model:\n",
    "\n",
    "- Assuming arbitrary parameters and a random network structure for each factor in the conditional probability distribution.\n",
    "- Assuming random variables with different cardinalities."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "5e046417",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAgMAAAGFCAYAAABg2vAPAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAA9hAAAPYQGoP6dpAABg00lEQVR4nO3dd3xN5x/A8c/NMoK0VotStZXYQSp27Vi1YxWhQ4ldqvaqPYrSxF6lRoJQe+9YSVC0fmLEbiQxIuv+/jjJIQMJNzl3fN+vV17JGfeeb8R97vee5/k+j06v1+sRQgghhMWy0joAIYQQQmhLkgEhhBDCwkkyIIQQQlg4SQaEEEIICyfJgBBCCGHhJBkQQgghLJwkA0IIIYSFs0nJSbGxsQQHB5M1a1Z0Ol1axySEEEIIA9Dr9YSHh5M3b16srF7/+T9FyUBwcDD58+c3WHBCCCGESD83b97kk08+ee3xFCUDWbNmVZ8sW7ZsholMCCGEEGkqLCyM/Pnzq+/jr5OiZCC+ayBbtmySDAghhBAm5m1d/DKAUAghhLBwkgwIIYQQFk6SASGEEMLCSTIghBBCWDhJBoQQQggLJ8mAEEIIYeEkGRBCCCEsnCQDQgghhIWTZEAIIYSwcJIMCCGEEBZOkgEhhBDCwkkyIIQQQlg4SQaEEEIICyfJgBBCCGHhUrSEsenQA+FAJGAHZAXevGyjEEKkDWmPhOkwg2QgEFgNnAD8gLBXjmUDKgFVADegdLpHJ4SwJNIeCdOk0+v1+redFBYWhoODA6GhoWTLli094koBX2AScAQlp4lBycQT0wHWQDRQDfgJaJxOMQohLIO0R8I4pfT92wTHDDxCyapdgWNx+6JJ/oVH3P7ouJ+PAU2AjsB/aRijEMIySHskzIOJJQP+wOfAurjt2FQ+Pv78tUBJIMBAcQkhLI+0R8J8mFAy4A9UR8nEY97zuWLinscFeQEKIVJP2iNhXkwkGXgE1AOe8v4vvHgxcc/3JXKLTgiRctIeCfNjIslAHwyTgScWn5H3MfDzCiHMl7RHwvyYQDLgC6zB8C+8eDEopUDb0uj5hRDmQ9ojYZ5MIBmYxNKlOnQ6EnzlygW1asHWrck/asAAKFtW+XnvXujeHUqUAHt7yJcPmjeH06fjz7ZGKQsSQog3mcSrzeacOUp7VPotUwa82h6Fh8OQIVC/vtKO6XQwevSrZ0t7JNKfkScDgSh1u0qZzpIlcOwYHD0Kv/8O1tbQtCls2ZL0kRs3QqtWys+//QbXr4OHB2zbBrNnw/37ULWqkigo2fhh4EJ6/FJCCJMU3x69rBpYvFj5fuECnDjx+ke+2h49eqS0Xy9eQIsWyZ0t7ZFIf0Y+A+FqlBCVutzSpaFSpZdHGzaEDz+ENWuUpCDeqVMQFPTyxTdvHuTOnfCZGzaEIkVg4kSoU4e466wGJqTZbyOEMGUJ2yM/Pzh/Hpo0AV9fWLQIqlRJ+qjE7dGnn0JIiHJH4OFD8PJK7lrSHon0ZeR3Bk7wcoKOpDJmBDs7sLVNuH/DBiheHEqVUrYTJwIAWbLA55/DzZvxe2LirieEEMlJ2B4tWqR8/+UX+OIL+OMPePYs6aMSt0fxXZ1vJu2RSF9GnAzoUeb2fikmBqKjISoKbt2Cfv3g6VNwc0v4yA0bXmbhrxMaCmfOvHyBKtc7xetnDhNCWK6E7dHz58odSScn5Y5l9+7KWIA//0z6yJS0R8lfT9ojkX6MOBkIJ+EiH0ofv62tcjcgf35YuBDmzoUGDV6ec/48/PPP2198vXsricTw4a/uDQOeGCh+IYT5SNgerV+vfKDo0UPZbtdOudsYf7cgXkrbo+RJeyTSjxGPGYhMsmf5cihZUvn54UPYtEl5U4+JgR9+UPZv2AAFC0KFCq9/5hEjYNUq+PVXqFjx7dcVQli6hO3CokWQKRO0b69sZ8kCbdoog5yvXoWiRZX9KWmPUnNdIdKKEd8ZsEuyp2RJZQBhpUrKAMCFC5XynCFD4PFj5Zz169+chY8ZA+PHw4QJLxOIt11XCGHpXrYL//wDBw8qAwf1eqXtefwYWrdWjsdXGMDb26PUXFeItGTEyUBWlPW/36xMGaX/7soVuHRJ+Xrdi2/MGKWed/Ro+Omn5M7IBmR554iFEObqZXu0eLGSBKxfr1QzxX81aaKcuWyZcrfybe3R20l7JNKPEXcT6IBKwN43nnXunPI9Vy7l1n/evMrYgsTGjVOSgJ9/hlGjXnc9p7jvQgjxKqU9io7ey7JlULhw8iWBW7fC9OmwfbvSNr2uPUrZ9aQ9EunHiJMBgCrAQeLLeQIDlWoCUCbu2LgRdu2Cli3hs8+UTP2rr5KW7UyfDiNHKl0LTZrA8eMJjysvVuu46wkhxEt6vZ4jR44QEnKfmBgIDobJk5UZUBMrXVoZ1LxoEfzvf8m3R6AkC0+fKhUIABcvKu0XQOPGkDmztEcifen0ev1ba1fCwsJwcHAgNDSUbNnefuvecAIBR5YuhW7dEh5xcFASgC5d4PvvlVLDIkVg376kL9JateDAgddf5eW/QCBQ6vUnCiEsxoMHD1i2bBleXl5cvnyZUqWUgYHbtintTa5cyT+uQwfljT06Ovn2CJRBhUFByT/+f/9Tjkt7JAwhpe/fRp4MgLLG9zFenQI0OVOmwLRpcOeOMk1x6lgDzsChd4pQCGEeYmNj2b17N56envj4+BAVFZXg+LFj1lSuHIuV1ZubTWmPhLEwo2TAF3BNp+s0TofrCCGMzc2bN1myZAmLFy8mKJmP7DVr1qRnz560bp2JDBneqzwghaQ9EoaR0vdvIx8zANAE6ACsI22WDbUG2iEvPCEsS1RUFFu3bsXLy4u//vqL2NiEdx9z587N119/TY8ePShWrNgrR6Q9EubHBJIBgF+BPcAjDPsCtAZyxD2/EMISXL16lUWLFrF06VLu3buX4JhOp6Nhw4b07NkTV1dXbBMvfALEt0exsQ+xsnpz92XqSHsktGMiyUAOYDfK+IGnGCYhsAbs4543uwGeTwhhrJ4/f87GjRvx9PTkQDKjiQsUKECPHj3o1q0b+fPnf8uz5SAwcBYFCriROTPYGKQVlfZIaMtEkgEAR5Q1vr/k/e8QxGfgu+OeVwhhjvz9/fH09GTlypU8jp+mNI6NjQ0tWrTA3d2dL7/8EusUjvQLCgriyy/7kysX7N4NOXPqsLZ+9wWF9HordDppj4S2TCgZAOWFcgnog7LWtzWpSwriz28PzEEycCHMT3h4OH/88Qeenp6cOnUqyfHixYvj7u5Oly5dyJ3c+uZvEBoaSpMmTbh37x737kGPHs74+BQA1pLa9ig6WrmrsH27Aw0aBGJt/ZpaRSHSgRFPR/w62YFVKKNtneP22fC6mbpiY5UljxXOcY9biSQCQpgPvV7P8ePH6dGjB3ny5KFXr14JEoGMGTPSpUsXDh48yKVLlxg0aFCqE4GoqCjatm3LhQsXAChatCjLlm3B2voPUtoeKfuVz2D+/vY0bgxNmoSwZIlPqmIRwtBMoLTwbS6g3CU4gbL+98tlRvX6bBw8+IzDh6Px8cnMoUP/kSFDBo3iFEIY2qNHj1ixYgVeXl7qm/SrypUrR8+ePXFzc+ODDz545+vo9Xq+++47Fi5cCED27Nk5fvw4ReOXJ1S9vj1S1hpwQplZ0I2DBx9Rs2ZNAHLlysXVq1dxcHB45xiFSE6K37/1KRAaGqoH9KGhoSk5XUOxer0+TK/XP4z7Hqvv2rWrHtADel9fX23DE0K8t5iYGP3u3bv17du319vZ2amv7/ivrFmz6r/99lu9n5+fwa45depU9fnt7Oz0Bw8eTMGjkrZHibVt21Z93oEDBxosXiHipfT92wzuDLzZ5s2bad68OQA9evTAK7nVRYQQRi84OJilS5eyaNEirl27luR4tWrVcHd3p02bNtjb2xvsuhs3bqR169bEN5UrV66kY8eOBnnuoKAgSpQoQUREBDY2NgQGBlK8eHGDPLcQkPL3bxMcM5A69erVUxsGHx8fYmLSYqIQIURaiI6OZsuWLTRv3pwCBQowfPjwBIlAjhw5GDBgABcvXuTw4cN8/fXXBk0ETp06RadOndREYPTo0QZLBAA+/fRThgwZAii/68CBAw323EKkhtnfGQBo06YN6+OWBNu/f7/aTyeEME7Xrl1j8eLFLFmyhODg4CTH69Wrh7u7O82bN0+zcUBBQUFUqVJFnZioU6dOLF++HF1yyxC+h6dPn1KiRAlu3boFwLZt22jUqJFBryEsl9wZeEXLli3Vnzdu3KhhJEKI13nx4gVr167lyy+/pHDhwkyYMCFBIpAvXz5GjBjBtWvX2LlzJ23btk2zRCA0NBRXV1c1EahRowZeXl4GTwQA7O3tmTJlirrdv3//JAskCZHWLOLOQGhoKLly5SIqKor8+fMTFBSUJi9qIUTqXbhwAS8vL1asWMGjR48SHLO2tqZp06a4u7vToEEDbAwz3d8bRUVF4erqys6dOwGlhPDYsWPkyJEjza6p1+upXr06R44cAWDmzJn069cvza4nLIcZrVpoGI0bN2b79u2A0g9YqVIljSMSwnI9ffqUtWvX4uXlxbFjx5IcL1y4MO7u7nTt2pU8efKkW1z6FJcQGt7p06dxcnJCr9fj4ODA1atXyZVLJiIS70e6CRKRrgIhtKXX6/Hz8+Obb74hT5489OjRI0EikCFDBtzc3Ni7dy9Xrlxh6NCh6ZoIAEyfPl1NBOzs7PD29k6XRACgYsWKdOvWDVDuZo4YMSJdrisEWNCdgfv37/Pxxx+j1+spUaIEly5d0jokISxCSEgIq1atwsvLi/Pnzyc5Xrp0aXr27EmnTp3Inl27mUE3bdpEq1at0qSEMKXu3r1LsWLFCA8PR6fTcebMGcqVK5euMQjzIncGEsmdOzcuLi4A/P3335IMCJGG9Ho9Bw4coHPnzuTNm5c+ffokSATs7e1xd3fn+PHj+Pv707dvX00TgVOnTtGxY8c0KyFMqY8//li9I6DX6+nXrx8p+LwmxHuzmGQA4KuvvlJ/3rRpk4aRCGGe7t27x5QpUyhevDi1atVi5cqVREREqMerVKmCp6cnd+7cwdPTkypVqmg+mDcoKIimTZvy/PlzQCkhHDlypGbx9O3blyJFigBw4MABNmzYoFkswnJYTDcBKC/6ggULAkr/nJ+fn7YBCWEGYmJi2LlzJ15eXmzevJno6OgExz/88EM6d+6Mu7s7jo7GtURvaGgoLi4uBAYGAlC9enV27dql+RomW7ZsoVmzZoAyMdGlS5fIlCmTpjEJ0yTdBMn49NNPqVChAqCM3A0KCtI4IiFMV1BQEKNHj+azzz6jcePGbNy4MUEiULt2bVavXk1wcDCzZ882ukQgfhXC+ESgaNGibNq0SfNEAMDV1ZX69esDyr/z9OnTNY5ImDuLSgYgYVeBt7e3doEIYYIiIyPZsGEDDRs25LPPPmPMmDHcvHlTPf7xxx8zdOhQrl69yt69e+nQoQMZM2bUMOLk6fV6+vTpo84lkD17dnx9fdN0LoHU0Ol0zJw5E2trawAmTZqkzlAoRFqw6GRASgyFSJnLly8zePBgPvnkE1q3bs2OHTvUgW1WVlY0adIEb29vbty4waRJk9Q+b2M1Y8YMzUoIU+rzzz+nd+/eADx79oyhQ4dqHJEwZxY1ZiBeiRIluHz5Mjqdjrt375I7d26tQxLC6Dx79oz169fj5eXFoUOHkhwvWLAgPXr04Ouvv+aTTz7RIMJ3k7iEcMWKFXTq1EnjqJIXEhJC0aJF1ZkZjxw5whdffKFxVMKUyJiBN4i/O6DX69m8ebPG0QhhXM6ePUvv3r3JmzcvXbt2TZAI2Nra0rZtW3bu3Mm///7Lzz//bFKJQHIlhMaaCIAy+HLs2LHqtoeHB7GxsRpGJMyVRd4ZOHXqFJUrVwagUaNGbNu2TeOIhNBWWFgYq1evxsvLi9OnTyc5XrJkSdzd3encubPJTpEbFBRE1apVuXv3LpB2qxAaWnR0NOXLl1cHOi5ZsoSvv/5a26CEyZC1Cd5Ar9dToEABbt26hZ2dHffv38fBwUHrsIRIV3q9nqNHj+Ll5cW6det49uxZguOZMmWiXbt29OzZE2dnZ6N/03wTYy0hTKk9e/bw5ZdfAsogzStXrpA1a1aNoxKmQLoJ3kCn06ldBZGRkXJnQFiUBw8eMGPGDEqVKoWLiwtLly5NkAhUrFiR3377jTt37rBkyRK++OILk04EjLmEMKXq1q2rrq9y9+5dJk6cqHFEwtxYZDIAsnCRsCyxsbHs2rWLdu3akS9fPgYOHJhgSm4HBwe+//57zpw5g5+fH99++61Z3C0z9hLC1Jg2bRp2dnaAUg3x77//ahyRMCcWmwy4uLiQM2dOALZv365ORSqEObl16xbjx4+ncOHC1K9fn3Xr1hEVFaUer169OsuXLyc4OJh58+ZRvnx5DaM1PFMoIUypQoUKMXDgQEC5ozlo0CCNIxLmxGKTARsbG3W6z6dPn7Jr1y6NIxLCMKKiovDx8cHV1ZVPP/2UESNGcP36dfV4rly5GDx4MH///TcHDx6kc+fOZM6cWbuA08imTZsYPHiwur1o0SKqV6+uYUTvb9iwYeqyzt7e3uzevVvjiIS5sNhkAGQCImFe/vnnH4YNG0aBAgVo0aIFvr6+ahmaTqejQYMGrF+/nlu3bqmLCZmrxCWEo0aNMuoSwpTKmjUrv/zyi7rdr1+/JGtBCPEuLLKaIF5ERAS5c+cmPDyc7Nmzc/fuXWxtbbUOS4gUi4iIYOPGjXh5ebFv374kx/Pnz0/37t3p1q0bn376qQYRpj9TLSFMqdjYWJydnTl58iQAc+fOVWcqFCIxKS1Mofbt27N27VoAdu/eTd26dTWOSIi3CwgIwMvLixUrVhASEpLgWHwXmLu7O/Xr11fnt7cEpl5CmFLHjx/H2dkZUCYmunr1qkkOihRpT0oLU0i6CoSpCA8Px8vLi6pVq1KmTBnmzJmTIBEoWrQokydP5ubNm2zYsIFGjRpZVCIQHR1Nu3btTLqEMKWqVq1K586dAWXK4tGjR2sbkDB5Fn9nIDw8nFy5cvHixQvy5s3LzZs3sbKy+BxJGAm9Xs/Jkyfx8vLijz/+4MmTJwmOZ8yYkdatW+Pu7k6NGjXM5lZ4aun1er777ju1ciB79uwcP37cZCsHUuL27dsUL16cp0+fYm1tzblz5yhdurTWYQkjI3cGUihr1qzUq1cPgODgYLUfTggt/ffff8yZM4eyZctStWpVvLy8EiQCZcqU4ddffyU4OJgVK1ZQs2ZNi00EIGEJoa2tLZs2bTLrRAAgX758/PTTTwDExMTQr18/UvDZTohkWXwyAAm7CjZt2qRhJMKSxcbGsm/fPjp27EjevHnx8PAgICBAPZ4lSxZ69erFqVOnOHfuHD/88AMffvihhhEbh8QlhIsXL6ZGjRoaRpR+BgwYQMGCBQFlymJZeE28K4vvJgB4+PAhH330EbGxsRQpUoQrV65Y9Kcskb7u3LnDsmXL8PLySnZWOWdnZ3r27EmbNm3IkiWLBhEar1OnTlGzZk110rBRo0ZZXP/5xo0badWqFaBMTHTx4kWzHCch3o10E6RCzpw5qVmzJqDUascPQBIirURHR+Pr60uLFi3Inz8/w4YNS5AI5MiRg/79+xMYGMjRo0fp1q2bJAKJ3Lhxg2bNmqmJQKdOnRg1apTGUaW/li1bUrt2bQCuXbvGrFmztA1ImCRJBuJIV4FID9evX2fEiBEULFgQV1dXfHx8iImJUY9/+eWXrFmzhtu3b6uLCYmkwsLCaNKkiTqXQPXq1fHy8rLIO3o6nY5Zs2apA5/Hjx/PnTt3NI5KmBpJBuK0aNFC/VlKDIUhvXjxgnXr1lG/fn0KFSrE+PHjuX37tno8b968DB8+nH///Zddu3bRvn17uc37BtHR0QlWISxSpIjZlhCmVJkyZfjmm28AePLkiTqwUIiUkjEDr6hSpYpaTfDvv/9SqFAhjSMSpuzSpUt4eXmxfPlyHj58mOCYtbU1TZo0wd3dnUaNGmFjY6NRlKbFEksIU+rhw4cULVqUx48fA3Dy5EmcnJy0DUpoTsYMvANZ1li8r6dPn7J06VJcXFz4/PPPmTFjRoJE4LPPPmPChAncuHEDHx8fmjZtKolAKlhiCWFK5cyZkzFjxqjbffv2lVJDkWJyZ+AVV65cURdv+eKLLzhy5IjGEQlTcfr0aby8vFi9ejVhYWEJjtnZ2fHVV1/h7u5O7dq1ZVKrd+Tt7c1XX32lvsGtWLHCLBYfMqSoqCjKli3LpUuXAFi5ciUdO3bUOCqhJVmb4B2VLl2aCxcuoNPpuH37trpcqBCJPX78mNWrV+Pl5cXZs2eTHP/888/p2bMnnTt3lnnj35Ofnx81atSw6BLClNqxYwcNGzYElPEoly9flkoUCybdBO8ovqtAr9fj4+OjcTTC2Oj1eg4dOkTXrl3JmzcvvXv3TpAIZM6cme7du3Ps2DECAwPp16+fJALv6caNGzRt2lRNBDp27GiRJYQp1aBBA5o2bQoos6q+uuSxEK8jdwYSOXv2LBUqVACgXr167Ny5U+OIhDG4f/8+y5cvx8vLi8uXLyc57uTkRM+ePWnXrp3Zv0bSU1hYGNWqVVMrB1xcXNi9e7dFVw6kxNWrVylVqhRRUVFkyJCBv//+W52pUFgWuTPwjsqVK6e+aPbt25dkeVhhOWJiYtixYwetW7cmX758DB48OEEi8MEHH9CnTx/OnTvHyZMn6dmzpyQCBpRcCaG3t7ckAilQtGhR+vXrByilra9O1yxEciQZSESn06ldBdHR0WzdulXjiER6u3HjBmPGjKFQoUI0bNiQDRs2EB0drR6vVasWK1euJDg4WF1MSBiWXq+nT58+7NixA1BKCH19faXLJRV+/vlncufODcD69evZv3+/tgEJoybJQDJenY1QSgwtQ1RUFBs3bqRx48YULFiQ0aNHc+PGDfX4Rx99xI8//siVK1fUxYQyZcqkYcTmbebMmSxYsAB4WUJYrFgxjaMyLdmyZWPixInqtoeHR4LZLoV4lYwZSEZMTAz58uXj3r17ZMqUiQcPHmBvb691WCINXLlyhUWLFrF06VLu37+f4JiVlRUNGzbE3d0dV1dXbG1tNYrSskgJoeHExMRQuXJlzpw5A8CCBQvUmQqFZZAxA+/B2tqa5s2bA/D8+XP1VqUwD8+fP2flypXUqlWL4sWLM2XKlASJQIECBRgzZgzXr1/H19eXli1bSiKQTvz8/HBzc1MTgVGjRkki8B6sra2ZPXu2uj18+HAZByWSJcnAa8jCRebn/Pnz9OnTh7x589K5c2cOHDigHrO1taV169b89ddfXLt2jZEjR5I/f34No7U8UkKYNlxcXGjfvj0Ajx49YuzYsRpHJIyRdBO8RmRkJLlz5yY0NBQHBwfu37+PnZ2d1mGJVAoLC+OPP/7A09MTPz+/JMeLFy+uTgwUP9hKpD8pIUxbN27coESJEjx//hwbGxv8/f0pWbKk1mGJdCDdBO/Jzs4OV1dXAEJDQ9m3b5/GEYmU0uv1HDt2jB49epAnTx6++eabBIlApkyZ6Nq1K4cOHeLSpUsMHDhQEgENySqEaa9AgQL8+OOPgPLv3b9/f1m3QCQgycAbSFeBaXn48CGzZs3C0dGRL774gsWLF/Ps2TP1ePny5Zk/fz7BwcHqYkI6nU7DiMXrSghz5sypcWTmZ/DgwWrX144dO9i2bZvGEQljIt0Eb/D06VNy5sxJREQEH330Ebdv38ba2lrrsMQrYmNj2bdvH56enmzatInIyMgEx7Nly0bHjh1xd3dXZ5YUxmPGjBkMHDgQUMZt7N69mxo1amgclflau3atOn6gWLFiBAQESPenmZNuAgOwt7dXF/y4d+8ex44d0zgiEe/27dtMmDCBIkWK8OWXX7J27doEiYCLiwtLly4lODiY+fPnSyJghLy9vRk0aJC6vXjxYkkE0ljbtm1xcXEBlLLauXPnahyRMBaSDLyFdBUYj+joaDZv3kyzZs0oUKAAP//8M//73//U4zlz5mTgwIFcvHhRXUxI5ocwTolLCEeOHCklhOlAp9Mxe/ZstXtszJgxSebXEJZJugneIiQkhNy5cxMdHU3BggW5du2a9DOns3///ZfFixezZMkS7ty5k+CYTqejXr16uLu707x5c7nlaQJu3LhBlSpVuHv3LgBubm6sXLlSXlfpqGfPnnh5eak///777xpHJNJKSt+/JRlIgfr167Nr1y4Azpw5Q/ny5TWOyPxFRETg7e2Nl5cXe/bsSXI8X758dO/ene7du8tqbCZESgiNw7179yhWrBhhYWHodDpOnz4t7ZqZkjEDBiRdBennwoUL9O/fn3z58tGhQ4cEiYC1tTUtWrTA19eXoKAgxo4dK4mACYmOjqZdu3ZSQmgEPvroI0aOHAkoFR0eHh5Samjh5M5ACty5c4d8+fKh1+spVaqU2pgJw3jy5Anr1q3D09OT48ePJzlepEgR3N3d6dq1Kx9//LEGEYr3pdfr+f7779XFh7Jnz86xY8dk8SENRUZGUrp0aa5evQoolQZt27bVOCphaNJNYGDVqlXj6NGjAFy+fFkasfek1+vx8/PD09OTNWvW8OTJkwTHM2TIQOvWrXF3d6dmzZrSn2zipITQOPn6+qqTqxUoUIBLly6ROXNmjaMShiTdBAYmXQWGERISwty5cylXrhyVK1fG09MzQSLg6OjInDlzCA4OVhcTkkTAtCUuIVy0aJEkAkaicePGavn0jRs3mDZtmsYRCa3InYEUunbtGoULFwagcuXKnDhxQuOITIder+fgwYN4enqyfv16Xrx4keB4lixZ6NChA+7u7jg5Ocmbvxnx8/OjRo0a6uJDI0eOZMyYMRpHJV516dIlypQpQ3R0NJkyZeLy5cuySJcZkTsDBlaoUCHKli0LwMmTJ7l165bGERm/u3fvMnnyZIoXL06tWrVYtWpVgkSgatWqeHl5cefOHX7//XcqV64siYAZSbwKoZubG6NHj9Y2KJFEyZIl+eGHHwBlee/4NQyEZZFkIBVe7Srw9vbWLhAjFhMTw7Zt2/jqq6/Inz8/Q4cOVQcogTJwzMPDA39/f3UxoSxZsmgYsUgLYWFhuLq6qnMJuLi4sHjxYkn2jNTIkSPV9SDWrFnD4cOHNY5IpDdJBlKhZcuW6s8bN27UMBLjExQUxKhRoyhYsCBNmjRh06ZNREdHq8fr1KnDmjVruH37trqYkDBP8SWEAQEBABQuXFhKCI3chx9+yPjx49VtDw8PYmNjNYxIpDcZM5AKer2eYsWK8c8//2Btbc3du3ctenW1yMhINm/ejKenJ7t27UpSp5wnTx66detG9+7d1fEWwrwlLiH88MMPOX78uFTfmICYmBgqVKiAv78/oAz07N69u8ZRifclYwbSgE6nU+8OxMTEsGXLFo0j0sbff//N4MGD+eSTT2jTpg07d+5UEwErKyuaNm2Kj48PN27cYMKECZIIWJCZM2eqiYCtrS3e3t6SCJgIa2trZs+erW4PGzaMsLAwDSMS6UmSgVSy1BLDZ8+esXz5cqpXr07JkiWZNm0aDx48UI9/9tlnjB8/nhs3bqiLCdnY2GgYsUhviUsIvby8pITQxNSqVYvWrVsDcP/+/QRdB8K8STdBKsXGxpI/f36Cg4Oxs7Nj9erVHD16lBIlSuDu7m52A6TOnj2Lp6cnq1atSvIpwc7OjpYtW+Lu7k6dOnWwspLc0lL5+flRs2ZNnj17BkgJoSm7fv06JUqU4MWLF9ja2nLhwgWKFi2qdVjiHaX4/VufAqGhoXpAHxoampLTzVp4eLi+fv36eiDJ18WLF7UOzyAeP36snz9/vr5ChQrJ/p4lS5bUz5gxQ//gwQOtQxVGICgoSP/xxx+r/z/c3Nz0sbGxWocl3sPw4cPVv2fTpk21Dke8h5S+f8udgVT4448/6N69u1o3ndilS5coUaJEOkdlGHq9nqNHj+Lp6cm6deuS/I6ZM2emXbt2uLu74+zsbHZ3QMS7CQsLw8XFRa0ccHFxYdeuXWTMmFHjyMT7ePLkCcWLFyc4OBiAv/76iwYNGmgclXgXKX3/lk7dVFiyZMlrEwEbGxuTHCj34MEDli9fjpeXF3///XeS45UqVcLd3Z0OHTpYdCIoknpdCaEkAqYvS5YsTJ48mc6dOwPQv39/zp8/j62trcaRibQinbyp0Ldv39f2ixcuXNhkXiixsbHs3LmTtm3bki9fPgYNGpQgEXBwcKB3796cPXuWU6dO8c0330giIBLQ6/X07duXv/76C1BKCLdt22bRpbbmxs3NjapVqwLKXc/ffvtN44hEWpJkIBWaNGnCihUrkr1FbgrdA7du3WLcuHEUKlSIBg0a8OeffxIVFaUer1GjBitWrODOnTvqYkJCJGfWrFnqm4OtrS2bNm2SEkIzY2VllaDUcNSoUTx8+FDDiERakmQgldzc3Fi8eHGS/cbaEEZFReHt7U2TJk349NNPGTlyJEFBQerx3LlzM2TIEC5fvsyBAwfo1KkTmTJl0jBiYey8vb3V5YhBKSGsWbOmhhGJtFK5cmW6du0KwOPHjxk5cqTGEYm0IsnAO/j6669ZuHBhgn329vZxP+mBMOBh3Pe3js9ME//88w9Dhw4lf/78tGzZkm3btqnTi+p0Oho1asSGDRu4efMmkydPNtpkRhiX06dP07FjR3WSqZEjR9KlSxeNoxJpadKkSer6IQsXLlRnKBRmxpClCZbGw8NDD+hLlUJ/9253vV5fR6/XZ9Pr9bzylS1u/zC9Xh+QpvE8f/5cv2rVKn2tWrWSLQnMnz+/fvTo0fqgoKA0jUOYJykhtFyTJk1S/+61a9eWv7sJkdLCdOFLePhPZM3qj1KYEUPydwJ0gDUQDVQDfgIaGyyKgIAAPD09WblyJSEhIQmO2djY0Lx5c9zd3alXrx7W1tYGu66wHIlLCKtVq8bu3bulcsBCREREUKpUKa5duwbAhg0bEszGKoxXSt+/JRl4J4+APsAalJ6W1KzuFX++G/ArkP21Z/7vf/8jV65cyS7xGx4ezh9//IGXlxcnT55McrxYsWK4u7vTpUsXPvroo1TEJ0RC0dHRNG3aVK0cKFy4MMePH5fKAQvj7e2trs1SsGBBLl26JMmgCZCFitKMP/A5sC5uO7XLfMafvxYoCQQke9bIkSMpVKgQ5cuX58mTJ4BSznX8+HHc3d3JkycPvXr1SpAIZMyYkc6dO3PgwAF1MSFJBMT70EsJoYjTvHlz6tatCyhTFs+cOVPjiIQhyZ2BVPEHqgNPUboE3pc1YA8cBhzVvXPmzMHDw0PdnjVrFqCM2g4MDEzyLOXKlaNnz564ubnxwQcfGCAuIRQzZ85kwIABgFJCuGvXLqkcsGCBgYGULVuW2NhY7O3tuXLlCnnz5tU6LPEG0k1gcI9Q7gg8wjCJQDxrIAdwCcjOunXraN++Pa/+WXQ6HYn/TFmzZqVjx464u7tToUIFmR5YGJyPjw8tW7ZU/+8tW7ZMKgcEP/zwA/PmzQOgS5cuLFu2TOOIxJtIMmBwbihdA4ZMBOJZA+3Yu7cHjRo1IjIy8rVnVqtWDXd3d9q0afNKOaMQhnX69Glq1KihrkI4YsQIxo4dq3FUwhg8evSIokWLqoOVjx8/TpUqVTSOSryOjBkwoKVL+6PTrUGni0Gng4wZ4eOPoXZtmDQJ7t9P/nH//Qe2trBpU9Jjz59DsWKg08G0aTHAambOrPfaRKB8+fJcuHCBw4cP8/XXX0siINLMjRs3cHV1VRMBNzc3WY5YqHLkyJEgMfTw8FDnMBGmS5KBFNkKwJIlcOwY7NoF8+ZBuXIweTKULAm7dyd9lI8P2NlBw4ZJj40YAU+fvtyOjoYhQ17/grp+/bpJLoQkTEtYWBiurq7cvXsXUO5ELVq0SLqhRALffvstpUqVAuDEiROsWrVK44jE+5Jk4K0CgX8AKF0aqlaF6tWhVSuYORP8/cHeHr76Cu7dS/jI9euhcWNIPLvvyZPw66/wyrTf2Ngoz/v558lHERISQmhoqOF+LSESSW4VQm9vbykfE0nY2NioA5sBfvzxR7XqSZgmWcL4rVbzprkEChSA6dOhbVtYuBDip+4OC1PuFiQeWxMZCd27Q+/eUKlSwmN6vTW+vu04ebI59+7dU7/+++8/GjZsSO7cuQ3+2wkByZcQ+vr6SgmheK0vv/yS5s2b4+Pjw507d5g0aRITJkzQOizxjiQZeKsTvG0ugcaNwdoaDh58uW/LFmU8QJMmCc8dO1bpHhg3Dh48SHhMp4ulYMF7FCzY1jChC5FCya1CWLx4cY2jEsZu2rRpbN++ncjISKZPn06PHj0oVKiQ1mGJdyDdBG+kB/zeepa9PeTMCcHBL/etXw/160PWrC/3nTsHU6bAggXKY5K/3im0WtxIWCYfHx9ZhVC8kyJFitC/f38AXrx4weDBgzWOSLwrSQbeKBxl5cG3e7VA8+lT2LFDGVcQLzpa6R5o1w4aNHjTM4UB0vcm0sfp06dxc3NT5xIYMWKEzCUgUmX48OF8/PHHAGzcuJG9e/dqHJF4F5IMvNHr6/1f9fQpPHoE8RNx+fpCVBQ0a/bynFmz4No1GDUKHj9WvsLi8oyICGU7Rp3CIGXXFeJ93Lx5k6ZNm0oJoXgvWbNmZdKkSep2v379iI6O1jAi8S4kGXgjuxSd5eurvJHXqqVsb9gAderAhx++PCcwEEJDoWhRZf+HH0LZssqxESOU7QB1mYKUXVeIdxUWFkaTJk24c+cOICWE4v106dKFSnEjouNXURWmRZKBN8oKvHnGxRs3YNAgcHCAb75RPuVv25awiwBg6FDYty/h15o1yrFvv1W2ixQh7npJVykUwlCkhFAYmpWVFXPmzFG3R4wYwX///adhRCK1pJrgjXRAJUDpAwsMVPr+o6OVWQcPHVImIrK2VmYZzJULvL3h2TNo0SLhM5UooXy96vp15XvhwvF3FXSAU9x3IQxPSghFWnF2dsbNzY3Vq1fz6NEjxowZw+xXJ1MRRk3uDLxVFeL/mbp1A2dnqFsXvvsOzp6FH3+Ev/9WpiYGpYugenV4tykBrOOuJ0TakBJCkZYmT55M5syZAZg3bx4XL17UOCKRUrJQ0VsF8urywm8SGakkAePGQZ8+73O9Uu/6YCFeS1YhFOlh3LhxjIybfa1evXrs2LFDxqJoSFYtNCgX4Bhvm3zo/VgDzsChNLyGsFSyCqFIL8+fP6dEiRLcuHEDgM2bN9O0aVONo7JcsmqhQQ0jbRMBUJZGHpbG1xCWKHEJYYcOHaSEUKSZTJkyMW3aNHV7wIABvHjxQsOIREpIMpAiTYAOKJ/e04I14AY0TqPnF5YquRLCxYsXy21bkaZat25NjRo1APjnn38SVBoI4yTJQIr9CuTA8AmBddzz/mrg5xWWTkoIhVZ0Oh2zZ89Wk85x48apy2IL4yTJQIrlAHYD9hgqIYiOBr0+c9zzZjfIcwoBSgmhh4eHlBAKzZQrV46ePXsCEB4ezvDhwzWOSLyJJAOp4ggcxhB3CKKjlSmM+/atQEzM54YITgjV7NmzmT9/PiAlhEI748ePx8HBAYAlS5bg5/f2hd+ENiQZSDVH4BLQLm47tUmBcv7GjbaULAlz5x5g2DAZOCgMx8fHhwEDBqjbsgqh0EquXLkYNWoU8PJuVQoK2IQGJBl4J9mBVYAvSjkgKJM5vm5Qlo6Xkz06A75kz76NsDAlMZg6dSpLlixJu3CFxZBVCIWx6d27t3pX6ujRo/zxxx8aRySSI/MMGMQFYDVwAjhFwmWPs6FMMVwFpWLg5YRCCxYs4LvvvgOUW7m7d+9WR+AKkVo3b96kSpUqauVAhw4dWLVqlVQOCM1t376dxo2VaqlPPvmEv//+G3t7e42jsgwyz0C6KgVMQBkI+BglGXgY9/1x3P4JJJ5Z8Ntvv6VP3FSFUVFRfPXVV1y7di29ghZmJDw8HFdXVykhFEapUaNGajJw69YtpkyZonFEIjG5M6Cx6OhoXF1d2bFjBwCff/45R48eVQfdCPE20dHRNGvWjO3btwNKCeHx48elckAYlcuXL1O6dGmio6PJmDEjly9fpkCBAlqHZfbkzoCJsLGxYe3atZQsWRKAixcv0r59e6KjozWOTJiC+EFZ8YmAlBAKY1W8eHH69u0LQEREBEOGDNE4IvEqSQaMgIODA1u2bCF7dmWugb/++otBgwZpHJUwBYlLCDdu3CglhMJojRgxgly5cgGwdu1aDh2StViMhSQDRqJw4cJs3LgRGxul6mD27NksXLhQ46iEMUuuhLBWrVraBSTEW3zwwQdMmDBB3fbw8CAmJkbDiEQ8SQaMSM2aNVmwYIG6/cMPP7B3714NIxLGSkoIhanq3r075cqVA+Ds2bNSVm0kJBkwMj169GDgwIGAMjCsdevWXLlyReOohDGRVQiFKbO2tmb27Nnq9k8//URoaKiGEQmQZMAoTZ48mSZNmgAQEhJC06ZNCQkJ0TgqYQykhFCYgxo1atC2bVsAHjx4wLhx4zSOSEgyYISsra1ZvXo1pUuXBuDKlSu0adOGqKgojSMTWopfhdDf3x9Qxpls2rRJViEUJmnKlCnq/93Zs2dz+fJljSOybJIMGKls2bKxZcsWdeTtnj17ZF5vC/a6EsL4/x9CmJpPP/1ULS+Mjo5Wu0eFNiQZMGIFCxZk06ZN2NnZAfDbb78xb948jaMSWpASQmGOhgwZwieffAKAr6+vmuyK9CfJgJGrVq0anp6e6raHhwc7d+7UMCKR3qSEUJgre3v7BFMT9+/fX7pDNSLJgAno0qULQ4cOBSA2NpY2bdpw6dIljaMS6SFxCeHPP/8sJYTCrLRv355q1aoBypTFcvdTG7I2gYmIjY2lVatWeHt7A1CoUCFOnjxJjhw5tA1MpBlZhVBYCj8/P5ycnABlRtarV6/KeBgDkbUJzIyVlRUrVqxQJ+u4du0arVq1IjIyUtvARJqQEkJhSSpVqkS3bt0ACA0NZcSIERpHZHkkGTAhWbJkYfPmzXz88ccAHDhwgO+//14qDMyMlBAKSzRx4kSyZMkCwO+//865c+e0DcjCSDJgYvLnz4+Pj4/6xrBo0SJmzpypcVTCUKSEUFiqjz/+WL0joNfr6devn3zQSUeSDJigypUrJ5jPe9CgQWzdulXDiIShSAmhsGQeHh4ULlwYUO58btiwQeOILIckAyaqffv2jBw5ElCy6A4dOhAYGKhxVOJ9JC4h9PT0lBJCYVEyZMjAjBkz1O1Bgwbx/PlzDSOyHJIMmLBRo0ap83s/efIEV1dX7t+/r3FU4l0kV0LYtWtXjaMSIv01bdqUevXqARAUFMT06dM1jsgySGmhiXv27Bk1a9bEz88PgC+++IK9e/eSIUMGjSMTKSUlhEIkdOHCBcqWLUtMTAyZM2fm8uXL6kyFInWktNBCZM6cGR8fH/LmzQvA0aNH6dWrlwy8MRFSQihEUqVKleL7778HlA888ZOuibQjyYAZyJs3L5s3byZTpkwALF++nMmTJ2sclXibxCWEhQoVkhJCIeKMHj2a7NmzA7Bq1SqOHj2qcUTmTZIBM1GxYkVWrFihbg8bNoxNmzZpGJF4k+RKCLdt2yYlhELEyZ49O+PGjVO3PTw8iI2N1TAi8ybJgBlp1aoV48ePV7c7derE2bNnNYxIvI6UEArxdr169aJ06dKAMmXx8uXLNY7IfEkyYGZ++uknOnbsCCh9bc2aNVP7o4Vx2Lx5s5QQCpECNjY2zJo1S90eOnQoYWFh2gVkxiQZMDM6nQ4vLy+qVq0KwK1bt2jRooXU6hqJ06dP06FDBykhFCKF6tatS8uWLQG4d+8eEydO1Dgi8ySlhWbq3r17ODk5cfPmTUCZpGj16tUySl1DiUsI5W8iRMpcu3aNkiVLEhkZiZ2dHRcvXlRnKhRvJqWFFu6jjz5iy5Yt2NvbA/DHH38kGE8g0lfiEsIvvviCJUuWSCIgRAoUKlSIgQMHAhAZGcmgQYM0jsj8SDJgxsqWLZvgk+fIkSNZt26dxlFZnuRKCL29vaWEUIhUGDZsGHny5AHA29ub3bt3axyReZFkwMw1a9YswZwDXbt25dSpUxpGZFniV1+TEkIh3k/WrFn55Zdf1O1+/foRHR2tYUTmRZIBCzBo0CC+/vprACIiImjevDm3b9/WNigLMWfOHObNmwdICaEQ76tTp05UrlwZUKYsXrhwocYRmQ9JBiyATqdjwYIFVK9eHYA7d+7QrFkznj59qnFk5m3z5s30799f3ZYSQiHej5WVFbNnz1a3R4wYwaNHjzSMyHxIMmAhMmTIwIYNG/jss88AOHPmDF27dpUZvdKIlBAKkTaqVq1K586dAQgJCWH06NHaBmQmJBmwILly5WLLli1kzZoVgA0bNjBq1CiNozI/N2/epGnTpjx79gxQSgjHjh2rcVRCmI9JkyaplVK//fYbgYGBGkdk+iQZsDClSpVi7dq1WFkpf/rx48ezatUqjaMyH1JCKETay5cvHz/99BMAMTEx9OvXT1ZqfU+SDFigRo0aMWPGDHW7R48eHDt2TMOIzEN0dDTt27eXEkIh0sGAAQMoWLAgAHv27GHz5s3aBmTiJBmwUH379qVXr14AvHjxghYtWhAUFKRxVKYrvoRw27ZtAHzwwQf4+vpKCaEQaSRjxoxMnz5d3R4wYAAvXrzQMCLTJsmAhdLpdMydO1cd3X7//n2aNWvGkydPtA3MRCUuIdy0aRMlSpTQOCohzFvLli3VNuzatWsJFjUSqSPJgAWztbVlw4YNFClSBAB/f386duxITEyMxpGZFikhFEIbOp2OWbNmJRgDJau0vhtJBixc9uzZ2bp1Kw4ODoDyxhY/MEe83ZkzZ6SEUAgNlS1bVu3yfPLkibRf70hWLRQA7Nq1i0aNGql3BZYsWaLOWiiSJ6sQCmEcHjx4QLFixXj8+DEAJ0+exMnJSdugjISsWihSpV69esyZM0fd7tWrF4cOHdIwIuMmJYRCGI9cuXIlmHyob9++UmqYSpIMCNX333/PDz/8AEBUVBQtW7bk2rVrGkdlfKSEUAjj8/3336uDdo8fP87q1as1jsi0SDIgEpg5cyb169cH4NGjRzRt2pSwsDCNozIeUkIohHGytbVNUE0wZMgQqY5KBUkGRAI2NjasXbtWXVnv4sWLtG/fXpYKjSMlhEIYrwYNGuDq6gpAcHBwguXbxZtJMiCS+OCDD9i6dSvZs2cHYPv27QwePFjjqLQnJYRCGL/p06dja2sLwNSpU7l+/bq2AZkISQZEsooUKcKGDRuwsbEBYNasWfz+++8aR6WdxCWEw4cPlxJCIYxQsWLF8PDwAJTZVeWDTMpIaaF4Iy8vL3r27AkoXQg7d+6kdu3aGkeVvpIrIVy1apU60YkQwriEhoZSrFgx7t+/D8C+ffss9i6elBYKg3B3d1dvjUdHR9OqVSuuXr2qcVTp53UlhJIICGG8HBwcmDhxorrt4eEhM6u+hbRo4q2mTp1K48aNAQgJCcHV1ZWQkBCNo0p7UkIohOn6+uuvqVChAqBMte7l5aVxRMZNkgHxVtbW1qxZs4ZSpUoBcOXKFdq2bUtUVJTGkaUdKSEUwrRZW1sze/ZsdXv48OEW8SHmXUkyIFIkW7ZsbNmyhZw5cwKwe/du+vXrp21QaejVEkIbGxs2btwoJYRCmBgXFxfat28PKPOmjB07VuOIjJckAyLFPvvsMzZt2oSdnR0A8+fPV98wzUlyJYSWNmhSCHMxefJkMmXKBMDcuXO5dOmSxhEZJ0kGRKq4uLgkKDH08PBg586dGkZkWMmVEMqCTUKYrgIFCvDjjz8Cyjig/v37y7oFyZBkQKRa165d1RdXTEwMbdu25e+//9Y4qvd369YtmjZtyrNnzwClhFBuKwph+gYPHkz+/PkB2LFjhzoWSLwkyYB4JxMnTqR58+aAUtPr6urKo0ePNI7q3cWXEAYHBwNSQiiEOcmcOTNTp05VtwcMGEBkZKSGERkfaenEO7GysmLlypWULVsWgH///ZdWrVqZ5AssvoTw/PnzgJQQCmGO2rZti4uLC6BURM2dO1fjiIyLJAPinWXJkoXNmzfz0UcfAXDgwAF++OEHk+qPkxJCISyDTqdj9uzZ6HQ6AMaMGaPOUCgkGRDvqUCBAnh7e5MhQwZAGXn/6jKixk5KCIWwHBUqVKBHjx6AMk3vzz//rHFExkPWJhAGsXr1ajp27AgoXQibN2+mSZMmGkf1Zlu2bKF58+bqnYwlS5ZI5YAQZu7evXsUK1aMsLAwdDodp0+fpnz58lqHlWZkbQKRrtzc3BgxYgQAsbGxdOjQgcDAQI2jer0zZ87Qvn17KSEUwsJ89NFHjBw5ElC6CT08PEyqazOtSDIgDGb06NG0bt0aUEbnN23a1Cj75BKXELZr105KCIWwIH369KFo0aIAHDp0iD///FPjiLQnyYAwGCsrK5YtW0bFihUBuH79Ol999RUvXrzQOLKXkishXLp0qZQQCmFB7OzsmDFjhro9ePBg9cOBpZIWUBhU5syZ8fHxIW/evAAcOXKEXr16GcVtOCkhFELEa9KkCQ0aNADgxo0bTJs2TeOItCXJgDC4fPny4ePjo84Hvnz5cqZMmaJxVNC/f38pIRRCAEqp4cyZM7G2tgbgl19+4ebNmxpHpR1JBkSaqFSpEsuWLVO3hw0bhre3t2bxzJkzR51kREoIhRAAJUuW5IcffgDg+fPn/Pjjjzx9+pR58+axYMECYmJiNI4w/UhpoUhT48aNU0fu2tvbc/jwYcqVK5euMUgJoRDidUJCQihatKg6nXquXLl48OABAL6+vjRu3FjL8N6blBYKo/Dzzz/ToUMHAJ4+fUqzZs24c+cOv/76KyVLlkwwiMcQYmNj6d+/P25ubly/fl1KCIUQb/Thhx/Ss2dPdTs+EQC4fPmyFiFpQu4MiDT3/PlzateuzYkTJ4CEmbeNjQ2PHz/G3t7eINfavHmzuoBSzpw5AXj48CGglBCuXr1aKgeEEIAyz8CQIUNeO3hw1KhRjB49On2DMjC5MyCMRqZMmfD29lYrDF7NvKOjozl+/LjBrrVnzx7154cPH6qJgJQQCiESu3379hurCB4/fpx+wWhMWkaRLh4/fqyO2k3s0KFDBrvO/v37k91ft25ddf0EIYQAyJMnDw0bNnzt8dDQ0HSMRluSDIg0Fxsbi6ur62vLdt6eDOiBMOBh3Pfke7YePXqEv79/ssfGjRvHggULUhqyEMICWFtbs3XrVqZPn57sfCO3bt1K5lEpa49MjSQDIs1ZWVkRFRX12uNHjx5N5ngg8BNQF/gAcAByxX3/IG7/T3HnKQ4ePPjGOAICAlIbuhDCzFlbWzNgwAD8/f2pXr16gmMXL16M+yn17ZGpkWRApIv9+/fTvXv3ZG/VR0REvHJ73xdwARyBqcA+lOz7VWFx+6fGnecCbEswr8Gr7Ozs6Ny5M+PHjzfAbyKEMEdFixZl//79TJ06FZ1OB8BXX2XgXdsjUyPVBCJdPXjwgIULFzJ//nzu3Lmj7h89ug+jRj0E1qDkqLGpeFbl/PXrM9Cr1wtCQpS9n3zyCd999x3u7u7kzp3bcL+EEMKsXb58lMjIb3F0DOBd2yNwA34FsqdBhCmX0vdvSQaEJiIjI1m0aBGjRo2iVKkYdu+2wto6BHj3Gb9iYnQ8fKinU6eP+PbbeTRv3hwbGxvDBS2EsAD+QD3gEe/THoE1kAPYjXLHQBuSDAgT4Q9UB57yfi88hV5vjU5nDxxGyxegEMIUGbY9UhICbdsjmWdAmIBHKBm4oV54oNPFxD3fl8B/BnlOIYQlMHx7pDyPabRHkgwIDfXh/W/FJScm7nn7GPh5hRDmy7LbI0kGhEZ8Wbp0DTpdDDodSb4GDUp49pw5kCMHREe/3BcYCG3aQK5ckCEDFCwI338ffzQGWI0pjuoVQqQ3X5TByzH4+0O3bvDZZ5AxI2TJAhUqwJQp8F+iD/ebN4ONDTx4AE+fQvv2ULw4ZM0K9vZQqhSMHw9Pnxp/eySjq4RGJgE6QM+SJZB4NeG4mYtVGzZA8+bKCw9g3z5o0gSqV4cFCyBnTrhxA86effVR1nHXMe1Vx4QQaW0SYIWnZyzff6+8oQ8eDJ9/DlFR4OentDPHjsGmTS8ftWED1KihfCB5/Bj0ehgwQEkkrKzg4EEYOxb274fdu427PZIBhEIDgYAjS5cqGfipU1Cp0uvPvndPSQ42b1YSgGfPoEgRJVvfskW5k/D265UyWPRCCHOitEfHjikfLurVA29v5W7jqyIj4a+/oFkzZTsqCnLnVj759+79+mf/8UflrsK//0KhQvHXS7/2SAYQCiO2mtTclNq0SblV9+WXyvaff8KdO0rm/vZEwCbuekIIkRylPZo4UWlPfv89aSIAYGf3MhEA2LMHQkOhZcs3P3uuXMp35a6m8bZHkgwIDZwAXnb+x8QoYwFe/XrVhg3g6vryBRo/63BMDLi4KC/SDz+EDh0gODjxtWLirieEEMk5QUxMNHv3QsWKkD9/yh61YQM4Oyft0tTrlTYsLEy5kzB9utI2FSgAxtweSTIg0pke8Euwp2pVsLVN+BWfEDx6pPS3tWr18vzbt5XvrVpBtWqwYwf88gvs2gU1ayrdCAmvdwpzWUxECGFISnv08KHSbnz2WcoeFROjdCW82i7FW7tWacMcHKBRI+Vr+fJXr2ec7ZEMIBTpLJzEc3svXw4lSyY8K36goI+P8sn/1VVGY+NmBm3XDiZPVn6uXRs+/hhatIDVq8Hd/dVnCwOeAFkN9lsIIcxB0vYoJQ4cgIcP4auvkh5r0EAZBxUergw4nDxZ+VCzaZMyqNBY2yNJBkQ6i0yyp2TJ1w8gXL9eyawzZ365L0cO5XuDBgnPbdBA6fM7cyZl1xVCWDqlXciZU2lj/ve/lD1q/XqlS6FgwaTHPvzwZXtWuzYULqyUHPr4vDq+wPjaI+kmEOnMLsVnhoYqg3QS34orU+bNj7NK9n91yq8rhLAUSrtgbQ1168Lp03Dr1psfERurfMpProsgOZUrK9+vXEl6XWMiyYBIZ1mBlJWnxpcNurom3N+ypbJ/+/aE+7dvVwbvVK2a+JmyAVneMV4hhPl62R4NG6a0Hz17KmWEiUVFKW3S0aNw927Kk4F9+5TvRYrE7zHO9ki6CUQ60wGVgL1vPXP9eqXmN2uirrUSJZS63vnzlWONGilZ988/Q/ny0LZt4us5xX0XQohXvWyPnJ3ht9+UWUwrVoTvvlNmEIyKUiYz+/13KF0aPv1U+V6sWMJnWrgQDh2C+vWVioSnT5XtX3+FL75QJk0z5vZIkgGhgSrAQV4tL0zs6VOlSuC335I/PmsWfPIJeHkpL7acOZV+uYkTlQGHL1nHXU8IIZLzsj3q2VO5rT9zpjLw7+5dpTKgWDFwc4MfflAShe7dkz6LoyNs3arcYXj4UBkEXbQo/PSTMiuhMijaeNsjmYFQaECZ8etN1q2Djh2V2QezZzfE9WQGQiFEct7eHsU7eRKqVAF/f+XN/92vZ3wzEEoyIDTiAhwDYtPsCnq9NTqdM3Aoza4hhDB9MTHO6HQnsLJKy/p/ayD92yOZjlgYuWGkZSIAoNPF8Pjxd2l6DSGEaTt27Bi9el1P40QAlNkHh6XxNd6dJANCI02ADijZsuFFR8OqVVC8eH+2bNmSJtcQQpiuqKgoRowYgYuLC4sX32X16qRToRuONeCGsa5YCJIMCE39CuTA0AlBbKwVISE6+vSB+/fv06xZM3r16sWTJ08Meh0hhGm6dOkSzs7OjB8/nti4KU2XL69EWrRHyvPlQGnvjJckA0JDOYDdgD2GewFaY2WVBSurvVSv/nKJMU9PT8qVK8exY8cMdB0hhKnR6/XMnTuXChUqcPr0aQBsbGwYP348W7cew8ZmH4Zuj5Tn2w2890joNCXJgNCYI3AYw2Tk8Rn4YXLkqIW3tzeenp7Y29sD8O+//+Li4sKIESOIiop6z2sJIUxJcHAwDRs2pE+fPkRERABQokQJjh8/zvDhw7GxsSGt2qOUVitoSZIBYQQcgUtAu7jt1L4I489vH/c8ygtPp9Ph7u7O+fPncXZ2BiA2Npbx48fj7OzM33///d6RCyGM359//knp0qXZuXOnuq9Pnz6cPn2aihUrJjo7bdojYyfJgDAS2YFVgC9K+Q0oc2K9bqYuHS/nzHKOe9xKkrsVV7hwYQ4ePMj48ePjsn84ffo05cuXZ+7cuaSgulYIYYJCQ0Pp3Lkzbdu2JSQkBIC8efOyY8cO5syZQ+ZXV0BLIO3aI2Ml8wwII3UBWA2cQFn/+9VlRrOhTOlZBWWEbson8Dh9+jSdOnVKcFegQYMGLF68mLx58xogbiGEMdi/fz9du3blxo0b6r42bdqwYMECsqd6JrO0aY/Sg0w6JMyIHmX970iU1b6y8D5zez979oyhQ4fy668vR/dmz56dhQsX0rp16/eMVQihpRcvXjB8+HBmzJih3vVzcHBg3rx5uLm5odO977oAhm2P0ppMOiTMiA5ldbEccd/f74WXOXNm5syZw44dO9S7Af/99x9t2rShS5cuhIaGvm/AQggN+Pv74+TkxPTp09VEoFatWvj7+9OxY0cDJAJg6PbIWEgyICxW/fr1CQgIoE2bNuq+FStWUKZMGQ4cOKBhZEKI1IiJiWHatGk4OTkREBAAgJ2dHdOmTWPPnj0UKFBA4wiNnyQDwqJlz56dtWvXsmLFCvUW2o0bN6hduzaDBw/mxYsXGkcohHiToKAg6taty+DBg4mMjASgTJky+Pn5MXDgQKys5G0uJeRfSVg8nU5Hp06dCAgIoFatWoAyOUn8Jw1/f39tAxRCJKHX65PcydPpdAwePJiTJ0/i+O7LClokSQaEiFOgQAH27NnDtGnTsLOzAyAgIAAnJyemTZumTlsqhNDWo0ePaNeuHV26dCEsTBnZX6BAAfbt28eUKVPIkCGDxhGaHkkGhHiFlZUVAwcOxM/PT/1kERkZyeDBg6lbty5BQUEaRyiEZduxYweOjo78+eef6r4uXbrg7+9PzZo1NYzMtEkyIEQyHB0dOXXqFIMHD1ZHIO/fv58yZcqwcuVKmahIiHT27Nkz+vTpQ8OGDblz5w6gjPn5888/WbZsGQ4ODhpHaNokGRDiNTJkyMCUKVPYt2+fOho5LCyMzp07065dOx49eqRxhEJYBj8/PypWrMjcuXPVfQ0aNCAgIEDmBjEQSQaEeIuaNWvi7+9P586d1X1//vknjo6OCeY6F0IYVnR0dJK1RDJlysTcuXPZvn27zBpqQJIMCJECDg4OLF++nHXr1qlTmd65c4cGDRrQp08fnj17pnGEQpiXf//9lxo1ajBixAiio6MBqFSpEmfOnKF3794GmkBIxJNkQIhUaNOmDQEBATRo0EDdN3fuXCpWrKiujy6EeHd6vR5PT0/Kli3LsWPHAGVg74gRIzh69CglSpTQOELzJMmAEKmUN29etm/fzty5c8mYMSMAf//9N1WrVmX8+PHqpxghROrcu3eP5s2b06tXL54+fQooq44eOXKEsWPHYmtrq3GE5kuSASHegU6no3fv3pw9e5ZKlSoBSv/miBEjqFGjBv/++6/GEQphWjZv3oyjoyNbtmxR9/Xq1Ytz585RtWpVDSOzDJIMCPEeSpQowdGjRxkxYoQ67emxY8coW7Ysnp6eUoIoxFuEh4fj7u5O8+bNefDgAQC5c+dmy5YtLFy4kCxZsmgcoWWQZECI92Rra8vYsWM5fPgwhQsXBuDp06f06tWL5s2bc+/ePY0jFMI4HT16lHLlyrFo0SJ1X7NmzQgICMDV1VXDyCyPJANCGIizszPnzp2jV69e6r4tW7bg6OjI5s2bNYxMCOMSGRnJzz//TPXq1bl27RoA9vb2eHl54e3tTe7cuTWO0PJIMiCEAWXJkoWFCxeyefNmtUF78OABzZs3p2fPnjx58kTjCIXQ1qVLl3B2dmbChAnqeh9ffPEF58+fp0ePHlIyqBFJBoRIA02bNiUgIIBmzZqp+7y8vChbtixHjx7VMDIhtBEbG8uvv/5KhQoVOHPmDAA2NjZMmDCBgwcPql1sQhuSDAiRRnLnzo23tzdeXl7Y29sDcO3aNapXr87PP/9MVFSUxhEKkT5u375No0aN6Nu3LxEREQCULFmSEydO8NNPP2Ftba1xhEKSASHSkE6no0ePHpw/f54vvvgCUD4hTZgwAWdnZy5duqRxhEKkrXXr1iWZurtv376cPn2aChUqaBiZeJUkA0Kkg8KFC3PgwAHGjx+PjY0NgNoYzp07V+07FcJcPH78WF3UKyQkBFAm7Nq5cyezZ88mU6ZMGkcoXiXJgBDpxMbGhuHDh3P8+HF1StWIiAj69OlDo0aNuH37tsYRCmEYry73Ha9t27YEBARQr149DSMTryPJgBDprGLFipw5c4Y+ffqo+3bu3ImjoyN//vmnhpEJ8X4iIiIYNGgQderU4ebNm4CyyNeqVav4448/1EW+hPGRZEAIDWTKlIk5c+awY8cOdRnWkJAQ2rZtS+fOnXn8+LG2AQqRSv7+/lSuXJnp06erM2/Wrl0bf39/3NzcpGTQyEkyIISG6tevT0BAAG3atFH3rVy5kjJlyrB//37tAhMihWJiYpg6dSpOTk4EBAQAYGdnx/Tp09m9ezcFChTQOEKREpIMCKGx7Nmzs3btWlauXImDgwMAN2/epE6dOgwaNEgtxRLC2AQFBVGnTh2GDBlCZGQkAGXKlOH06dMMGDBAXa9DGD/5SwlhBHQ6HR07dsTf359atWoByrru06dPp3Llyvj7+2sboBCv0Ov1LF++nDJlynDw4EFA+T88ZMgQTp48SenSpTWOUKSWJANCGJECBQqwZ88epk+fjp2dHQABAQE4OTkxdepUYmJiNI5QWLpHjx7Rtm1bunbtSlhYGACffvop+/fvZ/LkyWTIkEHjCMW7kGRACCNjZWXFgAED8PPzo0yZMoCysMuQIUOoW7cuQUFBGkcoLNVff/2Fo6Mj69evV/d17dqV8+fPU6NGDQ0jE+9LkgEhjJSjoyMnT55kyJAh6kjsAwcOUKZMGZYvX66O2BYirT179owffviBRo0acefOHUAZ67J+/XqWLl2qjnURpkuSASGMWIYMGZg8eTL79u1TR2WHhYXRtWtX2rZty6NHjzSOUJg7Pz8/KlSowLx589R9DRs2JDAwkFatWmkYmTAkSQaEMAE1a9bE39+fLl26qPvWr1+Po6MjO3bs0DAyYa6io6MZN24czs7OXL58GVDmx5g3bx7btm0jT548GkcoDEmSASFMhIODA8uWLePPP/9UZ3K7c+cODRs25IcffuDZs2caRyjMxT///EP16tUZOXIk0dHRADg5OXH27Fm+//57mUDIDEkyIISJad26NYGBgTRs2FDdN2/ePCpUqICfn5+GkQlTp9fr+f333ylbtizHjx8HwNrampEjR3LkyBGKFy+ucYQirUgyIIQJypMnD9u2bWPu3Lnq6m+XL1/G2dmZcePGqZ/mhEipe/fu0axZM7755hv1LlORIkU4fPgwY8aMwdbWVuMIRVqSZEAIE6XT6ejduzdnzpyhUqVKgNLPO3LkSKpXr84///yjcYTCVPj4+ODo6MjWrVvVfd988w1nz56latWqGkYm0oskA0KYuBIlSnD06FFGjBihTv96/PhxypUrh6enp5QgitcKDw/H3d2dFi1a8ODBAwBy587Nli1bWLBgAVmyZNE4QpFeJBkQwgzY2toyduxYjhw5QuHChQF4+vQpvXr1olmzZty7d0/jCIWxOXr0KOXKlWPRokXqvubNmxMYGIirq6uGkQktSDIghBmpWrUq586do1evXuq+rVu34ujoyObNmzWMTBiLyMhIhg8fTvXq1bl27RoAWbJkYdGiRWzatIlcuXJpHKHQgiQDQpiZLFmysHDhQrZs2ULu3LkBePDgAc2bN8fd3Z3w8HCNIxRauXTpEs7OzkycOJHY2FgAqlWrxvnz5+nevbuUDFowSQaEMFOurq4EBATQrFkzdd+iRYsoV64cR48e1TAykd5iY2OZM2cOFSpU4MyZMwDY2NgwceJEDhw4QKFChTSOUGhNkgEhzFju3Lnx9vZm0aJF6mCwa9euUb16dYYPH66uQS/M1+3bt2nYsCEeHh5EREQAULJkSU6cOMGwYcOwtrbWOEJhDCQZEMLM6XQ6unfvzvnz5/niiy8A5ZPixIkTcXZ25tKlSxpHKNLK2rVrcXR0ZNeuXeo+Dw8PTp8+TYUKFTSMTBgbSQaEsBCFChXi4MGDTJw4ERsbGwDOnDlDhQoVmDNnjtqHLEzf48eP6dSpE+3btyckJASAvHnzsnPnTmbNmqVOVCVEPEkGhLAg1tbWDBs2jBMnTlCyZEkAIiIi8PDwoGHDhty+fVvjCMX72rdvH2XKlGHVqlXqvnbt2hEQEEC9evU0jEwYM0kGhLBAFSpU4PTp0/Tt21fdt2vXLhwdHVm3bp2GkYl3FRERwcCBA6lTpw43b94ElMWtVq1axR9//KEubiVEciQZEMJCZcqUidmzZ7Nz507y5s0LQEhICO3ataNTp048fvxY2wBFip0/fx4nJydmzJih7qtTpw4BAQG4ublpGJkwFZIMCGHh6tWrR0BAAO3atVP3rVq1ijJlyrBv3z4NIxNvExMTw9SpU6lcuTKBgYEAZMiQgRkzZrBr1y7y58+vcYTCVEgyIIQge/bsrFmzhlWrVuHg4ADAzZs3qVu3LoMGDVJL0oTxuH79OnXq1GHIkCFqiWjZsmXx8/Ojf//+6joVQqSE/G8RQgBKCaKbmxv+/v7Url0bUNa3nz59Ok5OTpw/f17jCAUof5Ply5dTpkwZDh48CCh/ux9//JETJ05QunRpjSMUpkiSASFEAgUKFGD37t3MmDGDDBkyABAYGEjlypWZOnUqMTExGkdouR4+fEibNm3o2rWrOq30p59+yoEDB/jll1/Uv5cQqSXJgBAiCSsrK/r374+fnx9ly5YFlAVuhgwZQp06dbh+/bq2AVqgv/76C0dHRzZs2KDu69q1K/7+/lSvXl3DyIQ5kGRACPFapUuX5sSJEwwZMkRdxObgwYOUKVOG5cuXo9frNY7Q/D179ozevXvTqFEj7t69C0COHDlYv349S5cuJVu2bBpHKMyBJANCiDfKkCEDkydPZv/+/Xz66acAhIeH07VrV9q0acPDhw81jtB8nTp1ivLlyzN//nx1X6NGjQgICKBVq1YaRibMjSQDQogUqVGjBv7+/nTt2lXdt2HDBhwdHfnrr780jMz8REdHM3bsWJydnbly5QqgzAsxf/58fH19yZMnj8YRCnMjyYAQIsWyZcvG0qVLWb9+PTly5ADg7t27NGrUiN69e/Ps2TONIzR9V69excXFhVGjRqmDNZ2cnDh79izfffed2l0jhCFJMiCESLVWrVoREBBAw4YN1X3z58+nfPnynDp1SsPITJder+f333+nXLlynDhxAlDWkhg1ahRHjhyhePHiGkcozJkkA0KId5InTx62bdvG/Pnz1VXwrly5grOzM2PHjiU6OlrjCE3HvXv3aNasGd988416d6VIkSIcOXKE0aNHY2trq3GEwtxJMiCEeGc6nY7vvvuOs2fP4uTkBChT5I4aNQoXFxeuXr2qcYTGz8fHh9KlS7N161Z137fffsu5c+eoUqWKhpEJSyLJgBDivRUvXpwjR44watQorK2tAThx4gTlypXj999/lxLEZISHh+Pu7k6LFi3UioyPPvqIrVu38ttvv2Fvb69xhMKSSDIghDAIW1tbRo8ezZEjRyhSpAig1Mh/8803NGvWjHv37mkcofE4cuQIZcuWZdGiReq+Fi1aEBAQQJMmTTSMTFgqSQaEEAZVpUoVzp07xzfffKPu27p1K6VLl8bHx0fDyLQXGRnJ8OHDqVGjBv/73/8AyJIlC4sXL2bjxo3kypVL4wiFpZJkQAhhcPb29ixYsICtW7fy0UcfAcq8+i1atMDd3V2dV9+SXLx4kapVqzJx4kRiY2MBqFatGv7+/nTr1k1KBoWmJBkQQqSZJk2aEBAQQIsWLdR9ixYtomzZshw5ckS7wNJRbGwsc+bMoWLFipw9exZQulQmTZrEgQMH+OyzzzSOUAhJBoQQaSxXrlxs3LiRxYsXkyVLFgD+97//UaNGDYYPH05kZKTGEaadW7du0aBBAzw8PIiIiADg888/58SJEwwdOlQdbCmE1iQZEEKkOZ1OR7du3Th//jzVqlUDlE/MEydOpGrVqly8eFHjCA1v7dq1ODo6snv3bnVfv3798PPzo3z58hpGJkRSkgwIIdJNoUKFOHDgAJMmTVIn0jl79iwVK1Zkzpw5al+6KXv8+DEdO3akffv2PH78GIB8+fKxa9cuZs6cqU7QJIQxkWRACJGurK2tGTp0KCdOnODzzz8HICIiAg8PDxo2bMjt27c1jvDd7d27F0dHR1avXq3ua9++PQEBAXz55ZcaRibEm0kyIITQRPny5fHz86Nfv37qvl27duHo6MjatWu1C+wdREREMGDAAOrWrcutW7cAcHBwYPXq1axZs4YPP/xQ4wiFeDNJBoQQmsmUKRMzZ85k165d5MuXD4CQkBDat29Px44d1dvsxuz8+fNUqlSJmTNnqvvq1KlDQEAAHTp00DAyIVJOkgEhhOa+/PJLAgICaN++vbpv9erVODo6snfvXnXf5cuXGTt2LEFBQeke4549e5g6dSpPnz4FlDUYpkyZgpOTExcuXAAgQ4YManKTP3/+dI9RiHel06dg0vCwsDAcHBwIDQ0lW7Zs6RGXEMJCrVmzhu+++47Q0FB1X//+/enbty9OTk48fPiQYsWKERgYmMrV/PRAOBAJ2AFZgZRN9HP+/HkqVqxITEwMTZo04ddff6Vr164cOnRIPads2bKsXLmS0qVLpyImIdJWSt+/JRkQQhidmzdv8vXXXye4K5AlSxaePHmibs+bN4/vv//+Lc8UCKwGTgB+QNgrx7IBlYAqgBvw+jfxhg0bsmPHDnU7Q4YMvHjxAlDKJn/88UdGjx5NhgwZUvT7CZFeJBkQQpi0+Jn7hg4dqr7xvip37tz8888/ZM2aNZlH+wKTgCOADRCDcmcgMR1gDUQD1YCfgMYJzti9ezf16tVLNsaCBQuyfPlyqlevnuLfS4j0lNL3bxkzIIQwSlZWVvTr14/169djZZW0qbp//z7Tpk1LtPcRyqd8V+BY3L5okk8EiNsfHffzMaAJ0BH4D1ASkiFDhiT7yJw5c3L69GlJBIRZkGRACGHUJk+e/NrJiCZPnkxwcHDclj/wObAubju1ExjFn78WKAkEsGzZMnU9gcQePnyYYAliIUyZJANCCKMVERHB8ePHX3v8xYsXtGvXDiURqI5yZyDmPa8aAzxCr3dh3rxv33jm/v373/NaQhgHSQaEEEYrY8aMzJs3j7Jly5I5c+Zkz3nw4G+gHvCU908E4sUAT/H1jSS5+YJy585NnTp1GDVqlIGuJ4S2ZAChEMIk6PV67ty5w5UrV7h69So7duzgn3/+Ydu2D8ib9zCGSwReio3VsX27Azt2dKZq1aoUK1aMokWL4uDgYPBrCZEWpJpACGEBfFm61JVu3V5/xr59UKvWy+05c2DMGLh3D27dgs8+S/5xa9bAyzmQfElcZSCEKUjp+7dNOsYkhBAGNgmlPFDPkiVQokTSM+LWQlJt2ADNm4PNK61fnz7g5pbwvKJF43+yjruOJAPCfEkyIIQwUYEo8wgoSpeGSpXe/Ih79+DwYUhcLVigAFSt+rpHxQCHgQtAqXeOVghjJgMIhRAmajWp/TyzaRNkyQKpX03YJu56QpgnSQaEECbqBC8nDIKYGIiOTvgVk2hM4YYN4OoKiWcN/uUXsLODzJnBxQU2b058rZi46wlhniQZEEKYID3KWgMvVa0KtrYJv15903/0CPbvh1atXu7LkAF69oTffoO9e8HLS0kgmjdXfk54vVO8fiZDIUybjBkQQpigcBIuOgTLl0PJkgnP0r2yKKGPj/Lpv2HDl/vy5IHff0/4mDZtoEoVGDoUvv761YGGYcATlNUOhTAvkgwIIUxQZJI9JUu+eQDh+vXQqJHSFfAmtrbQrp2SDFy9mjjBSHpdIcyBdBMIIUyQXarODg2FPXsSdhG8SfzsK0nXR0rddYUwFZIMCCFMUFYg5ROgbdmidBm4ur793KgoWLsWcuaEIkVePZINyJLKOIUwDdJNIIQwQTqgErBX3RMYqFQQJFa4sNJFUK8eZE3U3T9ggPLmX60afPwx3LwJv/4K587BkiVgbf3q9ZzivgthfiQZEEKYqCrAQeLLC183JfGsWbBjh1IxkFjp0rBwIaxeDWFhSrJQubJyfv36r55pHXc9IcyTrE0ghDBRgYDjW89atw46dlRmH8ye/X2vJzMQCtOS0vdvGTMghDBRpYFqvK0Za9tW6Qp490TAGnBBEgFhziQZEEKYsGFAbBpfIybuOkKYL0kGhBAmrAnQAeXTe1qwBtyQFQuFuZNkQAhh4n4FcmD4hMA67nl/NfDzCmF8JBkQQpi4HMBuwB7DJQTWcc+3G3ivUYdCmARJBoQQZsAROIxh7hDE3xE4TEqqFYQwB5IMCCHMhCNwCWgXt53apCD+/PZxzyOJgLAckgwIIcxIdmAV4As4x+2z4fUzB+p4Ofeac9zjViJdA8LSyAyEQggz1Dju6wKwGjgBnCLhssfZUKYYroJSMSDzCAjLJcmAEMKMlQImxP2sB56gLENsh7LokKw1IARIMiCEsBg6lNUOhRCJyZgBIYQQwsJJMiCEEEJYOEkGhBBCCAsnyYAQQghh4SQZEEIIISycJANCCCGEhZNkQAghhLBwkgwIIYQQFk6SASGEEMLCSTIghBBCWDhJBoQQQggLJ8mAEEIIYeEkGRBCCCEsnCQDQgghhIVL0RLGer0egLCwsDQNRgghhBCGE/++Hf8+/jopSgbCw8MByJ8//3uGJYQQQoj0Fh4ejoODw2uP6/RvSxeA2NhYgoODyZo1KzqdzqABCiGEECJt6PV6wsPDyZs3L1ZWrx8ZkKJkQAghhBDmSwYQCiGEEBZOkgEhhBDCwkkyIIQQQlg4SQaEEEIICyfJgBBCCGHhJBkQQgghLJwkA0IIIYSF+z+5W9Xwr0Rf6gAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<Figure size 100x100 with 0 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "let_to_num = {\"A\":1, \"B\":2, \"C\":3, \"D\":4, \"E\":5, \"F\":6}\n",
    "randVar = ['A', 'B', 'C', 'D', 'E', 'F']\n",
    "values  = [ 1,   2,   3,   4,   5,   6]\n",
    "\n",
    "ltn = dict( zip(randVar, values) )\n",
    "\n",
    "# graphical model construction\n",
    "G = nx.DiGraph()\n",
    "G.add_edges_from([('B/2', 'D/4')])\n",
    "G.add_edges_from([('A/1', 'D/4')])\n",
    "G.add_edges_from([('A/1', 'C/3')])\n",
    "G.add_edges_from([('D/4', 'F/6')])\n",
    "G.add_edges_from([('A/1', 'F/6')])\n",
    "G.add_edges_from([('C/3', 'E/5')])\n",
    "G.add_edges_from([('F/6', 'E/5')])\n",
    "\n",
    "#pos = nx.spring_layout(G, seed=3) # positions for all nodes, seed for reproducibility\n",
    "pos = nx.nx_agraph.graphviz_layout(G, prog=\"dot\")\n",
    "\n",
    "# nodes\n",
    "nx.draw_networkx_nodes(G, pos, node_color= 'yellow', node_size=400)\n",
    "\n",
    "# edges\n",
    "nx.draw_networkx_edges(G, pos, width=2)\n",
    "labels = nx.get_edge_attributes(G,'weight')\n",
    "nx.draw_networkx_edge_labels(G,pos,edge_labels=labels, font_size=12)\n",
    "\n",
    "# node labels\n",
    "nx.draw_networkx_labels(G, pos, font_size=12, font_family=\"sans-serif\")\n",
    "\n",
    "plt.figure(figsize=(1, 1))\n",
    "plt.show()\n",
    "\n",
    "class Struct:\n",
    "    def __init__(self, var=None, card=None, val=None):\n",
    "        self.var = var\n",
    "        self.card = card\n",
    "        self.val = val\n",
    "\n",
    "# input probab tables\n",
    "probab_A = Struct([let_to_num['A']], [2], [0.4, 0.6])\n",
    "probab_B = Struct([let_to_num['B']], [3], [0.3, 0.4, 0.3])\n",
    "probab_C = Struct([let_to_num['C'], let_to_num['A']], [3, 2], [0.3, 0.25, 0.45, 0.05, 0.2, 0.75])\n",
    "probab_D = Struct([let_to_num['D'], let_to_num['A'], let_to_num['B']], [2, 2, 3], [0.05, 0.95, 0.22, 0.78, 0.93, 0.07, 0.53, 0.47, 0.4, 0.59, 0.3, 0.7])\n",
    "probab_E = Struct([let_to_num['E'], let_to_num['C'], let_to_num['F']], [2, 3, 2], [0.02, 0.98, 0.3, 0.7, 0.13, 0.87, 0.32, 0.68, 0.92, 0.08, 0.58, 0.42])\n",
    "probab_F = Struct([let_to_num['F'], let_to_num['A'], let_to_num['D']], [2, 2, 2], [0.95, 0.05, 0.43, 0.57, 0.1, 0.9, 0.28, 0.72])\n",
    "probab_tables = [probab_A, probab_B, probab_C, probab_D, probab_E, probab_F]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "52559756",
   "metadata": {},
   "outputs": [],
   "source": [
    "def IndexToAssignment(I, D):\n",
    "    \"\"\"\n",
    "    I: list of indices to be converted.\n",
    "    D: list of caradinalities.\n",
    "    \n",
    "    returns: np.array() of assignments.\n",
    "    \"\"\"\n",
    "    D = np.array(D).reshape( (1, len(D)) ) # ensure D is a row vector\n",
    "    I = np.array( I ).reshape( (len(I), 1) )\n",
    "    \n",
    "    den = np.matlib.repmat( np.cumprod( np.concatenate( (np.array([1]), D[0][:-1]) ) )  , I.shape[0], 1)\n",
    "    num = np.matlib.repmat( I-1, 1,  D.shape[1])\n",
    "    AA = np.floor( num / den )\n",
    "    BB = np.matlib.repmat( D, len(I), 1 )\n",
    "    result = np.mod( AA, BB ) + 1     # whether or not you want to add 1 is a modeling option, ofcourse.\n",
    "    return np.array(result, dtype=int)\n",
    "\"\"\"\n",
    "# to check whether the function is working as expected or not.\n",
    "D = [2, 2]\n",
    "I = list( range(1, np.prod(D) + 1) )\n",
    "A = IndexToAssignment(I, D)\n",
    "print(A, AssignmentToIndex(A, D))\n",
    "\"\"\"\n",
    "def AssignmentToIndex(A, D):\n",
    "    \"\"\"\n",
    "    A: np.array() of assignments \n",
    "    D: list of caradinalities.\n",
    "    \n",
    "    returns: np.array() of indices.\n",
    "    \"\"\"\n",
    "    D = np.array(D).reshape( (1, len(D)) ) # ensure D is a row vector\n",
    "    \n",
    "    if 1 in A.shape:\n",
    "        x = np.cumprod( np.concatenate( (np.array([1]), D[0][:-1]) ) )\n",
    "        x = np.expand_dims(x, axis=0)\n",
    "        I = np.dot( x, (A - 1).T) + 1\n",
    "    else:\n",
    "        I = np.sum( np.matlib.repmat( np.cumprod( np.concatenate( (np.array([1]), D[0][:-1]) ) ), A.shape[0], 1 ) * (A - 1) , axis=1) + 1\n",
    "        \n",
    "    return I\n",
    "\n",
    "def product_factor(factor1, factor2):  \n",
    "    \"\"\"\n",
    "    Computes the product of two factors (fundamental computation to solve many downstream tasks).\n",
    "    \n",
    "    input:  factor1: Struct \n",
    "            factor2: Struct\n",
    "           \n",
    "    output: factor3: Struct: another factor that computes their product.\n",
    "    \"\"\"\n",
    "    # determining variables in the new factor.\n",
    "    factor3 = Struct()\n",
    "    factor3.var = sorted( list( set(factor1.var + factor2.var) ) ) # a sorted non-repeated union\n",
    "    #print(\"factor3 variables: \", factor3.var)\n",
    "    \n",
    "    # determining the cardinality of the new factor.\n",
    "    _, mapfactor1 = ismember(factor1.var, factor3.var)\n",
    "    _, mapfactor2 = ismember(factor2.var, factor3.var)\n",
    "    factor3.card = [0] * len(factor3.var)\n",
    "    #print(\"mapfactor1: \", mapfactor1) # factor 1 me jo elements hai vo factor 3 mei kaha milengey?\n",
    "    #print(\"mapfactor2: \", mapfactor2)\n",
    "    \n",
    "    for index, element in enumerate(mapfactor1):\n",
    "        factor3.card[element] = factor1.card[index]\n",
    "    for index, element in enumerate(mapfactor2):\n",
    "        factor3.card[element] = factor2.card[index]\n",
    "    \n",
    "    # determining the probability values associated with the new factor.\n",
    "    totalNumberofValues = math.prod(factor3.card)\n",
    "    factor3.val = [-1] * totalNumberofValues\n",
    "    assignments = IndexToAssignment( list( range(1, totalNumberofValues+1) ),  factor3.card ) \n",
    "    #print(assignments)\n",
    "    \n",
    "    indxfactor1 = AssignmentToIndex( assignments[:, mapfactor1], factor1.card ).ravel()  # acessing factor1.val elements\n",
    "    indxfactor2 = AssignmentToIndex( assignments[:, mapfactor2], factor2.card ).ravel()  # acessing factor2.val elements\n",
    "    #print(indxfactor1, indxfactor2)\n",
    "    #print(indxfactor1.ravel(), indxfactor2.ravel())\n",
    "    \n",
    "    \n",
    "    for index in range(totalNumberofValues):\n",
    "        factor3.val[index] = factor1.val[ indxfactor1[index]-1 ] * factor2.val[ indxfactor2[index]-1 ]\n",
    "        \n",
    "    return factor3\n",
    "\n",
    "def compute_joint_distribution(probab_tables):\n",
    "    \"\"\"\n",
    "    Computes the joint distribution given a list of probability factors as input.\n",
    "    \n",
    "    input:   probab_tables: List of Struct's\n",
    "    \n",
    "    output:  joint:         Struct\n",
    "    \"\"\"\n",
    "    joint = probab_tables[0]\n",
    "    for index, probab_table in enumerate(probab_tables):\n",
    "        if index >0:\n",
    "            joint = product_factor(joint, probab_table)\n",
    "            \n",
    "    return joint\n",
    "\n",
    "def sum_marginalization(factor1, variableOut):\n",
    "    \"\"\"\n",
    "    Marginalizes out variables in variableOut from factor1.\n",
    "    \n",
    "    input:  factor1: Struct  \n",
    "            variableOut: List: A list of variables to be marginalized out from factor1.var\n",
    "            \n",
    "    output: Struct\n",
    "    \"\"\"\n",
    "    # determining variables in the new factor.\n",
    "    factor3 = Struct()\n",
    "    factor3.var = sorted( list( set(factor1.var) - set(variableOut) ) )\n",
    "    \n",
    "    # determining the cardinality of the new factor.\n",
    "    mapfactor3 = []\n",
    "    for randomVar in factor3.var:\n",
    "        for index in range(len(factor1.var)):\n",
    "            if factor1.var[index] == randomVar:\n",
    "                mapfactor3.append(index)\n",
    "                break\n",
    "    factor3.card = [0] * len(factor3.var)\n",
    "    \n",
    "    for idx in range(len(factor3.card)):\n",
    "        factor3.card[idx] = factor1.card[ mapfactor3[idx] ]\n",
    "    \n",
    "    # yeh pata karo ki factor3 me jo variables hai vo factor1 me kaha milengey.\n",
    "    #_, mapfactor1 = ismember(factor1.var, factor3.var)\n",
    "    #for idx in mapB:\n",
    "    #    B.card.append(A.card[idx])\n",
    "    \n",
    "    totalNumberofValuesfactor1 = math.prod(factor1.card)\n",
    "    assignments = IndexToAssignment( list( range(1, totalNumberofValuesfactor1+1) ), factor1.card )\n",
    "    indxfactor3 = AssignmentToIndex( assignments[:, mapfactor3], factor3.card ).ravel()\n",
    "   \n",
    "    # determining the probability values associated with the new factor.\n",
    "    totalNumberofValuesfactor3 = math.prod(factor3.card)\n",
    "    factor3.val = [0] * totalNumberofValuesfactor3\n",
    "    \n",
    "    for i in range(1, totalNumberofValuesfactor3+1):\n",
    "        for j in range(totalNumberofValuesfactor1):\n",
    "            if indxfactor3[j]==i:\n",
    "                factor3.val[i-1] = factor3.val[i-1] + factor1.val[j]\n",
    "                \n",
    "    return factor3\n",
    "\n",
    "def observe_evidence(probab_tables, evidence):\n",
    "    \"\"\"\n",
    "    Incorporates evidence in the \n",
    "    \n",
    "    input:  probab_tables: List of Struct's\n",
    "            evidence:      np.array([])\n",
    "    \n",
    "    output: probab_tables: List of Struct's: the function modifies the input factor to incorporate evidence.\n",
    "    \n",
    "    \"\"\"\n",
    "    if evidence.size!=0:\n",
    "        for rowIndex, element in enumerate(evidence):\n",
    "            #print(rowIndex, element)\n",
    "            variable = element[0] # which RV is evidence?\n",
    "            value    = element[1] # the value it takes?\n",
    "            #print(variable, value)\n",
    "\n",
    "            for probab_table in probab_tables:\n",
    "\n",
    "                if variable in probab_table.var:\n",
    "                    totalNumberofValues = math.prod(probab_table.card)\n",
    "                    assignments = IndexToAssignment( list( range(1, totalNumberofValues+1) ),  probab_table.card )\n",
    "                    #print(assignments)\n",
    "                    for ii, rowAssignment in enumerate(assignments):\n",
    "                        for jj, randVar in enumerate(probab_table.var):\n",
    "                            if randVar==variable:\n",
    "                                if assignments[ii, jj]!=value:\n",
    "                                    length = assignments[ii, :].shape[0]\n",
    "                                    #print( assignments[ii, :].reshape( (1, length)).shape )\n",
    "                                    #print(assignments[ii, :], probab_table.card)\n",
    "                                    index_to_be_assigned = AssignmentToIndex(assignments[ii, :].reshape( (1, length) ), probab_table.card)\n",
    "                                    index_to_be_assigned = index_to_be_assigned.item()\n",
    "                                    #print(index_to_be_assigned)\n",
    "                                    probab_table.val[ index_to_be_assigned-1 ] = 0\n",
    "                       \n",
    "    return probab_tables\n",
    "\n",
    "def compute_marginal(variables_in_marginal, probab_tables, evidence):\n",
    "    \"\"\"\n",
    "    input:  variables_in_marginal: list of variables that you want in the marginal\n",
    "                    probab_tables: List of Struct's\n",
    "                         evidence: np.array([])\n",
    "    \n",
    "    output:        marginalFactor: Struct\n",
    "    \n",
    "    \"\"\"\n",
    "    probab_tables_with_evidence = observe_evidence(probab_tables, evidence)\n",
    "    joint_with_evidence = compute_joint_distribution(probab_tables_with_evidence)\n",
    "    \n",
    "    variableOut = sorted( list( set(joint_with_evidence.var) - set(variables_in_marginal) ) )\n",
    "    marginalFactor = sum_marginalization(joint_with_evidence, variableOut)\n",
    "    #print(marginalFactor.var, marginalFactor.card, marginalFactor.val)\n",
    "    \n",
    "    marginalFactor.val = list( np.array(marginalFactor.val)/sum(marginalFactor.val) )\n",
    "    return marginalFactor"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "692f58d0",
   "metadata": {},
   "source": [
    "## 2. [Inference] Posterior Probability Inference: Variable Elimination.\n",
    "\n",
    "The goal is to compute the posterior distribution of any of the random variable given any other in order to confirm that the above solver is correct. Few of them that I compute are the following:\n",
    "$$ \n",
    "\\begin{align}\n",
    "\\mathbb{P}(\\; C \\;\\vert \\; A = 1)\\ ,\n",
    "\\end{align}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "36ca1333",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1, 3] [2, 3] [0.30000000000000004, 0.0, 0.25, 0.0, 0.44999999999999996, 0.0]\n"
     ]
    }
   ],
   "source": [
    "evidence = np.array([[let_to_num['A'], 1]])\n",
    "marginalFactor = compute_marginal([let_to_num['C'], let_to_num['A']], probab_tables, evidence)\n",
    "print(marginalFactor.var, marginalFactor.card, marginalFactor.val)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7570241c",
   "metadata": {},
   "source": [
    "Another one to compute is \n",
    "$$ \n",
    "\\begin{align}\n",
    "\\mathbb{P}(\\; E \\;\\vert \\; C=3, F=1)\\ ,\n",
    "\\end{align}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "40668a7f",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[3, 5, 6] [3, 2, 2] [0.0031934804413239716, 0.03991850551654965, 0.03113643430290873, 0.15648054162487465, 0.09314317953861584, 0.20837459879638917, 0.04490431293881645, 0.10758324974924774, 0.1220836008024072, 0.09542166499498496, 0.00935506519558676, 0.08840536609829487]\n"
     ]
    }
   ],
   "source": [
    "evidence = np.array([])\n",
    "marginalFactor = compute_marginal([let_to_num['E'], let_to_num['C'], let_to_num['F']], probab_tables, evidence)\n",
    "print(marginalFactor.var, marginalFactor.card, marginalFactor.val)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "3a00334a",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[3, 5, 6] [3, 2, 2] [0.0, 0.0, 0.13000000000000003, 0.0, 0.0, 0.87, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]\n"
     ]
    }
   ],
   "source": [
    "evidence = np.array([[let_to_num['C'], 3], \n",
    "                     [let_to_num['F'], 1]])\n",
    "marginalFactor = compute_marginal([let_to_num['E'], let_to_num['C'], let_to_num['F']], probab_tables, evidence)\n",
    "print(marginalFactor.var, marginalFactor.card, marginalFactor.val)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fc3de7ae",
   "metadata": {},
   "source": [
    "### 2.1 In order to interpret outcome of probability values from the output of the above factor, it is essential to use the following handy tool to see how to read off the .val factor fields."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "fc984643",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[1 1 1]\n",
      " [2 1 1]\n",
      " [3 1 1]\n",
      " [1 2 1]\n",
      " [2 2 1]\n",
      " [3 2 1]\n",
      " [1 1 2]\n",
      " [2 1 2]\n",
      " [3 1 2]\n",
      " [1 2 2]\n",
      " [2 2 2]\n",
      " [3 2 2]] [ 1  2  3  4  5  6  7  8  9 10 11 12]\n"
     ]
    }
   ],
   "source": [
    "D = [3, 2, 2]\n",
    "I = list( range(1, np.prod(D) + 1) )\n",
    "A = IndexToAssignment(I, D)\n",
    "print(A, AssignmentToIndex(A, D))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d1a7161f",
   "metadata": {},
   "source": [
    "## 3. [Learning] Parameter Learning of Each CPD:\n",
    "- In order to learn the parameters of the above graphical model, the first step is to generate a dataset (samples faithful to this particular graphical model).\n",
    "- Once a dataset is generated, the goal is to implement a learning algorithm that can learn all the parameters from that dataset of the graphical model.\n",
    "- Finally, the goal is to figure out a metric that quantifies the quality of learned parameters (one naive way is to naively compare each CPD parameter of the ground truth and learned but this get's out of hand quickly).\n",
    "\n",
    "### 3.1 This step has already been done in the other notebook of binary RV. Refer to that notebook for more details."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "eae3785c",
   "metadata": {},
   "source": [
    "### If all RVs Binary: 4 *effective* data points for linear regression.\n",
    "- C = 0, A = 0, P = 0 (Okay)\n",
    "- C = 1, A = 0, P = 0 (Okay)\n",
    "- C = 0, A = 1, P = 0 \n",
    "- C = 1, A = 1, P = 0 (Okay)\n",
    "- C = 0, A = 0, P = 1 \n",
    "- C = 1, A = 0, P = 1 \n",
    "- C = 0, A = 1, P = 1\n",
    "- C = 1, A = 1, P = 1 (Okay)\n",
    "\n",
    "### If all RVs but C$\\in \\{0, 1, 2\\}$ are Binary: 7 *effective* data points for linear regression when $|C| = 3$.\n",
    "- C = 0, A = 0, P = 0 (Okay)\n",
    "- C = 1, A = 0, P = 0 (Okay)\n",
    "- C = 2, A = 0, P = 0 (Okay)\n",
    "- C = 0, A = 1, P = 0 \n",
    "- C = 1, A = 1, P = 0 (Okay)\n",
    "- C = 2, A = 1, P = 0 (Okay)\n",
    "- C = 0, A = 0, P = 1 \n",
    "- C = 1, A = 0, P = 1 \n",
    "- C = 2, A = 0, P = 1 \n",
    "- C = 0, A = 1, P = 1\n",
    "- C = 1, A = 1, P = 1 (Okay)\n",
    "- C = 2, A = 1, P = 1 (Okay)\n",
    "\n",
    "### If all RVs but C$\\in \\{0, 1, 2\\}$,  A$\\in \\{0, 1, 2\\}$ are Binary: 9 *effective* data points for linear regression when $|C| = 3$ and  $|A| = 3$.\n",
    "- C = 0, A = 0, P = 0 (Okay)\n",
    "- C = 1, A = 0, P = 0 (Okay)\n",
    "- C = 2, A = 0, P = 0 (Okay)\n",
    "- C = 0, A = 1, P = 0 \n",
    "- C = 1, A = 1, P = 0 (Okay)\n",
    "- C = 2, A = 1, P = 0 (Okay)\n",
    "- C = 0, A = 2, P = 0 \n",
    "- C = 1, A = 2, P = 0 \n",
    "- C = 2, A = 2, P = 0 (Okay)\n",
    "- C = 0, A = 0, P = 1 \n",
    "- C = 1, A = 0, P = 1 \n",
    "- C = 2, A = 0, P = 1 \n",
    "- C = 0, A = 1, P = 1\n",
    "- C = 1, A = 1, P = 1 (Okay)\n",
    "- C = 2, A = 1, P = 1 (Okay)\n",
    "- C = 0, A = 2, P = 1\n",
    "- C = 1, A = 2, P = 1 \n",
    "- C = 2, A = 2, P = 1 (Okay)\n",
    "\n",
    "### If all RVs but C$\\in \\{0, 1, 2, 3\\}$,  A$\\in \\{0, 1, 2\\}$ are Binary: 14 *effective* data points for linear regression when $|C| = 4$ and  $|A| = 3$.\n",
    "- C = 0, A = 0, P = 0 (Okay)\n",
    "- C = 1, A = 0, P = 0 (Okay)\n",
    "- C = 2, A = 0, P = 0 (Okay)\n",
    "- C = 3, A = 0, P = 0 (Okay)\n",
    "- C = 0, A = 1, P = 0 \n",
    "- C = 1, A = 1, P = 0 (Okay)\n",
    "- C = 2, A = 1, P = 0 (Okay)\n",
    "- C = 3, A = 1, P = 0 (Okay)\n",
    "- C = 0, A = 2, P = 0 \n",
    "- C = 1, A = 2, P = 0 \n",
    "- C = 2, A = 2, P = 0 (Okay)\n",
    "- C = 3, A = 2, P = 0 (Okay)\n",
    "- C = 0, A = 0, P = 1 \n",
    "- C = 1, A = 0, P = 1 \n",
    "- C = 2, A = 0, P = 1 \n",
    "- C = 3, A = 0, P = 1 \n",
    "- C = 0, A = 1, P = 1\n",
    "- C = 1, A = 1, P = 1 (Okay)\n",
    "- C = 2, A = 1, P = 1 (Okay)\n",
    "- C = 3, A = 1, P = 1 (Okay)\n",
    "- C = 0, A = 2, P = 1\n",
    "- C = 1, A = 2, P = 1 \n",
    "- C = 2, A = 2, P = 1 (Okay)\n",
    "- C = 3, A = 2, P = 1 (Okay)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6343c6a1",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}