Python/Learning & Inference in Linear Gaussian Bayesian Networks.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "3f822df8",
   "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": "feb8d586",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import networkx as nx        # for graphical model visualization\n",
    "from numpy.linalg import inv\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy import linalg"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a4750924",
   "metadata": {},
   "source": [
    "## 1. [Modeling] Computing the Joint Distribution from the Conditional Probability Distribution (CPD) Parameters: Linear Gaussian Bayesian Networks.\n",
    "\n",
    "The Conditional Probability Distriibution (CPD) for each factor in a linear Guassian bayesian network can be parameterized via the equation:\n",
    "$$\n",
    "\\begin{align}\n",
    "    \\mathbb{P}_{X_{\\scriptsize \\mbox{child}}|U_{\\scriptsize \\mbox{parents}}} = \\mathcal{N}\\left(b + \\sum_{i \\in U_{parents}} w_{i} \\cdot u_{i} ; \\sigma^{2}\\right)\\ .\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "- i.e., each random variable is associated with the following parameters:\n",
    "    - 1 $\\times$ conditional variance parameter\n",
    "    - 1 $\\times$ conditional bias parameter\n",
    "    - numParents $\\times$ weight parameters"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "df43e04a",
   "metadata": {},
   "source": [
    "- In order to compute the joint distribution of the entire graphical model, one needs to compute mean $\\mu$ and covariance $\\Sigma$. \n",
    "- $\\mu$ and $\\Sigma$ can directly computed as follows:\n",
    "    - Let $A$ represent the lower triangular adjacency matrix of the graphical model (function of weight edges).\n",
    "    - Let $S$ denote the diagonal conditional variance matrix.\n",
    "    - Then, $\\Sigma$ can be computed via the formula:\n",
    "        - $\\Sigma = (I - A)^{-1} \\cdot S \\cdot [inv(I - A)^{-1}]^{\\top}$\n",
    "    - $\\mu$ can be computed hierarchically starting from the root node. Ofcourse, $\\mu$ is some function of $A$."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "695eadb7",
   "metadata": {},
   "source": [
    "- Before we fit this Gaussian model to real Amazon data, let's make sure that the learning algorithm developed actually works. To ensure that, I will create a random Gaussian BN with arbitrary parameters, generate data according to this model and then learn these parameters via MLE. To confirm whether the learning algorithm works or not, we can match the parameters :)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4209590a",
   "metadata": {},
   "source": [
    "### 1.1 Graphical Model Construction.\n",
    "\n",
    "- Assuming the following parameters for the conditional Gaussian graphical model (arbitrary parameters):\n",
    "\n",
    "|RV |Cond. Bias $b_{A}$| Weight $w_{AB}$ |Cond. Variance $\\sigma^{2}_{A}$|\n",
    "| :---|  ----:  |  ----:  |  ----:   |\n",
    "|A  |0.4                   |        N/A                   |  0.2         | \n",
    "\n",
    "\n",
    "|RV |Cond. Bias $b_{B \\vert A}$| Weight $w_{AB}$ | Cond. Variance $\\sigma^{2}_{B \\vert A}$|\n",
    "| :----|  ----:  |  ----:  |  ----:                           |\n",
    "|B     |0        |  0.3    |  0.4                             | \n",
    "\n",
    "\n",
    "|RV |Cond. Bias $b_{C \\vert A}$| Weight $w_{AC}$ | Cond. Variance $\\sigma^{2}_{C \\vert A}$|\n",
    "| :----|  ----:  |  ----:  |  ----:   |\n",
    "|C  |0.2                 | 2               | 0.3 |\n",
    "\n",
    "\n",
    "|RV |Cond. Bias $b_{D \\vert B,C}$| Weights $w_{BD}$, $w_{CD}$ | Cond. Variance $\\sigma^{2}_{D \\vert B,C}$|\n",
    "| :----|  ----:  |  ----: |  ----:  |\n",
    "|D  |0                     |$[1, 2]$                   |  0.3                               | \n",
    "\n",
    "\n",
    "|RV |Cond. Bias $b_{E \\vert C}$| Weight $w_{CE}$ | Cond. Variance $\\sigma^{2}_{E \\vert C}$|\n",
    "| :----|  ----:  |  ----:  |  ----:   |\n",
    "|E  |12                  |   1             | 0.2                              | "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "3c687bd8",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAgQAAAGFCAYAAACCBut2AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAA9hAAAPYQGoP6dpAABS3UlEQVR4nO3dd1hUZ9oG8HuoSgQs2NeuUVxb1Ji4pBiN0iKiggmWjQ01dgURRFCigKCkiFGMBQ1qEgE7MphslERNNFG/tcSNJRhbBBtFASlzvj+OqMcZFISZM+X+XReXznPemfOgzDs3pyoEQRBAREREJs1M7gaIiIhIfgwERERExEBAREREDAREREQEBgIiIiICAwERERGBgYCIiIgAWFRkkEqlwvXr12FrawuFQqHtnoiIiKgaCIKAvLw8NGnSBGZmz94GUKFAcP36dTRr1qxamiMiIiLdunLlCv7xj388c0yFAoGtre2jF7Szs6t6Z0RERKR1ubm5aNas2aPP8WepUCAo201gZ2fHQEBERGRgKrK7nwcVEhEREQMBERERMRAQERERGAiIiIgIDAREREQEBgIiIiICAwERERGBgYCIiIjAQEBERERgICAiIiIwEBAREREYCIiIiAgMBERERIQK3u2QiMh4FQM4CeAYgOMA/gbwAIA1gMYAugPoAaALAEuZeiTSPgYCIjJRfwFYDSAOwN2HNUuIAQFPPF798O91AEwCMBFACx31SKQ73GVARCYmB4AvgFYAovE4DADSMPD047sPx7d6+PxcLfZIpHsMBERkQvYB6AAgHoAAoLSSzy99+Lz4h6+zr1q7I5ITAwERmYgVAJwBZKHyQeBppQAyH77eF1V8LSL9wEBARCbgCwDTHv5dVU2vWfY6U8FQQMaAgYCIjNw+iB/a2jQV3H1Aho6BgIiMWA6AD6H9qc4MwGjwQEMyZAwERGTE/AHcRPXtJiiPCuKxCX5aXg+R9jAQEJGRugRgHap+AGFFlT5c3186Wh9R9WIgICIj9SU0TXHLlwMKBdCp07OfPXs20LXr48f37gEzZwJNmgA1agDdugHffPP0s8werpfI8DAQEJERKoZ4BUL1rQPr14t/njkDHDlS/its2wYMHfr48ZAhwMaNwIIFQGoq8OqrgI8PsGXLk88qBbAK6hc4ItJ/DAREZIROQnoFQtFvvwH//S/g7i4+XrdO87N//RX466/HgWDvXuC774CVK4GJE4F33gHWrAH69wfmzAFKJbnjLoBT1fi9EOkGAwERGaFjGqtlAWDJEuBf/xI3+efnq49LTgbatwf++U/x8fbtQK1agLe3dNyYMcD165q2NGheP5E+YyAgIiN0HE/fmbCgAPj6a3FTf6dOwNixQF4ekJio/uzkZOnugtOnAUdHwOKp28F16fJ4+WOWYCAgQ8RAQERG6G88vR8/KQnIyQHGjRMfv/+++Fv/07sN/vtf4MIFaSC4fRuoW1d9LWW127efrBYDuFHF/ol0j4GAiIzQA7XKunVAzZrABx+Ij8t2Afz0E3D+/ONxyclAy5ZA9+7S5ysU5a9NfVnhizRNJCsGAiIyQtaSRxcuAD/+KB5MKAhAdrb45eUlLi878wAQtyQ8uXUAAOrVe3orgOjOHfFP9a0HNarQO5E8GAiIyAg1xpPHEKxfLwaBpCSgTp3HX2VnG2zcKJ4pcPas+PV0IOjcWayXlEjrpx6eTCC9poElgEbV/P0QaR8DAREZoe4oO4agtFT8wG/TBti/X/3Lzw/4+2/x2gLJyeKFh15/XfpqgweLFyZKTpbWN24Ux7/22pPVYgA9tPi9EWmHxfOHEBEZmscfyKmp4qmBUVFAnz7qIzt1AlasEI8xyMgQL0D09DEBrq7iNQc++gjIzQXathXPWFAqgU2bAHPz8tdPZCgYCIjICHUBUAfAXaxbB1hZidcM0MTBQdwCkJQk7hL47DPN47ZtA4KDgdBQ8diBDh3EUFB2kOJjdQB0rq5vhEhnFIIgCM8blJubC3t7e+Tk5MDOzk4XfRERvbArV67gxAlXuLmdUbt2QHmio4Fly8TdB+q/8VeUOYC5AMJf9AWIqlVlPr95DAERGY28vDwEBwfj5ZdfxrRpZ2BWiRkuIADIyqpKGADE2yBPqMoLEMmGgYCIDF5JSQlWr16Ntm3bIiIiAoWFhbh8WTwu4OkzA7THHMA4AC10tUKiasVAQEQGTalUolu3bpg0aRKysrIky/z8gJs3FVCpnnFVoWphBqABgBgtr4dIexgIiMggnTp1Cs7OznB1dcWZM2fUlpubm2PUqMmoWfMbmJk991CpKlIB2ACAx1iR4eJZBkRkUG7cuIGQkBCsX78eKpVK45j33nsP0dHRcHR0fFi5CWCqFrv6AsAALb4+kfYxEBCRQcjPz0dMTAyioqJw//59jWO6deuGmJgY9O3b96klUx7+ORXihlHNQaJyyl7nCwCTq+H1iOTFQEBEek2lUiEhIQHBwcG4du2axjFNmjRBeHg4Ro0aBfNyTxOYAqAdgNEAsgCUVqErc4jHDGwAtwyQseAxBESkt/bv34+ePXti9OjRGsOAjY0NwsLCcO7cOYwePfoZYaDMAABnAYwBoID4wV4Z5g+fNwbA/8AwQMaEgYCI9M4ff/wBDw8P9O3bFydOnFBbrlAoMG7cOFy4cAGhoaF46aWXKvHq9gDWAMiAeBGhOo+WFBcrUHapNkEASkufDAx1Ho7PePh8HkBIxoW7DIhIb9y8eRNhYWGIi4tDaanmTfrvvvsuli1bhq5du1ZxbS0gXlFwIYBTAI4hJSUEQCZq1AAKC4GXX/4XOnYcBfHeBJ3x5B0UiYwNAwERya6wsBDLly9HeHg4cnNzNY5xdHTEsmXL4OrqCsXTdx+qEkuId0fsjmXLNuLQocxHS+LiRqBjR99qXBeR/mIgICLZCIKAb7/9FkFBQbh06ZLGMfXr18fHH3+M8ePHw6KiNyYgokrju4uIZHH48GH4+fnhl19+0bjc2toas2fPRmBgIG+qRqQDDAREpFN//vknAgMDkZiYWO6YESNGICIiAs2bN9dhZ0SmjYGAiHTi7t27CA8PR2xsLIqKijSOefPNNxETE4NXX31Vx90REQMBEWlVcXExVq1ahbCwMNy5c0fjmLZt2yIqKgqDBw+u5gMGiaiiGAiISCsEQcDOnTsREBCA8+fPaxxTp04dhIaGYvLkybCystJxh0T0JAYCIqp2x44dg5+fH9LT0zUut7S0xNSpUzF//nzUrVtXx90RkSYMBERUba5cuYLg4GAkJCSUO2bo0KGIiopCmzZtdNgZET0PAwERVVleXh6ioqIQExODwsJCjWN69eqFmJgYvPHGGzrujogqgoGAiF5YSUkJ1q9fj9DQUGRmZmoc06JFC0RGRuL999+HmRlvn0KkrxgIiOiFKJVK+Pv748yZMxqX29nZYd68eZgxYwZq1Kih4+6IqLIYCIioUk6dOgV/f3/s27dP43Jzc3NMnDgRCxcuRP369XXcHRG9KAYCIqqQGzduICQkBOvXr4dKpdI45r333kN0dDQcHR113B0RVRUDARE9U35+PmJiYhAVFYX79+9rHNO1a1fExMSgX79+Ou6OiKoLAwERaaRSqZCQkIDg4GBcu3ZN45jGjRsjPDwc//73v2Fubq7jDomoOjEQEJGa/fv3w8/PDydOnNC43MbGBgEBAfD398dLL72k4+6ISBsYCIjokT/++ANz5szB7t27NS5XKBQYM2YMFi1ahCZNmui4OyLSJgYCIsKtW7cQFhaGuLg4lJSUaBzz7rvvYtmyZejatauOuyMiXWAgIDJhhYWFiI2NRXh4OHJycjSOcXR0xLJly+Dq6so7ERIZMQYCIhMkCAK2bt2KwMBAXLp0SeOY+vXr4+OPP8b48eNhYcGpgsjY8V1OZGJ+/vlnzJ49G7/88ovG5dbW1pg1axaCgoJgZ2en4+6ISC4MBEQm4s8//0RgYCASExPLHTN8+HBERESgRYsWOuyMiPQBAwGRkbt79y7Cw8MRGxuLoqIijWPeeOMNxMTEoFevXjrujoj0BQMBkZEqLi7GqlWrEBYWhjt37mgc06ZNG0RHR2Pw4ME8YJDIxDEQEBkZQRCwc+dOBAQE4Pz58xrH1KlTB6GhoZg8eTKsrKx03CER6SMGAiIjcuzYMfj5+SE9PV3jcktLS0ydOhXz589H3bp1ddwdEekzBgIiI3DlyhUEBwcjISGh3DFDhgxBVFQU2rZtq8POiMhQMBAQGbC8vDxERUUhJiYGhYWFGse8+uqriImJwZtvvqnj7ojIkDAQEBmgkpISrF+/HqGhocjMzNQ4pnnz5oiMjMQHH3wAMzMzHXdIRIaGgYDIwCiVSvj7++PMmTMal9va2mLevHmYMWMGatasqePuiMhQMRAQGYhTp07B398f+/bt07jc3NwcEyZMwMKFC9GgQQMdd0dEho6BgEjP3bhxAyEhIVi/fj1UKpXGMe7u7oiOjkbHjh113B0RGQsGAiI9lZ+fj08++QRLlizB/fv3NY7p0qULYmJi8O677+q4OyIyNgwERHpGpVJh06ZNmDdvHq5du6ZxTOPGjbF48WJ8+OGHMDc313GHRGSMGAiI9MiBAwfg5+eH48ePa1xuY2ODOXPmwN/fH7Vq1dJxd0RkzBgIiPTAH3/8gYCAAOzatUvjcoVCgdGjR2PRokVo2rSpjrsjIlPAQEAko1u3biEsLAxxcXEoKSnROKZfv35YtmwZunXrptvmiMikMBAQyaCwsBCxsbEIDw9HTk6OxjGOjo5YunQp3NzceCdCItI6BgIiHRIEAVu3bkVgYCAuXbqkcUz9+vURFhYGX19fWFjwLUpEusHZhkhHfv75Z8yePRu//PKLxuXW1taYOXMmgoKCYG9vr+PuiMjUMRAQadmff/6JwMBAJCYmljtm+PDhiIiIQIsWLXTYGRHRYwwERFqSnZ2NxYsXIzY2FkVFRRrHODk54ZNPPkGvXr103B0RkRQDAVE1Ky4uRlxcHMLCwnD79m2NY9q0aYOoqCgMGTKEBwwSkV5gICCqJoIgYNeuXQgICMC5c+c0jqlTpw5CQkIwZcoUWFlZ6bhDIqLyMRAQVYNjx47Bz88P6enpGpdbWlpiypQpCAkJQd26dXXcHRHR8zEQEFXBlStXEBwcjISEhHLHDBkyBFFRUWjbtq0OOyMiqhwGAqIXkJeXh6ioKMTExKCwsFDjmFdffRUxMTF48803ddwdEVHlMRAQVYIgCNi0aRPmzJmDzMxMjWOaNWuGyMhI+Pj4wMzMTMcdUlW0atUK2dnZjx7XqVNHvmaIdIyBgKgSSktLUVhYqDEM2NraYt68eZgxYwZq1qwpQ3dUVc/a9UNk7PjrC5m0X375BYcOHcK5c+dQUFAAAFCpVOWOt7CwwLhx49CpU6dHNXNzc3z00Ue4cOECAgMDGQaIyCAxEJBJOnToEHr37g1vb28MHjwYPXv2xMKFCwHguZv5S0tL8cknnwAA3N3dcfLkSaxcuRINGjTQdttERFrDQEAmJS8vDzNmzICHhwc6duyI+Ph4fPvtt/jXv/6FFStWICUl5bmvYWlpif79++P48ePYs2cPOnbsqIPOiYi0i8cQkMkQBAHp6ek4ceIEoqOj4ePjA2tra5ibm8PW1hZOTk4VvmqgIAh45ZVXtNwxEZHuMBCQyVAoFGjRogX8/f3Rv39/1KxZ89HxArdv30aHDh3QuXPnR+MFQSg3IPByw0RkbBSCIAjPG5Sbmwt7e3vk5OTAzs5OF30R6czRo0cxatQoqFQqtG/fHm+88QbGjh2LBg0aPDMUEBHpu8p8fvMYAjJJKpUKKpUKQUFB6N27NxwdHTFy5EgAwMKFC/HBBx8A4JYAU/POO+8gJibmmWNu3bqFYcOGYf/+/Trqikg3uMuAjEZeXh7y8/NRv379R2cKlPcbftnyV199FQcPHkSPHj0e3WwoPDwcISEh2L9/P9555x1uJTAhP/74I9LT05Gfn4+QkBCNYxwcHHDnzh1s2rQJ77zzjo47JNIebiEgo+Dv74+uXbtiwIAB8PT0xLZt2wCU/xt+2Z6yIUOGoHfv3rCysnpUe/vtt6FQKHDx4sVnvgYZJ1dXVyxYsABhYWHljvH29sbhw4d12BWR9nELARm0mzdv4sMPP8Tly5cxe/Zs5OTk4Ntvv4WXlxc2btyI999//9GH/ZMf7Jo+5BUKBQRBwMGDB2FlZcWzCExUaGgoXn311UeBYMGCBWpjOnbsiCtXrui6NSKtYiAgg3bs2DEcOXIEmzdvhouLCwBg2LBhCA4OxuzZs1G7dm0MHDiwQr/lFxQU4IcffsDmzZvx0UcfoUuXLtpun/TUwoULYWFhgQULFqCkpASLFi2SLL958yZq164tT3NEWsJAQAYtNTUV9vb2cHJyelRr27YtYmJi8PrrryM2NhYdO3ZEmzZtyj0W4LvvvkNWVhZ2796Nbdu2wcfHB6GhobC0tNTlt0J6Zv78+bC0tMS8efOQkZGBlStXws7ODvfv38eyZcvQs2dPuVskqlYMBGTQmjRpgr///hvm5uYAxLMHzMzM0KxZM3z88cfw9fXFf/7zH7Ru3RoKhQKlpaWPxgJAZmYmIiMjcfv2bTRu3BgHDx5Er1695Pp2SM/MnTsXdnZ2mDlzJlJSUtC9e3f8/vvvuH37Nn755Re52yOqVjyokAxat27dYGZmhvj4eLVlH374If75z39i27ZtyM3NBSDeiOj06dOP9v82bNgQy5cvx+bNm6FUKhkGSM1HH32EI0eO4N1338X169fRuXNn7Nu3D927d5e7NaJqxQsTkUF78OABOnfujObNmyMxMRF16tR5tJUAAD799FMEBQXhr7/+QsOGDZGeno533nkHb775JtLT02XunohIu3hhIjIZ1tbWCAsLw7Fjx7Bx40YA0jMIWrVqBUEQcPz4cQDi8QXu7u544403ZOmXiEhfMRCQwXv//ffRtm1bxMXF4cyZM5JAkJubi5KSEjRr1gwA0LRpUyQmJiI8PFyudomI9BIDARmEkpISfPnllzh8+DCKi4sly8zMzLBu3Trk5+dj+PDhOHfuHHJycnDlyhVs374dnp6eaN269aPxNWrU0HX7RER6j2cZkN5LS0uDv78/Tp8+jW7duuHYsWNqY7p06YIvv/wSkyZNQs+ePdGtWzdcv34dxcXF+Oqrr2BjYyND50REhoMHFZLeOn36NPz9/ZGWliapx8fHY8SIERqvE3Dx4kX88MMPOH/+PJo1a4Zp06bpql0iIr1Tmc9vBgLSOzdu3EBoaCjWrVsHlUqltrxJkybIyMh4dDMioupy7do1PHjw4NFjBwcHznlk0HiWARmk/Px8LF68GO3atcOaNWs0hoEuXbpg48aNDAOkFe+//z7atGnz6Ovrr7+WuyUineExBCQ7lUqFTZs2ITg4GFevXtU4pnHjxli8eDE+/PBDyZUGiYioejAQkKwOHDgAPz+/R9cJeJqNjQ3mzJkDf39/1KpVS8fdERGZDgYCksW5c+cQEBCAnTt3alyuUCgwevRoLFq0CE2bNtVxd0REpoeBgHTq1q1b+Pjjj7Fq1SqUlJRoHNOvXz8sW7YM3bp1021zREQmjIGAdOLBgweIjY3F4sWLkZOTo3FMhw4dsGzZMri5uWm8TTEREWkPAwFplSAISExMRGBgIDIyMjSOcXBwQFhYGHx9fTVeW4CIiLSPgYC05ueff4afnx9+/vlnjcutra0xc+ZMBAUFwd7eXsfdERHRkxgIqNplZGQgMDAQW7duLXeMj48PIiIi0LJlS901RkRE5WIgoGqTnZ2N8PBwLF++HEVFRRrHODk5ISYmBq+99pqOuyMiomdhIKAqKy4uRlxcHMLCwnD79m2NY9q0aYOoqCgMGTKEBwwSEekhBgJ6YYIgYNeuXQgICMC5c+c0jqlTpw5CQkIwefJkWFtb67hDIiKqKAYCeiHHjx+Hn58fDhw4oHG5paUlpkyZgpCQENStW1e3zRERUaUxEFClXL16FcHBwUhISEB5N8ocMmQIoqKi0LZtWx13R0REL4qBgCokLy8P0dHRiImJQUFBgcYxr776KmJiYvDmm2/quDsiIqoqBgJ6ptLSUqxfvx4hISHIzMzUOKZZs2aIjIyEj48PzMx4R20iIkPEQEDlSktLg7+/P06fPq1xua2tLYKCgjBz5kzUrFlTx90REVF1YiAgNadPn4a/vz/S0tI0LjczM8OECRMQFhaGBg0a6Lg7IiLSBgYCeuTGjRsIDQ3FunXroFKpNI5xc3PD0qVL0bFjRx13R0RE2sRAQMjPz8enn36KJUuW4N69exrHdOnSBTExMXj33Xd13B0REekCA4EJU6lU2Lx5M+bNm4erV69qHNO4cWMsXrwYH374IczNzXXcIRER6QoDgYlKT0+Hn58fjh07pnF5zZo1MWfOHMyZMwe1atXScXdERKRrDAQm5ty5cwgICMDOnTs1LlcoFBg9ejQWLVqEpk2b6rg7IiKSCwOBibh9+zbCwsKwatUqlJSUaBzTr18/LFu2DN26ddNtc0REJDsGAiP34MEDxMbGYvHixcjJydE4pkOHDli2bBnc3Nx4J0IiIhPFQGCkBEFAYmIiAgMDkZGRoXGMg4MDwsLC4OvrC0tLSx13SERE+oSBwAj9/PPP8PPzw88//6xxubW1NWbOnImgoCDY29vruDsiItJHDARGJCMjA4GBgdi6dWu5Y3x8fBAREYGWLVvqrjEiItJ7DARGIDs7GxEREfj8889RVFSkcYyTkxNiYmLw2muv6bg7IiIyBAwEBqy4uBirV6/GwoULcfv2bY1j2rRpg6ioKAwZMoQHDBIRUbkYCAyQIAjYvXs3AgIC8Mcff2gcU7t2bYSGhmLy5MmwtrbWcYdERGRoGAgMzPHjx+Hn54cDBw5oXG5hYYEpU6YgJCQE9erV021zRERksBgIDMTVq1cRHByMhIQECIKgcczgwYMRFRWFdu3a6bg7IiIydAwEeu7evXuIiopCTEwMCgoKNI7p2bMnYmJi8NZbb+m4OyIiMhYMBHqqtLQU8fHxmD9/PjIzMzWOadasGSIjI+Hj4wMzMzMdd0hERMaEgUAP7du3D/7+/jh16pTG5ba2tggKCsLMmTNRs2ZNHXdHRETGiIFAj5w5cwb+/v5QKpUal5uZmWHChAkICwtDgwYNdNwdEREZMwYCPZCZmYnQ0FCsXbsWKpVK4xg3NzcsXboUHTt21HF3RERkChgIZFRQUIBPPvkES5Yswb179zSO6dKlC2JiYvDuu+/quDsiIjIlDAQyUKlU2Lx5M+bNm4erV69qHNO4cWMsXrwYH374IczNzXXcIRERmRoGAh1LT0+Hn58fjh07pnF5zZo1MWfOHMyZMwe1atXScXdEpiszMxOFhYWS2q1bt3D//n289NJLMnVFpDs8V01Hzp07B09PT/Tp00djGFAoFBg9ejTOnz+PsLAwhgEiHcnOzsZbb72FRo0aqb0358+fj9q1ayM8PFym7oh0h4FAy27fvo3p06fjn//8J3bu3KlxTN++fXH8+HHEx8ejadOmOu6QyLQlJCTgp59+Knd5SUkJQkNDy72BGJGx4C4DLXnw4AFiY2OxePFi5OTkaBzToUMHLF26FO7u7rwTIZFMLC0tnztGEIQKjSMyZNxCUM0EQUBiYiIcHR0xZ84cjWHAwcEBX3zxBU6ePIn33nuPYYBIRt7e3rC1tX3mmGHDhsHOzk5HHRHJg4GgGv3yyy944403MGzYMGRkZKgtt7a2xty5c3HhwgVMnjyZv3EQ6YF69eph+vTp5S5XKBQICQnRYUdE8mAgqAYZGRn44IMP0Lt3bxw+fFjjGB8fH/zvf//DkiVLYG9vr+MOiehZZs2aVe5WgmHDhuGf//ynjjsi0j0GgirIzs5GQEAAOnTogG+//VbjGCcnJ/zyyy/YsmULWrZsqdsGiahCyttKwK0DZEoYCF5AcXExVqxYgbZt22Lp0qUoKipSG9OmTRskJSXhp59+wmuvvSZDl0RUGbNmzVK7CFivXr24dYBMhhGcZVAM4CSAYwCOA/gbwAMA1gAaA+gOoAeALgAqvs++oKAABQUFqFu37qOaIAjYvXs3AgIC8Mcff2h8Xu3atREaGorJkyfD2tr6xb4lItIhcQ6pV+8YUlKaobDwEqytgQcPgF69GgH4Ei8yhxAZGoUgCMLzBuXm5sLe3h45OTl6dKTtXwBWA4gDcPdhzRLimxsaHtcBMAnARAAtnvnK33zzDUaPHo3i4mIEBgYiPDwcJ06cgJ+fH/bv36/xORYWFpgyZQpCQkJQr169F/6uiEhX1OcQQRDnDIUCEARAoXixOYRIX1Tm89sAA0EOAH8A6yDu8SitxHPNAagAjAMQA0D9e/n999/RvXt3PHjw4FFtwIAB+O6771DeP9XgwYMRFRWFdu3aVaIXIpKHducQIn1Smc9vAzuGYB+ADgDiAQio3BsZD8cLD5/f4eHrPVZcXIx///vfkjAAAPv27dMYBnr27In09HRs27aNYYDIIGh3DiEyZAYUCFYAcAaQhcq/iZ9WCiDz4et98agaGRlZ7k2HntSsWTNs2rQJR44cwVtvvVXFXohIN7Q/hxAZMgM5qPALANMe/l1VTa9Z9jpTAQDHj/fGxx9//Mxn1KpVC/PmzcPMmTNRs2bNauqDiLRP+3MIMKWaXpdIHgYQCPbh8RtOW6Zi8eKXUFr67N8apk6diqCgIC33QkTVSzdzCNAOwAAtr4dIe/R8l0EOgA+h7TZVKgVWrLiP51zOHHv37tVqH0RU3XQzh4ivPxpArpbXQ6Q9eh4I/AHcRPVt4tPMzExAgwbAsmXPHscrDRIZGt3MIeLrZwHw0/J6iLRHjwPBJYinBZViwwZAoZB+1a8P9OkD7Nmj+dmzZwNdu4p/z8sDAgKAAQPE5ykUwMKF0vEWFsD48UDz5o9rNWrUgIODA1q0aAEvLy/ExcVV9zdJRFpzCWVzSJmTJ4ExY4BWrYAaNYBatYDu3YHoaODOHemzd+0S54WbN4Fz5wB/f6BHD6B2baBuXcDJCUhKevIZpQ/X95eWvy8i7dDjQPAlnm4vPh74+Wfg8GHgyy8Bc3Ng4EBg9271Z2/bBgwdKv799m1x/IMHgKdn+WtUKMxw+vQ0ZGdno7i4GAUFBbh58yYuXbqExMRENG7cuNq+OyLSNukcsmaN+IH+66/AnDmAUgls3w54ewNxccC4cdJnJycDb70l/hKxbx+QkiLOKYmJwObNQLt24nOlxyKbPVwvkQESKiAnJ0cAIOTk5FRkeDUoEgShjiAIEAQBQnw8BADCr7/iUU0QIOTnQ7C2huDjI60fPSqOP31afKxSiV+CAOHmTXHZggXS5zz+qvNw/URkuKRzyOHDEMzNIbi4QCgsVH/fP3gAYefOx4+LiiDUrg1hxYrH80bZHPLkl7s7BBubp1+Tcwjpj8p8fuvpFoKTeHw54vLVqAFYWQGWT11ePDkZaN8eKLsnSdluhoq5C+BUJXolIv0jnUMiIsQ54MsvAU23GLGyAjw8Hj/+z3+AnBxg8GDxsYOD5jmkVy8gP//p3Q2cQ8gw6Wkg0HxxoNJSoKQEKC4Grl4FZs4E7t8Hhg+XjktOfry7oDrXT0SG4vF7uLQU+OEHcXdBs2YVe3ZyMtC7N9CkybPH7d8v7lJo0KD89RMZCj29DsFxqN+oCHj9dekoa2tgxQrA2flx7b//BS5cqEogsIT4ZvZ90RcgItk9nkNu3RJ/i2/VqmLPLC0FduwAnnfJkbVrgQMHgM8/F49neoxzCBkmPQ0Ef+PpMAAAX30FODqKf791SzwgaMoU8Q089eF1R5KTgZYtxSOHX0wxgBsv+mQi0gua55CKSE8X55chQ8ofk5oqzj1eXsC0aU8v5RxChklPA8EDjVVHR6Bnz8ePXVyAv/4STykcOVI8HSgpqaq7CwCgsKovQESyejyHODgANjZARkbFnpmUJO5eKO+yI2lpYljo318820Dz8UmcQ8jw6OkxBBqO+ilHly5AQYF4nvDZs+JX1QNBjaq+ABHJ6vEcYm4O9OsHHDsmHnv0LCqVuOWxvDkkLU08dfntt8WtkVZW5b0S5xAyPHoaCBpD3A/3fP/3f+Kf9euLb9AmTdSPNagcSwCNqvICRCQ76RwSFAQIAuDrCxQVqY8uLhavZ3L4MHDjhuZAsG+fGAbeeEM8xkDT2QoiziFkmPR0l0F3AKvVqqdPi2cZAOLFhrZtA777Tjw1qFUrcVPfkCGaN+GlpopnJOTliY9///3xVcbc3MRNiqJiAD2q+fshIt2SziG9ewOrVgGTJ4u7Az76SDwtubgYOHFCPB2xUyegRQvxz5dflr7awYNiGGjUCJg37/EvImU6dgTs7MoecQ4hw6SngUDzm2nMmMd/t7cXQ8Ann4hv8osXxTMMPvtM8yt+9JF4vEGZxETxCxD3LUr3F/LNTGTY1N/Dvr7idQM+/RSIihK3BFhaih/+w4eLByb36AGMHav+at9/L+6avHQJ6NtXffn+/eKl1J+1fiJ9pxAEQXjeoNzcXNjb2yMnJwd2j2OwFhUDaIiKXJyoTHS0eHOiv/9++hSgyqoDIBMV3WVBRPqo8nPI0aPAa6+J9zvo3Lkq6+YcQvqjMp/fenoMgSWASQAq/skeEABkZVU1DJgD+Ah8IxMZusrPIb16iccZVC0McA4hw6WngQAAJkL7tyx9mgrABB2vk4i0g3MIUWXocSBoAWAcKpPwq8b84fpa6Gh9RKRdnEOIKkOPAwEAxABoAO23afZwPTFaXg8R6ZY4h2j/Pm6cQ8jw6XkgsAOwAdrf7Kd6uB5dHDBJRLpjh9LSdVAoOIcQPY+eBwIAGABghVbXsGCBA7Kyuml1HUSke4IgYOrUXZgyRbvrmT7dHPv380BCMmwGEAgAYAoeh4Lqabm0VPxz8mTg449vwd3dHffu3auW1yYi/RAeHo64uDisXIlHoaDsvV9VT84hsbGl8PT0xH//+9/qeXEiGRhIIADEUJAG8dziqh0kJAhmuHPHAgMGiFcvA4DffvsNXl5eKC5+sTukEZF+Wb9+PUJCQh49XrkScHe3QElJPVR9DjHHvXs2kjkkNzcXrq6u+OvJK6ARGRADCgSAuPvgLIAxABSo/JvaHIACCsVYqFS/4/z5lpKlaWlpGD9+PCpwrSYi0mMpKSmYMEF6+p9CocCYMV/D2voiqj6HjEGtWldhZye96cHff/8NFxcX3L59uwrdE8nDwAIBANgDWAMgA8BciFcFK/P0PrwnH9d5OD4DwBo0bNgOaWlpqFevnuQZX331FYKDg6u9ayLSjaNHj2LYsGEofWrfwOeffw4vLy9U1xxibl4HmzZtwptvvil5xv/+9z94eHigoKCgWr4fIl3R00sXV0YxgFMAjj38ugHxXuQ1IN5xrMfDr87QdPWwI0eO4J133lF788bGxmLq1Kla7ZyIqte5c+fg5OSEW7duSepz587FkiVLynlW1eaQu3fv4s0338SZM2ck9UGDBiEpKQkWFnp6yxgyCZX5/DaCQFB1e/bsgaenp+Q3CoVCgcTERAwt78boRKRXbty4gX/961/IyMiQ1EeNGoWNGzdCoek2qNXkypUr6N27N65duyapT5w4EatWrdLquomexQjuZaBb7733HuLi4iQ1QRAwYsQI/PTTTzJ1RUQVlZeXB3d3d7UwMGDAAKxbt07rH8jNmjWDUqlE7dq1JfXVq1dj8eLFWl03UXVhIHho/PjxCAsLk9QePHgADw8PtU2BRKQ/ioqKMHToUBw/flxS79GjB5KSkmBpqZvrA3Tq1Ak7d+6EtbW1pB4aGop169bppAeiqmAgeEJISIjakcnZ2dlwcXHBlStXZOqKiMqjUqkwbtw4fPfdd5J669atkZKSAltbW53289Zbb2Hz5s1qWyQmTpyIlJQUnfZCVFkMBE9QKBT44osv4OHhIalfvXoVrq6uuHu34vdWJyLtCwoKwqZNmyQ1BwcHKJVKNGzYUJaehg4diuXLl0tqpaWl8Pb2xpEjR2TpiagiGAieYmFhga+//hq9e/eW1M+cOQNPT08UFhbK1BkRPWn58uWIjo6W1GxsbJCSkoJ27drJ1JVo6tSpCAwMlNQKCgrg7u6Oc+fOydQV0bMxEGhgY2OD3bt3o0OHDpL6jz/+iJEjR6qd30xEurV161bMnDlTUjM3N0dSUhJ69eolT1NPiYiIwL///W9J7fbt23B2dsaNGzdk6oqofAwE5ahXrx6USiUaN24sqScnJ2PmzJm8miGRTA4cOIBRo0apvQfXrl0LV1dXmbpSp1AosHbtWjg7O0vqly5dgpubG3Jzc2XqjEgzBoJnaNGiBVJTU9XO3VyxYgWioqJk6orIdJ06dQqenp4oKiqS1BcvXozRo0fL09QzWFpaIikpCT169JDUT5w4gaFDh6p9H0RyYiB4jq5du2L79u1qpy4FBQXhq6++kqkrItNz+fJluLi4ICcnR1KfPHky5s2bJ1NXz1erVi2kpKSgTZs2kvr333+PsWPHQqVSydQZkRQDQQX07dtX44f/uHHjoFQqZeiIyLTcuXMHLi4uuH79uqQ+ePBgLF++XO+vBNiwYUMolUrUr19fUt+8ebPawYdEcmEgqKAPPvgAn3zyiaRWUlICLy8v/PbbbzJ1RWT8CgoK4OHhgbNnz0rqTk5O2Lx5M8zNq3YrY11p27YtUlJSYGNjI6kvXboUn3/+uUxdET3GQFAJs2bNgp+fn6R2//59uLu74+LFizJ1RWS8SktLMXz4cBw6dEhSd3R0xK5du1CzZk2ZOnsxr776KpKSktRCzKxZs7B161aZuiISMRBUUnR0NIYPHy6pZWVlwdnZGVlZWTJ1RWR8BEHAtGnTsGPHDkm9adOmUCqVqFu3rjyNVZGrqyvWrl0rqQmCgFGjRmH//v0ydUXEQFBpZmZmiI+PR79+/ST1ixcvwt3dHffu3ZOpMyLjEhERgVWrVklq9vb2SE1NRfPmzWXqqnqMHj0a4eHhklpRURE8PT1x8uRJmboiU8dA8AKsrKywbds2dOvWTVL/7bff4O3tjeLiYnkaIzIS8fHxmD9/vqRmZWWFHTt2oHPnzjJ1Vb2CgoIwefJkSS03Nxeurq64fPmyTF2RKWMgeEF2dnbYu3cvWrZsKakrlUr4+vrywkVEL2jv3r3w9fWV1BQKBTZt2oQ+ffrI05QWKBQKLF++HEOGDJHUr1+/DhcXF9y5c0emzshUMRBUQePGjaFUKlGvXj1JfePGjQgODpapKyLDdfToUXh7e6tdHvyzzz6Dt7e3TF1pj7m5OTZt2oQ33nhDUj979iwGDhyIgoICmTojU8RAUEXt27fHnj171I52joyMxIoVK2TqisjwnD9/Hu7u7sjPz5fUAwICMH36dJm60r6aNWti165d6Nixo6R++PBhDB8+nPdOIZ1hIKgGr7/+Or799luYmUn/OadPn47k5GSZuiIyHDdu3ICzszNu3bolqY8cORKRkZEydaU7derUgVKpRNOmTSX1HTt2YOrUqdwFSTrBQFBNBg4ciNWrV0tqgiBgxIgR+Omnn2Tqikj/5eXlwd3dHRkZGZL6gAEDsG7dOrWgbayaNWsGpVIJe3t7ST0uLk7tjAQibTCNd5qOjB8/HgsXLpTUHjx4AA8PD5w5c0aepoj0WFFREYYOHYrjx49L6t27d0dSUhKsrKxk6kwenTp1ws6dO9W+75CQEKxbt06mrshUMBBUs9DQULUjpLOzs+Hi4oIrV67I1BWR/lGpVBg3bhy+++47Sb1169bYu3cvbG1tZepMXm+//TY2b96sdn+GiRMnIiUlRaauyBQwEFQzhUKBlStXwsPDQ1K/evUqXF1dcffuXZk6I9IvQUFB2LRpk6Tm4OAApVKJhg0bytSVfvDy8lK7v0FpaSm8vb1x5MgRmboiY8dAoAUWFhb4+uuv0bt3b0n9zJkz8PT0RGFhoUydEemH5cuXIzo6WlKzsbFBSkoK2rVrJ1NX+mXatGmYO3eupFZQUAB3d3ecO3dOpq7ImDEQaImNjQ12796NDh06SOo//vgjRo4cyVOJyGRt3boVM2fOlNTMzc2RlJSEXr16ydOUnoqMjMSoUaMktdu3b8PZ2Rk3btyQqSsyVgwEWlSvXj0olUo0btxYUk9OTsbMmTN5KhGZnAMHDmDUqFFqP/tr166Fq6urTF3pL4VCgXXr1mHAgAGS+qVLl+Dm5obc3FyZOiNjxECgZS1atEBqairs7Owk9RUrViAqKkqmroh079SpU/D09ERRUZGkvnjxYowePVqepgyApaUlkpKS0KNHD0n9xIkTGDp0qNq/J9GLYiDQga5du2L79u2wtLSU1IOCgvDVV1/J1BWR7ly+fBkuLi7IycmR1CdPnox58+bJ1JXhsLW1RUpKClq3bi2pf//99xg7dixUKpVMnZExYSDQkb59+2r88B83bhyUSqUMHRHpxp07d+Di4oLr169L6oMHD8by5cvVTq8jzRo2bAilUgkHBwdJffPmzQgMDJSpKzImDAQ69MEHH+CTTz6R1EpKSuDl5YXffvtNpq6ItKegoAAeHh44e/aspO7k5ITNmzfD3Nxcps4MU7t27ZCSkgIbGxtJfenSpWqnKRJVFgOBjs2aNQt+fn6S2v379+Hu7o6LFy/K1BVR9SstLcXw4cNx6NAhSd3R0RG7du1SuyEYVUyvXr2QlJSkFqZmzZqFrVu3ytQVGQMGAhlER0dj+PDhklpWVhacnZ2RlZUlU1dE1UcQBEybNg07duyQ1Js2bQqlUom6devK05iRcHV1xdq1ayU1QRAwatQo7N+/X6auyNAxEMjAzMwM8fHx6Nevn6R+8eJFuLu74969ezJ1RlQ9IiIisGrVKknN3t4eqampaN68uUxdGZfRo0dj8eLFklpRURE8PT1x8uRJmboiQ8ZAIBMrKyts27YN3bp1k9R/++03eHt7o7i4WJ7GiKooPj4e8+fPl9SsrKywY8cOdO7cWaaujNO8efMwefJkSS03Nxeurq64fPmyTF2RoWIgkJGdnR327t2Lli1bSupKpRK+vr68cBEZnL1796rd3EuhUGDTpk3o06ePPE0ZMYVCgeXLl2Pw4MGS+vXr1+Hi4oI7d+7I1BkZIgYCmTVu3BhKpRL16tWT1Ddu3Kj2WxaRPjt69Ci8vb3VLsv92WefwdvbW6aujJ+5uTk2b94MJycnSf3s2bMYOHAgCgoKZOqMDA0DgR5o37499uzZo3bUdUREBFauXClTV0QVd/78ebi7uyM/P19SDwgIwPTp02XqynTUrFkTu3btgqOjo6R++PBhDB8+nPdOoQphINATr7/+Or799luYmUn/S6ZOnYpt27bJ1BXR82VmZsLZ2Rm3bt2S1EeOHInIyEiZujI9devWhVKpRNOmTSX1HTt2YOrUqdwFSc/FQKBHBg4ciNWrV0tqgiBg+PDh+Omnn2Tqiqh8eXl5cHNzQ0ZGhqQ+YMAArFu3Ti3gknY1b94cqampsLe3l9Tj4uIQHh4uU1dkKPhu1TPjx4/HwoULJbUHDx7Aw8MDZ86ckacpIg2Kiorg5eWF48ePS+rdu3dHUlISrKysZOrMtHXu3Bk7duxQ+/cPCQnB+vXrZeqKDAEDgR4KDQ1VO1I7OzsbLi4uuHr1qkxdET0mCALGjx+Pffv2SeqtWrVCSkoKbG1tZeqMAKBPnz7YtGmT2n0iJkyYgJSUFJm6In3HQKCHFAoFVq5cCQ8PD0n96tWrcHFxwd27d2XqjEgUFBSEhIQESc3BwQFpaWlo1KiRTF3Rk7y9vfHZZ59JaqWlpRg2bBiOHj0qT1Ok1xgI9JSFhQW+/vpr9O7dW1I/c+YMPD09UVhYKFNnZOpiY2MRFRUlqdnY2CAlJQXt2rWTqSvSZPr06QgICJDU8vPz4e7ujnPnzsnUFekrBgI9ZmNjg927d6N9+/aS+o8//oiRI0fyVCLSucTERMyYMUNSMzc3R2JiInr16iVTV/QskZGRGDVqlKR269YtuLi44MaNGzJ1RfqIgUDP1atXD0qlUm0zbHJyMmbOnMlTiUhnDhw4gJEjR6r9zK1ZswZubm4ydUXPY2ZmhnXr1mHAgAGSekZGBtzd3ZGXlydTZ6RvGAgMQMuWLZGamqp2oNaKFSvUNt0SacOpU6fg6emJoqIiSX3RokUYM2aMTF1RRVlaWiIpKQk9evSQ1I8fP46hQ4eq/b+SaWIgMBDdunXDjh07YGlpKakHBQXhq6++kqkrMgWXL1+Gi4sLcnJyJPVJkyYhODhYpq6osmxtbZGSkoLWrVtL6t999x3GjRsHlUolU2ekLxgIDEjfvn01fviPGzcOSqVSho7I2N25cwcuLi64fv26pO7p6YkVK1aondZG+q1hw4ZQKpVwcHCQ1Ddt2oSgoCCZuiJ9wUBgYD744APExMRIaiUlJfDy8sJ///tfmboiY1R2QayzZ89K6k5OTtiyZQvMzc1l6oyqol27dkhJSYGNjY2kHh0djdjYWJm6In3AQGCAZs+ejdmzZ0tqTZs2hYODAw8ypGohCAIsLS3RoUMHSd3R0RG7du1SuxEXGZZevXohKSlJEupsbGzQpk0bziEmjIHAQC1duhQ+Pj4AgNdeew2//PILGjRowE24VC0UCgUUCgXWrl2L0NBQAGLoVCqVqFu3rszdUXVwdXXF2rVrAYgXlUpPT4ezszPnEBNmIXcD9GLMzMwQHx+P1q1bY/78+bCwsICFRcX+O3/44QcUFhbyVDF6prIPhrCwMLRq1Qo9evRA8+bNZe6KqtPo0aNx//59uLq6olmzZpXaDcR5xPgohApsH8rNzYW9vT1ycnJgZ2eni76oggRBgCAIFb6rXG5uLjp16gSFQoGJEydi3rx5Wu6QjIEgCPzN0UgJggCVSlWpMMB5xHBU5vObuwwMnEKhqHAYEAQBdnZ22LVrFwYMGICFCxdi7ty5Wu6Q9E1xcTFu3LiBv//+Gw8ePKjQcxgGjJdCoahUGOA8YrwYCEyIQqFASUkJunXr9ujGSVu3buVtlU3I/v37MWLECHTs2BE9evSAt7c3Ll++LHdbZEA4jxgvBgITY2FhAaVSiZEjR6Jjx45YvHgxb0hjIg4dOoT3338fCoUC/v7+ePfdd3Ho0CFMmjRJ7aJDRM/CecQ48aBCE5OQkICxY8eiZ8+eCA0NhaurKwDuIzZ2165dw+jRozFixAgsWLAAtWvXBgDMmjULcXFxyMrKgr29vbxNksHgPGKcuIXAhMTExODDDz9Ev379EB0dzTexiVCpVIiPj0ejRo0wYcIE1K5d+9G55lOnTkVpaSkOHToEADwHnZ6L84jxYiAwQrdv38bnn38uudxsQEAA5syZg2HDhiEiIgJvvvkmAL6JTUFJSQkuXryIdu3awdHREcDjgwTr168PW1vbR8cR8GeBynAeMT3cZWCEzp8/j1mzZuH27duYMWMG5syZgw0bNmDChAmYOXPmo6vP8U1sGqysrBAREYHCwkIA4hYDMzMzlJaW4qWXXsI//vEP3L17V/KckydPokuXLnK0S3qC84jpYSAwQq+//jrWr1+PsWPHYsuWLbh69SqCgoIwYcIEtGjRAgDfxKamcePGj/5edpqqmZkZFAoFbG1tHwWCst0Hffr0wVdffYWRI0fK0i/Jj/OI6WEgMFKjR49Gbm4uZs6cicGDB2PixImSq8zxTUxlPwOWlpa4desWAODUqVMYMWIERo0axTBAnEdMDAOBEZs+fTquXLmC1atX4/jx42jUqBGsrKzkbov0RNmug3r16uHmzZu4fPkyfHx80LFjR2zcuFHu9khPcB4xHbx0sQkYMGAAsrOz8f333/P/jx4pCwSTJk1Ceno6rK2toVKpcPLkSblbIz3EecQwVebzm1sITMC+ffvw119/lfvDUPbBQMapvP28Zf/n9erVwx9//IEOHTrg+PHjum6PDMTz5hEyfPwUMBFlBwE9rbi4GHfv3kVKSoqOOyJduHr1KtLT06FSqcodM2DAALzyyis4cOAAatSoocPuyNBomkdUKhVUKhXWr1//zJ8z0n8MBCasuLgY165dw2uvvQZPT0+kpaXJ3RJVozt37sDZ2Rn9+/d/5jEBb7/9Ng4ePIgGDRrosDsyBiUlJSguLsaQIUMwbtw4BAUFyd0SVQEDgYlSqVQ4deoUevXqhYsXL6KkpARDhw7FsWPH5G6NqkFBQQEGDRqE33//HSUlJRg7diwWL15c7viaNWvqsDsyBoIgID8/H3379sXOnTsBANHR0Vi+fLnMndGLYiAwUWZmZkhMTMTNmzcf1e7fvw83NzdcvHhRxs6oqkpLSzFixAgcPHhQUt+yZQvu378vU1dkbBQKBU6cOIEjR45I6jNnzsTWrVtl6oqqgoHAhIWHh8PHx0dSy8rKgouLC7KysmTqiqpCEARMnz4d27dvl9SbNGkCpVKJl156SabOyBi9/fbbWLNmjaQmCAJGjRqFAwcOyNMUvTAGAhNmZmaG+Ph49O3bV1K/cOEC3N3dce/ePZk6oxcVGRmJlStXSmp2dnZQKpWSC8oQVZcxY8Zg0aJFklpRURE8PT1x6tQpmbqiF8FAYOKsra2xfft2dO3aVVL/7bff4O3tjeLiYpk6o8rasGEDgoODJTUrKyvs3LkTnTt3lqkrMgXBwcGYNGmSpJaTkwMXF5dHN84i/cdAQLCzs0NqaqraKUVKpRK+vr68Ja4B2Lt3L8aPHy+pKRQKJCQkoE+fPvI0RSZDoVBgxYoV8PT0lNSvX78OFxcX3LlzR57GqFIYCAiAePObtLQ01K1bV1LfuHEj5s+fL1NXVBFHjx6Ft7c3SktLJfVPP/0Uw4YNk6krMjXm5ubYsmULnJycJPWzZ8/Cw8MDBQUFMnVGFcVAQI+0b98ee/bsUTsFLSIiQm2/NOmH8+fPw93dHfn5+ZL6nDlzMGPGDJm6IlNVs2ZN7Nq1C46OjpL6oUOHMHz4cLXQSvqFgYAkevfujW+++UbtUsZTp07Ftm3bZOqKNMnMzISzs/OjOxWWGTFiBJYsWSJTV2Tq6tatC6VSiSZNmkjqO3bswLRp07gLUo8xEJAaDw8PxMXFSWqCIGD48OH46aefZOqKnpSXlwc3NzdkZGRI6v3798f69et5bwqSVfPmzaFUKtXue7Bq1SpERETI1BU9D2cN0sjX1xcLFiyQ1B48eAAPDw+cOXNGpq4IEE/p8vLyUrsRUffu3ZGcnMxb05Je6Ny5M3bu3Kn28zh//nzEx8fL1BU9CwMBlWvBggXw9fWV1LKzs+Hi4oKrV6/K1JVpEwQB48ePx759+yT1Vq1aISUlBba2tjJ1RqSuT58+2LRpk9rdNn19fbF3716ZuqLyMBBQuRQKBVauXAkPDw9J/erVq3B1dUV2drY8jZmwoKAgJCQkSGoODg5IS0tDo0aNZOqKqHze3t747LPPJLXS0lJ4e3vj6NGj8jRFGjEQ0DNZWFjg66+/Ru/evSX106dPY9CgQSgsLJSpM9MTGxuLqKgoSc3GxgYpKSlo166dTF0RPd/06dMREBAgqeXn58Pd3R3nz5+XqSt6GgMBPZeNjQ12796N9u3bS+o//vgjRo0axVOJdCAxMVHtNEJzc3MkJiaiV69eMnVFVHGRkZEYOXKkpHbr1i04OzsjMzNTpq7oSQwEVCH16tWDUqlU2yydlJSEWbNm8VQiLUpPT8fIkSPV/o3XrFkDNzc3mboiqhwzMzOsW7cOAwYMkNQzMjLg5uaGvLw8mTqjMgwEVGEtW7ZEamqq2oFrsbGxiI6Olqkr43bq1CkMGjQIRUVFkvqiRYswZswYmboiejFWVlZISkpC9+7dJfXjx4/Dy8tL7eecdIuBgCqlW7du2LFjBywtLSX1wMBAtYPdqGquXLkCV1dX5OTkSOqTJk1Su4kRkaGwtbVFSkoKWrVqJanv27cP48eP59ZGGTEQUKX17dsXX331lVp97NixSEtLk6Ej43Pnzh24uLjg2rVrkrqnpydWrFihdhoXkSFp1KgR0tLS4ODgIKknJCQgKChIpq6IgYBeyAcffICYmBhJraSkBEOHDsWxY8dk6so4FBQUYNCgQfj9998ldScnJ2zZsgXm5uYydUZUfdq1a4eUlBTY2NhI6lFRUYiNjZWpK9PGQEAvbPbs2Zg9e7akdv/+fbi5ueHixYsydWXYSktLMWLECBw8eFBSd3R0xK5du9RuPEVkyHr16oXExES1kDtjxgwkJibK1JXpYiCgKlm6dCl8fHwktaysLLi4uCArK0umrgyTIAiYPn06tm/fLqk3adIESqVS7dbURMbAzc0Na9askdQEQcDIkSORnp4uU1emiYGAqsTMzAzx8fHo27evpH7hwgW89957uHfvnkydGZ7IyEi120zb2dlBqVSiefPmMnVFpH1jxozBokWLJLWioiIMGjQIp06dkqkr08NAQFVmbW2N7du3o2vXrpL6r7/+imHDhqG4uFimzgzHhg0b1M4csLKyws6dO9G5c2eZuiLSneDgYEyaNElSy8nJgYuLCy5fvixTV6aFgYCqhZ2dHVJTU9GiRQtJPTU1FRMmTOCpRM+QmpqK8ePHS2oKhQIJCQno06ePPE0R6ZhCocCKFSvg6ekpqV+/fh0uLi64c+eOPI2ZEAYCqjaNGzdGWlqa2r7uDRs2ICQkRKau9Nuvv/4KLy8vtcs/f/rppxg2bJhMXRHJw9zcHFu2bIGTk5OkfvbsWXh4eKCgoECmzkwDAwFVq/bt22PPnj1qR8OHh4dj1apVMnWlny5cuAB3d3fk5+dL6nPmzFG7bwGRqahZsyZ27doFR0dHSf3QoUMYPnw4752iRQwEVO169+6Nb775BmZm0h+vKVOmqB1Bb6oyMzPh7OyMmzdvSuojRozAkiVLZOqKSD/UrVsXSqUSTZo0kdR37NiBadOmcRekljAQkFZ4eHggLi5OUhMEAT4+Pmrn2Juae/fuwd3dHX/++aek3r9/f6xfv14tSBGZoubNm0OpVMLOzk5SX7VqFSIiImTqyrhx5iGt8fX1xYIFCyS1Bw8eYODAgWpX4TMVxcXF8PLyUrua4yuvvILk5GRYWVnJ1BmR/uncuTN27typ9r6YP38+4uPjZerKeDEQkFYtWLAAvr6+klp2djZcXFxw9epVmbqShyAIGD9+vNr9Hlq1aoW9e/eq3UWSiIA+ffogISFB7f4dvr6+2Lt3r0xdGScGAtIqhUKBlStXYuDAgZJ62Z38srOz5WlMBvPmzVO7KZSDgwOUSiUaNWokU1dE+m/YsGH49NNPJbXS0lJ4e3vj6NGjMnVlfBgISOssLCzwzTff4PXXX5fUT58+DU9PTxQWFsrUme7ExsaqHSxoY2ODPXv24OWXX5apKyLDMWPGDMyZM0dSy8/Ph7u7O86fPy9TV8aFgYB0wsbGBrt371b78EtPT8eoUaOM+lSixMREtdMIzc3NsXXrVrz22msydUVkeJYsWYIRI0ZIardu3YKzszMyMzNl6sp4MBCQzjg4OCAtLU1t83hSUhJmzZpllKcSpaenY+TIkWrf25dffgl3d3eZuiIyTGZmZli/fj369+8vqWdkZMDNzQ15eXkydWYcGAhIp1q2bInU1FS1A+hiY2MRHR0tU1facerUKQwaNAhFRUWS+scff4yxY8fK1BWRYbOyskJycjJeeeUVSf348ePw8vJSe79RxTEQkM5169YN27dvh6WlpaQeGBiIhIQEmbqqXmUHTebk5EjqEydOxPz582Xqisg42NraYu/evWjVqpWkvm/fPowfP94otzbqAgMByaJfv37YuHGjWn3s2LFqp+UZmjt37sDFxQXXrl2T1D09PfHFF1+onT5FRJXXqFEjKJVKODg4SOoJCQkICgqSqSvDxkBAsvHx8UFMTIykVlJSgqFDh6pduMdQFBQUYNCgQWoXXnJycsKWLVtgbm4uU2dExufll19GSkoKbGxsJPWoqCjExsbK1JXhYiAgWc2ePRuzZ8+W1O7fvw83NzdcvHhRpq5eTGlpKUaMGKF2aWZHR0fs2rVL7YZPRFR1vXr1QmJiolrYnjFjBhITE2XqyjAxEJDsli5dCh8fH0ktKysLLi4uyMrKkqmryhEEAdOnT1e7eVOTJk2gVCrVbglNRNXHzc0Na9askdQEQcDIkSORnp4uU1eGh4GAZGdmZob4+Hj07dtXUr9w4QLee+893Lt3T6bOKi4yMhIrV66U1Ozs7KBUKtG8eXOZuiIyHWPGjMGiRYsktaKiIgwaNAinTp2SqSvDwkBAesHa2hrbt29H165dJfVff/0Vw4YNQ3FxsUydPd+GDRsQHBwsqVlZWWHnzp3o3LmzTF0RmZ7g4GBMmjRJUsvJyYGrqyuuXLkiU1eGg4GA9IadnR1SU1PRokULST01NRUTJkzQy1OJUlNTMX78eElNoVAgISEBffr0kacpIhOlUCiwYsUKeHp6SurXrl2Di4sL7ty5I09jBoKBgPRK48aNkZaWprbPfcOGDQgJCZGpK81+/fVXeHl5qV12+dNPP8WwYcNk6orItJmbm2PLli1wcnKS1H///XcMGjQIBQUFMnWm/xgISO+0b98ee/bsUTsqPzw8HKtWrZKpK6kLFy7A3d0d+fn5kvqcOXPU7ltARLpVs2ZN7Nq1C46OjpL6wYMHMWLECKO+d0pVMBCQXurduze++eYbmJlJf0SnTJmidiS/rmVmZsLZ2Rk3b96U1EeMGKF2R0MikkfdunWhVCrRpEkTSX379u2YPn26Xu6ClBsDAektDw8PxMXFSWqCIMDHxwcHDx5EVlYWwsPD8fnnn6v9pl6dUlNTMW/ePBw+fBj37t2Du7s7/vzzT8mY/v37Y/369WoBhojk07x5cyiVStjZ2UnqK1euRGRkJIqKirBmzRosXLgQly9flqlL/aEQKhCTcnNzYW9vj5ycHLV/WCJtW7hwIcLCwiS1WrVqoUaNGrh16xYAYOjQoUhKSqr2de/YsQODBw9+9Pjll1/GuXPnJGNeeeUVpKenq92wiYj0w4EDB+Ds7Kx246P27dvjjz/+ACBeCvncuXNG9z6uzOc3f50hvbdgwQL4+vpKavfu3XsUBgAgOTlZKxcxevqYhafDQKtWrbB3716jm0SIjEmfPn2QkJCgdh+RsjAAADdu3MCuXbt03ZpeYSAgvadQKLBy5UoMHDjwmeO+//77al1vQUEBfvzxx3KX16tXD0qlEo0aNarW9RJR9Rs2bBg+/fTTZ47Zt2+fjrrRTwwEZBAsLCzw9ttvP3NMdb+ZDx48iMLCwnKXN2/eHM2aNavWdRKR9vTu3RvW1tblLv/uu+9M+mBDC7kbIKqIjIwMzJ0795ljyt7M6rcXLgZwEsAxAMcB/A3gAQBrAI0BdAfQA0AXAJaS13uWEydO4PPPP0dgYGClvhciksfo0aPx4MGDcpf//fffOHPmDDp16vTUkhebQwwNAwEZhBs3bjz33OHr168/9Wb+C8BqAHEA7j6sWUJ8c+OJx6sf/r0OgEkAJgJoUaHTG69evVrRb4GIZHbt2rXnjtm9e3e1zSGGhrsMyCC89tpr6Nev33PHifdAzwHgC6AVgGg8fiMD0jfy04/vPhzfCvfu+SAz88Iz11WvXj1MmDDh+c0TkV6YN2/ec8esXr0a1TGHiM/PrUK3usdAQAbBzMwMaWlp2Lp1K/r371/uOAuLHwB0ABAPQABQ2SuSlQIQYGPzLf73P0DTqlq3bo2IiAj8/vvv6NKlSyVfn4jkMnfuXBw5cgRjxoyBjY2NxjGdOl1Hdcwh4vM7ADCcAxV5HQIySH/++SfWrVuHlStXIjs7GwAwZQqwYgUg5lxVlddRWgqYm4uvu2qVAl5eXpgwYQL69u3LCxARGbicnBxs2bIFsbGxOHv2LIDqn0Mev84KAFOq4fUqrzKf3wwEZNCKi4vxxRdfQKFYiRkzzmttPbm5EbCzC9La6xORPARBwH/+8x/88cd0TJlyVotrkicUMBCQidkHwFkH60kDMEAH6yEi3TLeOYRXKiQTkgPgQ2j/R9kMwGgY2kFCRPQ8nEPKMBCQgfMHcBPVs7/vWVQAsgD4aXk9RKRbnEPKMBCQAbsEYB02bCiFQoFHXzVqAI0aAe+8A0RGAuXd4uDOHcDSEii73EBwMPDKK0DduuJrtG4NTJgA/PVX2TNKAayDeG4yERm+SxDf06XYsAGSeeTprwMHpM9cvhyoVw8oKVF/1e+/f/y8J265An2fQ3hhIjJgX0LMtOJpQfHxQIcOQHGxGAIOHgSiooBly4BvvwXefVf67J07ASsrwMVFfJydDfj4AI6OgK0t8PvvwOLFwK5dwJkz4ptfXN+XAMJ19U0SkdZI5xDg8TzytI4dpY+Tk4FBgwCLpz5F790DfH2BJk2A69c1rVN/5xAGAjJQxRCvHvb4jdypE9Cz5+MRQ4cCs2YBb7wBDBkCnD8PNGz4eHlSEuDmBtSsKT7+4gvpGvr0AVq1Esfs3AmMHYuH61sFYCEM+RKlRKQ+hwDq84gmmZniLxwBAerLAgOBOnUAd3fxFwp1+juHcJcBGaiTkF49TLPmzYGYGCAvD1i9+nE9N1fcrDd06LOfX7+++Kf0t4C7AE5Vsl8i0i8Vm0M02b4dqFVLfavjTz8BX34JrF0rXsOkfPo5hzAQkIE6VuGRbm7im/PJOxnv3i3u33N3Vx9fUgIUFAAnTgAzZwIvvyxuYXjR9RORPtL8Hi4tFeeAJ7+evo1KcjLw3nvAkzdOLCgAxo0T54zu3V98/XJiICADdRwV3dz20kuAg4N0f15SEjBggHiswJNu3BAPNLSxEd/UJSXA/v3ibwOPWUIf38xEVBma55DXXxfngCe/nvzgv31bPMDw6a2LISFicAgLq8i69XMO4TEEZKD+hvpNRsr35OW37t8H0tKAVavUxzk4AL/+Cjx4AJw9C0RHi2crHDgANG5cNqoYwI0Xb52I9IDmOeSrr8QDi5/05B3Vnz4YGQCOHgU++wxQKh8fk/Rs+jmHMBCQgSr/nuZPu39fTPWdO4uPU1LEMxE8PNTHWlg8PqDIyUl807dqBSxZAnz++ZMjC1+4cyLSB5rnEEfHZx9UmJQEuLqKWxHLjB0r7lbs2VM8WwkACh9OEbm54haGp7dG6uMcwkBABsr6+UMeSkkRN+X16SM+Tk4G+vYVjwR+nn/8Qzx96Ny5p5fUqPD6iUgfVXwOKZOTA/znP8CGDdL6mTPiV2Ki+nPatAG6dgX+7/+eXqJ/cwgDARmoxhD3wz17t8Hly4C/P2BvD0ycKKb2vXvFMw8q4sIF4OrVp7cmWAJo9GJtE5GeqNgc8qSyg5Hfe09a379ffeyGDcDGjcCOHUDTpk8v1c85hIGADFR3AKslldOnHx8VnJUlngIUHy+eYbB9u3gK4Y4dQH4+4OkpfbWTJ8VrFnh5iVcoNDMDTp0CPv1UvCCRv/+To4sB9NDqd0dE2qY+hwCP55GntWkj7i7o319983/Z1scnlV3Z0MlJPDZJSj/nEAYCMlDqb6YxY8Q/rayA2rXFfYFz5wLjxz++nkByMvDmm0CDBtLnNmwo7hqIiQH+/lucEP7xD/E3gXnzgGbNnr9+IjIkmt/DZfPI0z77rPyDkatz/XLi7Y/JQBUDaIjKXFikqEgMAosWAdOmVWXddQBkQt+uMkZElVG5OWTrVmDECPEqhXXrVnXduptDePtjMgGWACYBeOblwCSsrMQjgKsWBswBfASGASJDV7k5ZNgw8eykqocB/Z1DGAjIgE2E9m9Z+jQVgAk6XicRaQfnkCcxEJABawFgHCqzlaBqzB+ur4WO1kdE2sU55EkMBGTgYgA0gPZ/lM0erqeC5ysSkYHgHFKGgYAMnB2ADdD+Zj/Vw/XwoFoi48I5pAwDARmBAQBWaHkdXzxcDxEZH84hAAMBGY0pePyGrq4f67LX+QLA5Gp6TSLST5xDGAjIiEwBkAbx3OKqHiRk/vB10mAIb2Qiqg6mPYcwEJCRGQDgLIAxABSo/Jva/OHzxgD4H/R9Ex8RVTfTnUMYCMgI2QNYAyADwFyIVwUr8/TFQJ58XOfh+IyHz9ffg3+ISJtMcw7hpYvJBBQDOAXg2MOvGxDvRV4D4h3Hejz86gx9vHoYEcnNcOeQynx+MxAQEREZKd7LgIiIiCqFgYCIiIgYCIiIiIiBgIiIiMBAQERERGAgICIiIjAQEBERERgIiIiICAwEREREBAYCIiIiAgMBERERgYGAiIiIwEBAREREACwqMqjshoi5ublabYaIiIiqT9nndgVubFyxQJCXlwcAaNasWRXaIiIiIjnk5eXB3t7+mWMUQgVig0qlwvXr12FrawuFQlFtDRIREZH2CIKAvLw8NGnSBGZmzz5KoEKBgIiIiIwbDyokIiIiBgIiIiJiICAiIiIwEBAREREYCIiIiAgMBERERAQGAiIiIgLw/8K/LqU2L3c7AAAAAElFTkSuQmCC\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": [
    "randVar = ['A', 'B', 'C', 'D', 'E']\n",
    "values = [0, 1, 2, 3, 4]\n",
    "ltn = dict(zip(randVar, values))\n",
    "\n",
    "# graphical model construction\n",
    "G = nx.DiGraph()\n",
    "G.add_edges_from([('A/0', 'B/1')], weight=0.3)\n",
    "G.add_edges_from([('A/0', 'C/2')], weight=2)\n",
    "G.add_edges_from([('B/1', 'D/3')], weight=1)\n",
    "G.add_edges_from([('C/2', 'D/3')], weight=2)\n",
    "G.add_edges_from([('C/2', 'E/4')], weight=1)\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=500)\n",
    "\n",
    "# edges\n",
    "nx.draw_networkx_edges(G, pos, width=3)\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()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "89aa9ff1",
   "metadata": {},
   "source": [
    "### 1.2 Covariance Matrix Construction from Conditional Variances."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "fde6ca80",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "A:  [[0.  0.  0.  0.  0. ]\n",
      " [0.3 0.  0.  0.  0. ]\n",
      " [2.  0.  0.  0.  0. ]\n",
      " [0.  1.  2.  0.  0. ]\n",
      " [0.  0.  1.  0.  0. ]]\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([[ 0.2, -0. , -0. , -0. , -0. ],\n",
       "       [-0. ,  0.4, -0. ,  0. ,  0. ],\n",
       "       [-0. , -0. ,  0.3, -0. , -0. ],\n",
       "       [-0. ,  0. , -0. ,  0.3, -0. ],\n",
       "       [-0. ,  0. , -0. , -0. ,  0.2]])"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# constructing S (conditional variance of each node given parents)\n",
    "condVariance = [0.2, 0.4, 0.3, 0.3, 0.2]\n",
    "S = np.diag(condVariance)\n",
    "\n",
    "# constructing A (function of weights)\n",
    "# scope for more automation based on the GSO of the graphical model, e.g., can look at it's adjacency matrix.\n",
    "# build the graphical relation first.\n",
    "# then extract the weight matrix from that graph to develop A.\n",
    "A = nx.adjacency_matrix(G)\n",
    "A = A.todense()\n",
    "A = A.transpose()\n",
    "print('A: ', A)\n",
    "\n",
    "I = np.identity(len(values))\n",
    "temp = inv(I - A)\n",
    "cov = temp*S*(temp.transpose())\n",
    "cov"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4f256c36",
   "metadata": {},
   "source": [
    "### 1.3 Mean Vector Construction from Conditional Biases and Edge Weights."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "11cd4e75",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.4 , 0.12, 1.  , 2.12, 2.  ])"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# constructing Mean vector from CPDs. \n",
    "condBias = [0.4, 0, 0.2, 0, 1]\n",
    "condBias = np.array(condBias)\n",
    "mu = np.array( [-1.1]*len(condBias) )\n",
    "\n",
    "for rowIndex, _ in enumerate(condBias):\n",
    "    # find parents of \"rowIndex\"\n",
    "    nonZeroIndices = np.where( A[rowIndex, :]>0 )[0]\n",
    "    parentWeights = np.array(A[rowIndex, nonZeroIndices])\n",
    "    if len(nonZeroIndices)>0:\n",
    "        # has parents\n",
    "        temp = np.dot(np.squeeze(parentWeights), mu[nonZeroIndices])\n",
    "        mu[rowIndex] = condBias[rowIndex] + temp.item()\n",
    "    else:\n",
    "        # no parent\n",
    "        mu[rowIndex] = condBias[rowIndex]   \n",
    "    #print(mu)\n",
    "    #print(\"\")\n",
    "    \n",
    "mu"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7a056d5e",
   "metadata": {},
   "source": [
    "### 1.4 Simulated Dataset Generation by sampling from the induced Joint Gaussian Distribution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "8af3f300",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(5, 1000)"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# draw samples from the joint distribution (happens to be Multi variate gaussian)\n",
    "numSamples = 1000\n",
    "data = np.random.multivariate_normal(mu, cov, numSamples).T\n",
    "data.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bc6abb7b",
   "metadata": {},
   "source": [
    "## 2. [Learning] Parameter Learning of Each CPD from Samples Generated Above:\n",
    "\n",
    "- Inferring the conditional weights, biases and variances from data can be casted as a MLE problem.\n",
    "- When the data is complete (no latent variables), the MLE decomposes into a set of local learning problems, one for each random variable.\n",
    "- Simplifying the parameter learning for each CPD results in the following three key equations (rendering one each when differentiating w.r.t condtional bias, weights and condtional variance):\n",
    "    $$\n",
    "    \\begin{align}\n",
    "    E_{D}[X_{child}] &= b + \\sum_{i \\in U_{parents}} w_{i} \\cdot E_{D}[u_{i}]\\ , \\\\\n",
    "    E_{D}[X_{child} \\cdot u_{i}] &= b \\cdot E_{D}[u_{i}] + \\sum_{j \\in U_{parents}} w_{j} \\cdot E_{D}[u_{j} \\cdot u_{i}]\\ ,\\\\\n",
    "    \\sigma^{2} &= \\mbox{Cov}_{D}[X; X] - \\sum_{i} \\sum_{j} w_{i}\\cdot w_{j} \\cdot \\mbox{Cov}_{D}[u_{i}; u_{j}]\\ .\n",
    "    \\end{align}\n",
    "    $$\n",
    "- The first two system of linear equations can learn the biases and weights asscociated with each random variables.\n",
    "- After which, they can be leveraged to learn conditional variances $\\sigma^{2}$ for each CPD separately."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "205ba4c5",
   "metadata": {},
   "source": [
    "### 2.1 Parameterizing the Conditional Biases, Weights and Variances."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "59f28d2b",
   "metadata": {},
   "outputs": [],
   "source": [
    "learnBias = [None]*len(values)\n",
    "learnVariance = [None]*len(values)\n",
    "learnWeights = {}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "679a0829",
   "metadata": {},
   "source": [
    "### 2.2 Computing Essential Singleton and Pairwise Sufficient Statistics for Plug-in Parameter Estimation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "ef383d71",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.43151081 0.12565589 0.98224119 2.12450314 2.02905367]\n",
      "[0.36901627 0.43463116 1.27121347 4.8039744  4.32786599]\n"
     ]
    }
   ],
   "source": [
    "expectedX = np.average(data, axis=1)           #Singleton expected value for each RV.\n",
    "print(expectedX)\n",
    "expectedXSquared = np.average(data**2, axis=1) #Singleton expected value for each RV^{2}.\n",
    "print(expectedXSquared)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "39035870",
   "metadata": {},
   "source": [
    "###  2.3 Formulating the System of Linear Equations for Each Random Variable:\n",
    "- Simple algorithm for computing the unknown parameter:\n",
    "    - if the RV has parents, then solve a system of linear equations by dynamically constructing the first two equations above.\n",
    "        - To find out the number parents for the current node, find the number of non-zero elements in the $A$ matrix assocaited with the current node row.\n",
    "    - else, compute the mean and variance only."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "65b2b561",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "currNode:  0\n",
      "\n",
      "currNode:  1\n",
      "parent:          0\n",
      "\n",
      "currNode:  2\n",
      "parent:          0\n",
      "\n",
      "currNode:  3\n",
      "parents:         1   2\n",
      "\n",
      "currNode:  4\n",
      "parent:          2\n",
      "\n",
      "learned Bias:      [0.4315108090450466, 0.1291064827329677, 0.9831834927706627, 2.208173925256592, 1.9851544643661425]\n",
      "learned Variance:  [0.18281468693128786, 0.4188300671313852, 0.3064148478235934, 0.2882649482144718, 0.2101951275647364]\n",
      "learned Weights:   {'W_(1, 0)': -0.007996545313692333, 'W_(2, 0)': -0.002183734874086333, 'W_(3, 1)': -0.006745460846322168, 'W_(3, 2)': -0.08432060995473152, 'W_(4, 2)': 0.04469290323335007}\n"
     ]
    }
   ],
   "source": [
    "for currNode in range(len(values)):\n",
    "    print(\"currNode: \", currNode)\n",
    "    numParents = len( np.nonzero(A[currNode, :])[0] )\n",
    "    if numParents:\n",
    "        if numParents==1:\n",
    "            parent = np.nonzero(A[currNode, :])[0].item()\n",
    "            print(\"parent:         \", parent)\n",
    "            left = np.array( [ [1, expectedX[parent]], \n",
    "                               [expectedX[parent], expectedXSquared[parent]] ] \n",
    "                           )\n",
    "            right = np.array( [ expectedX[currNode], \n",
    "                                np.dot(data[currNode, :], data[parent,:])/numSamples]\n",
    "                            )\n",
    "            variables = linalg.solve(left, right)\n",
    "            learnBias[currNode] = variables[0]\n",
    "            #print(f\"W_{currNode, parent}\")\n",
    "            learnWeights[f\"W_{currNode, parent}\"] = variables[1]\n",
    "            #print(learnBias, learnWeights)\n",
    "            learnVariance[currNode] = expectedXSquared[currNode] - (expectedX[currNode])**2 - (learnWeights[f\"W_{currNode, parent}\"]**2)*(expectedXSquared[parent] - (expectedX[parent])**2)\n",
    "            #learnVariance\n",
    "        elif numParents==2:\n",
    "            parent1, parent2 = np.nonzero(A[currNode, :])[0]\n",
    "            print(\"parents:        \", parent1, \" \",  parent2)\n",
    "            left = np.array( [ [1, expectedX[parent1], expectedX[parent2]], \n",
    "                               [expectedX[parent1], expectedXSquared[parent1], np.dot(data[parent1, :], data[parent2,:])/numSamples], \n",
    "                               [expectedX[parent2], np.dot(data[parent2, :], data[parent1,:])/numSamples, expectedXSquared[parent2]] ] )\n",
    "            right = np.array( [expectedX[currNode], \n",
    "                               np.dot(data[currNode, :], data[parent1,:])/numSamples, \n",
    "                               np.dot(data[currNode, :], data[parent2,:])/numSamples ] )\n",
    "            variables = linalg.solve(left, right)\n",
    "            learnBias[currNode] = variables[0]\n",
    "            learnWeights[f\"W_{currNode, parent1}\"] = variables[1]\n",
    "            learnWeights[f\"W_{currNode, parent2}\"] = variables[2]\n",
    "            #print(learnBias, learnWeights)\n",
    "            learnVariance[currNode] = expectedXSquared[currNode] - (expectedX[currNode])**2 - 2*((learnWeights[f\"W_{currNode, parent1}\"]*learnWeights[f\"W_{currNode, parent2}\"])*(np.dot(data[parent1, :], data[parent2,:])/numSamples - expectedX[parent1]*expectedX[parent2])) -  (learnWeights[f\"W_{currNode, parent1}\"]**2)*(expectedXSquared[parent1] - (expectedX[parent1])**2) -  (learnWeights[f\"W_{currNode, parent2}\"]**2)*(expectedXSquared[parent2] - (expectedX[parent2])**2)\n",
    "            #learnVariance\n",
    "    else:\n",
    "        # no parents\n",
    "        learnBias[currNode] = expectedX[currNode]\n",
    "        learnVariance[currNode] = expectedXSquared[currNode] - (expectedX[currNode])**2\n",
    "    print(\"\")\n",
    "        \n",
    "print(\"learned Bias:     \", learnBias)\n",
    "print(\"learned Variance: \", learnVariance)\n",
    "print(\"learned Weights:  \", learnWeights)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0eebee87",
   "metadata": {},
   "source": [
    "## 3. [Inference] Computing Posteriors: Inference in Graphical models\n",
    "\n",
    "In order to compute any posterior in this model, we can simply leverage the property that all the RVs jointly induce a multivariate Gaussian density. Hence, any posterior ditribution of the form $\\mathbb{P}(X_{1} \\vert X_{2} = a)$ can be expressed as:\n",
    "$$\n",
    "\\begin{align}\n",
    "\\mathbb{P}(X_{1} \\vert X_{2} = a) &\\sim \\mathcal{N}(\\bar \\mu, \\bar \\Sigma), \\\\\n",
    "\\bar \\mu & = \\mu_{1} + \\Sigma_{12} \\cdot \\Sigma^{-1}_{22} \\cdot (a - \\mu_{2})\\ ,\\\\\n",
    "\\bar \\Sigma & = \\Sigma_{11} - \\Sigma_{12}\\cdot \\Sigma^{-1}_{22}\\cdot \\Sigma_{21}\\ .\n",
    "\\end{align}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "cd01d2d8",
   "metadata": {},
   "outputs": [],
   "source": [
    "def moment_estimation(weights, condVariance, condBias):\n",
    "    \"\"\"\n",
    "    Estimates the Mean and Covariance matrix of the joint Gaussian distribution given CPD parameters.\n",
    "    \"\"\"\n",
    "    G = nx.DiGraph()\n",
    "    G.add_edges_from([('A/0', 'B/1')], weight=weights[0])\n",
    "    G.add_edges_from([('A/0', 'C/2')], weight=weights[1])\n",
    "    G.add_edges_from([('B/1', 'D/3')], weight=weights[2])\n",
    "    G.add_edges_from([('C/2', 'D/3')], weight=weights[3])\n",
    "    G.add_edges_from([('C/2', 'E/4')], weight=weights[4])\n",
    "\n",
    "    S = np.diag(condVariance)\n",
    "    # constructing A (function of weights)\n",
    "    # scope for more automation based on the GSO of the graphical model, e.g., can look at it's adjacency matrix.\n",
    "    # build the graphical relation first.\n",
    "    # then extract the weight matrix from that graph to develop A.\n",
    "    A = nx.adjacency_matrix(G)\n",
    "    A = A.todense()\n",
    "    A = A.transpose()\n",
    "    #print('A: ', A)\n",
    "    I = np.identity(len(values))\n",
    "    temp = inv(I - A)\n",
    "    cov = temp*S*(temp.transpose())\n",
    "\n",
    "    condBias = np.array(condBias)\n",
    "    mu = np.array( [-1.1]*len(condBias) )\n",
    "\n",
    "    for rowIndex, _ in enumerate(condBias):\n",
    "        # find parents of \"rowIndex\"\n",
    "        nonZeroIndices = np.where( A[rowIndex, :]>0 )[0]\n",
    "        parentWeights = np.array(A[rowIndex, nonZeroIndices])\n",
    "        if len(nonZeroIndices)>0:\n",
    "            # has parents\n",
    "            temp = np.dot(np.squeeze(parentWeights), mu[nonZeroIndices])\n",
    "            mu[rowIndex] = condBias[rowIndex] + temp.item()\n",
    "        else:\n",
    "            # no parent\n",
    "            mu[rowIndex] = condBias[rowIndex]   \n",
    "            \n",
    "    return mu, cov"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "e724f069",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Mean:  [0.4  0.12 1.   2.12 2.  ]\n",
      "Covariance:  [[ 0.2 -0.  -0.  -0.  -0. ]\n",
      " [-0.   0.4 -0.   0.   0. ]\n",
      " [-0.  -0.   0.3 -0.  -0. ]\n",
      " [-0.   0.  -0.   0.3 -0. ]\n",
      " [-0.   0.  -0.  -0.   0.2]]\n"
     ]
    }
   ],
   "source": [
    "# True Moments\n",
    "condBias = [0.4, 0, 0.2, 0, 1]\n",
    "condVariance = [0.2, 0.4, 0.3, 0.3, 0.2]\n",
    "weights = [0.3, 2, 1, 2, 1]\n",
    "mu, cov = moment_estimation(weights, condVariance, condBias)\n",
    "print(\"Mean: \", mu)\n",
    "print(\"Covariance: \", cov)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "1a8efc46",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Mean:  [0.43151081 0.12910648 0.98318349 2.20817393 2.02909579]\n",
      "Covariance:  [[ 0.18281469 -0.         -0.          0.         -0.        ]\n",
      " [-0.          0.41883007  0.         -0.          0.        ]\n",
      " [-0.          0.          0.30641485 -0.          0.        ]\n",
      " [ 0.         -0.         -0.          0.28826495  0.        ]\n",
      " [-0.          0.          0.          0.          0.21019513]]\n"
     ]
    }
   ],
   "source": [
    "# Learned Moments\n",
    "learnWeights_list = [learnWeights['W_(1, 0)'], learnWeights['W_(2, 0)'], learnWeights['W_(3, 1)'], learnWeights['W_(3, 2)'], learnWeights['W_(4, 2)']]\n",
    "muhat, covhat = moment_estimation(learnWeights_list, learnVariance, learnBias)\n",
    "print(\"Mean: \", muhat)\n",
    "print(\"Covariance: \", covhat)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "90e065b1",
   "metadata": {},
   "source": [
    "Compute the posterior distribution:\n",
    "$$ \n",
    "\\begin{align}\n",
    "\\mathbb{P}(A \\vert D=1, E=2) = \\int_{B} \\int_{C} \\mathbb{P}(A, B, C \\vert D=1, E=2)\\ ,\\\\\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "Once the above conditional distribution $\\mathbb{P}(A, B, C \\vert D=1, E=2)$ is determined, the required posterior can be determined instantly. Leveraging the earlier mentioned property of joint Gaussians, we can compute the conditional distribution by assigning $X_{1}$ and $X_{2}$ as follows:\n",
    "$$\n",
    "\\begin{align}\n",
    "X_{1} \\in \\mathbb{R}^{3} & \\quad X_{1} = [A, B, C]^{\\top}\\ ,\\\\\n",
    "X_{2} \\in \\mathbb{R}^{2} & \\quad X_{2} = [D, E]^{\\top} \\quad a = [1, 2]^{\\top}\\ .\n",
    "\\end{align}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "592fdf26",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Conditional Mean:  [[0.43151081]\n",
      " [0.12910648]\n",
      " [0.98318349]]\n",
      "Conditional Covariance:  [[ 0.18281469 -0.         -0.        ]\n",
      " [-0.          0.41883007  0.        ]\n",
      " [-0.          0.          0.30641485]]\n"
     ]
    }
   ],
   "source": [
    "a = np.array([1, 2]) # evidence/observation\n",
    "muBar = muhat[:3].reshape((3,1)) + np.dot(np.dot(covhat[:3, 3:], np.linalg.inv(covhat[3:, 3:])), (a - muhat[3:]).reshape((2,1)))\n",
    "sigmaBar = covhat[:3, :3] - np.dot(np.dot(covhat[:3, 3:], np.linalg.inv(covhat[3:, 3:])), covhat[3:, :3])\n",
    "\n",
    "print(\"Conditional Mean: \", muBar)\n",
    "print(\"Conditional Covariance: \", sigmaBar)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "25d85d5c",
   "metadata": {},
   "source": [
    "Hence, the posterior $\\mathbb{P}(A \\vert D=1, E=2)$ is a Gaussian denssity with Mean: muBar[0] and Variance: sigmaBar[0]."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "98f76887",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Posterior Mean:  0.4315108090450466\n",
      "Posterior Covariance:  0.18281468693128786\n"
     ]
    }
   ],
   "source": [
    "print(\"Posterior Mean: \", muBar[0].item())\n",
    "print(\"Posterior Covariance: \", sigmaBar[0, 0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f4568199",
   "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
}