examples/BootstrapWithImputation.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "8d062890-891b-4008-9f53-259aa2c6d03a",
   "metadata": {},
   "source": [
    "# Bootstrap with imputation\n",
    "\n",
    "This notebook explains how to use `bootstrap_with_imputation`. `bootstrap_with_imputation` performs causal discovery with multiple imputations on data with missing values.\n",
    "\n",
    "Hereinafter, NaN or nan represents missing values."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d0bfe5f4-6474-4a37-a04d-9a5a6a52cfc8",
   "metadata": {},
   "source": [
    "## Import and settings"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "5ead87c6-e0a7-4838-af1b-ed695937065f",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['1.24.4', '1.9.0']\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "\n",
    "import lingam\n",
    "from lingam.utils import make_dot\n",
    "from lingam.tools import bootstrap_with_imputation\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "print([np.__version__, lingam.__version__])\n",
    "\n",
    "np.set_printoptions(precision=3, suppress=True)\n",
    "np.random.seed(0)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a68b6a12-2d4b-4b9f-8d22-57dc7cb7c218",
   "metadata": {},
   "source": [
    "## Test data\n",
    "\n",
    "Test data is generated according to the following adjacency matrix."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "c8b47f1c-657b-4fb1-a05b-e1235dca9672",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "image/svg+xml": [
       "<?xml version=\"1.0\" encoding=\"UTF-8\" standalone=\"no\"?>\n",
       "<!DOCTYPE svg PUBLIC \"-//W3C//DTD SVG 1.1//EN\"\n",
       " \"http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd\">\n",
       "<!-- Generated by graphviz version 2.43.0 (0)\n",
       " -->\n",
       "<!-- Title: %3 Pages: 1 -->\n",
       "<svg width=\"267pt\" height=\"218pt\"\n",
       " viewBox=\"0.00 0.00 267.00 218.00\" xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\">\n",
       "<g id=\"graph0\" class=\"graph\" transform=\"scale(1 1) rotate(0) translate(4 214)\">\n",
       "<title>%3</title>\n",
       "<polygon fill=\"white\" stroke=\"transparent\" points=\"-4,4 -4,-214 263,-214 263,4 -4,4\"/>\n",
       "<!-- x0 -->\n",
       "<g id=\"node1\" class=\"node\">\n",
       "<title>x0</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"158\" cy=\"-105\" rx=\"27\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"158\" y=\"-101.3\" font-family=\"Times,serif\" font-size=\"14.00\">x0</text>\n",
       "</g>\n",
       "<!-- x1 -->\n",
       "<g id=\"node2\" class=\"node\">\n",
       "<title>x1</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"139\" cy=\"-18\" rx=\"27\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"139\" y=\"-14.3\" font-family=\"Times,serif\" font-size=\"14.00\">x1</text>\n",
       "</g>\n",
       "<!-- x0&#45;&gt;x1 -->\n",
       "<g id=\"edge2\" class=\"edge\">\n",
       "<title>x0&#45;&gt;x1</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M158.52,-86.64C158.42,-76.91 157.65,-64.63 155,-54 154.23,-50.92 153.2,-47.78 152.03,-44.71\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"155.13,-43.07 147.95,-35.28 148.7,-45.85 155.13,-43.07\"/>\n",
       "<text text-anchor=\"middle\" x=\"173\" y=\"-57.8\" font-family=\"Times,serif\" font-size=\"14.00\">0.56</text>\n",
       "</g>\n",
       "<!-- x4 -->\n",
       "<g id=\"node5\" class=\"node\">\n",
       "<title>x4</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"27\" cy=\"-18\" rx=\"27\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"27\" y=\"-14.3\" font-family=\"Times,serif\" font-size=\"14.00\">x4</text>\n",
       "</g>\n",
       "<!-- x0&#45;&gt;x4 -->\n",
       "<g id=\"edge5\" class=\"edge\">\n",
       "<title>x0&#45;&gt;x4</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M142.61,-89.91C128.57,-77.22 108.73,-59.71 100,-54 87.13,-45.58 71.92,-37.96 58.73,-32\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"60.1,-28.78 49.54,-27.96 57.29,-35.19 60.1,-28.78\"/>\n",
       "<text text-anchor=\"middle\" x=\"134\" y=\"-57.8\" font-family=\"Times,serif\" font-size=\"14.00\">0.99</text>\n",
       "</g>\n",
       "<!-- x5 -->\n",
       "<g id=\"node6\" class=\"node\">\n",
       "<title>x5</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"232\" cy=\"-18\" rx=\"27\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"232\" y=\"-14.3\" font-family=\"Times,serif\" font-size=\"14.00\">x5</text>\n",
       "</g>\n",
       "<!-- x0&#45;&gt;x5 -->\n",
       "<g id=\"edge7\" class=\"edge\">\n",
       "<title>x0&#45;&gt;x5</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M173.06,-89.74C179.4,-83.54 186.74,-76.08 193,-69 200.33,-60.71 207.91,-51.22 214.42,-42.76\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"217.27,-44.8 220.53,-34.71 211.69,-40.56 217.27,-44.8\"/>\n",
       "<text text-anchor=\"middle\" x=\"221\" y=\"-57.8\" font-family=\"Times,serif\" font-size=\"14.00\">0.90</text>\n",
       "</g>\n",
       "<!-- x2 -->\n",
       "<g id=\"node3\" class=\"node\">\n",
       "<title>x2</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"67\" cy=\"-105\" rx=\"27\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"67\" y=\"-101.3\" font-family=\"Times,serif\" font-size=\"14.00\">x2</text>\n",
       "</g>\n",
       "<!-- x2&#45;&gt;x1 -->\n",
       "<g id=\"edge3\" class=\"edge\">\n",
       "<title>x2&#45;&gt;x1</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M63.15,-86.98C61.78,-76.66 61.85,-63.67 68,-54 76.07,-41.31 90.22,-33.07 103.77,-27.8\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"105.12,-31.04 113.43,-24.47 102.84,-24.42 105.12,-31.04\"/>\n",
       "<text text-anchor=\"middle\" x=\"84\" y=\"-57.8\" font-family=\"Times,serif\" font-size=\"14.00\">0.38</text>\n",
       "</g>\n",
       "<!-- x2&#45;&gt;x4 -->\n",
       "<g id=\"edge6\" class=\"edge\">\n",
       "<title>x2&#45;&gt;x4</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M43.16,-96.44C30.29,-91 15.54,-82.25 8,-69 3.41,-60.94 5.29,-51.56 9.3,-43.11\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"12.34,-44.84 14.27,-34.43 6.27,-41.37 12.34,-44.84\"/>\n",
       "<text text-anchor=\"middle\" x=\"26.5\" y=\"-57.8\" font-family=\"Times,serif\" font-size=\"14.00\">&#45;0.12</text>\n",
       "</g>\n",
       "<!-- x3 -->\n",
       "<g id=\"node4\" class=\"node\">\n",
       "<title>x3</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"102\" cy=\"-192\" rx=\"27\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"102\" y=\"-188.3\" font-family=\"Times,serif\" font-size=\"14.00\">x3</text>\n",
       "</g>\n",
       "<!-- x3&#45;&gt;x0 -->\n",
       "<g id=\"edge1\" class=\"edge\">\n",
       "<title>x3&#45;&gt;x0</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M112.53,-175.01C120.82,-162.43 132.52,-144.68 141.97,-130.33\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"144.94,-132.19 147.51,-121.92 139.09,-128.34 144.94,-132.19\"/>\n",
       "<text text-anchor=\"middle\" x=\"149\" y=\"-144.8\" font-family=\"Times,serif\" font-size=\"14.00\">0.90</text>\n",
       "</g>\n",
       "<!-- x3&#45;&gt;x2 -->\n",
       "<g id=\"edge4\" class=\"edge\">\n",
       "<title>x3&#45;&gt;x2</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M90.85,-175.54C86.99,-169.64 82.9,-162.72 80,-156 76.87,-148.75 74.36,-140.58 72.4,-132.98\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"75.78,-132.03 70.09,-123.1 68.96,-133.63 75.78,-132.03\"/>\n",
       "<text text-anchor=\"middle\" x=\"96\" y=\"-144.8\" font-family=\"Times,serif\" font-size=\"14.00\">0.90</text>\n",
       "</g>\n",
       "</g>\n",
       "</svg>\n"
      ],
      "text/plain": [
       "<graphviz.graphs.Digraph at 0x7fbfed27ccd0>"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "m = np.array([\n",
    "    [ 0.000,  0.000,  0.000,  0.895,  0.000,  0.000],\n",
    "    [ 0.565,  0.000,  0.377,  0.000,  0.000,  0.000],\n",
    "    [ 0.000,  0.000,  0.000,  0.895,  0.000,  0.000],\n",
    "    [ 0.000,  0.000,  0.000,  0.000,  0.000,  0.000],\n",
    "    [ 0.991,  0.000, -0.124,  0.000,  0.000,  0.000],\n",
    "    [ 0.895,  0.000,  0.000,  0.000,  0.000,  0.000]\n",
    "])\n",
    "\n",
    "make_dot(m)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "80b33192-08dc-49a4-b7e1-cc43fbde47fa",
   "metadata": {},
   "source": [
    "The variance of all variables in this data is approximately 1.0 and the mean is approximately 0.0, the variance of error terms excluding an exogenous variable is approximately 0.2."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "0023e24b-793a-400c-9440-5de52db5fe27",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1000, 6)"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sample_size = 1000\n",
    "\n",
    "error_vars = [0.2, 0.2, 0.2, 1.0, 0.2, 0.2]\n",
    "params = [0.5 * np.sqrt(12 * v) for v in error_vars]\n",
    "\n",
    "generate_error = lambda p: np.random.uniform(-p, p, size=sample_size)\n",
    "e = np.array([generate_error(p) for p in params])\n",
    "\n",
    "X = np.linalg.pinv(np.eye(len(m)) - m) @ e\n",
    "X = X.T\n",
    "\n",
    "X.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "c411ecd2-7004-4ca0-97ef-d0f690283660",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-0.057, -0.03 , -0.042, -0.056, -0.069, -0.061])"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X.mean(axis=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "4ae5dfbc-a141-41fe-b3d5-ea7a3c75a8f0",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.968, 0.989, 1.021, 0.982, 0.979, 0.973])"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X.var(axis=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "943fa72d-e1df-40a5-82d1-7d18410d253d",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.203, 0.215, 0.197, 0.982, 0.191, 0.202])"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "e.T.var(axis=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d5868a34-8a7a-413e-b02c-f788cb477314",
   "metadata": {},
   "source": [
    "Some data for x5 will be replaced by NaN using the MCAR method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "69539cca-8809-41bd-b025-26a64b7190c7",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "X_mcar = X.copy()\n",
    "\n",
    "prop_missing = [0, 0, 0, 0, 0, 0.1]\n",
    "\n",
    "for i, prop in enumerate(prop_missing):\n",
    "    mask = np.random.uniform(0, 1, size=len(X_mcar))\n",
    "    X_mcar[mask < prop, i] = np.nan"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "44c61343-e877-4fdd-b132-2b59b69d17f5",
   "metadata": {},
   "source": [
    "The proportaion of missing is the following:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "9469fdf1-ba8b-4047-8975-cbe161c4251a",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.   , 0.   , 0.   , 0.   , 0.   , 0.098])"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.isnan(X_mcar).sum(axis=0) / sample_size"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ac869c10-fa9a-48c7-8efe-fbbf8c2fe66f",
   "metadata": {},
   "source": [
    "## Causal discovery with missing data\n",
    "\n",
    "`bootstrap_with_imputation` discovers causality in given data with NaNs.\n",
    "\n",
    "In the following settings, `bootstrap_with_imputation` creates 30 bootstrap samples and repeat the imputation 5 times for each bootstrap sample and performs causal discovery assuming a common causal structure on the 5 imputed data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "bd34846f-fc46-4d19-8bd9-ca9f00a5c284",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "n_sampling = 30\n",
    "n_repeats = 5\n",
    "causal_orders, adj_matrices_list, resampled_indices, imputation_results = bootstrap_with_imputation(X_mcar, n_sampling, n_repeats)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "55fe01d0-0343-4360-bef4-378f4ecab0b9",
   "metadata": {},
   "source": [
    "`causal_orders` stores the causal order of the 30 bootstrap samples."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "58613e65-e7a7-4092-9874-012b918576e7",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(30, 6)"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "causal_orders.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "335d8ee1-f43b-4ad3-987f-c17f6167d5ae",
   "metadata": {},
   "source": [
    "`adj_matrices_list` is the list of adjacency matrices estimated from the results of 5 repeated imputations for each of the 30 bootstrap samples."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "75bdd668-62d4-45d9-a616-d50e36ee7438",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(30, 5, 6, 6)"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "adj_matrices_list.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "87bfce37-740b-4ac3-9cac-d87ece5932e4",
   "metadata": {},
   "source": [
    "`resampled_indices` contains lists of the original indices of the 30 bootstrap samples."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "053bd42f-7337-4538-b0da-44d4ea585ad3",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(30, 1000)"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "resampled_indices.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9948ecc0-28b3-4ef3-b80c-a03e19e53b36",
   "metadata": {},
   "source": [
    "`imputation_results` stores the result of imputations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "a1c10046-e376-4f77-842e-033dfe68a6ae",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(30, 5, 1000, 6)"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "imputation_results.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b0ad2c57-0278-4b5a-b1bc-29a3bfc47100",
   "metadata": {},
   "source": [
    "## Checking results\n",
    "\n",
    "### Comparing obtained results with results for data with no missing\n",
    "\n",
    "#### Histogram\n",
    "\n",
    "The results for data with no missing data are as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "c145158a-5b1d-4417-9052-15e3469deabf",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "model = lingam.DirectLiNGAM()\n",
    "bs_result = model.bootstrap(X, n_sampling)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "31375796-2194-4b88-b0ea-194432d055d0",
   "metadata": {},
   "source": [
    "The adjacency matirces is stored in `adjacency_matrices_`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "323762b5-bfd5-476e-bc8f-b8aa01e59964",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(30, 6, 6)"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "bs_result.adjacency_matrices_.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "864d08a1-a91e-450a-918e-c961ad01456b",
   "metadata": {},
   "source": [
    "The distributions of the elements of the estimated adjacency matrix of the two bootstrap results are compared as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "515b8965-9169-4050-ab52-c2511a46bfbc",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABAMAAAH/CAYAAAAi1hTrAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuNSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/xnp5ZAAAACXBIWXMAAAsTAAALEwEAmpwYAABCkklEQVR4nO3df9Asd10n+vcnHJMcFSUBQSSYk6Sg3KC7IlnW1dLKnrgHhBXQUGu8rCYKt1S2lto/rNpDsVZtcXfL4K29u1pYl6KICm4JKJYazbrekMNRr7UBAwQCupBzDvGaQxQ2CmoS4gLf+8f0ITNPnjnPPPOzn+nXq2rq6enp7vn0tz/9nX4+091TrbUAAAAAw3HBpgMAAAAA1ksxAAAAAAZGMQAAAAAGRjEAAAAABkYxAAAAAAZGMQAAAAAG5tA8M11wwQXt8OHDy45lq1z+uc8lSf704os3HEn/PPzww621tpZClFydjXydTr7SZ+P7rlxdD/3lcsjXfpHX060zVxP5Ogv5Ot1+83WuYsDhw4fz0EMPzTPrcFx77ejvyZObjKKXquqRdb2XXJ2RfJ1KvtJrY/uuXF0T/eVSyNeekddTrTNXE/k6E/k61X7z1WUCAAAAMDBznRnAbO4882BuOH7bxLhbjr1m4vl1R0+vMySYame+7szVRL7Cqh3Z4zMjsR+O29leSXLfzS/e93LuOHHV48ZpZ4bGcSsHiePW5XBmAAAAAAyMYgAAAAAMjGIAAAAADIxiAAAAAAyMYgAAAAAMjGIAAAAADIxiAAAAAAyMYgAAAAAMjGIAAAAADMyhTQcAANB3R47fliR5x5kHkyQ3HL8t99384k2GBAALcWYAAAAADIxiAAAAAAyMYgAAAAAMjGIAAAAADIxiAAAAAAyMXxMAALbGHSeumnh+3dHTG4oEAPrNmQEAAAAwMIoBAAAAMDCKAQAAADAwigEAAAAwMG4gCAAwBzcrBOAgc2YAAAAADIxiAAAAAAyMywQAgEE5cvy2iee3HNtQIACwQc4MAAAAgIFRDAAAAICBUQwAAACAgVEMAAAAgIFRDAAAAICBUQwAAACAgVEMAAAAgIE5tOkAAAB2uuPEVY8bd93R0xuIBAC2kzMDAAAAYGAUAwAAAGBgFAMAAABgYBQDAAAAYGAUAwAAAGBgFAMAAABgYBQDAAAAYGAUAwAAAGBgFAMAAABgYBQDAAAAYGAObToAAGB73XHiqseNu+7o6Q1EAgCMUwwAADbuyPHbJp7fcmxDgQDAQLhMAAAAAAZGMQAAAAAGRjEAAAAABsY9AwAA1sQNFQHoC8UAAIAV2dYbI+5cryS57+YXbyASAOblMgEAAAAYGMUAAAAAGBiXCQD0nNNxAQBYNmcGAAAAwMBUa23/M1V9Mckju7x0KMnnFw1qi2iPSefa43BrbS2FqPPk6ng8jGiPSX3L10X1ffuKb34HLVf73JaL2tZ1W8V69SVft3WbzUNbTFp735r4P2sftMekufJ1rmLA1IVV3dVau2ZpCzzgtMekvrVH3+LZNO0xadvao+/rI7759Tm23Ry0ePdjW9dtW9cr2e512y9tMalv7dG3eDZNe0yatz1cJgAAAAADoxgAAAAAA7PsYsCbl7y8g057TOpbe/Qtnk3THpO2rT36vj7im1+fY9vNQYt3P7Z13bZ1vZLtXrf90haT+tYefYtn07THpLnaY6n3DAAAAAD6z2UCAAAAMDALFQOq6tKqur2q7u3+XjJlui9U1d3d49ZF3rOPquqFVfWxqjpVVcd3ef2iqnpn9/p7q+rIBsJcmxna46aq+vRYTrxqTXENPl/l6qS+5up+7COvb+ymubeqbuzGPXFs3e6uqv9ZVf+5e20p675IfN34k902OhfHU7vxC+fqgm335VV1W1X9j6r6aFXdPDb9Qm23yH5aVa/txn+sql4w6zKXZYbYv76q3lNVH6yqD1fVi7rxR6rqkbE2e9OqYpzHDOt1eVXd0a3Tyaq6bOy1XXO7DxZcrwP3WbmPff7Ardt+LNLHbKMZ2sNx6wbJ18esJFdba3M/kvx0kuPd8PEkb5gy3d8u8j59fiR5QpLTSa5McmGSDyW5esc0r07ypm74hiTv3HTcG26Pm5K8cQOxDTpf5epc7bGRXN3neuyZ10kuTXKm+3tJN3zJLtO9P8l3LnPdF40vyckk1+wyz8K5ukhsSb48yT/pprkwyR8k+e5F226R/TTJ1d30FyW5olvOE2ZZ5pJycZbY35zkx8fiva8bPpLkI5vYh5a0Xr+a5MZu+GiSX9ortzf9WGS9uucH7rNyln3+oK7bkre7Y4HJaebu0xeMTb7K1/22xb5zddHLBF6a5K3d8FuTvGzB5R1Ez09yqrV2prX2d0nekVG7jBtvp3clua6qao0xrtMs7bEpQ89XuTqpz7m6H7Pk9QuS3N5a+8vW2l8luT3JC8cnqKpnJ3lqRv/U9i6+PZY7b67OHVtr7eHW2nuSpMufDyS5bJf592uR/fSlSd7RWnu0tfaJJKe65a0r12d5n5bkq7rhr07yyRXEsWyzrNfVSU50w+8Ze32e3F6XRdbroBr6cUDiWGCnPh8LyFf5Om4lubpoMeBprbUHuuE/T/K0KdNdXFV3VdWdVfWyBd+zb56R5M/Gnt/fjdt1mtba55N8NsmT1xLd+s3SHklyfXfa4buq6pnrCW3w+SpXJ/U5V/djlryeZV3PVdPH7yq7jHVfRny/0J3u9pNjH/DLyNWltF1VPSnJ9yS5Y2z0vG23yH46bd5Zc31Rs7zPv0vyL6rq/iT/Ncm/Gnvtiu7ygd+rqu9YQXzzmmW9PpTk+7rh703yxKo63zbpg0XWKzmYn5VDPw5IHAvs1OdjAfkqX8etJFcP7TVBVb07ydfu8tLrxp+01lpVTftpgstba2er6sokJ6rqntba6b3em631W0ne3lp7tKp+NKNq3tFlLFi+smQry9X9WFJe7+WGJD849nzmdV9xfK/o9scnJvm1Lsa3zTrzqtuuqg4leXuSn22tnelG9yJveuoHkvxia+0/VtU/TvJLVfWNSR5I8vWttQer6nlJfqOqntNa++uNRju7n0jyxqq6KcnvJzmb5AsbjWg5zrdevfysdBzACjhu5aDYd67uWQxorX3XtNeq6i+q6umttQeq6ulJPjVlGWe7v2eq6mSS52Z0zcM2OJtkvOpyWTdut2nu7w4cvzrJg+sJb+32bI/W2vi6vyWja6KWQr6el1ydtNFc3Y8l5PXZJNeOPb8so2vxzy3jHyQ51Fp7/9h7zrzuq4xvbH/8m6r65YxOk3tbZszVVbddRtfA39ta+89j77lI3iyyn55v3r2WuQyzxP7KdKfJt9b+e1VdnOQprbVPJXm0G//+qjqd5NlJ7lpBnPs1S1/xyXTfoFfVVya5vrX2maraK382ae716l7r5Wel44A9ORaY5Li13+TrY1aSq4teJnBrknN3xr0xyW/unKCqLqmqi7rhpyT59iR/vOD79skfJXlWVV1RVRdm9O3azjt5jrfTy5Oc2HEq7jbZsz26Du2clyT5kzXFNvR8lauT+pyr+7FnXif53STHuvy+JMmxbtw5P5DRt9tfssR1nzu+qjrU7Yepqi9L8s+SfGSX5c6bqwu1XVX9+4wOOv71+AwLtt0i++mtSW6o0Z2Vr0jyrCTvm3GZyzDL+/x/Sa5Lkqr6e0kuTvLpqvqaqnpCN/7KLvYz6YdZ+oqnVNW5Y6rXJvn5bnivfW+T5l6vA/xZOfTjgMSxwE59PhaQr/J13GpytS12V8MnZ3SN5L1J3p3k0m78NUne0g1/W5J7Mrru7J4kr1zkPfv4SPKiJB/PqAr3um7c65O8pBu+OKM78p7K6MDsyk3HvOH2+KkkH+1y4j1JvmFNcQ0+X+XqwcjVfa7DnnndPf+RbrueSvLDO5ZxZue6LWvdF4kvyVdk9AsHH+5i+ZkkT1hWri4Y22UZ3QzvT5Lc3T1etYy2W2Q/zejU0tNJPpbu1w2mLXND+9TVSf6wa5u7kxzrxl/ftdndGd2M8Xs2vW/tc71e3uXRxzP6NuaiWfa9TT/mXa8c0M/KWfb5g7puS97ujgV6cCwgX2fePoPJ11XkanUzAgAAAAOx6GUCAAAAwAGjGAAAAAADoxgAAAAAA6MYAAAAAAOjGAAAAAADoxgAAAAAA6MYAAAAAAOjGAAAAAADoxgAAAAAA6MYAAAAAAOjGAAAAAADoxgAAAAAA3NonpkuuOCCdvjw4WXHsnUu/9znkiR/evHFG46kXx5++OHWWltLIUquzkauTidfN09+zkau0nfj+7J87Rf97HTrzNVk93y1fZjVfvN1rmLA4cOH89BDD80z67Bce+3o78mTm4yid6rqkXW9l1ydkVydSr72gPyciVyl98b2ZfnaM/rZqdaZq8mUfLV9mNF+83WuYsC2uePEVRPPrzt6+nHTHDl+28Tz+25+8UpjgnW788yDuWEsz2859prHTbPbvgGLmKdv3TnPrPPttLPvT+Q4zGq3/XDn54b96WDZeRyQ2KZ9YvuwCu4ZAAAAAAOjGAAAAAADoxgAAAAAA6MYAAAAAAOjGAAAAAADoxgAAAAAA6MYAAAAAANzaNMBAMA5d5y46nHjrsszNxAJAMB2c2YAAAAADIxiAAAAAAyMYgAAAAAMjHsGrNnO62GvO3p6Q5EAAAAwVM4MAAAAgIFRDAAAAICBUQwAAACAgVEMAAAAgIEZ3A0Ejxy/7XHjbjm2gUAAAABgQ5wZAAAAAAOjGAAAAAADoxgAAAAAA6MYAAAAAAOjGAAAAAADM7hfEwBg+9xx4qqJ59cdPb0V7wUAsCpbVQzY5MHgqt8PAAAAlsVlAgAAADAwigEAAAAwMFt1mQAALNOR47c9btwtxzYQCADAkjkzAAAAAAZGMQAAAAAGRjEAAAAABsY9AwBYiXX+3CsAAPvjzAAAAAAYGGcGrNCdZx5MktwwdjfqWe5CvfPu1ffd/OKlxgUAAMCwKQYAMEFBEgBg+7lMAAAAAAZGMQAAAAAGRjEAAAAABubA3jNg5zWtyWw35wMAAIChO7DFAIAhu+PEVRPPrzt6ekORMIud2yuxzQCAzXKZAAAAAAyMMwMAGISdl5e5tAwAGDLFgAPA6aUAAAAsk8sEAAAAYGAUAwAAAGBgFAMAAABgYBQDAAAAYGDcQBBgQHbekHSWm5G6iSkAwPZZeTFgngNPlkPbw3bY+ZN4yWw/izfvfPNY53sdBPP+jOHO+e67+cXLCgkAYEIvzgxw8HOw7HbQb5sBAAAcHO4ZAAAAAAOjGAAAAAADoxgAAAAAA9OLewYAANvBfYAAYHlWeVP4aq3tf6aqLyZ5ZGlRTDqU5PMrWvYy9Dm+PseWPBbf4dbaWs5K2SNX+95e66Y9JvUtXzdtaPlxENd3m3O1r9tDXPvXt761z221btpi0tpzNTlvvvZx+/Qtpr7Fk6w/pn3l61zFgFWqqrtaa9dsOo5p+hxfn2NL+hdf3+LZNO0xSXtMGlp7DG19+66v20Nc+9e32PoWzyZpi0l9a4++xZP0L6a+xZP0M6Zx7hkAAAAAA6MYAAAAAAPTx2LAmzcdwB76HF+fY0v6F1/f4tk07TFJe0waWnsMbX37rq/bQ1z717fY+hbPJmmLSX1rj77Fk/Qvpr7Fk/Qzpi/p3T0DAAAAgNXq45kBAAAAwAopBgAAAMDArLUYUFUvrKqPVdWpqjq+y+s3VdWnq+ru7vGqsddurKp7u8eNPYvtC2Pjb112bLPE103zz6vqj6vqo1X1y2PjN9p2e8S28rYbe69Lq+r2rh1ur6pLpky3tpjWbYY8v6iq3tm9/t6qOrKBMNdmkf1+m8zQDj9WVfd0bfD/VtXV3fgjVfXIWPu8af3R788s/VU33fVV1arqmrFxr+3m+1hVvWA9EQ/HPvro/1ZVn6mq394x/her6hNj+fjNPYnriq4/PdX1rxeuOa5djwGq6mSXy+fa66kLxjP358u69i3HASOOBSZt6lhgFfvMrJ9xy4ynqv5pVb2/RscJ76+qo2PzLNTPLBDT1OOTqnpeF+upqvrZqqo1xPOKsVjurqovVvcZtey+eN9aa2t5JHlCktNJrkxyYZIPJbl6xzQ3JXnjLvNemuRM9/eSbviSPsTWvfa3PWi7ZyX54Ll2SfLUHrXdrrGto+12xPHTSY53w8eTvGET23NTjxm31auTvKkbviHJOzcd94bbY+p+vy2PGdvhq8aGX5Lkv3XDR5J8ZNPrsMx17aZ7YpLfT3Jnkmu6cVd301+U5IpuOU/Y9Dpt02MfffR1Sb4nyW/vGP+LSV7ew7h+JckN3fCbkvz4uuI63zFAkpPn8nsJscz9+bLOfWsf23IrjwMW3Vbb+JixPW7Kko8FVrHPzLLMFcXz3CRf1w1/Y5KzY/PM3c8sGNORTDk+SfK+JN+apJL8TpLvXnU8O6b5piSnl9FGy3is88yA5yc51Vo701r7uyTvSPLSGed9QZLbW2t/2Vr7qyS3J3lhT2Jbh1ni+9+T/FzXPmmtfaob34e2mxbbur00yVu74bcmedmG4tiUWbbVeBu9K8l1+6mYHjB93+/XZc92aK399djTr0hyUO88O+s2/z+SvCHJ58bGvTTJO1prj7bWPpHkVLc8lmemPrq1dkeSv1lTTMkCcXX959GM+tPzzr+iuFZ9DHDOIp8v69y3hn4ckDgW2GlTxwKr2GcWWZe542mtfbC19slu/EeTHK6qi2Z835XENG2BVfX0jL7guLON/hN/W2bvB5YVzw908/bCOosBz0jyZ2PP7+/G7XR9VX24qt5VVc/c57ybiC1JLq6qu6rqzqp62RLj2k98z07y7Kr6wy6OF+5j3k3Flqy+7cY9rbX2QDf850meNmW6dca0TrNsqy9N01r7fJLPJnnyWqJbv0X3+20xUztU1b+sqtMZfbP2mrGXrqiqD1bV71XVd6w21IXtua5V9S1Jntlau22/87KwWfvo8/kP3f76n5Z0MLpoXE9O8pmuP02WmzezxLVX3v5Cd1rqTy74z94iny/r3LeGfhyQOBbYaVPHAqvYZxbZl5aVF9cn+UBr7dGxcfP2M4vGtNvxyTO65ZxvmauK55zvT/L2HeOW1Rfv26F1vtkMfivJ21trj1bVj2ZUWTm6xzzrcr7YLm+tna2qK5OcqKp7Wmun1xzfoYxOx782yWVJfr+qvmnNMUyza2yttc9kyW1XVe9O8rW7vPS68SettVZV077d7MP2pB/63CetVWvt55L8XFX9b0n+bZIbkzyQ5Otbaw9W1fOS/EZVPWfHmQQHRlVdkOT/yuiUUFZgSX30NK/N6B+8CzP6Xed/k+T1PYhrbiuO6xXdZ90Tk/xakh/M6FuyA81xACvgWGAGVfWcjM6qOzY2elP9zK7HJ2t43/Oqqn+U5OHW2kfGRm+0L15nMeBskvFK2mXduC9prT049vQtGX0DdW7ea3fMe7InsaW1drb7e6aqTmZ07cwyPzT2jC+j6tR7W2v/K8knqurjGf0DvvG2O09sf7Tstmutfde016rqL6rq6a21B7rThHa9XGEN23NTZtlW56a5v6oOJfnqJA9mOy2032+RWfJi3DuS/N9J0lX+H+2G39+dOfDsJHetJtSF7bWuT8zoeseTXWH+a5PcWlUvmWFeZrCMPvo8yz73je+jVfULSX6iB3E9mORJVXWo+5ZoX3mzhLimHgOMfdb9TY1u7Pv8zH8Ausjny1L3LccBe3IsMGlTxwKr2mfm3ZcWyouquizJryf5ofHC2YL9zNwxdZcA7HZ8crZbzvmWufR4xl6/ITvOClhyX7xv67xM4I+SPKtGd9W9MKPGmLhDa9cxn/OSJH/SDf9ukmNVdUmN7vx6rBu38di6mC7qhp+S5NuT/PESY5spviS/ke4Dv4vj2RndKGjjbTcttjW13bhbM/o2M93f39w5wQZiWqdZttV4G708yYmuQ91Gi/RJ22SWdnjW2NMXJ7m3G/81VfWEbvjKjIp8Z9YS9XzOu66ttc+21p7SWjvSWjuS0Q0EX9Jau6ub7oYa3Sn4iozW9X3rX4WttmcffT7n9tcaVXJeluQj551hDXF1/ed7MupP9z3/EuLa9Rigqg51n3Gpqi9L8s+yWHst8vmyzn1r6McBiWOBnTZ1LLCKfWaWZS49nqp6UpLbMro55x+em3gJ/cwiMe16fNIVjf+6qr61+6z4oczeJy+079To7MN/nrH7BaygL96/tsa7FSZ5UZKPZ1RdfV037vUZHWwlyU9ldOOJD2X04fkNY/P+SEY3yDiV5If7EluSb0tyTzf+niSv3FDbVUant/5xF8cNPWq7XWNbV9uNxfnkJHdk9I/Mu5Nc2o2/JslbNhHTuh8zbKuLk/xqlyvvS3LlpmPecHtM7ZO26TFDO/xM1w53d+3wnG789WPjP5Dkeza9Louu645pT2bsDr8ZnWZ8OsnHMuPdhz32tW327KO753+Q5NNJHsnozLMXdONPdP32R5L8lyRf2ZO4ruz601Nd/3rRmuN63DFARjcCfX+SD3f78M9kwTv4z9CPTP18Wde+NUubZcuPAxbdVtv4mKE9VnIssIp9ZrdlrjqejC4dfCijY4Fzj6cuo59ZIKapxyfd/v6RbplvTFJr2mbXJrlzx/KW3hfv91FdIAAAAMBArPMyAQAAAKAHFAMAAABgYBQDAAAAYGAUAwAAAGBgFAMAAABgYBQDAAAAYGAUAwAAAGBgFAMAAABgYBQDAAAAYGAUAwAAAGBgFAMAAABgYA7NM9MFF1zQDh8+vOxYtsrln/tckuRPL754w5H0z8MPP9xaa2spRMnV2cjX6eRrv8jV6bY1V23z7bSt+XpQ2c+mW2euJvJ1FvJ1uv3m61zFgMOHD+ehhx6aZ9bhuPba0d+TJzcZRS9V1SPrei+5OiP5OpV87Rm5OtXW5qptvpW2Nl8PKvvZVOvM1US+zkS+TrXffHWZAAAAAAzMXGcGMJs7zzyYG47fNjHulmOvmXh+3dHT6wwJptqZrztzNZGv9IO+lWW748RVjxu3rBw6siNX77v5xUtZLiybvpWDxHHrcjgzAAAAAAZGMQAAAAAGRjEAAAAABkYxAAAAAAZGMQAAAAAGRjEAAAAABkYxAAAAAAZGMQAAAAAGRjEAAAAABkYxAAAAAAZGMQAAAAAGRjEAAAAABkYxAAAAAAZGMQAAAAAGRjEAAAAABkYxAAAAAAZGMQAAAAAG5tCmAwAA6Lsjx2973Lj7bn7xBiIBgOVwZgAAAAAMjDMDAADmcMeJqyaeX3f09IYiAYD9c2YAAAAADIxiAAAAAAyMYgAAAAAMjHsGAAAHkjv8A8D8nBkAAAAAA6MYAAAAAAPjMgEAgDXZ+XOEiZ8kBGAznBkAAAAAA6MYAAAAAAPjMgEAYFB2/grBLcc2FMgUfiUBgHVwZgAAAAAMjGIAAAAADIxiAAAAAAyMewYAAFtj50/3Dfln+/yMIQDn48wAAAAAGBjFAAAAABgYxQAAAAAYGMUAAAAAGBjFAAAAABgYvyYAAKyMO9ovh19JAGDZnBkAAAAAA6MYAAAAAAOjGAAAAAADoxgAAAAAA6MYAAAAAAOjGAAAAAAD46cFAYDeOfdTet/ymU8mST5w4io/p7cm8/4cpJ8/BDhYFAMAALbAkeO3TTy/5dj+55l1PgAOPpcJAAAAwMAoBgAAAMDAKAYAAADAwCgGAAAAwMAoBgAAAMDAVGtt/zNVfTHJI7u8dCjJ5xcNaotoj0nn2uNwa20thajz5Op4PIxoj0l9ylfbZpL2mNSnXN20PudGn2NL1h9fX/K179tlnbTFpLX3rYljgX3QHpPmyte5igFTF1Z1V2vtmqUt8IDTHpP61h59i2fTtMekPrVHn2LpA+0xSXs8ps9t0efYkv7HtypDXe/daItJfWuPvsWzadpj0rzt4TIBAAAAGBjFAAAAABiYZRcD3rzk5R102mNS39qjb/FsmvaY1Kf26FMsfaA9JmmPx/S5LfocW9L/+FZlqOu9G20xqW/t0bd4Nk17TJqrPZZ6zwAAAACg/1wmAAAAAAOzUDGgqi6tqtur6t7u7yVTpvtCVd3dPW5d5D37qKpeWFUfq6pTVXV8l9cvqqp3dq+/t6qObCDMtZmhPW6qqk+P5cSr1hTX4PNVrk6Sq/0mXyf1NV/XYYZ1v7yq7qiqD1fVyaq6bOy1G7t96d6qurFP8VXVN1fVf6+qj3avfX9fYht7/auq6v6qeuOyY1s3feuIvnVSX/tW+ToiXx+zklxtrc39SPLTSY53w8eTvGHKdH+7yPv0+ZHkCUlOJ7kyyYVJPpTk6h3TvDrJm7rhG5K8c9Nxb7g9bkryxg3ENuh8latztYdc7ff2ka+T02wkX3uy7r+a5MZu+GiSX+qGL01ypvt7STd8SY/ie3aSZ3XDX5fkgSRP6kNsY6//TJJf3obc0rfqW+dsD8cC/d4+g8jXVeXqopcJvDTJW7vhtyZ52YLLO4ien+RUa+1Ma+3vkrwjo3YZN95O70pyXVXVGmNcp1naY1OGnq9ydZJc7Tf5OqnP+bpqs6z71UlOdMPvGXv9BUlub639ZWvtr5LcnuSFfYmvtfbx1tq93fAnk3wqydf0IbYkqarnJXlakv9niTFtkr5V37pTn/tW+Spfx60kVxctBjyttfZAN/znGX1g7Obiqrqrqu6sqpct+J5984wkfzb2/P5u3K7TtNY+n+SzSZ68lujWb5b2SJLru1MS31VVz1xPaIPPV7k6Sa72m3yd1Od8XbVZ1v1DSb6vG/7eJE+sqifPOO8m4/uSqnp+Rt/2nO5DbFV1QZL/mOQnlhjPpulb9a079blvla/yddxKcvXQXhNU1buTfO0uL71u/ElrrVXVtJ8muLy1draqrkxyoqruaa0t88OOg+W3kry9tfZoVf1oRtW8o8tYsHxlyeQqB8nK8vUA+Ikkb6yqm5L8fpKzSb6w0YgmnTe+qnp6kl/K6HT9L/Yktlcn+a+ttfsP0pds+lZWwLEAB8W+c3XPYkBr7bumvVZVf1FVT2+tPdB9kH1qyjLOdn/PVNXJJM/Ncivfm3Q2yXjV5bJu3G7T3F9Vh5J8dZIH1xPe2u3ZHq218XV/S0bXRC2FfD0vuTpJrvabfJ200XzdsFnW/ZPpvt2uqq9Mcn1r7TNVdTbJtTvmPdmX+LrnX5XktiSva63d2ZfYquofJ/mOqnp1kq9McmFV/W1r7XE3reoTfeue9K2THAv0m3x9zEpyddHLBG5Ncu7OvDcm+c2dE1TVJVV1UTf8lCTfnuSPF3zfPvmjJM+qqiuq6sKMblyx806e4+308iQnWneXhy20Z3t0Hdo5L0nyJ2uKbej5KlcnydV+k6+T+pyvqzbLuj+lO609SV6b5Oe74d9NcqzbXy5Jcqwb14v4uul/PcnbWmvvWnJcC8XWWntFa+3rW2tHMjp74G19LwTMQN+qb92pz32rfJWv41aTq9PuLDjLI6PrMe5Icm+Sdye5tBt/TZK3dMPfluSejK5JuyfJKxd5zz4+krwoycczqsK9rhv3+iQv6YYvzuhuvaeSvC/JlZuOecPt8VNJPtrlxHuSfMOa4hp8vspVuXqQHvJ13+2xkXztybq/vNtfPp7RtyEXjc37I12OnEryw32KL8m/SPK/ktw99vjmPsS2Yxk3ZTt+TUDfOltO6Ft70LfK15m3z2DydRW5Wt2MAAAAwEAsepkAAAAAcMAoBgAAAMDAKAYAAADAwCgGAAAAwMAoBgAAAMDAKAYAAADAwCgGAAAAwMAoBgAAAMDAKAYAAADAwCgGAAAAwMAoBgAAAMDAKAYAAADAwByaZ6YLLrigHT58eNmxbJXLP/e5JMmfXnzxhiPpn4cffri11tZSiJKrs5Gv08nXfpGr08nV/pGv08nXfpGr060zVxP5Ogv5Ot1+83WuYsDhw4fz0EMPzTPrcFx77ejvyZObjKKXquqRdb2XXJ2RfJ1KvvaMXJ1KrvaQfJ1KvvaMXJ1qnbmayNeZyNep9puvLhMAAACAgZnrzABmc+eZB3PD8dsmxt1y7DUTz687enqdIcFUO/N1Z64m8pV+0LdykOhbOSj0rRwk+tblcGYAAAAADIxiAAAAAAyMYgAAAAAMjGIAAAAADIxiAAAAAAyMYgAAAAAMjGIAAAAADIxiAAAAAAyMYgAAAAAMjGIAAAAADIxiAAAAAAyMYgAAAAAMjGIAAAAADIxiAAAAAAyMYgAAAAAMjGIAAAAADIxiAAAAAAyMYgAAAAAMjGIAAAAADIxiAAAAAAyMYgAAAAAMjGIAAAAADIxiAAAAAAyMYgAAAAAMjGIAAAAADIxiAAAAAAyMYgAAAAAMjGIAAAAADIxiAAAAAAyMYgAAAAAMjGIAAAAADIxiAAAAAAyMYgAAAAAMjGIAAAAADIxiAAAAAAyMYgAAAAAMjGIAAAAADIxiAAAAAAyMYgAAAAAMjGIAAAAADIxiAAAAAAyMYgAAAAAMjGIAAAAADIxiAAAAAAyMYgAAAAAMjGIAAAAADIxiAAAAAAyMYgAAAAAMjGIAAAAADIxiAAAAAAyMYgAAAAAMjGIAAAAADIxiAAAAAAxMtdb2P1PVF5M8sstLh5J8ftGgtoj2mHSuPQ631tZSiDpPro7Hw4j2mNSnfLVtJmmPSX3K1fF4GNEek+Rrf2mLSWvP1cSxwD5oj0lz5etcxYCpC6u6q7V2zdIWeMBpj0l9a4++xbNp2mNSn9qjT7H0gfaY1Lf26Fs8m6Y9JvWtPfoWzyZpi0l9a4++xbNp2mPSvO3hMgEAAAAYGMUAAAAAGJhlFwPevOTlHXTaY1Lf2qNv8Wya9pjUp/boUyx9oD0m9a09+hbPpmmPSX1rj77Fs0naYlLf2qNv8Wya9pg0V3ss9Z4BAAAAQP+5TAAAAAAGZqFiQFVdWlW3V9W93d9Lpkz3haq6u3vcush79lFVvbCqPlZVp6rq+C6vX1RV7+xef29VHdlAmGszQ3vcVFWfHsuJV60prsHnq1ydJFf7Tb5Okq/9Jl8fI1f7Ta5Okq/9Jl8fs5Jcba3N/Ujy00mOd8PHk7xhynR/u8j79PmR5AlJTie5MsmFST6U5Ood07w6yZu64RuSvHPTcW+4PW5K8sYNxDbofJWrc7WHXO339pGvk9PI135vn0Hkq1zt90OuztUe8rXf22cQ+bqqXF30MoGXJnlrN/zWJC9bcHkH0fOTnGqtnWmt/V2Sd2TULuPG2+ldSa6rqlpjjOs0S3tsytDzVa5Okqv9Jl8nydd+k6+Pkav9Jlcnydd+k6+PWUmuLloMeFpr7YFu+M+TPG3KdBdX1V1VdWdVvWzB9+ybZyT5s7Hn93fjdp2mtfb5JJ9N8uS1RLd+s7RHklxfVR+uqndV1TPXE9rg81WuTpKr/SZfJ8nXfpOvj5Gr/SZXJ8nXfpOvj1lJrh7aa4KqeneSr93lpdeNP2mttaqa9tMEl7fWzlbVlUlOVNU9rbXTe703W+u3kry9tfZoVf1oRtW8o8tYsHxlyeQqB4l85aCQqxwk8pWDYt+5umcxoLX2XdNeq6q/qKqnt9YeqKqnJ/nUlGWc7f6eqaqTSZ6b0TUP2+BskvGqy2XduN2mub+qDiX56iQPrie8tduzPVpr4+v+loyuiVoK+XpecnWSXO03+TpJvvabfH2MXO03uTpJvvabfH3MSnJ10csEbk1yYzd8Y5Lf3DlBVV1SVRd1w09J8u1J/njB9+2TP0ryrKq6oqouzOjGFTvv5DneTi9PcqJ1d3nYQnu2R9ehnfOSJH+yptiGnq9ydZJc7Tf5Okm+9pt8fYxc7Te5Okm+9pt8fcxqcnU/dxvc+cjoeow7ktyb5N1JLu3GX5PkLd3wtyW5J6M7Ht6T5JWLvGcfH0lelOTjGVXhXteNe32Sl3TDFyf51SSnkrwvyZWbjnnD7fFTST7a5cR7knzDmuIafL7KVbl6kB7ydd/tIV/7vX0Gk69ytd8Pubrv9pCv/d4+g8nXVeRqdTMCAAAAA7HoZQIAAADAAaMYAAAAAAOjGAAAAAADoxgAAAAAA6MYAAAAAAOjGAAAAAADoxgAAAAAA6MYAAAAAAOjGAAAAAADoxgAAAAAA6MYAAAAAANzaJ6ZLrjggnb48OFlx7J1Lv/c55Ikf3rxxRuOpF8efvjh1lpbSyFKrjKP8X1Xvm6Xbe6X+5Cr29y+LFcf8pXH2HenW2euJrvnq+3DrPabr3MVAw4fPpyHHnponlmH5dprR39PntxkFL1TVY+s673kKnMZ23fl65bZ4n65F7m6xe3LcvUiX3mMfXeqdeZqMiVfbR9mtN98dZkAAAAADMxcZwYwvztOXDXx/Lqjp5ey3CPHb3vcuPtufvFSlg0H2c5945Zjr3ncNMvaD5k0S3+3c/uss9/aGV8iF5bhzjMP5oax7Wqfg4Nh576bPH7/te9uju3DKjgzAAAAAAZGMQAAAAAGxmUCAEzY5Kn7AACshzMDAAAAYGAUAwAAAGBgFAMAAABgYBQDAAAAYGAUAwAAAGBg/JoAAAfeHSeumnh+3dHTG4oEAOBgcGYAAAAADIxiAAAAAAyMYgAAAAAMjGIAAAAADIxiAAAAAAyMYgAAAAAMjJ8W7KEjx2+beH7fzS/eUCQAAABsI2cGAAAAwMA4MwCA3rjjxFWPG3fd0dMbiAQAYLs5MwAAAAAGRjEAAAAABkYxAAAAAAZGMQAAAAAGRjEAAAAABsavCQCwsCPHb3vcuFuObSAQAABm4swAAAAAGBjFAAAAABgYxQAAAAAYGMUAAAAAGBjFAAAAABgYxQAAAAAYGMUAAAAAGJhDmw6Avd1x4qrHjbvu6Ol9zzfLPAB9d+T4bY8bd8ux9b3XfTe/eDVvBgCwRs4MAAAAgIFRDAAAAICBUQwAAACAgXHPgDnNex0/wDLMe08Q9xJZ3Drb0GcNALAqigErdOeZB5MkN4zdgGpVN7kCAACAWblMAAAAAAbGmQEDt/Nns/xkFvTPvD+lt86f4AMA4GBxZgAAAAAMjDMDWAo3JQMAADg4nBkAAAAAAzO4MwN2u4bWdfL74zpkAACAg82ZAQAAADAwgzszAADWbed9VRL3VgEANsuZAQAAADAwzgwAYBCOHL8t7zjzYJLkhuO39fJeJzvvydLHGAGA7bDUYoCb8wEAAED/HdgzA9ZdePBtDQAAANuil8UAN1oCAACA1XEDQQAAABiYXp4ZAJvgjBQAmM1ul2vecuw1E899hgL0m2IAAAAAbNi674tXrbX9z1T1xSSPnGeSQ0k+P29QKyKm2awjpsOttbVcojJDri6qj9vwnD7Hlhyc+LYpX/uk79t/mda1rn3J1SFt28T6zku+9o+2mLT244DkvPna9+3T5/j6HFuy3Pj2la9zFQP2XGjVXa21a5a+4AWIaTZ9jKnP+txefY4tEd/QDal9h7SuifXddtu2vtu2PovQFpP61h59i2enPsfX59iSzcbnBoIAAAAwMIoBAAAAMDCrKga8eUXLXYSYZtPHmPqsz+3V59gS8Q3dkNp3SOuaWN9tt23ru23rswhtMalv7dG3eHbqc3x9ji3ZYHwruWcAAAAA0F8uEwAAAICB2VcxoKpeWFUfq6pTVXV8l9cvr6o7qurDVXWyqi4be+3Gqrq3e9y4jOCXENMXquru7nHrEmP6+ar6VFV9ZMrrVVU/28X84ar6lrHXVtVOi8S0knbqs6q6tKpu77bD7VV1yZTpHre9quqJY+11d1X9z6r6z91rN1XVp8dee9U6Y+vGn+z2mXMxPLUbf1FVvbPLgfdW1ZH9xrZofFX15VV1W1X9j6r6aFXdPDb93G03Qz8xdd2r6rXd+I9V1QtmXSb7yoX/VlWfqarf3jH+F6vqE2Pb/JvXEvgclrCuV3S5d6rLxQvXE/n+rGJf6rN517eqjlTVI2O5+6a1Bz+HGdb3O6vqA1X1+ap6+Y7XVnL8skz72E+3+rhnkf14G83QHgsfu835vhs7Nlmg7/unVfX+qrqn+3t0bJ5djz/XHN/UvrmqntfFfapG/xPVmmN7RU3+//DF6o57ltl2j9Nam+mR5AlJTie5MsmFST6U5Ood0/xqkhu74aNJfqkbvjTJme7vJd3wJbO+9ypi6p7/7aIxTInrO5N8S5KPTHn9RUl+J0kl+dYk711lOy0S0yrbqc+PJD+d5Hg3fDzJG3aZZqbtleT9Sb6zG74pyRs3GVuSk0mu2WWeVyd5Uzd8Q5J3rju+JF+e5J9001yY5A+SfPcibTdjP7Hruie5upv+oiRXdMt5wizL9JgtF7rXrkvyPUl+e8f4X0zy8k2vx5rW9VeS3NANvynJj296nXaJfen70qbXaYXreyRTPm/7+phxfY8k+ftJ3ja+b07r0ze9Trus46z76dYe9yyS19v4mLE9bsqCx27L3A7T+tNZlrmG2J6b5Ou64W9McnZsnpPZ5fhzzfFN7ZuTvC+j/4Eqo/+Jvnudse2Y5puSnF522+322M+ZAc9Pcqq1dqa19ndJ3pHkpTumuTrJiW74PWOvvyDJ7a21v2yt/VWS25O8cB/vvYqYVqa19vtJ/vI8k7w0ydvayJ1JnlRVT8/q2mmRmIbqpUne2g2/NcnLdplmz+1VVc9O8tSM/qntVWx7LPddSa6bsyo6d3yttYdba+9Jkm6f/kCSy3aZfz9m6SemrftLk7yjtfZoa+0TSU51y5tlmcyWC2mt3ZHkb9YU06rMva5drh3NKPfOO/+GrWJf6rNF1vcg2nN9W2v3tdY+nOSLO+Zd2fHLks20n265oeX1Xjb1ed7nY5O5Y2utfbC19slu/EeTHK6qi+aIYSXxTVtg9z/PV7XW7myj/77flvn6h2XF9gPdvCu3n2LAM5L82djz+7tx4z6U5Pu64e9N8sSqevKM885jkZiS5OKququq7qyqly0hnllNi3tV7bRITMnm2mmTntZae6Ab/vMkT9tlmlm217mKXxsbd32NLsV4V1U9c0Ox/UJ3mtFPjnVAX5qntfb5JJ9N8uTs31LarqqelNE3qHeMjZ6n7WbZTtPWvY/76kEySy7s5T902/w/reCAYpkWWdcnJ/lMl3tJf/NpFftSny2yvklyRVV9sKp+r6q+Y9XBLsEi2+igbN9Z99NtPu5ZNK+3zay5u+ix2zzvu6ljk2XlyPVJPtBae3Rs3G7Hn+uOb7e++Rndcs63zHXEds73J3n7jnHLaLvHObSsBXV+Iskbq+qmJL+f5GySLyz5PfbrfDFd3lo7W1VXJjlRVfe01k5vKM4+28p2qqp3J/naXV563fiT1lqrqrbLdLO4IckPjj3/rSRvb609WlU/mlFl8OjOmVYc2yu67fnEJL/Wxfe2/Sxg1W1XVYcy6gR/trV2phs9U9uxXivOhddmdMB+YUY/u/Nvkrx+njiXYU19BtvhgSRf31p7sKqel+Q3quo5rbW/3nRg225J++lWHvcwN8cf+1RVz0nyhiTHxkYvfPy5BLv2zWuO4byq6h8lebi1Nn6ft5W13X6KAWeTjFfCLuvGfUl3Wsj3JUlVfWWS61trn6mqs0mu3THvyTniXVpM3Wtnu79nqupkRte5rKOznxb3qtppkZg22U4r1Vr7rmmvVdVfVNXTW2sPdKcOfWqXyc67varqHyQ51Fp7/9h7Pjg2/Vsyun5xrbGNbc+/qapfzuiUprflsRy4v/tn/KuTjMe7lvg6b05yb2vtP4+950xtN+W9zttPZPq6n2/evZY5CEvIhfMt+9w3eI9W1S9kVNzdmBWu64MZXZp1qPuWoK/5tKp9qa/mXt/ubLBHk6S19v6qOp3k2UnuWnnU81tkG23y+GXCMvbTbT3u6SyyH2+jWf6XmPf4Y6H3zeaOTRbKkRrdrP3Xk/zQeBHtPMefa4vvPH3z2UxelrqRtuvckB1nBSyx7R5nP5cJ/FGSZ9XojscXdoFO3GG1qp5SVeeW+dokP98N/26SY1V1SY3u3HqsG7eouWPqYrno3DRJvj3JHy8hplncmuSHauRbk3y2O+hdVTvNHdOG22mTbk1yYzd8Y5Lf3GWavbbXD2THzlyT92F4SZI/WWdsVXWo246pqi9L8s+SnKs8ji/35UlO7Li8YeXxdXH9+4w6xn89PsMCbbdnP5Hp635rkhtqdOfXK5I8K6MbzMyyTGbLhanObfOqqoyu3dv111B6Yu517XLtPRnl3r7nX6NV7Et9Nvf6VtXXVNUTkqT7dvlZGd1Ur88W6dc2efyyH3vupwM47llkP95Gs/wvsYxjt32/bzZ3bLJI3/ekJLdldKPOPzw38R7Hn+uMb9e+ufs/7K+r6lu7Y44fynyfwwvtX93/rP88Y/cLWHLbPV7b3x0SX5Tk4xlVR1/XjXt9kpd0wy9Pcm83zVuSXDQ2749kdIOLU0l+eD/vu4qYknxbknsyuqfAPUleucSY3p7RaSj/K6NrRV6Z5MeS/Ni5bZ3k57qY78nY3SFX2E5zxbTKdurzI6Nrd+7ocufdSS7txl+T5C2zbK+MDvy+Yce4n8rohiofyujg/xvWGVuSr8jo1w0+3MXxM+nu6J3k4ox+feNURh8qV6677TKqoLaMPmjv7h6vWrTtsnc/MXXdMzq19HSSj2XszrK7LdNj7lz4gySfTvJIRv3TC7rxJ7p+5yNJ/kuSr9z0Oq1wXa/scu9Ul4sXrXsdZlzPpe9LfX7Mu74ZXSv70a4P+0CS79n0uixpff9hl7cPZfRN1kfH5l3J8cuS12/P/TQDOO5ZZD/exscM7bHwsduyt8O0/nS3Za4ztiT/tusf7h57PDXnOf5cc3xT++auH/hIt8w3JqkNbNdrk9y5Y3lLbbudj+reBAAAABiI/VwmAAAAAGwBxQAAAAAYGMUAAAAAGBjFAAAAABgYxQAAAAAYGMUAAAAAGBjFAAAAABgYxQAAAAAYGMUAAAAAGBjFAAAAABgYxQAAAAAYGMUAAAAAGJhD88x0wQUXtMOHDy87FpJc/rnPJUn+9OKLNxzJ6jz88MOttbaWQtS0XB1CO7McfchX9maflqt9JC+n60O+2j7MYp25muhfZ2HfnW6/+TpXMeDw4cN56KGH5pmVvVx77ejvyZObjGKlquqRdb3X1FwdQDuzHL3IV/Zmn5arfSQvp+pFvto+zGCduZroX2di351qv/nqMgEAAAAYmLnODGA2R47f9rhxtxx7zcTz646eXlc47HDnmQdzw9g22rltEtsH+mi3vvWDn3lvkuQDJ65KYt+lP3zW9JvtAweTfXc5nBkAAAAAA6MYAAAAAAOjGAAAAAADoxgAAAAAA6MYAAAAAAOjGAAAAAAD46cFe2jnT2Xcd/OLNxgNAAAA28aZAQAAADAwigEAAAAwMIoBAAAAMDCKAQAAADAwigEAAAAwMIoBAAAAMDCKAQAAADAwigEAAAAwMIoBAAAAMDCKAQAAADAwigEAAAAwMIoBAAAAMDCKAQAAADAwigEAAAAwMIc2HQCrc8eJqyaeX3f09IYiAQAAoE+cGQAAAAADoxgAAAAAA6MYAAAAAAOjGAAAAAAD4waCB8DOGwEmbgYIAADA/JwZAAAAAAOjGAAAAAAD4zKBLXHk+G2PG3fLsQ0EAgAAQO8pBszJdfwAAAAcVC4TAAAAgIFRDAAAAICBUQwAAACAgVEMAAAAgIFRDAAAAICBUQwAAACAgfHTgjM6cvy2iee3HNtQIAAAALAgZwYAAADAwCgGAAAAwMAoBgAAAMDAKAYAAADAwCgGAAAAwMAoBgAAAMDAKAYAAADAwCgGAAAAwMAc2nQA8zpy/LbHjbvv5hdvIBIAAAA4WHpRDNj5j71/6vtNIQYAAOBgc5kAAAAADIxiAAAAAAyMYgAAAAAMTC/uGQCwKttwTxL36QAAYNmWWgw4CAesByFGgGXQ3x0821C8AliVeT/X9K2wu2qt7X+mqi8meWQJ738oyeeXsJxlEtNsFonpcGttLZeo7JGr29auq9THuNYVU1/ydVP6uO2nOUixJsuPty+5etC2w6ppj0nn2kO+9o+2mLT2XE3Om6+2zyTtMWmufJ2rGLAsVXVXa+2ajQWwCzHNpo8x7Vcf16GPMSX9jKuPMW2jg9TOBynW5ODFO6ttXa95aY9JfWuPvsWzSdpiUt/ao2/xbJr2mDRve7iBIAAAAAyMYgAAAAAMzKaLAW/e8PvvRkyz6WNM+9XHdehjTEk/4+pjTNvoILXzQYo1OXjxzmpb12te2mNS39qjb/FskraY1Lf26Fs8m6Y9Js3VHhu9ZwAAAACwfps+MwAAAABYs5UVA6rqhVX1sao6VVXHd3n966vqPVX1war6cFW9aOy113bzfayqXrDpmKrqSFU9UlV3d483rTGmy6vqji6ek1V12dhrN1bVvd3jxp7E9IWxdrp1WTEtIeaLquqd3evvraojY6/tmm97LXOVcVXVP62q91fVPd3fo2PznOyWea6dn7qmmKbuB1X1vC7WU1X1s1VV+4lpwbheMRbT3VX1xar65u61hdpqKKrq0qq6vetLbq+qS6ZMt2ufs452XsU+viqr2Mf6bB/5s5HPh3VZJEe3zQxtcVNVfXosH161gRh602esw9D6pb3M0B7fWVUfqKrPV9XLd7y2kuPvXWLQt0bfOm4lfWtrbemPJE9IcjrJlUkuTPKhJFfvmObNSX68G746yX1jwx9KclGSK7rlPGHDMR1J8pENtdOvJrmxGz6a5Je64UuTnOn+XtINX7LJmLrnf7uKnFpCzK9O8qZu+IYk7zxfvs2yzBXH9dwkX9cNf2OSs2PznExyzQbaaup+kOR9Sb41SSX5nSTfva64dkzzTUlOL6OthvRI8tNJjnfDx5O8YZdppvY5q27nVezjPY116j7W58cs+dO9tvbPhzW2wVL6sG14zNgWNyV5Yx+3x7r7jB5tk63ql5bQHkeS/P0kb0vy8rHxKzn+nhKnvlXfut+22HffuqozA56f5FRr7Uxr7e+SvCPJS3dM05J8VTf81Uk+2Q2/NMk7WmuPttY+keRUt7xNxrQqs8R0dZIT3fB7xl5/QZLbW2t/2Vr7qyS3J3nhhmPalFlifmmSt3bD70pyXVVVpufbLMtcWVyttQ+21s7l30eTHK6qi/b5/kuNadoCq+rpSb6qtXZnG/VEb0vysg3F9QPdvOzPeNu+Nbtvv1X1ObNYxT7ex1gPqlnyZ9sNcbtPs4zPz3XE0Jc+Yx3k56Q926O1dl9r7cNJvrhj3nV+Fupb5e64lfStqyoGPCPJn409v78bN+7fJfkXVXV/kv+a5F/tY951x5QkV9To8oHfq6rvWEI8s8b0oSTf1w1/b5InVtWTZ5x33TElycVVdVdV3VlVL1tCPLOYJeYvTdNa+3ySzyY5Xzsuo30XiWvc9Uk+0Fp7dGzcL3Sn//zkPju8RWPabT94Rrec8y1z1XGd8/1J3r5j3LxtNSRPa6090A3/eZKn7TLNXttole28in18VVaxj/XdLPmTbObzYV2W1Ydtg1n3uetrdLnhu6rqmRuIoS99xjoMsV86n0W28TrzQ9+qbx23kr710LKim8MPJPnF1tp/rKp/nOSXquobNxjP+WJ6IMnXt9YerKrnJfmNqnpOa+2v1xDTTyR5Y1XdlOT3k5xN8oU1vO/5nC+my1trZ6vqyiQnquqe1trpDcV54FXVc5K8IcmxsdGv6Nr4iUl+LckPZvRt/Krtuh+s4X1nUlX/KMnDrbWPjI3eVFv1TlW9O8nX7vLS68aftNZaVbV9Ll47L8cmP2vOa0n54/OBc34rydtba49W1Y9m9K3e0T3mYTN62y9tA30rS7bvvnVVZwacTTJeibisGzfulUl+JUlaa/89ycVJnjLjvGuNqTtV7MFu/Pszul7j2euIqbX2ydba97XWnpuuY2itfWbG9Vl3TGmtne3+nsnoOuLnLiGmvczSFl+apqoOZXQZyIPnmXcZ7btIXKnRjRl/PckPjXfqY238N0l+Ofs7fXHumM6zH5ztlnO+Za4srrHXb8iOswIWbKut0lr7rtbaN+7y+M0kf9Fd7nHuso9P7bKIqdtoDe28in18VVaxj23cEvJnU58P67KMPmxbzHIc8eDY2W5vSfK8dceQ/vQZ67CV/dICFtnGS80Pfeue9K2PWU3f2lZzg4NDGd1Q44o8doOD5+yY5neS3NQN/72Mrs+vJM/J5I1bzmQ5NxBcJKavORdDRjdtOJvk0jXF9JQkF3TD/yHJ67vhS5N8IqObl1zSDW86pkuSXDQ2zb3Z5033Vhjzv8zkzUV+pRveNd9mWeaK43pSN/337bLMp3TDX5bRtVE/tqaYpu4HefwNBF+0rrbqnl/QxXPlstpqSI8k/2cmb1L007tMs2ufs452XsU+vsK2XMk+1ufHjPmzkc+HNbbBQn3YNj1mbIunjw1/b5I7+7I91t1n9GibbFW/tGh7jE37i3n8DQSXfvw95b31rfrW/bbFvvvWVQb8oiQfz6iC+Lpu3OuTvKQbvjrJH3YrcneSY2Pzvq6b72PZ513JVxFTRtdsf7Qb94Ek37PGmF7e7dgfz6jCc9HYvD+S0Y1tTiX54U3HlOTbktzTtd89SV65jp1jxpgvzuhXEE5l9I/r+D+Nu+bbbstcV1xJ/m2Sh7qcO/d4apKvSPL+JB/ucvJnss+DlAVimrofJLkmyUe6Zb4xSa15G16bHR3eMtpqKI+Mrq27o9uv353HijzXJHnL2HSP63PW1c6r2MdX2J5L38f6/Jglf7LBz4c1tsPcObptjxna4qe6XP9QRjci/oY+bY919xk92SZb1S8toT3+YUbXZD+U0bfMHx2bdyXH37vEqG+dbVvpWxfoW6ubEQAAABiIVd0zAAAAAOgpxQAAAAAYGMUAAAAAGBjFAAAAABgYxQAAAAAYGMUAAAAAGBjFAAAAABgYxQAAAAAYmP8f4j23zBqHG34AAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 1036.8x518.4 with 72 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "n_features = X.shape[1]\n",
    "\n",
    "fig, axes = plt.subplots(n_features, n_features, figsize=(n_features * 2.4, n_features * 1.2))\n",
    "\n",
    "for i in range(n_features):\n",
    "    for j in range(n_features):\n",
    "        result_missing = adj_matrices_list[:, :, i, j].flatten()\n",
    "        result_no_missing = bs_result.adjacency_matrices_[:, i, j].flatten()\n",
    "        \n",
    "        xrange = (\n",
    "            min(*result_missing, *result_no_missing),\n",
    "            max(*result_missing, *result_no_missing)\n",
    "        )\n",
    "        \n",
    "        hist, edges = np.histogram(result_missing, range=xrange)\n",
    "        hist2, edges2 = np.histogram(result_no_missing, range=xrange)\n",
    "        \n",
    "        width = (edges[1] - edges[0]) / 3\n",
    "        \n",
    "        ax = axes[i, j]\n",
    "        ax.bar(edges[:-1], hist, width=width, align=\"edge\", color=\"tab:blue\")\n",
    "        ax.set_yticks([])\n",
    "        \n",
    "        ax = axes[i, j].twinx()\n",
    "        ax.bar(edges2[:-1] + width, hist2, width=width, align=\"edge\", color=\"tab:olive\")\n",
    "        ax.set_yticks([])\n",
    "        \n",
    "        ax.axvline(m[i, j], color=\"red\")\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ad2a0415-b19e-4131-be24-b140052ef920",
   "metadata": {},
   "source": [
    "Red vertical lines indicate true values. The [i, j] of this scatter matrix indicates the distribution of the [i, j] of the adjacency matrix. Blue histograms are results with missing data, orange histograms are results with no missing data. \n",
    "\n",
    "#### Median value of each element of the matrix\n",
    "\n",
    "Median values of each element of the adjacency matrix estimated on the data with missing values are as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "5874ec01-c068-4e9b-a609-b14fb8127c86",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.   ,  0.   ,  0.   ,  0.882,  0.   ,  0.   ],\n",
       "       [ 0.555,  0.   ,  0.381,  0.   ,  0.   ,  0.   ],\n",
       "       [ 0.   ,  0.   ,  0.   ,  0.918,  0.   ,  0.   ],\n",
       "       [ 0.   ,  0.   ,  0.   ,  0.   ,  0.   ,  0.   ],\n",
       "       [ 1.025,  0.   , -0.147,  0.   ,  0.   ,  0.   ],\n",
       "       [ 0.892,  0.   ,  0.   ,  0.   ,  0.   ,  0.   ]])"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.median(adj_matrices_list, axis=(0, 1))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "439b3319-6c9f-4c7e-8649-5a0dc63cbc7d",
   "metadata": {},
   "source": [
    "The average value of each element of the adjacency matrix estimated on the data with no missing values is as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "2dde5539-ddc6-47c9-a7eb-b8c4681d5b1e",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.   ,  0.   ,  0.   ,  0.882,  0.   ,  0.   ],\n",
       "       [ 0.546,  0.   ,  0.383,  0.   ,  0.   ,  0.   ],\n",
       "       [ 0.   ,  0.   ,  0.   ,  0.912,  0.   ,  0.   ],\n",
       "       [ 0.   ,  0.   ,  0.   ,  0.   ,  0.   ,  0.   ],\n",
       "       [ 1.002,  0.   , -0.136,  0.   ,  0.   ,  0.   ],\n",
       "       [ 0.89 ,  0.   ,  0.   ,  0.   ,  0.   ,  0.   ]])"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.median(bs_result.adjacency_matrices_, axis=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "79397436-7886-4c40-aaef-2fea9b46ac07",
   "metadata": {},
   "source": [
    "Comparing each result with the adjacency matrix used in the data generation, it is clear that they are all able to estimate the presence of true edges."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "54127232-2cbb-4030-af98-6c66a1c7c7ea",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.   ,  0.   ,  0.   ,  0.895,  0.   ,  0.   ],\n",
       "       [ 0.565,  0.   ,  0.377,  0.   ,  0.   ,  0.   ],\n",
       "       [ 0.   ,  0.   ,  0.   ,  0.895,  0.   ,  0.   ],\n",
       "       [ 0.   ,  0.   ,  0.   ,  0.   ,  0.   ,  0.   ],\n",
       "       [ 0.991,  0.   , -0.124,  0.   ,  0.   ,  0.   ],\n",
       "       [ 0.895,  0.   ,  0.   ,  0.   ,  0.   ,  0.   ]])"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "m"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d4a21bb3-1598-41be-a04a-38d77b441156",
   "metadata": {
    "tags": []
   },
   "source": [
    "## Getting bootstrap samples and imputed bootstrap samples\n",
    "\n",
    "To get bootstrap samples, do the following:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "06d4501c-5073-40c8-b1fa-c34b51d80728",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(30, 1000, 6)"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "bootstrap_samples = np.array([X_mcar[indices] for indices in resampled_indices])\n",
    "bootstrap_samples.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5a2ca9e6-dfb7-4e58-9fbe-2e665a0f5385",
   "metadata": {},
   "source": [
    "`bootstrap_samples` contains NaNs. For example, the number of NaNs that the first bootstrap sample has is as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "fb665ef2-433b-459b-9ae4-6891ea8d8d48",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0,  0,  0,  0,  0, 85])"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.isnan(bootstrap_samples[0]).sum(axis=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "38a38d46-59e7-4075-b361-356efcd8bc77",
   "metadata": {},
   "source": [
    "`imputation_results` stores the result of the imputation. When X_mcar[i, j] is NaN, the complementary values are stored in the imputation_results[i, j], and when X_mcar[i, j] is non-NaN, NaN is stored in the imputation_results[i, j].\n",
    "\n",
    "The first 8 rows of the first imputation result for the first bootstrap sample is as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "df79359d-4b8e-4904-a985-92788a147c05",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[  nan,   nan,   nan,   nan,   nan,   nan],\n",
       "       [  nan,   nan,   nan,   nan,   nan,   nan],\n",
       "       [  nan,   nan,   nan,   nan,   nan,   nan],\n",
       "       [  nan,   nan,   nan,   nan,   nan,   nan],\n",
       "       [  nan,   nan,   nan,   nan,   nan,   nan],\n",
       "       [  nan,   nan,   nan,   nan,   nan,   nan],\n",
       "       [  nan,   nan,   nan,   nan,   nan, 0.188],\n",
       "       [  nan,   nan,   nan,   nan,   nan,   nan]])"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "imputation_results[0, 0, :8]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c137119a-d52f-4eab-a71b-433cf3567e2f",
   "metadata": {
    "tags": []
   },
   "source": [
    "The first 8 rows of each of the first bootstrap_sample is as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "06002377-094e-4d27-96a8-18d9b0be3810",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-0.193, -0.072,  1.018,  0.328, -0.863, -0.622],\n",
       "       [-1.972, -1.887, -1.687, -1.489, -1.216, -1.622],\n",
       "       [-0.444, -1.18 , -0.49 ,  0.269, -0.483, -0.011],\n",
       "       [ 0.531,  0.996,  1.201,  0.968,  1.035,  0.708],\n",
       "       [-2.202, -1.915, -0.988, -1.704, -1.42 , -1.531],\n",
       "       [ 1.6  ,  1.756,  0.591,  1.127,  1.494,  0.743],\n",
       "       [ 0.441, -0.17 , -0.933, -0.244,  1.117,    nan],\n",
       "       [ 0.746,  1.629,  1.329,  0.976,  0.302,  0.756]])"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "bootstrap_samples[0, :8]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7cad4949-ee08-4853-9a6b-23ec45eff77a",
   "metadata": {
    "tags": []
   },
   "source": [
    "`bootstrap_samples` is the data before the imputation is performed, so the missing values are still NaN. \n",
    "\n",
    "To obtain the imputed bootstrap_samples, do the following:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "686af4e3-27da-40d6-8f48-29081aaeb06a",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "imputed_bootstrap_samples = []\n",
    "\n",
    "for i in range(n_sampling):\n",
    "    imputeds = []\n",
    "    for j in range(n_repeats):\n",
    "        pos = ~np.isnan(imputation_results[i, j])\n",
    "        \n",
    "        imputed = bootstrap_samples[i].copy()\n",
    "        imputed[pos] = imputation_results[i, j, pos]\n",
    "        \n",
    "        imputeds.append(imputed)\n",
    "    imputed_bootstrap_samples.append(imputeds)\n",
    "imputed_bootstrap_samples = np.array(imputed_bootstrap_samples)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0a65d2dc-b59b-438f-a7ed-eee43829b3cd",
   "metadata": {},
   "source": [
    "It is confirmed that NaN was imputed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "8f5b0640-ef04-405b-80b9-006d49c200e3",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-0.193, -0.072,  1.018,  0.328, -0.863, -0.622],\n",
       "       [-1.972, -1.887, -1.687, -1.489, -1.216, -1.622],\n",
       "       [-0.444, -1.18 , -0.49 ,  0.269, -0.483, -0.011],\n",
       "       [ 0.531,  0.996,  1.201,  0.968,  1.035,  0.708],\n",
       "       [-2.202, -1.915, -0.988, -1.704, -1.42 , -1.531],\n",
       "       [ 1.6  ,  1.756,  0.591,  1.127,  1.494,  0.743],\n",
       "       [ 0.441, -0.17 , -0.933, -0.244,  1.117,  0.188],\n",
       "       [ 0.746,  1.629,  1.329,  0.976,  0.302,  0.756]])"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "imputed_bootstrap_samples[0, 0, :8]"
   ]
  }
 ],
 "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.10"
  },
  "widgets": {
   "application/vnd.jupyter.widget-state+json": {
    "state": {},
    "version_major": 2,
    "version_minor": 0
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}