{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a href=\"https://colab.research.google.com/github/gtbook/gtsam-examples/blob/main/GaussianMRFExample.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Gaussian Markov Random Fields\n",
    "\n",
    "A Markov random field or MRF is an *undirected* graphical model where the energy is given by the som of unary and multi-variable potentials $\\psi$:\n",
    "\n",
    "$$\n",
    "E(X) = \\sum_j \\psi_j(x_j) + \\sum_c \\psi_c(X_c)\n",
    "$$\n",
    "\n",
    "where $j$ ranges over the variable nodes $x_j$ and $c$ indexes over cliques in the graph of variables. We use $X_c$ to denote all variables associated with clique $c$. A typical case is a pairwise Gaussian MRF, where\n",
    "\n",
    "- the $\\psi_j(x_j)$ are quadratic potentials $\\|A_j x_j-b_j\\|^2_{R_j}$ associated with measurements $b_j$.\n",
    "- the $\\psi_c(X_c)=\\psi_c(x_{j_1},x_{j_2})$ are pairwise potentials $\\|A_{j_1} x_{j_1}+A_{j_2} x_{j_2}-b_c\\|^2_{Q_c}$ over pairwise cliques $X_c=\\{x_{j_1},x_{j_2}\\}$, typically enforcing smoothness."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It will not surprise anyone, that Gaussian MRF's can easily and advantageously be represented using Gaussian factor graphs. Indeed, we can rewrite the energy $E(X)$ as\n",
    "\n",
    "$$\n",
    "E(X) = \\sum_i \\phi_i(X_i)\n",
    "$$\n",
    "\n",
    "where the factors $\\phi_i(X_i)$ encompass both unary and binary factors, and of course we can add arbitrary arity factors as well.\n",
    "\n",
    "Below we create a small \"image \"example, and show how we can solve it the same tools as the ones we used in the Kalman smoothing examples."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Note: you may need to restart the kernel to use updated packages.\n"
     ]
    }
   ],
   "source": [
    "%pip -q install gtbook  # also installs latest gtsam pre-release"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import plotly.express as px\n",
    "try:\n",
    "    import google.colab\n",
    "except:\n",
    "    import plotly.io as pio\n",
    "    pio.renderers.default = \"png\"\n",
    "    \n",
    "import gtsam\n",
    "from gtbook.display import show\n",
    "from gtbook.gaussian import sample_bayes_net\n",
    "from gtsam import noiseModel"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## A Small \"Image\" MRF\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We first define the image dimensions and set up some keys referring to rows by a letter and columns by an integer:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "M, N = 3, 4  # try playing with this\n",
    "row_symbols = [chr(ord('a')+row) for row in range(M)]\n",
    "keys = {(row, col): gtsam.symbol(row_symbols[row], col+1)\n",
    "       for row in range(M) for col in range(N)}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We create random data, and define a noise model for the data with the right  standard deviation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "data_model = isotropic dim=1 sigma=0.5\n",
      "\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAArwAAAH0CAYAAADfWf7fAAAgAElEQVR4Xu2dd5xV1dW/FwwwVAELRmMnihJLTMRETSwxir3FErBhQ6yogCBIBBUVaSqKSFEERUSsKAqKokRjiRqjaSbGEhsgAmJB6u9zb35OMpfhZc8+Z9+19ubhP52zyn6+O+/75ORyp86qVatWCX8gAAEIQAACEIAABCCQKIE6CG+iyXIsCEAAAhCAAAQgAIEiAYSXiwABCEAAAhCAAAQgkDQBhDfpeDkcBCAAAQhAAAIQgADCyx2AAAQgAAEIQAACEEiaAMKbdLwcDgIQgAAEIAABCEAA4eUOQAACEIAABCAAAQgkTQDhTTpeDgcBCEAAAhCAAAQggPByByAAAQhAAAIQgAAEkiaA8CYdL4eDAAQgAAEIQAACEEB4uQMQgAAEIAABCEAAAkkTQHiTjpfDQQACEIAABCAAAQggvNwBCEAAAhCAAAQgAIGkCSC8ScfL4SAAAQhAAAIQgAAEEF7uAAQgAAEIQAACEIBA0gQQ3qTj5XAQgAAEIAABCEAAAggvdwACEIAABCAAAQhAIGkCCG/S8XI4CEAAAhCAAAQgAAGElzsAAQhAAAIQgAAEIJA0AYQ36Xg5HAQgAAEIQAACEIAAwssdgAAEIAABCEAAAhBImgDCm3S8HA4CEIAABCAAAQhAAOHlDkAAAhCAAAQgAAEIJE0A4U06Xg4HAQhAAAIQgAAEIIDwcgcgAAEIQAACEIAABJImgPAmHS+HgwAEIAABCEAAAhBAeLkDEIAABCAAAQhAAAJJE0B4k46Xw0EAAhCAAAQgAAEIILzcAQhAAAIQgAAEIACBpAkgvEnHy+EgAAEIQAACEIAABBBe7gAEIAABCEAAAhCAQNIEEN6k4+VwEIAABCAAAQhAAAIIL3cAAhCAAAQgAAEIQCBpAghv0vFyOAhAAAIQgAAEIAABhJc7AAEIQAACEIAABCCQNAGEN+l4ORwEIAABCEAAAhCAAMLLHYAABCAAAQhAAAIQSJoAwpt0vBwOAhCAAAQgAAEIQADh5Q5AAAIQgAAEIAABCCRNAOFNOl4OBwEIQAACEIAABCCA8HIHIAABCEAAAhCAAASSJoDwJh0vh4MABCAAAQhAAAIQQHi5AxCAAAQgAAEIQAACSRNAeJOOl8NBAAIQgAAEIAABCCC83AEIQAACEIAABCAAgaQJILxJx8vhIAABCEAAAhCAAAQQXu4ABCAAAQhAAAIQgEDSBBDepOPlcBCAAAQgAAEIQAACCC93AAIQgAAEIAABCEAgaQIIb9LxcjgIQAACEIAABCAAAYSXOwABCEAAAhCAAAQgkDQBhDfpeDkcBCAAAQhAAAIQgADCyx2AAAQgAAEIQAACEEiaAMKbdLwcDgIQgAAEIAABCEAA4c3hDixdukwWLPpSWm3YQurUqZNDR1pAAAIQgAAEIAABCORFAOHNQHLVqlVy6/hH5JY7Hix2Wb9FM7n5motkl7ata+w6c/ZrcmHfm1b72WszRktlg/oZNqEUAhCAAAQgAAEIQGBNBBDeDHfj9bf+ISedP0AmDO8tO22/jdw09gF5bObv5al7h0rduqu/6X1q9qty2TWjZcro/tWmbvH9VrwZzpADpRCAAAQgAAEIQOD/IoDwZrgfQ0ZOlr/+830ZM7hHscvczxbKfsdeVBTaHbbdcrXOBeHtP2SczH5oeIaplEIAAhCAAAQgAAEI1IYAwlsbWiXPdr/yVmnZvKn06Xpy1U9+uG8nGXHtxbLPHrvUKLxd+w6XI9vvJZWVDWS3XdpI+33bSb2KigxbUAoBCEAAAhCAAAQgwBveWt6Bjz/9TB6b+eIaq0769YHSqGED6dxjsLRpvYV063J81bPtDu4i/bp3kkP3/9lq9W/+7V2ZPutlad6siXw8Z75MfuQZ6Xj0/tWEuZar8jgEIAABCEAAAhCAwFoI8Ia3BkDvfzhHJj389BrRXXD60dK4UUMpvOEt/EW13hee5PSGt7ThA9Oek77X3y5vzBxbfMv7zilduLAJExg7aZ+ET8fRevaZBISECTRo8XXCp+Nojbo+CYTECSC8GQIufIb37+98IKMGdS92WdtneEtHzX7pTenSc4i8On2UNKxsgPBmyCKGUoQ3hpT8d0R4/dnFUInwxpCS/44Irz+7WCoR3gxJ/fdbGvrITjtsIzeOmSLTZr5Y9S0N4yY/IYWvIit8i0Phz8QHZ0qb1ptL2+22kkWLv5QeV46U+vUq5PZhPYs/5w1vhjAiKEV4Iwgpw4oIbwZ4EZQivBGElGFFhDcDvEhKEd4MQRW+h/fmOx6UkeMfKXYpfMxh1KBusuuO2xb/edCISTJ56ix55fGRxX8eettkGXvPtKqJO7dtLYP6dpHNNtkI4c2QQyylCG8sSfntifD6cYulCuGNJSm/PRFeP24xVSG8OaS15Nul8vmCL+R7rTao8ft3/3dE4dl58xdKsyaNpUXzptWm84Y3hzAMt0B4DYeTw2oIbw4QDbdAeA2Hk8NqCK87xJUrV0nhhV9FRV33IgNPIrwGQvhuBYTXUBgBVkF4A0A11BLhNRRGgFUQ3gBQDbVEeN3CKIhuvyHjig/3736aW5GRpxBeI0EU1kB4DYURYBWENwBUQy0RXkNhBFgF4Q0A1VBLhHftYRS+VvXqGybI5wsXy7GH7YPwrh0ZT6yJAMKb9t1AeNPOF+FNO1+EN+18Ed615/v1N9/KF19+JcNG3Vf8Zine8K6dGU+sgQDCm/bVQHjTzhfhTTtfhDftfBFe93yvHDZeVqxYgfC6I+PJUgIIb9p3AuFNO1+EN+18Ed6080V43fNFeN1Z8SRveNfJO4Dwph07wpt2vghv2vlaFt5F326jAr955b9qnIvwqsSR1lDe8KaVZ+lpEN6080V4084X4U07X9PCu0RJeBsivGnfesXTIbyK8MswGuEtA2TFEQivIvwyjEZ4ywBZcYRp4f1aSXgbVxfeFStWysqVK+XqGyfI8uUrpF+3TlJRUbHW3z+gGGu10XwtmZUk+FoyQ0mEWQXhDcPVSleE10oSYfZAeMNwtdLVsvB+8ZWO8K7XpLrwTn7kGek/9M5qkV116elyzCF7W4nx/9wD4TUUE294DYURYBWENwBUQy0RXkNhBFgF4Q0A1VBL08L7pZLwNq35Iw2GYqvVKghvrXCFfRjhDctXuzvCq51A2PkIb1i+2t0RXu0Ews43LbyLlYS3GcIb9tatw90R3rTDR3jTzhfhTTtfhDftfC0L7+JFOsLbrDnCm/atVzwdwqsIvwyjEd4yQFYcgfAqwi/DaIS3DJAVR5gW3oVKwtsC4VW8kmmPRnjTzhfhTTtfhDftfBHetPM1LbwLWqvAb9byHZW5oYbyGd5QZD36Irwe0CIqQXgjCstjVYTXA1pEJQhvRGF5rGpZeL/8XEd4m66P8HpcJUpcCCC8LpTifQbhjTc7l80RXhdK8T6D8MabncvmpoX3MyXh3RDhdbk7PONBAOH1gBZRCcIbUVgeqyK8HtAiKkF4IwrLY1XTwjtPSXg3Qng9rhIlLgQQXhdK8T6D8MabncvmCK8LpXifQXjjzc5lc8vC+9VcHeFt0grhdbk7PONBAOH1gBZRCcIbUVgeqyK8HtAiKkF4IwrLY1XTwjtHSXg3Rng9rhIlLgQQXhdK8T6D8MabncvmCK8LpXifQXjjzc5lc9PC+4mS8G6C8LrcHZ7xIIDwekCLqAThjSgsj1URXg9oEZUgvBGF5bGqZeH9+uMfeJwoe0njTf+ZvYmhDnwtmaEwEF5DYQRYBeENANVQS4TXUBgBVkF4A0A11NK08H6kJLzfR3gNXdG0VkF408qz9DQIb9r5Irxp54vwpp2vZeH95kMd4W20GcKb9q1XPB3Cqwi/DKMR3jJAVhyB8CrCL8NohLcMkBVHmBbeD5SEdwuEV/FKpj0a4U07X4Q37XwR3rTzRXjTzte08L6vJLxbIrxp33rF0yG8ivDLMBrhLQNkxREIryL8MoxGeMsAWXGEZeFd8p6O8DbcCuFVvJJpj0Z4084X4U07X4Q37XwR3rTzNS287yoJ79YIb9q3XvF0CK8i/DKMRnjLAFlxBMKrCL8MoxHeMkBWHGFaeP+lJLzbILyKVzLt0Qhv2vkivGnni/CmnS/Cm3a+loX3239uqwK/8gf/UJkbaijfwxuKrEdfhNcDWkQlCG9EYXmsivB6QIuoBOGNKCyPVU0L7z+UhHdbhNfjKlHiQgDhdaEU7zMIb7zZuWyO8LpQivcZhDfe7Fw2Ny28bysJ73YIr8vd4RkPAgivB7SIShDeiMLyWBXh9YAWUQnCG1FYHqtaFt6lf9cR3gZtEF6Pq0SJCwGE14VSvM8gvPFm57I5wutCKd5nEN54s3PZ3LTw/lVJeHdAeF3uDs94EEB4PaBFVILwRhSWx6oIrwe0iEoQ3ojC8ljVtPD+RUl42yK8HleJEhcCCK8LpXifQXjjzc5lc4TXhVK8zyC88Wbnsrll4V32Zx3hrf9DhNfl7vCMBwGE1wNaRCUIb0RheayK8HpAi6gE4Y0oLI9VTQvvW9t5nCh7Sf0d387exFAHvpbMUBgIr6EwAqyC8AaAaqglwmsojACrILwBoBpqaVp4/6QkvDsjvIauaFqrILxp5Vl6GoQ37XwR3rTzRXjTztey8C5/Q0d46+2C8KZ96xVPh/Aqwi/DaIS3DJAVRyC8ivDLMBrhLQNkxRGmhfePSsL7I4RX8UqmPRrhTTtfhDftfBHetPNFeNPO17Twvq4kvLsivGnfesXTIbyK8MswGuEtA2TFEQivIvwyjEZ4ywBZcYRl4V3xqo7wVvwE4VW8kmmPRnjTzhfhTTtfhDftfBHetPM1Lbx/UBLe3RDetG+94ukQXkX4ZRiN8JYBsuIIhFcRfhlGI7xlgKw4wrTwvqIkvO0QXsUrmfZohDftfBHetPNFeNPOF+FNO1/TwvtyGxX4Fbv/XWVuqKF8D28osh59EV4PaBGVILwRheWxKsLrAS2iEoQ3orA8VjUtvC8pCe9PEV6Pq0SJCwGE14VSvM8gvPFm57I5wutCKd5nEN54s3PZ3LTw/n57lyPk/kzFHn/LvadmQ97watIvmY3wGgojwCoIbwCohloivIbCCLAKwhsAqqGWpoX3BSXh3RPhNXRF01oF4U0rz9LTILxp54vwpp0vwpt2vqaF9/kdVOBX7PVXlbmhhvKGNxRZj74Irwe0iEoQ3ojC8lgV4fWAFlEJwhtRWB6rmhbe37X1OFH2koqf/yV7E0MdEF6FMBZ/+bUsX7FCWjZvVm06wqsQRhlHIrxlhK0wCuFVgF7GkQhvGWErjDItvLOVhPcXCK/CVUxj5NffLJGeV98mTz//evFAO7dtLcOvvlA2XL958Z8R3jRyXtMpEN6080V4084X4U07X9PC+9wPVeBX7P1nlbmhhvKGNxTZGvqOmfiY3Dd1lkwY3kcaNWwg5/QaJltvsYlcdenpCG8Zc9AahfBqkS/PXIS3PJy1piC8WuTLM9e08D6rJLz7ILzluX0JTjn2rCuk/b7t5KwTDyuebvqsl+WSfiPkrWfukDp16vCGN8HM//dICG/aASO8aeeL8Kadr2nhfWZHFfgV+72lMjfUUN7whiJbQ992B3eRq3ueUZTewp+/vP2eHNe5n7ww9RZp3qwJwlvGLDRGIbwa1Ms3E+EtH2uNSQivBvXyzTQtvE8rCe8vEd7y3cCEJq1atUp23O80GXHtxbLPHrsUT/bOex/JEZ36yFP3DpFNNt4A4U0o75qOgvCmHTDCm3a+CG/a+ZoW3pk7qcCv2P9NlbmhhvKGNxTZNbzhHdDrTDlwn914w1tG7lZGIbxWkgizB8IbhquVrgivlSTC7GFaeJ9SEt5fIbxhbts60LXwGd6D9ttdzux4aPG0fIZ3HQj9f46I8KadN8Kbdr4Ib9r5mhbeJ3dWgV9xwJ9U5oYayhveUGRr6Dv67kdlyqPPFr+loXGjSunScyjf0lBG/tqjEF7tBMLOR3jD8tXujvBqJxB2vmnhnaEkvAcivGFvXcLdv/p6iXS/8lZ57sU3iqfcsc3WMnxAV2m1YYviP/M9vAmHLyIIb9r5Irxp54vwpp2vaeGd/p+/91PuPxXt/+MqqfzhDa9CkosWfyXLli2v+oUT362A8CqEUcaRCG8ZYSuMQngVoJdxJMJbRtgKo0wL7+NKwnswwqtwFdeNkQhv2jkjvGnni/CmnS/Cm3a+poV32o9U4Fcc8keVuaGG8oY3FFmPvgivB7SIShDeiMLyWBXh9YAWUQnCG1FYHquaFt7HlIT3UITX4ypR4kIA4XWhFO8zCG+82blsjvC6UIr3GYQ33uxcNrcsvMsf3dXlCLk/U++w13PvqdmQN7ya9EtmI7yGwgiwCsIbAKqhlgivoTACrILwBoBqqKVp4Z2qJLyHI7yGrmhaqyC8aeVZehqEN+18Ed6080V4087XtPA+8mMV+PWOeE1lbqihvOENRdajL8LrAS2iEoQ3orA8VkV4PaBFVILwRhSWx6qmhffhn3icKHtJvSNfzd7EUAeE11AYCK+hMAKsgvAGgGqoJcJrKIwAqyC8AaAaamlaeB9UEt6jEV5DVzStVRDetPIsPQ3Cm3a+CG/a+SK8aedrWngf2E0Ffr1j/qAyN9RQ3vCGIuvRF+H1gBZRCcIbUVgeqyK8HtAiKkF4IwrLY1XTwnu/kvD+GuH1uEqUuBBAeF0oxfsMwhtvdi6bI7wulOJ9BuGNNzuXzU0L75R2LkfI/Zl6x76Se0/Nhrzh1aRfMhvhNRRGgFUQ3gBQDbVEeA2FEWAVhDcAVEMtTQvvfUrCexzCa+iKprUKwptWnqWnQXjTzhfhTTtfhDftfE0L7+TdVeDXO/5llbmhhvKGNxRZj74Irwe0iEoQ3ojC8lgV4fWAFlEJwhtRWB6rmhbee5WE9wSE1+MqUeJCAOF1oRTvMwhvvNm5bI7wulCK9xmEN97sXDY3LbyTfupyhNyfqfebl3LvqdmQN7ya9EtmI7yGwgiwCsIbAKqhlgivoTACrILwBoBqqKVp4Z2oJLwdEV5DVzStVRDetPIsPQ3Cm3a+CG/a+SK8aedrWnjv/pkK/HonvqgyN9RQ3vCGIuvRF+H1gBZRCcIbUVgeqyK8HtAiKkF4IwrLY1XTwnuXkvCehPB6XCVKXAggvC6U4n0G4Y03O5fNEV4XSvE+g/DGm53L5qaFd8IeLkfI/Zl6J/8+956aDXnDq0m/ZDbCayiMAKsgvAGgGmqJ8BoKI8AqCG8AqIZamhbe8UrCewrCa+iKprUKwptWnqWnQXjTzhfhTTtfhDftfE0L7517qsCvd+oLKnNDDeUNbyiyHn0RXg9oEZUgvBGF5bEqwusBLaIShDeisDxWNS2845SEtxPC63GVKHEhgPC6UIr3GYQ33uxcNkd4XSjF+wzCG292LptbFt5ld+zlcoTcn6l/2vO599RsyBteTfolsxFeQ2EEWAXhDQDVUEuE11AYAVZBeANANdTStPDe/nMVUvVP/53K3FBDEd5QZD36Irwe0CIqQXgjCstjVYTXA1pEJQhvRGF5rGpZeJeO1RHeBmcgvB5XiRIXAgivC6V4n0F4483OZXOE14VSvM8gvPFm57K5aeEd8wuXI+T+TIMzZ9fY87PPF0njRg2lcaPK3GeGbMgb3pB0a9kb4a0lsMgeR3gjC6yW6yK8tQQW2eMIb2SB1XJd08I7eu9aniafxxuc9Vy1Rh98NEe69Bwq7384p/jvjzlkb/ntJadK/XoVNQ4ceMs9Mv6+6dV+tuuO28pdN/fJZ8FadkF4awks5OMIb0i6+r0RXv0MQm6A8Iakq98b4dXPIOQGloX321E6wlvZubrwdu4xWJo2aSQDep0ln86dL8ef3V9+e/EpcviBNX+LxHU3T5R/fzxXLj23Q1V0lZX15XsbrR8yyjX2RnhVsNc8FOE1FEaAVRDeAFANtUR4DYURYBWENwBUQy1NC+9t+6iQqjz72aq5ixZ/JXsefl7x7WzhLW3hz4AbJ8incz+X4QO61rhfQXgXfvGlXNe7s8r+pUMRXhMx/GcJhNdQGAFWQXgDQDXUEuE1FEaAVRDeAFANtTQtvCP3VSFV2WVW1dx33vtIjujUR2bdf4NstEGL4r+fMGWGPDz9eZkyuv8ahXfGs6/Iz37cVlo2bya//PmP5Sc7b6dylsJQhFcN/eqDEV5DYQRYBeENANVQS4TXUBgBVkF4A0A11NKy8C65VUd4G57zX+F9/a1/yEnnD5AXpt4izZs1KSY3eeosGTn+YXn6vmE1Jjl1xgvy3oefSmWD+vLW39+VmbNfk6H9zpX2++6ukjzCq4K95qEIr6EwAqyC8AaAaqglwmsojACrILwBoBpqaVp4b9lPhVTD855Z7Q3vsw/cKBuu39zpDW/p0r2uGSULFy2WkQO7qZwH4VXBjvAawl62VRDesqFWGYTwqmAv21CEt2yoVQZZFt5vbv6lDpPzn66aW9NneK8aNl7mfrZgjZ/hLV36htFT5NU/vS0ThvdWOQ/Cq4Id4TWEvWyrILxlQ60yCOFVwV62oQhv2VCrDDItvMOVhPeC/wpvIZQzuw+S9Zo2kQG9zlztWxoWf/m1nHbxQDmjwyFy8C9/Wsxw2Kj75IgD95QtNvue/P2dD+S0iwbKmR0PlbNPPlwlY4RXBTvCawh72VZBeMuGWmUQwquCvWxDEd6yoVYZZFp4b9pfh8mFM6vNffeDT4rfw/vhJ/OK//6og34u/bp1kvr168miL76SPY84Ty6/6GTpcNR/9j3h7P7Fz+5+96fwfN+LT5GGlQ1UzoPwqmBHeA1hL9sqCG/ZUKsMQnhVsJdtKMJbNtQqgywL79c3/kqFSeOuT9U4d868BcXv423SuOFa9yq8+V2waLFstEFLadRQR3S/WxLhXWtc5XuAv7RWPtYakxBeDerlm4nwlo+1xiSEV4N6+WaaFt4blIT3opqFt3yp5DsJ4c2XZ6ZuCG8mfOaLEV7zEWVaEOHNhM98McJrPqJMC5oW3mEHZDqbb3Hji5/0LTVZh/AaigXhNRRGgFUQ3gBQDbVEeA2FEWAVhDcAVEMtLQvvV0MPVCHV5JIZKnNDDUV4Q5H16IvwekCLqAThjSgsj1URXg9oEZUgvBGF5bGqaeEdoiS83RBej6tEiQsBhNeFUrzPILzxZueyOcLrQineZxDeeLNz2dyy8H45uL3LEXJ/pmn36bn31GzIG15N+iWzEV5DYQRYBeENANVQS4TXUBgBVkF4A0A11NK08A46SIVU0x5PqMwNNRThDUXWoy/C6wEtohKEN6KwPFZFeD2gRVSC8EYUlseqpoX3eiXhvRTh9bhKlLgQQHhdKMX7DMIbb3YumyO8LpTifQbhjTc7l80tC+/igQe7HCH3Z5r1fDz3npoNecOrSb9kNsJrKIwAqyC8AaAaaonwGgojwCoIbwCohlqaFt7rDlEh1azXNJW5oYYivKHIevRFeD2gRVSC8EYUlseqCK8HtIhKEN6IwvJY1bTwXqskvJchvB5XiRIXAgivC6V4n0F4483OZXOE14VSvM8gvPFm57K5ZeH94ppDXY6Q+zPr9X4s956aDXnDq0m/ZDbCayiMAKsgvAGgGmqJ8BoKI8AqCG8AqIZamhbeAYepkFqvz6Mqc0MNRXhDkfXoi/B6QIuoBOGNKCyPVRFeD2gRlSC8EYXlsapp4b1aSXgvR3g9rhIlLgQQXhdK8T6D8MabncvmCK8LpXifQXjjzc5lc8vCu+iqw12OkPszzftOzb2nZkPe8OZAf+nSZbJg0ZfSasMWUqdOHe+OCK83uigKEd4oYvJeEuH1RhdFIcIbRUzeS5oW3iuP8D5XlsLmv30kS7m5WoQ3QySrVq2SW8c/Irfc8WCxy/otmsnN11wku7RtXWPXmbNfkwv73rTaz16bMVoqG9QXhDdDGBGUIrwRhJRhRYQ3A7wIShHeCELKsKJl4V3YX0d4W1yB8Ga4UmmVvv7WP+Sk8wfIhOG9Zaftt5Gbxj4gj838vTx171CpW3f1N71PzX5VLrtmtEwZ3b8aiC2+36r4ZhjhTet+lJ4G4U07X4Q37XwR3rTzNS28/Y5Ugd+i38Mqc0MN5Q1vBrJDRk6Wv/7zfRkzuEexy9zPFsp+x15UFNodtt1ytc4F4e0/ZJzMfmh4jVMR3gxhRFCK8EYQUoYVEd4M8CIoRXgjCCnDiqaF94qjMpzMv7RF/4f8iw1WIrwZQul+5a3SsnlT6dP15KouP9y3k4y49mLZZ49dahTern2Hy5Ht95LKygay2y5tpP2+7aReRUXxWYQ3QxgRlCK8EYSUYUWENwO8CEoR3ghCyrCiZeFd8Fsd4W15JcKb4UqlVdq5x2Bp03oL6dbl+KqDtTu4i/Tr3kkO3f9nqx32zb+9K9NnvSzNmzWRj+fMl8mPPCMdj96/SpgR3rTuR+lpEN6080V4084X4U07X9PC2/doFfgtr/rP309K5Q9veDMkWXjDW/iLar0vPMnpDW/pqAemPSd9r79d3pg5tviWd8OGgzJsQ6l1Aoct29T6iuyXgUC3I17NUE2pdQIfvsd/fq1nlGW/g1/vnqU8aO2Cy48J2n9NzVte/YDK3FBDEd4MZAuf4f37Ox/IqEH/+Q/K2j7DWzpq9ktvSpeeQ+TV6aOkYWUDhDdDFjGUIrwxpOS/I8Lrzy6GSoQ3hpT8d7QsvJ/30RHe9QcgvP43KrHK/35LQx/ZaYdt5MYxU2TazBervqVh3OQnpPBVZIVvcSj8mfjgTGnTenNpu91Wsmjxl9LjypFSv16F3D6sZ/HnvOFN7IKUHAfhTTtfhDftfBHetPM1Lby9f60Cf/1r7leZG2oob3gzkC18D+/NdzwoI8f/57vqGjdqKKMGdZNdd9y2+M+DRkySyVNnySuPjyz+89DbJsvYe6ZVTdy5bWsZ1LeLbLbJRghvhhxiKbbnxfUAACAASURBVEV4Y0nKb0+E149bLFUIbyxJ+e1pWXjnX3as36EyVm1w7ZSMHWyVI7w55LHk26Xy+YIv5HutNqjx+3f/d0Th2XnzF0qzJo2lRfOm1abzhjeHMAy3QHgNh5PDaghvDhANt0B4DYeTw2qmhbeXkvBeh/DmcLVoURMBhDfte4Hwpp0vwpt2vghv2vlaFt7Peh2nAn/D6+5TmRtqKG94Q5H16IvwekCLqAThjSgsj1URXg9oEZUgvBGF5bGqaeHt+d+vPvU4mnfJhgMne9daLER4DaWC8BoKI8AqCG8AqIZaIryGwgiwCsIbAKqhlpaFd96lOsK70fUIr6ErmtYqCG9aeZaeBuFNO1+EN+18Ed608zUtvD1OUIG/0aB7VeaGGsob3lBkPfoivB7QIipBeCMKy2NVhNcDWkQlCG9EYXmsall453bXEd5WgxFej6tEiQsBhNeFUrzPILzxZueyOcLrQineZxDeeLNz2dy08Hb7jcsRcn+m1ZBJuffUbMgbXk36JbMRXkNhBFgF4Q0A1VBLhNdQGAFWQXgDQDXU0rLwzrmkgwqpjYfeozI31FCENxRZj74Irwe0iEoQ3ojC8lgV4fWAFlEJwhtRWB6rWhbeTy/WEd7vDUN4Pa4SJS4EEF4XSvE+g/DGm53L5givC6V4n0F4483OZXPTwntRR5cj5P7M926YmHtPzYa84dWkXzIb4TUURoBVEN4AUA21RHgNhRFgFYQ3AFRDLS0L7yddT1QhtcmNd6vMDTUU4Q1F1qMvwusBLaIShDeisDxWRXg9oEVUgvBGFJbHqgjv6tAQXo+LRIkbAYTXjVOsTyG8sSbntjfC68Yp1qcQ3liTc9vbsvB+fOFJbofI+alNb7or54667XjDq8u/2nSE11AYAVZBeANANdQS4TUURoBVEN4AUA21NC28F5ysQmrT4RNU5oYaivCGIuvRF+H1gBZRCcIbUVgeqyK8HtAiKkF4IwrLY1XLwvvR+TrC+/2bEV6Pq0SJCwGE14VSvM8gvPFm57I5wutCKd5nEN54s3PZ3LTwnneKyxFyf+b7t4zPvadmQ97watIvmY3wGgojwCoIbwCohloivIbCCLAKwhsAqqGWloX3w3N1hHezEQivoSua1ioIb1p5lp4G4U07X4Q37XwR3rTztSy8/z7nVBX4m996p8rcUEN5wxuKrEdfhNcDWkQlCG9EYXmsivB6QIuoBOGNKCyPVU0Lb5dOHifKXrL5yHHZmxjqgPAaCgPhNRRGgFUQ3gBQDbVEeA2FEWAVhDcAVEMtLQvvB2frCO8WtyG8hq5oWqsgvGnlWXoahDftfBHetPNFeNPO17Twdj5NBf4Wo+5QmRtqKG94Q5H16IvwekCLqAThjSgsj1URXg9oEZUgvBGF5bGqZeF9v/PpHifKXrLlqNuzNzHUAeE1FAbCayiMAKsgvAGgGmqJ8BoKI8AqCG8AqIZamhbes5SEdzTCa+iKprUKwptWnqWnQXjTzhfhTTtfhDftfC0L73tnnqECf6sxY1XmhhrKG95QZD36Irwe0CIqQXgjCstjVYTXA1pEJQhvRGF5rGpaeM9QEt6xCK/HVaLEhQDC60Ip3mcQ3nizc9kc4XWhFO8zCG+82blsbll43z39TJcj5P7M1rePyb2nZkPe8GrSL5mN8BoKI8AqCG8AqIZaIryGwgiwCsIbAKqhlpaF91+nnaVCaps7RqvMDTUU4Q1F1qMvwusBLaIShDeisDxWRXg9oEVUgvBGFJbHqqaFt5OS8I5DeD2uEiUuBBBeF0rxPoPwxpudy+YIrwuleJ9BeOPNzmVzy8L7zqmdXY6Q+zOt7xyVe0/Nhrzh1aRfMhvhNRRGgFUQ3gBQDbVEeA2FEWAVhDcAVEMtTQvvKWerkGo9/jaVuaGGIryhyHr0RXg9oEVUgvBGFJbHqgivB7SIShDeiMLyWNWy8P7zZB3h/cEEhNfjKlHiQgDhdaEU7zMIb7zZuWyO8LpQivcZhDfe7Fw2ty28XVyOkPszP5gwMveemg15w6tJv2Q2wmsojACrILwBoBpqifAaCiPAKghvAKiGWloW3n+cdI4KqW3vulVlbqihCG8osh59EV4PaBGVILwRheWxKsLrAS2iEoQ3orA8VjUtvCcqCe/dCK/HVaLEhQDC60Ip3mcQ3nizc9kc4XWhFO8zCG+82blsbll43+54rssRcn9mu4kjcu+p2ZA3vJr0S2YjvIbCCLAKwhsAqqGWCK+hMAKsgvAGgGqopWXh/XsHHeFtcw/Ca+iKprUKwptWnqWnQXjTzhfhTTtfhDftfE0L72/OU4HfZtItKnNDDeUNbyiyHn0RXg9oEZUgvBGF5bEqwusBLaIShDeisDxWtSy8fzvhfI8TZS/Z/t6bszcx1AHhNRQGwmsojACrILwBoBpqifAaCiPAKghvAKiGWpoW3uOVhHcywmvoiqa1CsKbVp6lp0F4084X4U07X4Q37XwtC+9fj7tABf4O9w1XmRtqKG94Q5H16IvwekCLqAThjSgsj1URXg9oEZUgvBGF5bGqZeH9y3EXepwoe0nb+27K3sRQB4TXUBgIr6EwAqyC8AaAaqglwmsojACrILwBoBpqaVl4/3xsVxVSP5xyo8rcUEMR3lBkPfoivB7QIipBeCMKy2NVhNcDWkQlCG9EYXmsalp4f60kvPcjvB5XiRIXAgivC6V4n0F4483OZXOE14VSvM8gvPFm57K5ZeF965iLXI6Q+zM7PnBD7j01G/KGV5N+yWyE11AYAVZBeANANdQS4TUURoBVEN4AUA21tCy8bx59sQqpnR4cpjI31FCENxRZj74Irwe0iEoQ3ojC8lgV4fWAFlEJwhtRWB6rmhbeo5SE9yGE1+MqUeJCAOF1oRTvMwhvvNm5bI7wulCK9xmEN97sXDa3LLx/OuoSlyPk/szODw3NvadmQ97watIvmY3wGgojwCoIbwCohloivIbCCLAKwhsAqqGWloX3jSO7qZDa5eEhKnNDDUV4Q5H16IvwekCLqAThjSgsj1URXg9oEZUgvBGF5bGqZeH94xHdPU6UveRHjwzO3sRQB4TXUBgIr6EwAqyC8AaAaqglwmsojACrILwBoBpqaVp4D1cS3qkIr6ErGucqq1atkhUrV0q9iopqB0B448zTdWuE15VUnM8hvHHm5ro1wutKKs7nLAvv64f1UIG666ODapz72eeLpHGjhtK4UaXKXr5DecPrSy5D3dQZL8iw0ffJ0/dV/xuQCG8GqBGUIrwRhJRhRYQ3A7wIShHeCELKsKJl4X3t0EsznMy/9MePXV+t+IOP5kiXnkPl/Q/nFP/9MYfsLb+95FSpX6/6yzv/iWErEd6wfFe7LGd1HywffjJPNt6oJcJbRvYWRiG8FlIItwPCG46thc4Ir4UUwu1gWngPURLeadWFt3OPwdK0SSMZ0Oss+XTufDn+7P7y24tPkcMP3DNcMDl2RnhzhLm2VstXrJDC/xTw9O9elzETH0V41wYssZ8jvIkFWnIchDftfBHetPO1LLyvHtJTBf5Ppg2smrto8Vey5+HnyV0395Fdd9y2+O8H3DhBPp37uQwfoPOrj2sLBeGtLbEcnn/86Zdk0K2TEN4cWMbUAuGNKa3a74rw1p5ZTBUIb0xp1X5Xy8L7h4N71f5AOVTs9vh1VV3eee8jOaJTH5l1/w2y0QYtiv9+wpQZ8vD052XK6P45TAvfAuENz3i1CQivAnQDIxFeAyEEXAHhDQjXQGuE10AIAVewLLyvHHRZwJOvuXW7J66t+uHrb/1DTjp/gLww9RZp3qxJ8d9PnjpLRo5/eLWXdyrLOgxFeB0g5f0Iwps30Tj6Ibxx5OS7JcLrSy6OOoQ3jpx8tzQtvO2VhHf6f4X3uze8zz5wo2y4fnPe8PpetHWtDuFd1xL/z3kR3rRzR3jTzhfhTTtfy8L78oG9VeDvPuOaqrk1fYb3qmHjZe5nC/gMr0o6xocWvn93+fIV8sQzLxe/lmz6xEFSp26dqu/j5WvJjAeYcT2ENyNA4+UIr/GAMq6H8GYEaLzcsvC+dEAfFXo/fXJAtblndh8k6zVtIgN6ncm3NKgkEtHQf777kRx5WvWLW/g6j+t6dy6eAuGNKEyPVRFeD2gRlSC8EYXlsSrC6wEtohLLwvviAZerkPzZk1dXm/vuB58Uv4e38NWqhT9HHfRz6detk9SvX09lv9oO5TO8tSUW8HmENyBcA60RXgMhBFwB4Q0I10BrhNdACAFXMC28v1IS3qeqC+93+OfMW1D8Pt4mjRsGTCT/1ghv/ky9OyK83uiiKER4o4jJe0mE1xtdFIUIbxQxeS9pWXh/v39f73NlKdxj5lVZys3VIryGIkF4DYURYBWENwBUQy0RXkNhBFgF4Q0A1VBLy8L7wi9/q0Jqz6evVJkbaijCG4qsR1+E1wNaRCUIb0RheayK8HpAi6gE4Y0oLI9VTQvvfkrC+wzC63GVKHEhgPC6UIr3GYQ33uxcNkd4XSjF+wzCG292LptbFt7n973C5Qi5P7PXrDh+g5rrwXnD60qqDM8hvGWArDgC4VWEX4bRCG8ZICuOQHgV4ZdhtGXh/d0+/cpAYPURP39WZ26owyK8och69EV4PaBFVILwRhSWx6oIrwe0iEoQ3ojC8ljVsvDO3kfnTesvntV5s+wRn1MJwuuEqTwPIbzl4aw1BeHVIl+euQhveThrTUF4tciXZ65p4d1bSXifQ3jLc/vWwSkIb9qhI7xp54vwpp0vwpt2vpaF97lf6Pzlsb1n6/xluVA3jTe8och69EV4PaBFVILwRhSWx6oIrwe0iEoQ3ojC8ljVsvA++3Od78Pd53c63//rEZ9TCcLrhKk8DyG85eGsNQXh1SJfnrkIb3k4a01BeLXIl2euaeHdS0l4n0d4y3P71sEpCG/aoSO8aeeL8KadL8Kbdr6WhXfWnjX/it/Qiez7gs6vNA51Lt7whiLr0Rfh9YAWUQnCG1FYHqsivB7QIipBeCMKy2NVy8L7zJ4DPE6UvWS/F/pkb2KoA8JrKAyE11AYAVZBeANANdQS4TUURoBVEN4AUA21tCy8T+9xjQqpX/6+t8rcUEMR3lBkPfoivB7QIipBeCMKy2NVhNcDWkQlCG9EYXmsalp4f6YkvC8ivB5XiRIXAgivC6V4n0F4483OZXOE14VSvM8gvPFm57K5ZeGd+dNrXY6Q+zP7v3RZ7j01G/KGV5N+yWyE11AYAVZBeANANdQS4TUURoBVEN4AUA21tCy8T+1+nQqpX73cS2VuqKEIbyiyHn0RXg9oEZUgvBGF5bEqwusBLaIShDeisDxWtSy8T7bTEd4DXkF4Pa4SJS4EEF4XSvE+g/DGm53L5givC6V4n0F4483OZXPLwjtjt4EuR8j9mQP/0DP3npoNecOrSb9kNsJrKIwAqyC8AaAaaonwGgojwCoIbwCohlpaFt7pSsLbHuE1dEMTWwXhTSzQkuMgvGnni/CmnS/Cm3a+poX3J9erwG//6qUqc0MN5Q1vKLIefRFeD2gRlSC8EYXlsSrC6wEtohKEN6KwPFa1LLxP/HiQx4mylxz0Wo/sTQx1QHgNhYHwGgojwCoIbwCohloivIbCCLAKwhsAqqGWloX38V11hPfg1xFeQ1c0rVUQ3rTyLD0Nwpt2vghv2vkivGnna1l4p/1osAr8Q/7YXWVuqKG84Q1F1qMvwusBLaIShDeisDxWRXg9oEVUgvBGFJbHqpaF97FdhnicKHvJoW90y97EUAeE11AYCK+hMAKsgvAGgGqoJcJrKIwAqyC8AaAaamlZeB9VEt7DEF5DNzSxVRDexAItOQ7Cm3a+CG/a+SK8aedrWXin7jxUBf7hf7pEZW6oobzhDUXWoy/C6wEtohKEN6KwPFZFeD2gRVSC8EYUlseqloX3kZ2GeZwoe8kRb16cvYmhDgivoTAQXkNhBFgF4Q0A1VBLhNdQGAFWQXgDQDXU0rTw7qgkvG8hvIauaFqrILxp5Vl6GoQ37XwR3rTzRXjTztey8D78wxtU4B/554tU5oYayhveUGQ9+iK8HtAiKkF4IwrLY1WE1wNaRCUIb0RheaxqWXgfaqsjvEf9BeH1uEqUuBBAeF0oxfsMwhtvdi6bI7wulOJ9BuGNNzuXzS0L74M73OhyhNyfOfqvXXPvqdmQN7ya9EtmI7yGwgiwCsIbAKqhlgivoTACrILwBoBqqKVl4X1gh5tUSB3z1wtV5oYaivCGIuvRF+H1gBZRCcIbUVgeqyK8HtAiKkF4IwrLY1XLwnv/9jrC++u/IbweV4kSFwIIrwuleJ9BeOPNzmVzhNeFUrzPILzxZueyuWXhndJmuMsRcn/m2L9fkHtPzYa84dWkXzIb4TUURoBVEN4AUA21RHgNhRFgFYQ3AFRDLS0L733b3axC6ri3z1eZG2oowhuKrEdfhNcDWkQlCG9EYXmsivB6QIuoBOGNKCyPVU0L77ZKwvsPhNfjKlHiQgDhdaEU7zMIb7zZuWyO8LpQivcZhDfe7Fw2tyy8k39wi8sRcn/m+H+el3tPzYa84dWkXzIb4TUURoBVEN4AUA21RHgNhRFgFYQ3AFRDLS0L772tR6iQOuGdc1XmhhqK8IYi69EX4fWAFlEJwhtRWB6rIrwe0CIqQXgjCstjVcvCO0lJeH+D8HrcJEqcCCC8TpiifQjhjTY6p8URXidM0T6E8EYbndPiloX3nm1udTpD3g91+Nc5ebdU7ccbXlX81YcjvIbCCLAKwhsAqqGWCK+hMAKsgvAGgGqopWXhnbi1jvB2fBfhNXRF01oF4U0rz9LTILxp54vwpp0vwpt2vpaF9+6tRqrAP/G9LipzQw3lDW8osh59EV4PaBGVILwRheWxKsLrAS2iEoQ3orA8VjUtvFve5nGi7CUnvn929iaGOiC8hsJAeA2FEWAVhDcAVEMtEV5DYQRYBeENANVQS8vCe9cWOsJ70gcIr6ErmtYqCG9aeZaeBuFNO1+EN+18Ed6087UsvBO2GKUC/+QPOqvMDTWUN7yhyHr0RXg9oEVUgvBGFJbHqgivB7SIShDeiMLyWNWy8I7ffLTHibKXnPLvs7I3MdQB4TUUBsJrKIwAqyC8AaAaaonwGgojwCoIbwCohlpaFt47N9MR3lM/RHgNXdG0VkF408qz9DQIb9r5Irxp54vwpp2vZeEd9/0xKvA7fXSmytxQQ3nDG4qsR1+E1wNaRCUIb0RheayK8HpAi6gE4Y0oLI9VLQvvHZvqCO9pHyO8Hlcp/ZJly1dI/XoVmQ6K8GbCZ74Y4TUfUaYFEd5M+MwXI7zmI8q0oGXhvX2TsZnO5lt8+idn+JaarOMNbw6xfPDRXDn4xEvlyUmDZdPvbbjGjjNnvyYX9r1ptZ+/NmO0VDaoLwhvDmEYboHwGg4nh9UQ3hwgGm6B8BoOJ4fVTAvv927P4YS1b3H6p6fXvshwBcKbMZwO514lf/rLO8UuaxPep2a/KpddM1qmjO5fbeoW328lderUQXgzZmG9HOG1nlC2/RDebPysVyO81hPKtp9l4R2rJLxnILzZLlVq1XM/Wyifzp0vBfF1Ed7+Q8bJ7IeG14iBN7yp3Y7q50F4084X4U07X4Q37XwtC++Yje9QgX/mnNNU5oYayhveHMjOmbdAfnncxU7C27XvcDmy/V5SWdlAdtuljbTft53Uq/jPZ38R3hzCMNwC4TUcTg6rIbw5QDTcAuE1HE4Oq1kW3lGtdIS381yEN4erlVYLV+F982/vyvRZL0vzZk3k4znzZfIjz0jHo/eXPl1PRnjTuhI1ngbhTTtkhDftfBHetPO1LLy3KQnv2Qhv2pfe53Suwlva+4Fpz0nf62+XN2aOLb7lfbl9H5/x1ERC4NkXt49kU9b0IXBu37t9yqiJhMCKJQ0i2ZQ1fQisd/kjPmVlqRm50biyzCkd0mVeJ5W5oYbykYYcyPoK7+yX3pQuPYfIq9NHScPKBghvDllYboHwWk4n+24Ib3aGljsgvJbTyb6bZeG9VUl4z0F4s1+slDoUvn+38JfWDup4qUy7a2Dxa8m++z7ecZOfkMJXkU0Y3rt45IkPzpQ2rTeXttttJYsWfyk9rhxZfPb2YT2LP+cNb0o3Y/WzILxp54vwpp0vwpt2vpaFd8SGOm94z/2MN7xp3/panq7dwV3k62+WVFWt36JZ1bcwDBoxSSZPnSWvPD6y+POht02WsfdMq3p257atZVDfLrLZJhshvLXkHuPjCG+MqbnvjPC6s4rxSYQ3xtTcd7YsvLcoCe95CK/7BeLJ1Qks+XapzJu/UJo1aSwtmjet9gBveNO+MQhv2vkivGnni/Cmna9l4b15gztV4J8//1SVuaGG8hneUGQ9+iK8HtAiKkF4IwrLY1WE1wNaRCUIb0RheaxqWXiHr68jvBd8jvB6XCVKXAggvC6U4n0G4Y03O5fNEV4XSvE+g/DGm53L5paF9yYl4b0Q4XW5OjzjQwDh9aEWTw3CG09WPpsivD7U4qlBeOPJymdTy8J7Q0udN7wXLeANr89dosaBAMLrACniRxDeiMNzWB3hdYAU8SMIb8ThOaxuWXiHtRzvcIL8H7l4wSn5N1XsyGd4FeGXjkZ4DYURYBWENwBUQy0RXkNhBFgF4Q0A1VBLy8I7tIWO8F6yEOE1dEXTWgXhTSvP0tMgvGnni/CmnS/Cm3a+loV3SHMd4e22COFN+9Yrng7hVYRfhtEIbxkgK45AeBXhl2E0wlsGyIojLAvvYCXh7Y7wKt7IxEcjvGkHjPCmnS/Cm3a+CG/a+VoW3kHrTVCB3+OLk1XmhhrKZ3hDkfXoi/B6QIuoBOGNKCyPVRFeD2gRlSC8EYXlsapl4b1eSXgvRXg9bhIlTgQQXidM0T6E8EYbndPiCK8TpmgfQnijjc5pccvCO7CZzhvenot5w+t0eXio9gQQ3tozi6kC4Y0prdrvivDWnllMFQhvTGnVflfLwnudkvD2Qnhrf5GocCOA8LpxivUphDfW5Nz2RnjdOMX6FMIba3Jue1sW3mua3uV2iJyf6v3lSTl31G3HZ3h1+VebjvAaCiPAKghvAKiGWiK8hsIIsArCGwCqoZaWhXdAEx3h7fMVwmvoiqa1CsKbVp6lp0F4084X4U07X4Q37XwtC+/VSsJ7OcKb9qXXPB3Cq0k//GyENzxjzQkIryb98LMR3vCMNSdYFt6rGt+tgqbv1yd6z125cpXMnb9ANly/udSrqPDuk2chH2nIk2bGXghvRoDGyxFe4wFlXA/hzQjQeDnCazygjOtZFt4rlYT3t57C++zv35DuV94qX3+zpJjKFZecKscfsd8aEzri1N7yzvsfV/v5eZ2OknM7HZUx1erlCG+uOLM1Q3iz8bNejfBaTyjbfghvNn7WqxFe6wll28+y8PZvpPOG94pvav+G95slS2Xvoy+U808/Wk485lcy64U/Ste+w2X6PYNks002qjGkgvAe+qs95KD9dq/6efNmTaRF86bZQi2pRnhzxZmtGcKbjZ/1aoTXekLZ9kN4s/GzXo3wWk8o236WhbdfQx3h7bek9sJbeLt77mXD5PUZo6VBg/rFUA45qWdRfk885oA1Cm+nEw6SYw7ZO1uIa6lGeIPirV1zhLd2vGJ7GuGNLbHa7Yvw1o5XbE8jvLElVrt9LQvvFQ0n1u4wOT3df0nHWneaPHWWjLv3cZl218Cq2gv63Chbbb6JdOty/BqFt0mTRtJ6y01l0403kMMO2EO2+P7GtZ69tgKEd22EyvhzhLeMsBVGIbwK0Ms4EuEtI2yFUQivAvQyjrQsvL+t1BHeK7+tLrxTZ7wgn877vMZU2m63lezVbkcZM/ExeeKZl2XK6P5VzxU+z9u0cSPp171TjbW33PGg1K2oK6tWiTz9u9fk/Q/nyP1j+ucuvQhvGf8DtbZRCO/aCMX9c4Q37vzWtj3CuzZCcf8c4Y07v7Vtb1l4L1cS3qtLhPfuB56SDz+ZVyPKH++0rRyw927i84b3fxsuW7Zc2nfsISf/+kA57TcHry22Wv0c4a0VrrAPI7xh+Wp3R3i1Ewg7H+ENy1e7O8KrnUDY+ZaFt08DnTe8A5bW/iMN332G949PjpH69esVQ2vfoYecctyBa/wMb2myJ5zdX/bZ80dy7qlH5ho6wpsrzmzNEN5s/KxXI7zWE8q2H8KbjZ/1aoTXekLZ9rMsvL0b3JPtcJ7V1yztUOvKr7/5VtodfLb0PK+DdKzhWxpe+ePfZOAt98iQK86VLTfbWD74aI48/fzrxW9o2KBlc5n+zMvSc8BtMv6m3vKTnber9fz/qwDhzRVntmYIbzZ+1qsRXusJZdsP4c3Gz3o1wms9oWz7WRbey+rrCO+1y2ovvIUUCgJb+Itq3/25/KKTpcNR+xf/8ZkXXpfze98oD4y9Stq03rwovJ0uuk7mzFtQ9XxBlk85rn22QGuoRnhzR+rfEOH1ZxdDJcIbQ0r+OyK8/uxiqER4Y0jJf0fLwturno7wXrfcT3gLKaxYsbL4F9xabdCi6qMNa0pn1apV8vnCxcVfVLHJxhsE+81sCK//fz5yr0R4c0dqqiHCayqO3JdBeHNHaqohwmsqjtyXsSy8PZWEd2AG4c09oBwaIrw5QMyrBcKbF0mbfRBem7nktRXCmxdJm30QXpu55LWVZeG9tGJSXsesVZ/rV/ymVs9bfxjhNZQQwmsojACrILwBoBpqifAaCiPAKghvAKiGWloW3h5KwjsI4TV0QxNbBeFNLNCS4yC8aeeL8KadL8Kbdr6Whbd7XZ03vINX8oY37VuveDqEVxF+GUYjvGWArDgC4VWEX4bRCG8ZICuOsCy8lygJ71CEV/FGJj4a4U07YIQ37XwR3rTzRXjTztey8F6sJLzDEN60L73m6RBeTfrhZyO84RlrTkB4NemHn43whmesOcGy8F6kvFyF2AAAGRxJREFUJLw3ILyaVzLt2Qhv2vkivGnni/CmnS/Cm3a+loW3q5Lw3ojwpn3pNU+H8GrSDz8b4Q3PWHMCwqtJP/xshDc8Y80JloX3QiXhvQnh1bySac9GeNPOF+FNO1+EN+18Ed6087UsvBfU1flNa8NX+v+mNYu3he/hNZQKwmsojACrILwBoBpqifAaCiPAKghvAKiGWloW3vOUhPcWhNfQDU1sFYQ3sUBLjoPwpp0vwpt2vghv2vlaFt5zlYR3BMKb9qXXPB3Cq0k//GyENzxjzQkIryb98LMR3vCMNSdYFt5zlIT3VoRX80qmPRvhTTtfhDftfBHetPNFeNPO17Lwdqk7UQX+yJUdVeaGGspneEOR9eiL8HpAi6gE4Y0oLI9VEV4PaBGVILwRheWxqmXhPVtJeG9DeD1uEiVOBBBeJ0zRPoTwRhud0+IIrxOmaB9CeKONzmlxy8LbWUl4RyG8TneHhzwIILwe0CIqQXgjCstjVYTXA1pEJQhvRGF5rGpZeM9UEt4xCK/HTaLEiQDC64Qp2ocQ3mijc1oc4XXCFO1DCG+00Tktbll4z6h7t9MZ8n5o7MoT826p2o/P8Krirz4c4TUURoBVEN4AUA21RHgNhRFgFYQ3AFRDLS0L7+lKwns7wmvohia2CsKbWKAlx0F4084X4U07X4Q37XwtC+9pSsJ7B8Kb9qXXPB3Cq0k//GyENzxjzQkIryb98LMR3vCMNSdYFt5OSsI7DuHVvJJpz0Z4084X4U07X4Q37XwR3rTztSy8p9a9SwX+nStPUpkbaiif4Q1F1qMvwusBLaIShDeisDxWRXg9oEVUgvBGFJbHqpaF92Ql4Z2A8HrcJEqcCCC8TpiifQjhjTY6p8URXidM0T6E8EYbndPiloX3pAqdN7x3reANr9Pl4aGaCSxbvkI+m79Q1m+5nlQ2qF/tIYQ37VuD8KadL8Kbdr4Ib9r5WhbeE5WE926EN+1LH/J0o+9+VG4YPaVqRPt928kVl3SS5us1Kf47hDckff3eCK9+BiE3QHhD0tXvjfDqZxByA8vC27FiQsijr7H3xBUnq8wNNZTP8IYiW0Pf+x6dJZtv2kp2afsD+ffHc+WMSwbKGR0OlU4nHITwljEHrVEIrxb58sxFeMvDWWsKwqtFvjxzLQtvByXhvQfhLc/lWxem9L3+dvnok3ly+7CeCO86EDjCm3bICG/a+SK8aedrWXh/oyS8kxDetC99uU5X+Cxv+w7d5dD995BuXY5HeMsFXnEOwqsIvwyjEd4yQFYcgfAqwi/DaMvCe7yS8E5GeMtw89aBEVcMvkOmzXxJHptwnbTasAXCuw5kjvCmHTLCm3a+CG/a+VoW3uOUhPc+hDftS1+O040Y95DcMu4hmTTyCtlp+62rRvKX1spBX28GwqvHvhyTEd5yUNabgfDqsS/HZMvCe2zF+HIgWG3GlBWnqMwNNZS/tBaKbA19V65cJUNG3iuTp86SO2/sJW2326raUwhvGcNQGIXwKkAv40iEt4ywFUYhvArQyzjSsvD+Wkl470d4y3gDExt1+cCx8uDjs2XkwG6yzZabVJ1u441aSr2KCr6WLLG8S4+D8KYdMMKbdr4Ib9r5WhbeY5SE9wGEN+1LH/J07Tv0kA8/mbfaiGl3DZQtN9sY4Q0J30BvhNdACAFXQHgDwjXQGuE1EELAFSwL79FKwvsgwhvwxq3jrflIQ9oXAOFNO1+EN+18Ed6087UsvEdW3KkC/+EVp6rMDTWUz/CGIuvRF+H1gBZRCcIbUVgeqyK8HtAiKkF4IwrLY1XLwnuEkvA+gvB63CRKnAggvE6Yon0I4Y02OqfFEV4nTNE+hPBGG53T4paF93Al4Z2K8DrdHR7yIIDwekCLqAThjSgsj1URXg9oEZUgvBGF5bGqZeE9TEl4H0V4PW4SJU4EEF4nTNE+hPBGG53T4givE6ZoH0J4o43OaXHLwntoxTinM+T90GMrOuXdUrUfn+FVxV99OMJrKIwAqyC8AaAaaonwGgojwCoIbwCohlpaFt5DlIR3GsJr6IYmtgrCm1igJcdBeNPOF+FNO1+EN+18LQvvwUrC+zjCm/al1zwdwqtJP/xshDc8Y80JCK8m/fCzEd7wjDUnWBbe9krCOx3h1bySac9GeNPOF+FNO1+EN+18Ed6087UsvAfWu0MF/ozlp6nMDTWUz/CGIuvRF+H1gBZRCcIbUVgeqyK8HtAiKkF4IwrLY1XLwnuAkvA+ifB63CRKnAggvE6Yon0I4Y02OqfFEV4nTNE+hPBGG53T4paFd/96tzudIe+HZi4/Pe+Wqv14w6uKv/pwhNdQGAFWQXgDQDXUEuE1FEaAVRDeAFANtbQsvL+sN1aF1NPLz1CZG2oowhuKrEdfhNcDWkQlCG9EYXmsivB6QIuoBOGNKCyPVS0L775KwjsL4fW4SZQ4EUB4nTBF+xDCG210TosjvE6Yon0I4Y02OqfFLQvvPvXGOJ0h74eeXX5m3i1V+/GGVxV/9eEIr6EwAqyC8AaAaqglwmsojACrILwBoBpqaVl4f6EkvLMRXkM3NLFVEN7EAi05DsKbdr4Ib9r5Irxp52tZeH9eb7QK/N8tP0tlbqihvOENRdajL8LrAS2iEoQ3orA8VkV4PaBFVILwRhSWx6qWhXfP+qM8TpS95IVlnbM3MdQB4TUUBsJrKIwAqyC8AaAaaonwGgojwCoIbwCohlpaFt49lIT39wivoRua2CoIb2KBlhwH4U07X4Q37XwR3rTztSy8P61/mwr8l5adrTI31FDe8IYi69EX4fWAFlEJwhtRWB6rIrwe0CIqQXgjCstjVcvC267+SI8TZS95ZVmX7E0MdUB4DYWB8BoKI8AqCG8AqIZaIryGwgiwCsIbAKqhlpaFdzcl4f0Dwmvohia2CsKbWKAlx0F4084X4U07X4Q37XwtC++P69+qAv+1ZeeozA01lDe8och69EV4PaBFVILwRhSWx6oIrwe0iEoQ3ojC8ljVsvDuWn+Ex4myl7y+7NzsTQx1QHgNhYHwGgojwCoIbwCohloivIbCCLAKwhsAqqGWloV3lwY6wvvGUoTX0BVNaxWEN608S0+D8KadL8Kbdr4Ib9r5WhbenRvcogL/T0vPU5kbaihveEOR9eiL8HpAi6gE4Y0oLI9VEV4PaBGVILwRheWxqmXh3VFJeN9CeD1uEiVOBBBeJ0zRPoTwRhud0+IIrxOmaB9CeKONzmlxy8L7wwY3O50h74f+vPT8vFuq9uMNryr+6sMRXkNhBFgF4Q0A1VBLhNdQGAFWQXgDQDXU0rLw7tBguAqpvy69QGVuqKEIbyiyHn0RXg9oEZUgvBGF5bEqwusBLaIShDeisDxWtSy82ysJ798QXo+bRIkTAYTXCVO0DyG80UbntDjC64Qp2ocQ3mijc1rcsvBu1+AmpzPk/dDbSy/Mu6VqP97wquKvPhzhNRRGgFUQ3gBQDbVEeA2FEWAVhDcAVEMtLQvvtg1uVCH1j6VdVeaGGorwhiLr0Rfh9YAWUQnCG1FYHqsivB7QIipBeCMKy2NVy8LbulJHeN/5FuH1uEqUuBBAeF0oxfsMwhtvdi6bI7wulOJ9BuGNNzuXzS0L7zaVN7gcIfdn/vXtRbn31GzIG15N+iWzEV5DYQRYBeENANVQS4TXUBgBVkF4A0A11NKy8G5VOUyF1HvfXqwyN9RQhDcUWY++CK8HtIhKEN6IwvJYFeH1gBZRCcIbUVgeq1oW3i2VhPd9hNfjJlHiRADhdcIU7UMIb7TROS2O8DphivYhhDfa6JwWtyy8m1cOdTpD3g/9+9tL8m6p2o83vKr4qw9HeA2FEWAVhDcAVEMtEV5DYQRYBeENANVQS8vCu5mS8H6I8Bq6oYmtgvAmFmjJcRDetPNFeNPOF+FNO1/Lwrtp5RAV+B9/201lbqihvOENRdajL8LrAS2iEoQ3orA8VkV4PaBFVILwRhSWx6qWhXeThoM9TpS95JMl3bM3MdQB4TUUBsJrKIwAqyC8AaAaaonwGgojwCoIbwCohlpaFt6NlYR3DsJr6IYmtgrCm1igJcdBeNPOF+FNO1+EN+18LQtvq4aDVODPXdIj09zlK1ZI3Tp1pW7dOpn65FXMG968SObQB+HNAaLhFgiv4XByWA3hzQGi4RYIr+FwcljNsvBu2PD6HE5Y+xafLbm09kX/v+KbJUvlhLP7SeeTDpfDDtjDu0+ehQhvnjQz9kJ4MwI0Xo7wGg8o43oIb0aAxssRXuMBZVzPsvBuoCS88z2Fd/DIe+WOSY8XExnY52yEN+PdTLIc4U0y1qpDIbxp54vwpp0vwpt2vpaFt2XDgSrwFyzp6TV34aIvZcnSpdLx3Kvkks7HI7xeFBMvQnjTDhjhTTtfhDftfBHetPO1LLwtGl6nAn/hkl6Z5rbv0EMuOP0YhDcTxUSLEd5Eg/3/x0J4084X4U07X4Q37XwtC+96SsL7RYnwTp3xgnw67/MaL0Lb7baSvdrtWO1nCG/a/5nJdDqENxM+88UIr/mIMi2I8GbCZ74Y4TUfUaYFLQtvpoPlWHz3A0/Jh5/Mq7Hjj3faVg7YezeEN0feSbdCeJOOVxDetPNFeNPOF+FNO1+EN/98ecObP9NkOiK8yURZ40EQ3rTzRXjTzhfhTTtfhDe/fAvfv7tq5So57JTLpMspR8hhv9pD6tevl98Az058LZknuBBlCG8IqnZ6Irx2sgixCcIbgqqdngivnSxCbILw5kf1kn4jZPqsl6s1fHT8tbL1FpvkN8SjE8LrAe1/S1auXCWfL/yi+N9emjdrkqkbwpsJn/lihNd8RJkWRHgz4TNfjPCajyjTgghvJnxRFCO8GWL6/R/+LBf2HS5ff7Ok2KXdj7aX7uecIDu22brGrjNnvyYX9r1ptZ+9NmO0VDaoLwhvhjAiKEV4Iwgpw4oIbwZ4EZQivBGElGFFhDcDvEhKEd4MQb342l9k3mcLZe89dpElS5bKlcPulMIb31uvu7jGrk/NflUuu2a0TBndv9rPt/h+K6lTpw7CmyGLGEoR3hhS8t8R4fVnF0MlwhtDSv47Irz+7GKpRHhzTKrwHXW9rhklb8wcK/UqKlbrXBDe/kPGyeyHhtc4lTe8OYZhsBXCazCUHFdCeHOEabAVwmswlBxXQnhzhGm0FcKbYzAF2f3nux+t9gb3uxEF4e3ad7gc2X4vqaxsILvt0kba79uuSo4R3hzDMNgK4TUYSo4rIbw5wjTYCuE1GEqOKyG8OcI02grhzSmY797ujhncQ/bY7Yc1dn3zb+8W/+Zi4S+3fTxnvkx+5BnpePT+0qfrycXnEd6cwjDaBuE1GkxOayG8OYE02gbhNRpMTmshvDmBNNwG4c0hnOdfeUs69xgsV1xyqhx/xH7OHR+Y9pz0vf72NX4EwrkRD0IAAhCAAAQgAAEIrJEAwpvxchTe2Ba+c+7qnmfI0Qf/olbdZr/0pnTpOURenT5KGlY2qFUtD0MAAhCAAAQgAAEIuBFAeN041fjUw9Ofl97XjpZe53eUX/78x1XPtGzeVBo3aijjJj8hha8imzC8d/FnEx+cKW1aby5tt9tKFi3+UnpcOVLq16uQ24f1zLAFpRCAAAQgAAEIQAAC/xcBhDfD/bhy2Hi59+GnV+vw3dveQSMmyeSps+SVx0cWnxl622QZe8+0qud3bttaBvXtIpttslGGLeIsXfzl11L49YMtmzeL8wBs/X8SWLVqlaxYubLGbysBXdwEli1fIZ/NXyjrt1yv+P3h/EmLQOH/Ln/2+aLir4ZttWFLqaiom9YBOc06SwDhLXP0S75dKvPmL5RmTRpLi+ZNyzxdf1zhl3T0vPo2efr514vLFKR/+NUXyobrN9dfjg1yI1D4S5zDRt8nT983LLeeNNInMPruR+WG0VOqFil8y8wVl3SS5utl+y2T+idjgwKBwgucwouc7/5svFFLuenqC9f4y5SgBoGYCCC8MaWVwK5jJj4m902dJROG95FGDRvIOb2GFX+/9lWXnp7A6TjCBx/NkbO6D5YPP5knhf9nifCmdSfue3SWbL5pK9ml7Q/k3x/PlTMuGShndDhUOp1wUFoHXUdPU/gvqoUXMT/ZuU3xf4Hr3n+ELF++go/draP3IbVjI7ypJWr8PMeedUXxu4fPOvGw4qbf/aW/t565o/jb5vgTN4Hv/ufQp3/3uoyZ+CjCG3eca92+8C0zH30yDyFaK6k4H+h+5a3F3x46tN+5cR6ArSHwPwQQXq5DWQm0O7hL8RstCtJb+POXt9+T4zr3kxem3lL8fmL+pEHg8adfkkG3TkJ404izxlMUPsvbvkN3OXT/PaRbl+MTPum6d7RHZjwvhf/S+va//i1D+50n2/9gi3UPAidOjgDCm1ykdg9U+ItMO+53moy49mLZZ49diou+895HckSnPvLUvUNkk403sLs8m9WKAMJbK1xRPnzF4Dtk2syX5LEJ10mrDVtEeQaWrplA4XPar/7pbZn72QK56tIzZPddtwcVBKIngPBGH2FcByi84R3Q60w5cJ/deMMbV3S12hbhrRWu6B4eMe4huWXcQzJp5BWy0/ZbR7c/C7sRuG3CVLnr/hky+6HhbgU8BQHDBBBew+GkuFrhM7wH7be7nNnx0OLx+AxviimLILxp5lr4POeQkfcWv27xzht7Fb9TnD/pEpjx7B/k4itu5reBphvxOnUyhHedilv/sIWvNZry6LPFb2lo3KhSuvQcyrc06MeS2waFj60U/lb3E8+8XPxasukTB0mdunX4Pt7cCOs2unzgWHnw8dkycmA32WbLTaqWKXwjR72KCt3lmJ6ZQOHN/V6771T8BUnzF3whhb+01qiyAX8pMTNZGlgggPBaSGEd2uGrr5cU/4/ocy++UTz1jm22luEDuvIZwETuwD/f/UiOPK1PtdMcfuCecl3vzomccN0+RvsOPYpfOVf6Z9pdA2XLzTZet+EkcPo+142Rh574XdVJdt1xW7muT+d18pcjJRAnRyghgPByJVQILFr8lSxbtpxfOKFCn6EQgAAEaiawdOkymTt/oTRt3Gid/OVI3It0CSC86WbLySAAAQhAAAIQgAAERATh5RpAAAIQgAAEIAABCCRNAOFNOl4OBwEIQAACEIAABCCA8HIHIAABCEAAAhCAAASSJoDwJh0vh4MABCAAAQhAAAIQQHi5AxCAAAQgAAEIQAACSRNAeJOOl8NBAAIQgAAEIAABCCC83AEIQAACEIAABCAAgaQJILxJx8vhIAABCEAAAhCAAAQQXu4ABCAAAQhAAAIQgEDSBBDepOPlcBCAAAQgAAEIQAACCC93AAIQgAAEIAABCEAgaQIIb9LxcjgIQAACEIAABCAAAYSXOwABCEAAAhCAAAQgkDQBhDfpeDkcBCAAAQhAAAIQgADCyx2AAAQgAAEIQAACEEiaAMKbdLwcDgIQgAAEIAABCEAA4eUOQAACEIAABCAAAQgkTQDhTTpeDgcBCEAAAhCAAAQggPByByAAAQhAAAIQgAAEkiaA8CYdL4eDAAQgAAEIQAACEEB4uQMQgAAEIAABCEAAAkkTQHiTjpfDQQACEIAABCAAAQggvNwBCEAAAhCAAAQgAIGkCSC8ScfL4SAAAQhAAAIQgAAEEF7uAAQgAAEIQAACEIBA0gQQ3qTj5XAQgAAEIAABCEAAAggvdwACEIAABCAAAQhAIGkCCG/S8XI4CEAAAhCAAAQgAAGElzsAAQhAAAIQgAAEIJA0AYQ36Xg5HAQgAAEIQAACEIAAwssdgAAEIAABCEAAAhBImgDCm3S8HA4CEIAABCAAAQhAAOHlDkAAAhCAAAQgAAEIJE0A4U06Xg4HAQhAAAIQgAAEIIDwcgcgAAEIQAACEIAABJImgPAmHS+HgwAEIAABCEAAAhBAeLkDEIAABCAAAQhAAAJJE0B4k46Xw0EAAhCAAAQgAAEIILzcAQhAAAIQgAAEIACBpAkgvEnHy+EgAAEIQAACEIAABBBe7gAEIAABCEAAAhCAQNIEEN6k4+VwEIAABCAAAQhAAAIIL3cAAhCAAAQgAAEIQCBpAghv0vFyOAhAAAIQgAAEIAABhJc7AAEIQAACEIAABCCQNAGEN+l4ORwEIAABCEAAAhCAAMLLHYAABCAAAQhAAAIQSJoAwpt0vBwOAhCAAAQgAAEIQADh5Q5AAAIQgAAEIAABCCRNAOFNOl4OBwEIQAACEIAABCCA8HIHIAABCEAAAhCAAASSJoDwJh0vh4MABCAAAQhAAAIQQHi5AxCAAAQgAAEIQAACSRNAeJOOl8NBAAIQgAAEIAABCCC83AEIQAACEIAABCAAgaQJILxJx8vhIAABCEAAAhCAAAQQXu4ABCAAAQhAAAIQgEDSBBDepOPlcBCAAAQgAAEIQAACCC93AAIQgAAEIAABCEAgaQIIb9LxcjgIQAACEIAABCAAAYSXOwABCEAAAhCAAAQgkDQBhDfpeDkcBCAAAQhAAAIQgADCyx2AAAQgAAEIQAACEEiaAMKbdLwcDgIQgAAEIAABCEAA4eUOQAACEIAABCAAAQgkTeD/AVEJOdj3ITvlAAAAAElFTkSuQmCC"
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "sigma = 0.5\n",
    "rng = np.random.default_rng(42)\n",
    "data = rng.normal(loc=0, scale=sigma, size=(M, N, 1))\n",
    "data_model = noiseModel.Isotropic.Sigmas([sigma])\n",
    "print(f\"data_model = {data_model}\")\n",
    "px.imshow(data[:,:,0], zmin=-1, zmax=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For smoothness, we use a different noise model:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "smoothness_model = isotropic dim=1 sigma=0.5\n",
      "\n"
     ]
    }
   ],
   "source": [
    "smoothness_sigma = 0.5\n",
    "smoothness_model = noiseModel.Isotropic.Sigmas([smoothness_sigma])\n",
    "print(f\"smoothness_model = {smoothness_model}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then we create the factor graph:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "I = np.eye(1, 1, dtype=float)\n",
    "zero = np.zeros((1,1))\n",
    "graph = gtsam.GaussianFactorGraph()\n",
    "for row in range(M):\n",
    "    for col in range(N):\n",
    "        # add data terms:\n",
    "        j = keys[(row,col)]\n",
    "        graph.add(j, I, np.array(data[row,col]), data_model)\n",
    "        # add smoothness terms:\n",
    "        if col>0:\n",
    "            j1 = keys[(row,col-1)]\n",
    "            graph.add(j, I, j1, -I, zero, smoothness_model)\n",
    "        if row>0:\n",
    "            j2 = keys[(row-1,col)]\n",
    "            graph.add(j, I, j2, -I, zero, smoothness_model)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can display the graph, although graphviz (our graph plotting engine) does not respect the \"image\" nature. Note that we omit the \"dot\" for binary factors, which both reveals the undirected MRF and is easier on the eye."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "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.50.0 (0)\n -->\n<!-- Pages: 1 -->\n<svg width=\"360pt\" height=\"326pt\"\n viewBox=\"0.00 0.00 360.00 325.80\" xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\">\n<g id=\"graph0\" class=\"graph\" transform=\"scale(0.88 0.88) rotate(0) translate(4 366.88)\">\n<polygon fill=\"white\" stroke=\"transparent\" points=\"-4,4 -4,-366.88 405.81,-366.88 405.81,4 -4,4\"/>\n<!-- vara1 -->\n<g id=\"node1\" class=\"node\">\n<title>vara1</title>\n<ellipse fill=\"none\" stroke=\"black\" cx=\"158.81\" cy=\"-240.47\" rx=\"27\" ry=\"18\"/>\n<text text-anchor=\"middle\" x=\"158.81\" y=\"-236.77\" font-family=\"Times,serif\" font-size=\"14.00\">a1</text>\n</g>\n<!-- factor0 -->\n<g id=\"node13\" class=\"node\">\n<title>factor0</title>\n<ellipse fill=\"black\" stroke=\"black\" cx=\"16.4\" cy=\"-300.09\" rx=\"1.8\" ry=\"1.8\"/>\n</g>\n<!-- vara1&#45;&#45;factor0 -->\n<g id=\"edge1\" class=\"edge\">\n<title>vara1&#45;&#45;factor0</title>\n<path fill=\"none\" stroke=\"black\" d=\"M135.93,-250.05C99.14,-265.46 29.74,-294.51 18.09,-299.39\"/>\n</g>\n<!-- vara2 -->\n<g id=\"node2\" class=\"node\">\n<title>vara2</title>\n<ellipse fill=\"none\" stroke=\"black\" cx=\"230.81\" cy=\"-240.47\" rx=\"27\" ry=\"18\"/>\n<text text-anchor=\"middle\" x=\"230.81\" y=\"-236.77\" font-family=\"Times,serif\" font-size=\"14.00\">a2</text>\n</g>\n<!-- vara2&#45;&#45;vara1 -->\n<g id=\"edge3\" class=\"edge\">\n<title>vara2&#45;&#45;vara1</title>\n<path fill=\"none\" stroke=\"black\" d=\"M203.59,-240.47C197.86,-240.47 191.81,-240.47 186.07,-240.47\"/>\n</g>\n<!-- factor1 -->\n<g id=\"node14\" class=\"node\">\n<title>factor1</title>\n<ellipse fill=\"black\" stroke=\"black\" cx=\"105.23\" cy=\"-326.2\" rx=\"1.8\" ry=\"1.8\"/>\n</g>\n<!-- vara2&#45;&#45;factor1 -->\n<g id=\"edge2\" class=\"edge\">\n<title>vara2&#45;&#45;factor1</title>\n<path fill=\"none\" stroke=\"black\" d=\"M211.74,-253.49C179.63,-275.41 117.23,-318 106.74,-325.16\"/>\n</g>\n<!-- vara3 -->\n<g id=\"node3\" class=\"node\">\n<title>vara3</title>\n<ellipse fill=\"none\" stroke=\"black\" cx=\"302.81\" cy=\"-240.47\" rx=\"27\" ry=\"18\"/>\n<text text-anchor=\"middle\" x=\"302.81\" y=\"-236.77\" font-family=\"Times,serif\" font-size=\"14.00\">a3</text>\n</g>\n<!-- vara3&#45;&#45;vara2 -->\n<g id=\"edge5\" class=\"edge\">\n<title>vara3&#45;&#45;vara2</title>\n<path fill=\"none\" stroke=\"black\" d=\"M275.59,-240.47C269.86,-240.47 263.81,-240.47 258.07,-240.47\"/>\n</g>\n<!-- factor3 -->\n<g id=\"node15\" class=\"node\">\n<title>factor3</title>\n<ellipse fill=\"black\" stroke=\"black\" cx=\"208.07\" cy=\"-345.52\" rx=\"1.8\" ry=\"1.8\"/>\n</g>\n<!-- vara3&#45;&#45;factor3 -->\n<g id=\"edge4\" class=\"edge\">\n<title>vara3&#45;&#45;factor3</title>\n<path fill=\"none\" stroke=\"black\" d=\"M288.83,-255.97C264.95,-282.45 218.02,-334.49 209.43,-344.01\"/>\n</g>\n<!-- vara4 -->\n<g id=\"node4\" class=\"node\">\n<title>vara4</title>\n<ellipse fill=\"none\" stroke=\"black\" cx=\"374.81\" cy=\"-240.47\" rx=\"27\" ry=\"18\"/>\n<text text-anchor=\"middle\" x=\"374.81\" y=\"-236.77\" font-family=\"Times,serif\" font-size=\"14.00\">a4</text>\n</g>\n<!-- vara4&#45;&#45;vara3 -->\n<g id=\"edge7\" class=\"edge\">\n<title>vara4&#45;&#45;vara3</title>\n<path fill=\"none\" stroke=\"black\" d=\"M347.59,-240.47C341.86,-240.47 335.81,-240.47 330.07,-240.47\"/>\n</g>\n<!-- factor5 -->\n<g id=\"node16\" class=\"node\">\n<title>factor5</title>\n<ellipse fill=\"black\" stroke=\"black\" cx=\"336.01\" cy=\"-361.08\" rx=\"1.8\" ry=\"1.8\"/>\n</g>\n<!-- vara4&#45;&#45;factor5 -->\n<g id=\"edge6\" class=\"edge\">\n<title>vara4&#45;&#45;factor5</title>\n<path fill=\"none\" stroke=\"black\" d=\"M369.09,-258.27C359.31,-288.67 340.09,-348.41 336.57,-359.34\"/>\n</g>\n<!-- varb1 -->\n<g id=\"node5\" class=\"node\">\n<title>varb1</title>\n<ellipse fill=\"none\" stroke=\"black\" cx=\"158.81\" cy=\"-168.47\" rx=\"27\" ry=\"18\"/>\n<text text-anchor=\"middle\" x=\"158.81\" y=\"-164.77\" font-family=\"Times,serif\" font-size=\"14.00\">b1</text>\n</g>\n<!-- varb1&#45;&#45;vara1 -->\n<g id=\"edge9\" class=\"edge\">\n<title>varb1&#45;&#45;vara1</title>\n<path fill=\"none\" stroke=\"black\" d=\"M158.81,-186.64C158.81,-197.47 158.81,-211.18 158.81,-222.06\"/>\n</g>\n<!-- factor7 -->\n<g id=\"node17\" class=\"node\">\n<title>factor7</title>\n<ellipse fill=\"black\" stroke=\"black\" cx=\"1.8\" cy=\"-183.58\" rx=\"1.8\" ry=\"1.8\"/>\n</g>\n<!-- varb1&#45;&#45;factor7 -->\n<g id=\"edge8\" class=\"edge\">\n<title>varb1&#45;&#45;factor7</title>\n<path fill=\"none\" stroke=\"black\" d=\"M131.82,-171.07C90.79,-175.02 16.13,-182.2 3.61,-183.4\"/>\n</g>\n<!-- varb2 -->\n<g id=\"node6\" class=\"node\">\n<title>varb2</title>\n<ellipse fill=\"none\" stroke=\"black\" cx=\"230.81\" cy=\"-168.47\" rx=\"27\" ry=\"18\"/>\n<text text-anchor=\"middle\" x=\"230.81\" y=\"-164.77\" font-family=\"Times,serif\" font-size=\"14.00\">b2</text>\n</g>\n<!-- varb2&#45;&#45;vara2 -->\n<g id=\"edge12\" class=\"edge\">\n<title>varb2&#45;&#45;vara2</title>\n<path fill=\"none\" stroke=\"black\" d=\"M230.81,-186.64C230.81,-197.47 230.81,-211.18 230.81,-222.06\"/>\n</g>\n<!-- varb2&#45;&#45;varb1 -->\n<g id=\"edge11\" class=\"edge\">\n<title>varb2&#45;&#45;varb1</title>\n<path fill=\"none\" stroke=\"black\" d=\"M203.59,-168.47C197.86,-168.47 191.81,-168.47 186.07,-168.47\"/>\n</g>\n<!-- factor9 -->\n<g id=\"node18\" class=\"node\">\n<title>factor9</title>\n<ellipse fill=\"black\" stroke=\"black\" cx=\"74.28\" cy=\"-194.08\" rx=\"1.8\" ry=\"1.8\"/>\n</g>\n<!-- varb2&#45;&#45;factor9 -->\n<g id=\"edge10\" class=\"edge\">\n<title>varb2&#45;&#45;factor9</title>\n<path fill=\"none\" stroke=\"black\" d=\"M204.26,-172.81C163.44,-179.49 88.65,-191.73 76.09,-193.78\"/>\n</g>\n<!-- varb3 -->\n<g id=\"node7\" class=\"node\">\n<title>varb3</title>\n<ellipse fill=\"none\" stroke=\"black\" cx=\"302.81\" cy=\"-168.47\" rx=\"27\" ry=\"18\"/>\n<text text-anchor=\"middle\" x=\"302.81\" y=\"-164.77\" font-family=\"Times,serif\" font-size=\"14.00\">b3</text>\n</g>\n<!-- varb3&#45;&#45;vara3 -->\n<g id=\"edge15\" class=\"edge\">\n<title>varb3&#45;&#45;vara3</title>\n<path fill=\"none\" stroke=\"black\" d=\"M302.81,-186.64C302.81,-197.47 302.81,-211.18 302.81,-222.06\"/>\n</g>\n<!-- varb3&#45;&#45;varb2 -->\n<g id=\"edge14\" class=\"edge\">\n<title>varb3&#45;&#45;varb2</title>\n<path fill=\"none\" stroke=\"black\" d=\"M275.59,-168.47C269.86,-168.47 263.81,-168.47 258.07,-168.47\"/>\n</g>\n<!-- factor12 -->\n<g id=\"node19\" class=\"node\">\n<title>factor12</title>\n<ellipse fill=\"black\" stroke=\"black\" cx=\"164.72\" cy=\"-154.4\" rx=\"1.8\" ry=\"1.8\"/>\n</g>\n<!-- varb3&#45;&#45;factor12 -->\n<g id=\"edge13\" class=\"edge\">\n<title>varb3&#45;&#45;factor12</title>\n<path fill=\"none\" stroke=\"black\" d=\"M275.86,-165.72C239.75,-162.04 178.82,-155.84 166.81,-154.61\"/>\n</g>\n<!-- varb4 -->\n<g id=\"node8\" class=\"node\">\n<title>varb4</title>\n<ellipse fill=\"none\" stroke=\"black\" cx=\"374.81\" cy=\"-168.47\" rx=\"27\" ry=\"18\"/>\n<text text-anchor=\"middle\" x=\"374.81\" y=\"-164.77\" font-family=\"Times,serif\" font-size=\"14.00\">b4</text>\n</g>\n<!-- varb4&#45;&#45;vara4 -->\n<g id=\"edge18\" class=\"edge\">\n<title>varb4&#45;&#45;vara4</title>\n<path fill=\"none\" stroke=\"black\" d=\"M374.81,-186.64C374.81,-197.47 374.81,-211.18 374.81,-222.06\"/>\n</g>\n<!-- varb4&#45;&#45;varb3 -->\n<g id=\"edge17\" class=\"edge\">\n<title>varb4&#45;&#45;varb3</title>\n<path fill=\"none\" stroke=\"black\" d=\"M347.59,-168.47C341.86,-168.47 335.81,-168.47 330.07,-168.47\"/>\n</g>\n<!-- factor15 -->\n<g id=\"node20\" class=\"node\">\n<title>factor15</title>\n<ellipse fill=\"black\" stroke=\"black\" cx=\"337.41\" cy=\"-277.88\" rx=\"1.8\" ry=\"1.8\"/>\n</g>\n<!-- varb4&#45;&#45;factor15 -->\n<g id=\"edge16\" class=\"edge\">\n<title>varb4&#45;&#45;factor15</title>\n<path fill=\"none\" stroke=\"black\" d=\"M368.8,-186.05C359.31,-213.84 341.54,-265.82 338.03,-276.09\"/>\n</g>\n<!-- varc1 -->\n<g id=\"node9\" class=\"node\">\n<title>varc1</title>\n<ellipse fill=\"none\" stroke=\"black\" cx=\"158.81\" cy=\"-96.47\" rx=\"27\" ry=\"18\"/>\n<text text-anchor=\"middle\" x=\"158.81\" y=\"-92.77\" font-family=\"Times,serif\" font-size=\"14.00\">c1</text>\n</g>\n<!-- varc1&#45;&#45;varb1 -->\n<g id=\"edge20\" class=\"edge\">\n<title>varc1&#45;&#45;varb1</title>\n<path fill=\"none\" stroke=\"black\" d=\"M158.81,-114.64C158.81,-125.47 158.81,-139.18 158.81,-150.06\"/>\n</g>\n<!-- factor18 -->\n<g id=\"node21\" class=\"node\">\n<title>factor18</title>\n<ellipse fill=\"black\" stroke=\"black\" cx=\"9.5\" cy=\"-64.55\" rx=\"1.8\" ry=\"1.8\"/>\n</g>\n<!-- varc1&#45;&#45;factor18 -->\n<g id=\"edge19\" class=\"edge\">\n<title>varc1&#45;&#45;factor18</title>\n<path fill=\"none\" stroke=\"black\" d=\"M132.8,-90.91C94.03,-82.62 24.26,-67.7 11.52,-64.98\"/>\n</g>\n<!-- varc2 -->\n<g id=\"node10\" class=\"node\">\n<title>varc2</title>\n<ellipse fill=\"none\" stroke=\"black\" cx=\"230.81\" cy=\"-96.47\" rx=\"27\" ry=\"18\"/>\n<text text-anchor=\"middle\" x=\"230.81\" y=\"-92.77\" font-family=\"Times,serif\" font-size=\"14.00\">c2</text>\n</g>\n<!-- varc2&#45;&#45;varb2 -->\n<g id=\"edge23\" class=\"edge\">\n<title>varc2&#45;&#45;varb2</title>\n<path fill=\"none\" stroke=\"black\" d=\"M230.81,-114.64C230.81,-125.47 230.81,-139.18 230.81,-150.06\"/>\n</g>\n<!-- varc2&#45;&#45;varc1 -->\n<g id=\"edge22\" class=\"edge\">\n<title>varc2&#45;&#45;varc1</title>\n<path fill=\"none\" stroke=\"black\" d=\"M203.59,-96.47C197.86,-96.47 191.81,-96.47 186.07,-96.47\"/>\n</g>\n<!-- factor20 -->\n<g id=\"node22\" class=\"node\">\n<title>factor20</title>\n<ellipse fill=\"black\" stroke=\"black\" cx=\"94.35\" cy=\"-35.1\" rx=\"1.8\" ry=\"1.8\"/>\n</g>\n<!-- varc2&#45;&#45;factor20 -->\n<g id=\"edge21\" class=\"edge\">\n<title>varc2&#45;&#45;factor20</title>\n<path fill=\"none\" stroke=\"black\" d=\"M208.28,-86.34C173.15,-70.54 108.12,-41.29 96.24,-35.95\"/>\n</g>\n<!-- varc3 -->\n<g id=\"node11\" class=\"node\">\n<title>varc3</title>\n<ellipse fill=\"none\" stroke=\"black\" cx=\"302.81\" cy=\"-96.47\" rx=\"27\" ry=\"18\"/>\n<text text-anchor=\"middle\" x=\"302.81\" y=\"-92.77\" font-family=\"Times,serif\" font-size=\"14.00\">c3</text>\n</g>\n<!-- varc3&#45;&#45;varb3 -->\n<g id=\"edge26\" class=\"edge\">\n<title>varc3&#45;&#45;varb3</title>\n<path fill=\"none\" stroke=\"black\" d=\"M302.81,-114.64C302.81,-125.47 302.81,-139.18 302.81,-150.06\"/>\n</g>\n<!-- varc3&#45;&#45;varc2 -->\n<g id=\"edge25\" class=\"edge\">\n<title>varc3&#45;&#45;varc2</title>\n<path fill=\"none\" stroke=\"black\" d=\"M275.59,-96.47C269.86,-96.47 263.81,-96.47 258.07,-96.47\"/>\n</g>\n<!-- factor23 -->\n<g id=\"node23\" class=\"node\">\n<title>factor23</title>\n<ellipse fill=\"black\" stroke=\"black\" cx=\"207.16\" cy=\"-4.64\" rx=\"1.8\" ry=\"1.8\"/>\n</g>\n<!-- varc3&#45;&#45;factor23 -->\n<g id=\"edge24\" class=\"edge\">\n<title>varc3&#45;&#45;factor23</title>\n<path fill=\"none\" stroke=\"black\" d=\"M287.02,-81.31C262.4,-57.67 216.82,-13.91 208.49,-5.91\"/>\n</g>\n<!-- varc4 -->\n<g id=\"node12\" class=\"node\">\n<title>varc4</title>\n<ellipse fill=\"none\" stroke=\"black\" cx=\"374.81\" cy=\"-96.47\" rx=\"27\" ry=\"18\"/>\n<text text-anchor=\"middle\" x=\"374.81\" y=\"-92.77\" font-family=\"Times,serif\" font-size=\"14.00\">c4</text>\n</g>\n<!-- varc4&#45;&#45;varb4 -->\n<g id=\"edge29\" class=\"edge\">\n<title>varc4&#45;&#45;varb4</title>\n<path fill=\"none\" stroke=\"black\" d=\"M374.81,-114.64C374.81,-125.47 374.81,-139.18 374.81,-150.06\"/>\n</g>\n<!-- varc4&#45;&#45;varc3 -->\n<g id=\"edge28\" class=\"edge\">\n<title>varc4&#45;&#45;varc3</title>\n<path fill=\"none\" stroke=\"black\" d=\"M347.59,-96.47C341.86,-96.47 335.81,-96.47 330.07,-96.47\"/>\n</g>\n<!-- factor26 -->\n<g id=\"node24\" class=\"node\">\n<title>factor26</title>\n<ellipse fill=\"black\" stroke=\"black\" cx=\"337.99\" cy=\"-1.8\" rx=\"1.8\" ry=\"1.8\"/>\n</g>\n<!-- varc4&#45;&#45;factor26 -->\n<g id=\"edge27\" class=\"edge\">\n<title>varc4&#45;&#45;factor26</title>\n<path fill=\"none\" stroke=\"black\" d=\"M367.89,-78.66C358.39,-54.24 342.09,-12.35 338.64,-3.48\"/>\n</g>\n</g>\n</svg>\n",
      "text/plain": [
       "<gtbook.display.show at 0x113dcffa0>"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "position_hints = {c:float(1-i) for i,c in enumerate(row_symbols)}\n",
    "show(graph, binary_edges=True, hints=position_hints)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e59dca51",
   "metadata": {},
   "source": [
    "## Inference using a Direct Solver\n",
    "\n",
    "As in Kalman smoothing, we can use multi-frontal elimination, using a \"nested dissection\" ordering:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "789bd863",
   "metadata": {},
   "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.50.0 (0)\n -->\n<!-- Title: G Pages: 1 -->\n<svg width=\"519pt\" height=\"188pt\"\n viewBox=\"0.00 0.00 519.14 188.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 184)\">\n<title>G</title>\n<polygon fill=\"white\" stroke=\"transparent\" points=\"-4,4 -4,-184 515.14,-184 515.14,4 -4,4\"/>\n<!-- 0 -->\n<g id=\"node1\" class=\"node\">\n<title>0</title>\n<ellipse fill=\"none\" stroke=\"black\" cx=\"205.39\" cy=\"-162\" rx=\"59.59\" ry=\"18\"/>\n<text text-anchor=\"middle\" x=\"205.39\" y=\"-158.3\" font-family=\"Times,serif\" font-size=\"14.00\">b1, a3, b2, c3</text>\n</g>\n<!-- 1 -->\n<g id=\"node2\" class=\"node\">\n<title>1</title>\n<ellipse fill=\"none\" stroke=\"black\" cx=\"62.39\" cy=\"-90\" rx=\"62.29\" ry=\"18\"/>\n<text text-anchor=\"middle\" x=\"62.39\" y=\"-86.3\" font-family=\"Times,serif\" font-size=\"14.00\">a1 : a3, b1, b2</text>\n</g>\n<!-- 0&#45;&gt;1 -->\n<g id=\"edge1\" class=\"edge\">\n<title>0&#45;&gt;1</title>\n<path fill=\"none\" stroke=\"black\" d=\"M175.42,-146.33C154.11,-135.9 125.28,-121.78 102.06,-110.42\"/>\n<polygon fill=\"black\" stroke=\"black\" points=\"103.35,-107.15 92.82,-105.9 100.27,-113.44 103.35,-107.15\"/>\n</g>\n<!-- 3 -->\n<g id=\"node4\" class=\"node\">\n<title>3</title>\n<ellipse fill=\"none\" stroke=\"black\" cx=\"205.39\" cy=\"-90\" rx=\"62.29\" ry=\"18\"/>\n<text text-anchor=\"middle\" x=\"205.39\" y=\"-86.3\" font-family=\"Times,serif\" font-size=\"14.00\">c1 : b1, b2, c3</text>\n</g>\n<!-- 0&#45;&gt;3 -->\n<g id=\"edge3\" class=\"edge\">\n<title>0&#45;&gt;3</title>\n<path fill=\"none\" stroke=\"black\" d=\"M205.39,-143.7C205.39,-135.98 205.39,-126.71 205.39,-118.11\"/>\n<polygon fill=\"black\" stroke=\"black\" points=\"208.89,-118.1 205.39,-108.1 201.89,-118.1 208.89,-118.1\"/>\n</g>\n<!-- 5 -->\n<g id=\"node6\" class=\"node\">\n<title>5</title>\n<ellipse fill=\"none\" stroke=\"black\" cx=\"361.39\" cy=\"-90\" rx=\"76.09\" ry=\"18\"/>\n<text text-anchor=\"middle\" x=\"361.39\" y=\"-86.3\" font-family=\"Times,serif\" font-size=\"14.00\">b3, b4 : a3, b2, c3</text>\n</g>\n<!-- 0&#45;&gt;5 -->\n<g id=\"edge5\" class=\"edge\">\n<title>0&#45;&gt;5</title>\n<path fill=\"none\" stroke=\"black\" d=\"M237.34,-146.67C260.76,-136.16 292.81,-121.77 318.47,-110.26\"/>\n<polygon fill=\"black\" stroke=\"black\" points=\"319.91,-113.45 327.6,-106.16 317.05,-107.06 319.91,-113.45\"/>\n</g>\n<!-- 2 -->\n<g id=\"node3\" class=\"node\">\n<title>2</title>\n<ellipse fill=\"none\" stroke=\"black\" cx=\"62.39\" cy=\"-18\" rx=\"61.99\" ry=\"18\"/>\n<text text-anchor=\"middle\" x=\"62.39\" y=\"-14.3\" font-family=\"Times,serif\" font-size=\"14.00\">a2 : a1, a3, b2</text>\n</g>\n<!-- 1&#45;&gt;2 -->\n<g id=\"edge2\" class=\"edge\">\n<title>1&#45;&gt;2</title>\n<path fill=\"none\" stroke=\"black\" d=\"M62.39,-71.7C62.39,-63.98 62.39,-54.71 62.39,-46.11\"/>\n<polygon fill=\"black\" stroke=\"black\" points=\"65.89,-46.1 62.39,-36.1 58.89,-46.1 65.89,-46.1\"/>\n</g>\n<!-- 4 -->\n<g id=\"node5\" class=\"node\">\n<title>4</title>\n<ellipse fill=\"none\" stroke=\"black\" cx=\"205.39\" cy=\"-18\" rx=\"61.99\" ry=\"18\"/>\n<text text-anchor=\"middle\" x=\"205.39\" y=\"-14.3\" font-family=\"Times,serif\" font-size=\"14.00\">c2 : b2, c1, c3</text>\n</g>\n<!-- 3&#45;&gt;4 -->\n<g id=\"edge4\" class=\"edge\">\n<title>3&#45;&gt;4</title>\n<path fill=\"none\" stroke=\"black\" d=\"M205.39,-71.7C205.39,-63.98 205.39,-54.71 205.39,-46.11\"/>\n<polygon fill=\"black\" stroke=\"black\" points=\"208.89,-46.1 205.39,-36.1 201.89,-46.1 208.89,-46.1\"/>\n</g>\n<!-- 6 -->\n<g id=\"node7\" class=\"node\">\n<title>6</title>\n<ellipse fill=\"none\" stroke=\"black\" cx=\"347.39\" cy=\"-18\" rx=\"48.99\" ry=\"18\"/>\n<text text-anchor=\"middle\" x=\"347.39\" y=\"-14.3\" font-family=\"Times,serif\" font-size=\"14.00\">a4 : a3, b4</text>\n</g>\n<!-- 5&#45;&gt;6 -->\n<g id=\"edge6\" class=\"edge\">\n<title>5&#45;&gt;6</title>\n<path fill=\"none\" stroke=\"black\" d=\"M357.93,-71.7C356.39,-63.98 354.54,-54.71 352.82,-46.11\"/>\n<polygon fill=\"black\" stroke=\"black\" points=\"356.21,-45.22 350.82,-36.1 349.34,-46.6 356.21,-45.22\"/>\n</g>\n<!-- 7 -->\n<g id=\"node8\" class=\"node\">\n<title>7</title>\n<ellipse fill=\"none\" stroke=\"black\" cx=\"462.39\" cy=\"-18\" rx=\"48.99\" ry=\"18\"/>\n<text text-anchor=\"middle\" x=\"462.39\" y=\"-14.3\" font-family=\"Times,serif\" font-size=\"14.00\">c4 : b4, c3</text>\n</g>\n<!-- 5&#45;&gt;7 -->\n<g id=\"edge7\" class=\"edge\">\n<title>5&#45;&gt;7</title>\n<path fill=\"none\" stroke=\"black\" d=\"M384.82,-72.76C398.85,-63.04 416.84,-50.57 431.96,-40.09\"/>\n<polygon fill=\"black\" stroke=\"black\" points=\"434.25,-42.76 440.48,-34.19 430.27,-37.01 434.25,-42.76\"/>\n</g>\n</g>\n</svg>\n",
      "text/plain": [
       "<gtbook.display.show at 0x1263344c0>"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "nested = gtsam.Ordering.MetisGaussianFactorGraph(graph)\n",
    "nested_bayes_tree = graph.eliminateMultifrontal(nested)\n",
    "show(nested_bayes_tree)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this case, the structure of the Bayes tree reveals that (for M=3,N=4) the image is split into three regions that are separately solved.\n",
    "\n",
    "Unfortunately, the wrapper does not support easy visualization of the \"divide and conquer\" approach of nested dissection, but here is an ascii image showing the first \"cut\", dividing up the image into regions handled by the left, middle, and right branches, respectively:\n",
    "\n",
    "```\n",
    "  1 2 3 4\n",
    "  _______\n",
    "a|l l * r|\n",
    "b|* * r r|\n",
    "c|m m * r|\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cd2d57b3",
   "metadata": {},
   "source": [
    "The MAP solution is obtained by back-substitution, and can be seen as a smoother version of the data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAArwAAAH0CAYAAADfWf7fAAAgAElEQVR4Xu2dd5hV1dWHF8zAUAUUMfaCihLUaMTE8lmj2FssARs2xAbSBEEUVBSkiQgiTQRFBawoChFEsRs1RpMYSyyxASIgFqR+z70+TpxhCGf2Ofuutfe8/KecVfb72/m+Nyd37lRbu3btWuEPBCAAAQhAAAIQgAAEIiVQDeGNNFmOBQEIQAACEIAABCCQJ4DwchEgAAEIQAACEIAABKImgPBGHS+HgwAEIAABCEAAAhBAeLkDEIAABCAAAQhAAAJRE0B4o46Xw0EAAhCAAAQgAAEIILzcAQhAAAIQgAAEIACBqAkgvFHHy+EgAAEIQAACEIAABBBe7gAEIAABCEAAAhCAQNQEEN6o4+VwEIAABCAAAQhAAAIIL3cAAhCAAAQgAAEIQCBqAghv1PFyOAhAAAIQgAAEIAABhJc7AAEIQAACEIAABCAQNQGEN+p4ORwEIAABCEAAAhCAAMLLHYAABCAAAQhAAAIQiJoAwht1vBwOAhCAAAQgAAEIQADh5Q5AAAIQgAAEIAABCERNAOGNOl4OBwEIQAACEIAABCCA8HIHIAABCEAAAhCAAASiJoDwRh0vh4MABCAAAQhAAAIQQHi5AxCAAAQgAAEIQAACURNAeKOOl8NBAAIQgAAEIAABCCC83AEIQAACEIAABCAAgagJILxRx8vhIAABCEAAAhCAAAQQXu4ABCAAAQhAAAIQgEDUBBDeqOPlcBCAAAQgAAEIQAACCC93AAIQgAAEIAABCEAgagIIb9TxcjgIQAACEIAABCAAAYSXOwABCEAAAhCAAAQgEDUBhDfqeDkcBCAAAQhAAAIQgADCyx2AAAQgAAEIQAACEIiaAMIbdbwcDgIQgAAEIAABCEAA4eUOQAACEIAABCAAAQhETQDhjTpeDgcBCEAAAhCAAAQggPByByAAAQhAAAIQgAAEoiaA8EYdL4eDAAQgAAEIQAACEEB4uQMQgAAEIAABCEAAAlETQHijjpfDQQACEIAABCAAAQggvNwBCEAAAhCAAAQgAIGoCSC8UcfL4SAAAQhAAAIQgAAEEF7uAAQgAAEIQAACEIBA1AQQ3qjj5XAQgAAEIAABCEAAAggvdwACEIAABCAAAQhAIGoCCG/U8XI4CEAAAhCAAAQgAAGElzsAAQhAAAIQgAAEIBA1AYQ36ng5HAQgAAEIQAACEIAAwssdgAAEIAABCEAAAhCImgDCG3W8HA4CEIAABCAAAQhAAOHlDkAAAhCAAAQgAAEIRE0A4Y06Xg4HAQhAAAIQgAAEIIDwcgcgAAEIQAACEIAABKImgPBGHS+HgwAEIAABCEAAAhBAeLkDEIAABCAAAQhAAAJRE0B4o46Xw0EAAhCAAAQgAAEIILwZ3IEVK1bK4qXfSpPGDaVatWoZdKQFBCAAAQhAAAIQgEBWBBDeFCTXrl0rt098VEbc+VC+y8YN68ttN14hezRvWmHX2fNelw69b13n716fNUZKatZIsQmlEIAABCAAAQhAAALrI4Dwprgbb7z9npx5WT+ZNLyn7LbLDnLruAfl8dkvylP3D5Hq1dd90/vUvNfkqhvHyLQxfctM3WbLJrwZTpEDpRCAAAQgAAEIQOB/EUB4U9yPwaOmyD/f/1jGDuqW77LgqyVyyClX5IV21522XadzTnj7Dp4g8x4enmIqpRCAAAQgAAEIQAAClSGA8FaGVrlnu153uzRqUE96dTyr9G9+fXBbGXlTJzlo3z0qFN6OvYfLCa32l5KSmrL3Hs2k1cEtpbioKMUWlEIAAhCAAAQgAAEI8Ia3knfg8y+/ksdnv7TeqjP/eITUrlVT2nUbJM2abiNd2p9W+mzLo9pLn65t5ZjDfr9O/VvvfCgz574iDerXlc/nL5Ipjz4tbU46rIwwV3JVHocABCAAAQhAAAIQ2AAB3vBWAOjjT+fLfY/MWS+6y887SerUriW5N7y5H1Tr2eHMRG94yzd8cMaz0vvm8fLm7HH5t7xP7DmICxsxgR9+KIn4dBxt193eA0LEBOpu/G3Ep+No29wxHgiRE0B4UwSc+wzvvz74REYP7JrvsqHP8JYfNe/lt6R998Hy2szRUqukJsKbIosQShHeEFJy3xHhdWcXQiXCG0JK7jsivO7sQqlEeFMk9d9vaeglu+26gwwbO01mzH6p9FsaJkx5UnJfRZb7Fofcn8kPzZZmTbeW5jtvJ0uXfSvdrhslNYqLZPzQ7vm/5w1vijACKEV4AwgpxYoIbwp4AZQivAGElGJFhDcFvEBKEd4UQeW+h/e2Ox+SURMfzXfJfcxh9MAusmeLnfL/PHDkfTJl+lx59YlR+X8ecscUGXfvjNKJuzdvKgN7t5etNt8U4U2RQyilCG8oSbntifC6cQulCuENJSm3PRFeN24hVSG8GaS1/McV8vXib+RXTTap8Pt3fzki9+zCRUukft060rBBvTLTecObQRiGWyC8hsPJYDWENwOIhlsgvIbDyWA1hDc5xDVr1kruhV9RUfXkRQaeRHgNhPDzCgivoTA8rILweoBqqCXCaygMD6sgvB6gGmqJ8CYLIye6fQZPyD/ct+u5yYqMPIXwGgkitwbCaygMD6sgvB6gGmqJ8BoKw8MqCK8HqIZaIrwbDiP3tao33DJJvl6yTE459iCEd8PIeGJ9BBDeuO8Gwht3vghv3PkivHHni/BuON/vf/hRvvn2Oxk6emr+m6V4w7thZjyxHgIIb9xXA+GNO1+EN+58Ed6480V4k+d73dCJsnr1aoQ3OTKeLE8A4Y37TiC8ceeL8MadL8Ibd74Ib/J8Ed7krHiSN7xV8g4gvHHHjvDGnS/CG3e+loV36Y87qMBvUPLvCucivCpxxDWUN7xx5Vn+NAhv3PkivHHni/DGna9p4V2uJLy1EN64b73i6RBeRfgFGI3wFgCy4giEVxF+AUYjvAWArDjCtPB+ryS8dcoK7+rVa2TNmjVyw7BJsmrVaunTpa0UFRVt8PcPKMZaZjRfS2YlCb6WzFASflZBeP1wtdIV4bWShJ89EF4/XK10tSy833ynI7wb1S0rvFMefVr6DrmrTGTXX3menHz0gVZi/J97ILyGYuINr6EwPKyC8HqAaqglwmsoDA+rILweoBpqaVp4v1US3noVf6TBUGyVWgXhrRQuvw8jvH75andHeLUT8Dsf4fXLV7s7wqudgN/5poV3mZLw1kd4/d66Ktwd4Y07fIQ37nwR3rjzRXjjztey8C5bqiO89RsgvHHfesXTIbyK8AswGuEtAGTFEQivIvwCjEZ4CwBZcYRp4V2iJLwNEV7FKxn3aIQ37nwR3rjzRXjjzhfhjTtf08K7uKkK/PqNPlCZ62son+H1RdahL8LrAC2gEoQ3oLAcVkV4HaAFVILwBhSWw6qWhffbr3WEt97GCK/DVaIkCQGENwmlcJ9BeMPNLsnmCG8SSuE+g/CGm12SzU0L71dKwtsY4U1yd3jGgQDC6wAtoBKEN6CwHFZFeB2gBVSC8AYUlsOqpoV3oZLwborwOlwlSpIQQHiTUAr3GYQ33OySbI7wJqEU7jMIb7jZJdncsvB+t0BHeOs2QXiT3B2ecSCA8DpAC6gE4Q0oLIdVEV4HaAGVILwBheWwqmnhna8kvJshvA5XiZIkBBDeJJTCfQbhDTe7JJsjvEkohfsMwhtudkk2Ny28XygJ7+YIb5K7wzMOBBBeB2gBlSC8AYXlsCrC6wAtoBKEN6CwHFa1LLzff76jw4nSl9TZ4v30TQx14GvJDIWB8BoKw8MqCK8HqIZaIryGwvCwCsLrAaqhlqaF9zMl4d0S4TV0ReNaBeGNK8/yp0F4484X4Y07X4Q37nwtC+8Pn+oIb+2tEN64b73i6RBeRfgFGI3wFgCy4giEVxF+AUYjvAWArDjCtPB+oiS82yC8ilcy7tEIb9z5Irxx54vwxp0vwht3vqaF92Ml4d0W4Y371iueDuFVhF+A0QhvASArjkB4FeEXYDTCWwDIiiMsC+/yj3SEt9Z2CK/ilYx7NMIbd74Ib9z5Irxx54vwxp2vaeH9UEl4t0d44771iqdDeBXhF2A0wlsAyIojEF5F+AUYjfAWALLiCNPC+28l4d0B4VW8knGPRnjjzhfhjTtfhDfufBHeuPO1LLw/vr+TCvySHd9TmetrKN/D64usQ1+E1wFaQCUIb0BhOayK8DpAC6gE4Q0oLIdVTQvve0rCuxPC63CVKElCAOFNQincZxDecLNLsjnCm4RSuM8gvOFml2Rz08L7rpLw7ozwJrk7PONAAOF1gBZQCcIbUFgOqyK8DtACKkF4AwrLYVXLwrviXzrCW7MZwutwlShJQgDhTUIp3GcQ3nCzS7I5wpuEUrjPILzhZpdkc9PC+08l4d0V4U1yd3jGgQDC6wAtoBKEN6CwHFZFeB2gBVSC8AYUlsOqpoX3H0rC2xzhdbhKlCQhgPAmoRTuMwhvuNkl2RzhTUIp3GcQ3nCzS7K5ZeFd+Xcd4a3xa4Q3yd3hGQcCCK8DtIBKEN6AwnJYFeF1gBZQCcIbUFgOq5oW3rd3djhR+pIaLd5N38RQB76WzFAYCK+hMDysgvB6gGqoJcJrKAwPqyC8HqAaamlaeP+mJLy7I7yGrmhcqyC8ceVZ/jQIb9z5Irxx54vwxp2vZeFd9aaO8BbvgfDGfesVT4fwKsIvwGiEtwCQFUcgvIrwCzAa4S0AZMURpoX3r0rC+xuEV/FKxj0a4Y07X4Q37nwR3rjzRXjjzte08L6hJLx7Irxx33rF0yG8ivALMBrhLQBkxREIryL8AoxGeAsAWXGEZeFd/ZqO8Bb9FuFVvJJxj0Z4484X4Y07X4Q37nwR3rjzNS28f1ES3r0R3rhvveLpEF5F+AUYjfAWALLiCIRXEX4BRiO8BYCsOMK08L6qJLwtEV7FKxn3aIQ37nwR3rjzRXjjzhfhjTtf08L7SjMV+EX7/Etlrq+hfA+vL7IOfRFeB2gBlSC8AYXlsCrC6wAtoBKEN6CwHFY1LbwvKwnv7xBeh6tESRICCG8SSuE+g/CGm12SzRHeJJTCfQbhDTe7JJubFt4Xd0lyhMyfKdr3ncx7ajbkDa8m/XKzEV5DYXhYBeH1ANVQS4TXUBgeVkF4PUA11NK08L6gJLz7IbyGrmhcqyC8ceVZ/jQIb9z5Irxx54vwxp2vaeF9flcV+EX7/1Nlrq+hvOH1RdahL8LrAC2gEoQ3oLAcVkV4HaAFVILwBhSWw6qmhfe55g4nSl9SdMA/0jcx1AHhVQhj2bffy6rVq6VRg/plpiO8CmEUcCTCW0DYCqMQXgXoBRyJ8BYQtsIo08I7T0l4/w/hVbiKcYz8/ofl0v2GO2TO82/kD7R786Yy/IYO0njjBvl/RnjjyHl9p0B4484X4Y07X4Q37nxNC++zv1aBX3Tg31Xm+hrKG15fZCvoO3by4zJ1+lyZNLyX1K5VUy7uMVS232Zzuf7K8xDeAuagNQrh1SJfmLkIb2E4a01BeLXIF2auaeF9Rkl4D0J4C3P7IpxyyoXXSquDW8qFZxybP93Mua9I5z4j5e2n75Rq1arxhjfCzH95JIQ37oAR3rjzRXjjzte08D7dQgV+0SFvq8z1NZQ3vL7IVtC35VHt5Ybu5+elN/fnH+9+JKe26yMvTB8hDerXRXgLmIXGKIRXg3rhZiK8hWOtMQnh1aBeuJmmhXeOkvAeivAW7gZGNGnt2rXS4pBzZeRNneSgfffIn+yDjz6T49v2kqfuHyybb7YJwhtR3hUdBeGNO2CEN+58Ed648zUtvLN3U4FfdNhbKnN9DeUNry+y63nD26/HBXLEQXvzhreA3K2MQnitJOFnD4TXD1crXRFeK0n42cO08D6lJLx/QHj93LYq0DX3Gd4jD9lHLmhzTP60fIa3CoT+iyMivHHnjfDGnS/CG3e+poX3z7urwC86/G8qc30N5Q2vL7IV9B1zz2My7bFn8t/SUKd2ibTvPoRvaSggf+1RCK92An7nI7x++Wp3R3i1E/A737TwzlIS3iMQXr+3LuLu332/XLped7s8+9Kb+VO2aLa9DO/XUZo0bpj/Z76HN+LwRQThjTtfhDfufBHeuPM1Lbwzf/q5n0L/KWr1k6vE8oc3vApJLl32naxcuar0F078vALCqxBGAUcivAWErTAK4VWAXsCRCG8BYSuMMi28TygJ71EIr8JVrBojEd64c0Z4484X4Y07X4Q37nxNC++M36jALzr6rypzfQ3lDa8vsg59EV4HaAGVILwBheWwKsLrAC2gEoQ3oLAcVjUtvI8rCe8xCK/DVaIkCQGENwmlcJ9BeMPNLsnmCG8SSuE+g/CGm12SzS0L76rH9kxyhMyfKT72jcx7ajbkDa8m/XKzEV5DYXhYBeH1ANVQS4TXUBgeVkF4PUA11NK08E5XEt7jEF5DVzSuVRDeuPIsfxqEN+58Ed6480V4487XtPA+upcK/OLjX1eZ62sob3h9kXXoi/A6QAuoBOENKCyHVRFeB2gBlSC8AYXlsKpp4X3ktw4nSl9SfMJr6ZsY6oDwGgoD4TUUhodVEF4PUA21RHgNheFhFYTXA1RDLU0L70NKwnsSwmvoisa1CsIbV57lT4Pwxp0vwht3vghv3PmaFt4H91aBX3zyX1Tm+hrKG15fZB36IrwO0AIqQXgDCsthVYTXAVpAJQhvQGE5rGpaeB9QEt4/IrwOV4mSJAQQ3iSUwn0G4Q03uySbI7xJKIX7DMIbbnZJNjctvNNaJjlC5s8Un/Jq5j01G/KGV5N+udkIr6EwPKyC8HqAaqglwmsoDA+rILweoBpqaVp4pyoJ76kIr6ErGtcqCG9ceZY/DcIbd74Ib9z5Irxx52taeKfsowK/+LRXVOb6GsobXl9kHfoivA7QAipBeAMKy2FVhNcBWkAlCG9AYTmsalp471cS3tMRXoerREkSAghvEkrhPoPwhptdks0R3iSUwn0G4Q03uySbmxbe+36X5AiZP1P8p5cz76nZkDe8mvTLzUZ4DYXhYRWE1wNUQy0RXkNheFgF4fUA1VBL08I7WUl42yC8hq5oXKsgvHHlWf40CG/c+SK8ceeL8Madr2nhvef3KvCLz3hJZa6vobzh9UXWoS/C6wAtoBKEN6CwHFZFeB2gBVSC8AYUlsOqpoX3biXhPRPhdbhKlCQhgPAmoRTuMwhvuNkl2RzhTUIp3GcQ3nCzS7K5aeGdtG+SI2T+TPFZL2beU7Mhb3g16ZebjfAaCsPDKgivB6iGWiK8hsLwsArC6wGqoZamhXeikvCejfAauqJxrYLwxpVn+dMgvHHni/DGnS/CG3e+poX3rv1U4Bef84LKXF9DecPri6xDX4TXAVpAJQhvQGE5rIrwOkALqAThDSgsh1VNC+8EJeFti/A6XCVKkhBAeJNQCvcZhDfc7JJsjvAmoRTuMwhvuNkl2dyy8K68c/8kR8j8mRrnPp95T82GvOHVpF9uNsJrKAwPqyC8HqAaaonwGgrDwyoIrweohlqaFt7xB6iQqnHecypzfQ1FeH2RdeiL8DpAC6gE4Q0oLIdVEV4HaAGVILwBheWwqmXhXTFOR3hrno/wOlwlSpIQQHiTUAr3GYQ33OySbI7wJqEU7jMIb7jZJdnctPCO/b8kR8j8mZoXzKuw51dfL5U6tWtJndolmc/02ZA3vD7pVrI3wltJYIE9jvAGFlgl10V4KwkssMcR3sACq+S6poV3zIGVPE02j9e88NkyjT75bL607z5EPv50fv7fn3z0gXJN53OkRnFRhQMHjLhXJk6dWebv9myxk9x9W69sFqxkF4S3ksB8Po7w+qSr3xvh1c/A5wYIr0+6+r0RXv0MfG5gWXh/HK0jvCXtygpvu26DpF7d2tKvx4Xy5YJFctpFfeWaTmfLcUdU/C0S/W+bLP/5fIFceUnr0uhKSmrIrzbd2GeU6+2N8Kpgr3gowmsoDA+rILweoBpqifAaCsPDKgivB6iGWpoW3jsOUiFVctEzpXOXLvtO9jvu0vzb2dxb2tyffsMmyZcLvpbh/TpWuF9OeJd8863079lOZf/yQxFeEzH8tATCaygMD6sgvB6gGmqJ8BoKw8MqCK8HqIZamhbeUQerkCppP7d07gcffSbHt+0lcx+4RTbdpGH+30+aNksemfm8TBvTd73CO+uZV+X3ezWXRg3qy6EH7CW/3X1nlbPkhiK8aujXHYzwGgrDwyoIrweohloivIbC8LAKwusBqqGWloV3+e06wlvr4v8K7xtvvydnXtZPXpg+QhrUr5tPbsr0uTJq4iMyZ+rQCpOcPusF+ejTL6WkZg15+18fyux5r8uQPpdIq4P3UUke4VXBXvFQhNdQGB5WQXg9QDXUEuE1FIaHVRBeD1ANtTQtvCMOUSFV69Kn13nD+8yDw6Txxg0SveEtv3SPG0fLkqXLZNSALirnQXhVsCO8hrAXbBWEt2CoVQYhvCrYCzYU4S0YapVBloX3h9sOVWFS+7I5pXMr+gzv9UMnyoKvFq/3M7zll75lzDR57W/vyqThPVXOg/CqYEd4DWEv2CoIb8FQqwxCeFWwF2wowlsw1CqDTAvvcCXhvfy/wpsL5YKuA2WjenWlX48L1vmWhmXffi/ndhog57c+Wo469Hf5DIeOnirHH7GfbLPVr+RfH3wi514xQC5oc4xcdNZxKhkjvCrYEV5D2Au2CsJbMNQqgxBeFewFG4rwFgy1yiDTwnvrYSpManeYXWbuh598kf8e3k+/WJj/9yceeYD06dJWatQolqXffCf7HX+pXH3FWdL6xJ/2Pf2ivvnP7v78J/d8705nS62SmirnQXhVsCO8hrAXbBWEt2CoVQYhvCrYCzYU4S0YapVBloX3+2F/UGFSp+NTFc6dv3Bx/vt469aptcG9cm9+Fy9dJptu0khq19IR3Z+XRHg3GFfhHuCH1grHWmMSwqtBvXAzEd7CsdaYhPBqUC/cTNPCe4uS8F5RsfAWLpVsJyG82fJM1Q3hTYXPfDHCaz6iVAsivKnwmS9GeM1HlGpB08I79PBUZ3MtrtPpz66lJusQXkOxILyGwvCwCsLrAaqhlgivoTA8rILweoBqqKVl4f1uyBEqpOp2nqUy19dQhNcXWYe+CK8DtIBKEN6AwnJYFeF1gBZQCcIbUFgOq5oW3sFKwtsF4XW4SpQkIYDwJqEU7jMIb7jZJdkc4U1CKdxnEN5ws0uyuWXh/XZQqyRHyPyZel1nZt5TsyFveDXpl5uN8BoKw8MqCK8HqIZaIryGwvCwCsLrAaqhlqaFd+CRKqTqdXtSZa6voQivL7IOfRFeB2gBlSC8AYXlsCrC6wAtoBKEN6CwHFY1Lbw3KwnvlQivw1WiJAkBhDcJpXCfQXjDzS7J5ghvEkrhPoPwhptdks0tC++yAUclOULmz9Tv/kTmPTUb8oZXk3652QivoTA8rILweoBqqCXCaygMD6sgvB6gGmppWnj7H61Cqn6PGSpzfQ1FeH2RdeiL8DpAC6gE4Q0oLIdVEV4HaAGVILwBheWwqmnhvUlJeK9CeB2uEiVJCCC8SSiF+wzCG252STZHeJNQCvcZhDfc7JJsbll4v7nxmCRHyPyZjXo+nnlPzYa84dWkX242wmsoDA+rILweoBpqifAaCsPDKgivB6iGWpoW3n7HqpDaqNdjKnN9DUV4fZF16IvwOkALqAThDSgsh1URXgdoAZUgvAGF5bCqaeG9QUl4r0Z4Ha4SJUkIILxJKIX7DMIbbnZJNkd4k1AK9xmEN9zskmxuWXiXXn9ckiNk/kyD3tMz76nZkDe8GdBfsWKlLF76rTRp3FCqVavm3BHhdUYXRCHCG0RMzksivM7ogihEeIOIyXlJ08J73fHO50pT2OCaR9OUm6tFeFNEsnbtWrl94qMy4s6H8l02blhfbrvxCtmjedMKu86e97p06H3rOn/3+qwxUlKzhiC8KcIIoBThDSCkFCsivCngBVCK8AYQUooVLQvvkr46wtvwWoQ3xZWKq/SNt9+TMy/rJ5OG95TddtlBbh33oDw++0V56v4hUr36um96n5r3mlx14xiZNqZvGRDbbNkk/2YY4Y3rfpQ/DcIbd74Ib9z5Irxx52taePucoAK/YZ9HVOb6Gsob3hRkB4+aIv98/2MZO6hbvsuCr5bIIadckRfaXXfadp3OOeHtO3iCzHt4eIVTEd4UYQRQivAGEFKKFRHeFPACKEV4AwgpxYqmhffaE1OczL20Yd+H3YsNViK8KULpet3t0qhBPenV8azSLr8+uK2MvKmTHLTvHhUKb8few+WEVvtLSUlN2XuPZtLq4JZSXFSUfxbhTRFGAKUIbwAhpVgR4U0BL4BShDeAkFKsaFl4F1+jI7yNrkN4U1ypuErbdRskzZpuI13an1Z6sJZHtZc+XdvKMYf9fp3DvvXOhzJz7ivSoH5d+Xz+Ipny6NPS5qTDSoUZ4Y3rfpQ/DcIbd74Ib9z5Irxx52taeHufpAK/0fU//XxSLH94w5siydwb3twPqvXscGaiN7zlRz0441npffN4eXP2uPxb3iENJ6bYhlLrBBrX/9H6iuyXgkDzFv9OUU2pdQIltfjPr/WM0uy320ND0pR7rV189cle+6+veaMbHlSZ62sowpuCbO4zvP/64BMZPbBrvsuGPsNbftS8l9+S9t0Hy2szR0utkpoIb4osQihFeENIyX1HhNedXQiVCG8IKbnvaFl4v+6lI7wb90N43W9UZJX//ZaGXrLbrjvIsLHTZMbsl0q/pWHClCcl91VkuW9xyP2Z/NBsadZ0a2m+83aydNm30u26UVKjuEjGD+2e/3ve8EZ2QcodB+GNO1+EN+58Ed648zUtvD3/qAJ/4xsfUJnrayhveFOQzX0P7213PiSjJv70XXV1ateS0QO7yJ4tdsr/88CR98mU6XPl1SdG/SS0d0yRcffOKJ24e/OmMrB3e9lq800R3hQ5hAl4cUgAACAASURBVFKK8IaSlNueCK8bt1CqEN5QknLb07LwLrrqFLdDpaza5KZpKTvYKkd4M8hj+Y8r5OvF38ivmmxS4ffv/nJE7tmFi5ZI/bp1pGGDemWm84Y3gzAMt0B4DYeTwWoIbwYQDbdAeA2Hk8FqpoW3h5Lw9kd4M7hatKiIAMIb971AeOPOF+GNO1+EN+58LQvvVz1OVYHfuP9Ulbm+hvKG1xdZh74IrwO0gEoQ3oDCclgV4XWAFlAJwhtQWA6rmhbe7v/96lOHozmXNB4wxbnWYiHCaygVhNdQGB5WQXg9QDXUEuE1FIaHVRBeD1ANtbQsvAuv1BHeTW9GeA1d0bhWQXjjyrP8aRDeuPNFeOPOF+GNO1/TwtvtdBX4mw68X2Wur6G84fVF1qEvwusALaAShDegsBxWRXgdoAVUgvAGFJbDqpaFd0FXHeFtMgjhdbhKlCQhgPAmoRTuMwhvuNkl2RzhTUIp3GcQ3nCzS7K5aeHt8qckR8j8mSaD78u8p2ZD3vBq0i83G+E1FIaHVRBeD1ANtUR4DYXhYRWE1wNUQy0tC+/8zq1VSG025F6Vub6GIry+yDr0RXgdoAVUgvAGFJbDqgivA7SAShDegMJyWNWy8H7ZSUd4fzUU4XW4SpQkIYDwJqEU7jMIb7jZJdkc4U1CKdxnEN5ws0uyuWnhvaJNkiNk/syvbpmceU/Nhrzh1aRfbjbCaygMD6sgvB6gGmqJ8BoKw8MqCK8HqIZaWhbeLzqeoUJq82H3qMz1NRTh9UXWoS/C6wAtoBKEN6CwHFZFeB2gBVSC8AYUlsOqCO+60BBeh4tESTICCG8yTqE+hfCGmlyyvRHeZJxCfQrhDTW5ZHtbFt7PO5yZ7BAZP7XFrXdn3FG3HW94dfmXmY7wGgrDwyoIrweohloivIbC8LAKwusBqqGWpoX38rNUSG0xfJLKXF9DEV5fZB36IrwO0AIqQXgDCsthVYTXAVpAJQhvQGE5rGpZeD+7TEd4t7wN4XW4SpQkIYDwJqEU7jMIb7jZJdkc4U1CKdxnEN5ws0uyuWnhvfTsJEfI/JktR0zMvKdmQ97watIvNxvhNRSGh1UQXg9QDbVEeA2F4WEVhNcDVEMtLQvvp5foCO9WIxFeQ1c0rlUQ3rjyLH8ahDfufBHeuPNFeOPO17Lw/ufic1Tgb337XSpzfQ3lDa8vsg59EV4HaAGVILwBheWwKsLrAC2gEoQ3oLAcVjUtvO3bOpwofcnWoyakb2KoA8JrKAyE11AYHlZBeD1ANdQS4TUUhodVEF4PUA21tCy8n1ykI7zb3IHwGrqica2C8MaVZ/nTILxx54vwxp0vwht3vqaFt925KvC3GX2nylxfQ3nD64usQ1+E1wFaQCUIb0BhOayK8DpAC6gE4Q0oLIdVLQvvx+3OczhR+pJtR49P38RQB4TXUBgIr6EwPKyC8HqAaqglwmsoDA+rILweoBpqaVp4L1QS3jEIr6ErGtcqCG9ceZY/DcIbd74Ib9z5Irxx52tZeD+64HwV+NuNHacy19dQ3vD6IuvQF+F1gBZQCcIbUFgOqyK8DtACKkF4AwrLYVXTwnu+kvCOQ3gdrhIlSQggvEkohfsMwhtudkk2R3iTUAr3GYQ33OySbG5ZeD8874IkR8j8me3Hj828p2ZD3vBq0i83G+E1FIaHVRBeD1ANtUR4DYXhYRWE1wNUQy0tC++/z71QhdQOd45RmetrKMLri6xDX4TXAVpAJQhvQGE5rIrwOkALqAThDSgsh1VNC29bJeGdgPA6XCVKkhBAeJNQCvcZhDfc7JJsjvAmoRTuMwhvuNkl2dyy8H5wTrskR8j8maZ3jc68p2ZD3vBq0i83G+E1FIaHVRBeD1ANtUR4DYXhYRWE1wNUQy1NC+/ZF6mQajrxDpW5voYivL7IOvRFeB2gBVSC8AYUlsOqCK8DtIBKEN6AwnJY1bLwvn+WjvDuOAnhdbhKlCQhgPAmoRTuMwhvuNkl2RzhTUIp3GcQ3nCzS7K5beFtn+QImT+z46RRmffUbMgbXk365WYjvIbC8LAKwusBqqGWCK+hMDysgvB6gGqopWXhfe/Mi1VI7XT37SpzfQ1FeH2RdeiL8DpAC6gE4Q0oLIdVEV4HaAGVILwBheWwqmnhPUNJeO9BeB2uEiVJCCC8SSiF+wzCG252STZHeJNQCvcZhDfc7JJsbll4321zSZIjZP7MzpNHZt5TsyFveDXpl5uN8BoKw8MqCK8HqIZaIryGwvCwCsLrAaqhlpaF91+tdYS32b0Ir6ErGtcqCG9ceZY/DcIbd74Ib9z5Irxx52taeP90qQr8ZveNUJnrayhveH2RdeiL8DpAC6gE4Q0oLIdVEV4HaAGVILwBheWwqmXhfef0yxxOlL5kl/tvS9/EUAeE11AYCK+hMDysgvB6gGqoJcJrKAwPqyC8HqAaamlaeE9TEt4pCK+hKxrXKghvXHmWPw3CG3e+CG/c+SK8cedrWXj/eerlKvB3nTpcZa6vobzh9UXWoS/C6wAtoBKEN6CwHFZFeB2gBVSC8AYUlsOqloX3H6d2cDhR+pLmU29N38RQB4TXUBgIr6EwPKyC8HqAaqglwmsoDA+rILweoBpqaVl4/35KRxVSv542TGWur6EIry+yDn0RXgdoAZUgvAGF5bAqwusALaAShDegsBxWNS28f1QS3gcQXoerREkSAghvEkrhPoPwhptdks0R3iSUwn0G4Q03uySbWxbet0++IskRMn+mxYO3ZN5TsyFveDXpl5uN8BoKw8MqCK8HqIZaIryGwvCwCsLrAaqhlpaF962TOqmQ2u2hoSpzfQ1FeH2RdeiL8DpAC6gE4Q0oLIdVEV4HaAGVILwBheWwqmnhPVFJeB9GeB2uEiVJCCC8SSiF+wzCG252STZHeJNQCvcZhDfc7JJsbll4/3Zi5yRHyPyZ3R8eknlPzYa84dWkX242wmsoDA+rILweoBpqifAaCsPDKgivB6iGWloW3jdP6KJCao9HBqvM9TUU4fVF1qEvwusALaAShDegsBxWRXgdoAVUgvAGFJbDqpaF96/Hd3U4UfqS3zw6KH0TQx0QXkNhILyGwvCwCsLrAaqhlgivoTA8rILweoBqqKVp4T1OSXinI7yGrmiYq6xdu1ZWr1kjxUVFZQ6A8IaZZ9KtEd6kpMJ8DuENM7ekWyO8SUmF+Zxl4X3j2G4qUPd8bGCFc7/6eqnUqV1L6tQuUdnLdShveF3JpaibPusFGTpmqsyZWvYnIBHeFFADKEV4AwgpxYoIbwp4AZQivAGElGJFy8L7+jFXpjiZe+lej99cpviTz+ZL++5D5ONP5+f//clHHyjXdD5HahSXfXnnPtFvJcLrl+86l+XCroPk0y8WymabNkJ4C8jewiiE10IK/nZAeP2xtdAZ4bWQgr8dTAvv0UrCO6Os8LbrNkjq1a0t/XpcKF8uWCSnXdRXrul0thx3xH7+gsmwM8KbIcwNtVq1erXk/qeAOc+9IWMnP4bwbghYZH+P8EYWaLnjILxx54vwxp2vZeF97ejuKvB/O2NA6dyly76T/Y67VO6+rZfs2WKn/L/vN2ySfLngaxneT+dXH1cWCsJbWWIZPP/EnJdl4O33IbwZsAypBcIbUlqV3xXhrTyzkCoQ3pDSqvyuloX3L0f1qPyBMqjY+4n+pV0++OgzOb5tL5n7wC2y6SYN8/9+0rRZ8sjM52XamL4ZTPPfAuH1z3idCQivAnQDIxFeAyF4XAHh9QjXQGuE10AIHlewLLyvHnmVx5Ovv3XLJ28q/cs33n5Pzrysn7wwfYQ0qF83/++nTJ8royY+ss7LO5VlEwxFeBNAyvoRhDdromH0Q3jDyMl1S4TXlVwYdQhvGDm5bmlaeFspCe/M/wrvz294n3lwmDTeuAFveF0vWlWrQ3irWuI/nRfhjTt3hDfufBHeuPO1LLyvHNFTBf4+s24snVvRZ3ivHzpRFny1mM/wqqRjfGju+3dXrVotTz79Sv5ryWZOHijVqlcr/T5evpbMeIAp10N4UwI0Xo7wGg8o5XoIb0qAxsstC+/Lh/dSofe7P/crM/eCrgNlo3p1pV+PC/iWBpVEAhr6/oefyQnnlr24ua/z6N+zXf4UCG9AYTqsivA6QAuoBOENKCyHVRFeB2gBlVgW3pcOv1qF5O//fEOZuR9+8kX+e3hzX62a+3PikQdIny5tpUaNYpX9KjuUz/BWlpjH5xFej3ANtEZ4DYTgcQWE1yNcA60RXgMheFzBtPD+QUl4nyorvD/jn79wcf77eOvWqeUxkexbI7zZM3XuiPA6owuiEOENIibnJRFeZ3RBFCK8QcTkvKRl4X3xsN7O50pTuO/s69OUm6tFeA1FgvAaCsPDKgivB6iGWiK8hsLwsArC6wGqoZaWhfeFQ69RIbXfnOtU5voaivD6IuvQF+F1gBZQCcIbUFgOqyK8DtACKkF4AwrLYVXTwnuIkvA+jfA6XCVKkhBAeJNQCvcZhDfc7JJsjvAmoRTuMwhvuNkl2dyy8D5/8LVJjpD5M/vPDeM3qCU9OG94k5IqwHMIbwEgK45AeBXhF2A0wlsAyIojEF5F+AUYbVl4nzuoTwEIrDvigGd05vo6LMLri6xDX4TXAVpAJQhvQGE5rIrwOkALqAThDSgsh1UtC++8g3TetP7fMzpvlh3iS1SC8CbCVJiHEN7CcNaagvBqkS/MXIS3MJy1piC8WuQLM9e08B6oJLzPIryFuX1VcArCG3foCG/c+SK8ceeL8Madr2Xhffb/dH547MB5Oj8s5+um8YbXF1mHvgivA7SAShDegMJyWBXhdYAWUAnCG1BYDqtaFt5nDtD5PtyDntP5/l+H+BKVILyJMBXmIYS3MJy1piC8WuQLMxfhLQxnrSkIrxb5wsw1Lbz7Kwnv8whvYW5fFZyC8MYdOsIbd74Ib9z5Irxx52tZeOfuV/Gv+PWdyMEv6PxKY1/n4g2vL7IOfRFeB2gBlSC8AYXlsCrC6wAtoBKEN6CwHFa1LLxP79fP4UTpSw55oVf6JoY6ILyGwkB4DYXhYRWE1wNUQy0RXkNheFgF4fUA1VBLy8I7Z98bVUgd+mJPlbm+hiK8vsg69EV4HaAFVILwBhSWw6oIrwO0gEoQ3oDCcljVtPD+Xkl4X0J4Ha4SJUkIILxJKIX7DMIbbnZJNkd4k1AK9xmEN9zskmxuWXhn/+6mJEfI/JnDXr4q856aDXnDq0m/3GyE11AYHlZBeD1ANdQS4TUUhodVEF4PUA21tCy8T+3TX4XUH17poTLX11CE1xdZh74IrwO0gEoQ3oDCclgV4XWAFlAJwhtQWA6rWhbeP7fUEd7DX0V4Ha4SJUkIILxJKIX7DMIbbnZJNkd4k1AK9xmEN9zskmxuWXhn7T0gyREyf+aIv3TPvKdmQ97watIvNxvhNRSGh1UQXg9QDbVEeA2F4WEVhNcDVEMtLQvvTCXhbYXwGrqhka2C8EYWaLnjILxx54vwxp0vwht3vqaF97c3q8Bv9dqVKnN9DeUNry+yDn0RXgdoAZUgvAGF5bAqwusALaAShDegsBxWtSy8T+410OFE6UuOfL1b+iaGOiC8hsJAeA2F4WEVhNcDVEMtEV5DYXhYBeH1ANVQS8vC+8SeOsJ71BsIr6ErGtcqCG9ceZY/DcIbd74Ib9z5Irxx52tZeGf8ZpAK/KP/2lVlrq+hvOH1RdahL8LrAC2gEoQ3oLAcVkV4HaAFVILwBhSWw6qWhffxPQY7nCh9yTFvdknfxFAHhNdQGAivoTA8rILweoBqqCXCaygMD6sgvB6gGmppWXgfUxLeYxFeQzc0slUQ3sgCLXcchDfufBHeuPNFeOPO17LwTt99iAr84/7WWWWur6G84fVF1qEvwusALaAShDegsBxWRXgdoAVUgvAGFJbDqpaF99HdhjqcKH3J8W91St/EUAeE11AYCK+hMDysgvB6gGqoJcJrKAwPqyC8HqAaamlaeFsoCe/bCK+hKxrXKghvXHmWPw3CG3e+CG/c+SK8cedrWXgf+fUtKvBP+PsVKnN9DeUNry+yDn0RXgdoAZUgvAGF5bAqwusALaAShDegsBxWtSy8DzfXEd4T/4HwOlwlSpIQQHiTUAr3GYQ33OySbI7wJqEU7jMIb7jZJdncsvA+tOuwJEfI/JmT/tkx856aDXnDq0m/3GyE11AYHlZBeD1ANdQS4TUUhodVEF4PUA21tCy8D+56qwqpk//ZQWWur6EIry+yDn0RXgdoAZUgvAGF5bAqwusALaAShDegsBxWtSy8D+yiI7x/fAfhdbhKlCQhgPAmoRTuMwhvuNkl2RzhTUIp3GcQ3nCzS7K5ZeGd1mx4kiNk/swp/7o8856aDXnDq0m/3GyE11AYHlZBeD1ANdQS4TUUhodVEF4PUA21tCy8U3e+TYXUqe9epjLX11CE1xdZh74IrwO0gEoQ3oDCclgV4XWAFlAJwhtQWA6rmhbenZSE9z2E1+EqUZKEAMKbhFK4zyC84WaXZHOENwmlcJ9BeMPNLsnmloV3yo4jkhwh82dOe//SzHtqNuQNryb9crMRXkNheFgF4fUA1VBLhNdQGB5WQXg9QDXU0rLw3t90pAqp0z+4RGWur6EIry+yDn0RXgdoAZUgvAGF5bAqwusALaAShDegsBxWtSy89ykJ758QXoebREkiAghvIkzBPoTwBhtdosUR3kSYgn0I4Q02ukSLWxbee3e4PdEZsn6o9b8vzrqlaj/e8KriLzsc4TUUhodVEF4PUA21RHgNheFhFYTXA1RDLS0L7+TtdYS3zYcIr6ErGtcqCG9ceZY/DcIbd74Ib9z5Irxx52tZeO/ZbpQK/DM+aq8y19dQ3vD6IuvQF+F1gBZQCcIbUFgOqyK8DtACKkF4AwrLYVXTwrvtHQ4nSl9yxscXpW9iqAPCaygMhNdQGB5WQXg9QDXUEuE1FIaHVRBeD1ANtbQsvHdvoyO8Z36C8Bq6onGtgvDGlWf50yC8ceeL8MadL8Ibd76WhXfSNqNV4J/1STuVub6G8obXF1mHvgivA7SAShDegMJyWBXhdYAWUAnCG1BYDqtaFt6JW49xOFH6krP/c2H6JoY6ILyGwkB4DYXhYRWE1wNUQy0RXkNheFgF4fUA1VBLy8J711Y6wnvOpwivoSsa1yoIb1x5lj8Nwht3vghv3PkivHHna1l4J2w5VgV+288uUJnrayhveH2RdeiL8DpAC6gE4Q0oLIdVEV4HaAGVILwBheWwqmXhvXMLHeE993OE1+EqxV+yctVqqVFclOqgCG8qfOaLEV7zEaVaEOFNhc98McJrPqJUC1oW3vGbj0t1Ntfi874437XUZB1veDOI5ZPPFshRZ1wpf75vkGzxq8br7Th73uvSofet6/z967PGSEnNGoLwZhCG4RYIr+FwMlgN4c0AouEWCK/hcDJYzbTw/mp8BiesfIvzvjyv8kWGKxDelOG0vuR6+ds/Psh32ZDwPjXvNbnqxjEybUzfMlO32bKJVKtWDeFNmYX1coTXekLp9kN40/GzXo3wWk8o3X6WhXeckvCej/Cmu1SxVS/4aol8uWCR5MQ3ifD2HTxB5j08vEIMvOGN7XaUPQ/CG3e+CG/c+SK8cedrWXjHbnanCvwL5p+rMtfXUN7wZkB2/sLFcuipnRIJb8few+WEVvtLSUlN2XuPZtLq4JZSXPTTZ38R3gzCMNwC4TUcTgarIbwZQDTcAuE1HE4Gq1kW3tFNdIS33QKEN4OrFVeLpML71jsfysy5r0iD+nXl8/mLZMqjT0ubkw6TXh3PQnjjuhIVngbhjTtkhDfufBHeuPO1LLx3KAnvRQhv3Jfe5XRJhbd87wdnPCu9bx4vb84el3/L+3DzYS7jqYEABAwQ2GW39w1swQq+CNSq94Ov1vQ1QGC7cTpf/ZXk6KM2nZDkscyfab+wbeY9NRvykYYM6LsK77yX35L23QfLazNHS62SmghvBlnQAgJaBBBeLfKFmYvwFoaz1hTLwnu7kvBejPBqXUebc3Pfv5v7obUj21wpM+4ekP9asp+/j3fClCcl91Vkk4b3zC8/+aHZ0qzp1tJ85+1k6bJvpdt1o/LPjh/aPf/3vOG1mTFbQSAJAYQ3CaVwn0F4w80uyeaWhXdkY503vJd8xRveJHenyjzT8qj28v0Py0vPu3HD+qXfwjBw5H0yZfpcefWJUfm/H3LHFBl374zSZ3dv3lQG9m4vW22+KcJbZW4MB42VAMIba7I/nQvhjTtfy8I7Qkl4L0V44770vk+3/McVsnDREqlft440bFCvzDje8PqmT38I+COA8Ppja6EzwmshBX87WBbe2za5y9/B/0fnyxadozLX11A+w+uLrENfhNcBGiUQMEIA4TUShKc1EF5PYI20tSy8wzfWEd7Lv0Z4jVzP+NZAeOPLlBNVHQIIb9xZI7xx52tZeG9VEt4OCG/cl17zdAivJn1mQyAdAYQ3HT/r1Qiv9YTS7WdZeG9ppPOG94rFvOFNd6uoXi8BhJfLAYFwCSC84WaXZHOENwmlcJ+xLLxDG01UAdtp8dkqc30N5TO8vsg69EV4HaBRAgEjBBBeI0F4WgPh9QTWSFvLwjukoY7wdl6C8Bq5nvGtgfDGlyknqjoEEN64s0Z4487XsvAObqAjvF2WIrxx33rF0yG8ivAZDYGUBBDelACNlyO8xgNKuZ5l4R2kJLxdEd6Ut4ry9RJAeLkcEAiXAMIbbnZJNkd4k1AK9xnLwjtwo0kqYLt9c5bKXF9D+QyvL7IOfRFeB2iUQMAIAYTXSBCe1kB4PYE10tay8N6sJLxXIrxGbmeEayC8EYbKkaoMAYQ37qgR3rjztSy8A+rrvOHtvow3vHHfesXTIbyK8BkNgZQEEN6UAI2XI7zGA0q5nmXh7a8kvD0Q3pS3ivL1EkB4uRwQCJcAwhtudkk2R3iTUAr3GcvCe2O9u1XA9vz2TJW5vobyGV5fZB36IrwO0CiBgBECCK+RIDytgfB6AmukrWXh7VdXR3h7fYfwGrme8a2B8MaXKSeqOgQQ3rizRnjjztey8N6gJLxXI7xxX3rN0yG8mvSZDYF0BBDedPysVyO81hNKt59l4b2+zj3pDudY3fv7MxwrRdasWSsLFi2Wxhs3kOKiIuc+WRbykYYsaabshfCmBEg5BBQJILyK8AswGuEtAGTFEZaF9zol4b3GUXifefFN6Xrd7fL9D8vziV7b+Rw57fhD1pvu8ef0lA8+/rzM31/a9kS5pO2Jmd4IhDdTnOmaIbzp+FENAU0CCK8mff+zEV7/jDUnWBbevrV13vBe+0Pl3/D+sHyFHHhSB7nsvJPkjJP/IHNf+Kt07D1cZt47ULbafNMKI84J7zF/2FeOPGSf0r9vUL+uNGxQL9MrgfBmijNdM4Q3HT+qIaBJAOHVpO9/NsLrn7HmBMvC26eWjvD2WV554c293b3kqqHyxqwxUrNmjXykR5/ZPS+/Z5x8+HqFt+3pR8rJRx/o9QogvF7xVq45wls5XjwNAUsEEF5LaWS/C8KbPVNLHS0L77W1Jqug6ru8TaXnTpk+Vybc/4TMuHtAae3lvYbJdltvLl3an7Ze4a1bt7Y03XYL2WKzTeTYw/eVbbbcrNKzN1SA8G6IUAH/HuEtIGxGQSBjAghvxkCNtUN4jQWS8TqWhfeaEh3hve7HssI7fdYL8uXCrysk33zn7WT/li1k7OTH5cmnX5FpY/qWPpf7PG+9OrWlT9e2FdaOuPMhqV5UXdauFZnz3Ovy8afz5YGxfTOXXoQ34//QpGmH8KahRy0EdAkgvLr8fU9HeH0T1u1vWXivVhLeG8oJ7z0PPiWffrGwwqD22m0nOfzAvcXlDe8vG65cuUpatekmZ/3xCDn3T0dleikQ3kxxpmuG8KbjRzUENAkgvJr0/c9GeP0z1pxgWXh71dR5w9tvReU/0vDzZ3j/+uexUqNGcT7SVq27ydmnHrHez/CWz/30i/rKQfv9Ri4554RMrwTCmynOdM0Q3nT8qIaAJgGEV5O+/9kIr3/GmhMsC2/PmveqoLlxRetKz/3+hx+l5VEXSfdLW0ubCr6l4dW/viMDRtwrg6+9RLbdajP55LP5Muf5N/Lf0LBJowYy8+lXpHu/O2TirT3lt7vvXOn5/6sA4c0UZ7pmCG86flRDQJMAwqtJ3/9shNc/Y80JloX3qho6wnvTysoLby7DnMDmflDt5z9XX3GWtD7xsPw/Pv3CG3JZz2Hy4LjrpVnTrfPC2/aK/jJ/4eLS53OyfPaprTK/Dghv5kjdGyK87uyohIA2AYRXOwG/8xFev3y1u1sW3h7FOsLbf5Wb8OayXL16Tf4H3Jps0rD0ow3ry3jt2rXy9ZJl+V9Usflmm3j7zWwIr/Z/yn4xH+E1FAarQKCSBBDeSgIL7HGEN7DAKrmuZeHtriS8A1IIbyXxF+RxhLcgmJMNQXiTceIpCFgkgPBaTCW7nRDe7Fha7GRZeK8suk8F2c2r/6Qy19dQhNcXWYe+CK8DNEogYIQAwmskCE9rILyewBppa1l4uykJ70CE18jtjHANhDfCUDlSlSGA8MYdNcIbd76WhbdrdZ03vIPW8IY37luveDqEVxE+oyGQkgDCmxKg8XKE13hAKdezLLydlYR3CMKb8lZRvl4CCC+XAwLhEkB4w80uyeYIbxJK4T5jWXg7KQnvUIQ33AttfXOE13pC7AeB9RNAeOO+HQhv3PlaFt4rlIT3FoQ37kuveTqEV5M+syGQjgDCm46f9WqE13pC6fazLLwdlYR3GMKb7lJRvX4CCC+3AwLhEkB4w80uyeYIbxJK4T5jWXg7KAnvrQhvuBfa+uYIr/WE2A8CfKShqt4BhDfu5C0L7+XVdX7T2vA17r9pzeJt4Xt4DaWC8BoKg1UgUEkCvOGtJLDAHkd4AwusMXgc0AAAF2pJREFUkutaFt5LlYR3BMJbyVvE44kJILyJUfEgBMwRQHjNRZLpQghvpjjNNbMsvJcoCe9IhNfcPY1mIYQ3mig5SBUkgPDGHTrCG3e+loX3YiXhvR3hjfvSa54O4dWkz2wIpCOA8KbjZ70a4bWeULr9LAtv++qT0x3OsXrUmjaOlTbL+AyvoVwQXkNhsAoEKkkA4a0ksMAeR3gDC6yS61oW3ouUhPcOhLeSt4jHExNAeBOj4kEImCOA8JqLJNOFEN5McZprZll42ykJ72iE19w9jWYhhDeaKDlIFSSA8MYdOsIbd76WhfcCJeEdi/DGfek1T4fwatJnNgTSEUB40/GzXo3wWk8o3X6Whff86vekO5xj9bg1ZzhW2izjM7yGckF4DYXBKhCoJAGEt5LAAnsc4Q0ssEqua1l4z1MS3vEIbyVvEY8nJoDwJkbFgxAwRwDhNRdJpgshvJniNNfMsvCeqyS8dyK85u5pNAshvNFEyUGqIAGEN+7QEd6487UsvG2VhHcCwhv3pdc8HcKrSZ/ZEEhHAOFNx896NcJrPaF0+1kW3nOq353ucI7Vd60507HSZhmf4TWUC8JrKAxWgUAlCSC8lQQW2OMIb2CBVXJdy8J7lpLwTkJ4K3mLeDwxAYQ3MSoehIA5AgivuUgyXQjhzRSnuWaWhffMIp03vHev5g2vuYsa0kIrV62WrxYtkY0bbSQlNWuUWR3hDSlJdoVAWQIIb9w3AuGNO1/LwnuGkvDeg/DGfel9nm7MPY/JLWOmlY5odXBLubZzW2mwUd38v0N4fdKnNwT8EkB4/fLV7o7waifgd75l4W1TNMnv4dfTffLqs1Tm+hrKZ3h9ka2g79TH5srWWzSRPZrvKP/5fIGc33mAnN/6GGl7+pEIbwFzYBQEfBBAeH1QtdMT4bWThY9NLAtvayXhvRfh9XHVqmbP3jePl8++WCjjh3ZHeKvmFeDUERFAeCMKs4KjILxx52tZeP+kJLz3IbxxX/pCnS73Wd5WrbvKMYftK13an4bwFgo8cyDgiQDC6wmskbYIr5EgPK1hWXhPUxLeKQivp9tWxdpeO+hOmTH7ZXl8Un9p0rghwlvF8ue48RFAeOPL9JcnQnjjztey8J6qJLxTEd64L30hTjdywsMyYsLDct+oa2W3XbYvHckPrRWCPjMg4IcAwuuHq5WuCK+VJPzsYVl4Tyma6OfQG+g6bfXZKnN9DeWH1nyRraDvmjVrZfCo+2XK9Lly17Ae0nzn7co8hfAWMAxGQSBjAghvxkCNtUN4jQWS8TqWhfePSsL7AMKb8S2rQu2uHjBOHnpinowa0EV22Hbz0pNvtmkjKS4q4mvJqtBd4KjxEUB448v0lydCeOPO17LwnqwkvA8ivHFfep+na9W6m3z6xcJ1Rsy4e4Bsu9VmCK9P+PSGgGcCCK9nwMrtEV7lADyPtyy8JykJ70MIr+dbV4Xb85GGKhw+Rw+eAMIbfIT/8wAIb9z5WhbeE4ruUoH/yOpzVOb6GspneH2RdeiL8DpAowQCRgggvEaC8LQGwusJrJG2loX3eCXhfRThNXI7I1wD4Y0wVI5UZQggvHFHjfDGna9l4T1OSXinI7xxX3rN0yG8mvSZDYF0BBDedPysVyO81hNKt59l4T1WSXgfQ3jTXSqq108A4eV2QCBcAghvuNkl2RzhTUIp3GcsC+8xRRNUwD6+uq3KXF9D+QyvL7IOfRFeB2iUQMAIAYTXSBCe1kB4PYE10tay8B6tJLwzEF4jtzPCNRDeCEPlSFWGAMIbd9QIb9z5Whbeo5SE9wmEN+5Lr3k6hFeTPrMhkI4AwpuOn/VqhNd6Qun2syy8rZSEdybCm+5SUb1+AggvtwMC4RJAeMPNLsnmCG8SSuE+Y1l4jyi+UwXsrFXnqsz1NZTP8Poi69AX4XWARgkEjBBAeI0E4WkNhNcTWCNtLQvv4UrC+2eE18jtjHANhDfCUDlSlSGA8MYdNcIbd76Whfew4vEq8GevOk9lrq+hvOH1RdahL8LrAI0SCBghgPAaCcLTGgivJ7BG2loW3kOLx6lQmrPqfJW5voYivL7IOvRFeB2gUQIBIwQQXiNBeFoD4fUE1khby8J7sJLwzkV4jdzOCNdAeCMMlSNVGQIIb9xRI7xx52tZeA8qHqsC/5lVF6jM9TWUN7y+yDr0RXgdoFECASMEEF4jQXhaA+H1BNZIW8vC+39KwjsP4TVyOyNcA+GNMFSOVGUIILxxR43wxp2vZeE9oHiMCvznVl2oMtfXUN7w+iLr0BfhdYBGCQSMEEB4jQThaQ2E1xNYI20tC+9+NUarUHphZTuVub6GIry+yDr0RXgdoFECASMEEF4jQXhaA+H1BNZIW8vCu6+S8L6I8Bq5nRGugfBGGCpHqjIEEN64o0Z4487XsvD+rsYdKvBfXnmRylxfQ3nD64usQ1+E1wEaJRAwQgDhNRKEpzUQXk9gjbS1LLwta4xSofTqyvYqc30NRXh9kXXoi/A6QKMEAkYIILxGgvC0BsLrCayRtpaFd28l4f0Lwmvkdka4BsIbYagcqcoQQHjjjhrhjTtfy8K7V43bVeC/vvJilbm+hvKG1xdZh74IrwM0SiBghADCayQIT2sgvJ7AGmlrWXj3rDFShdIbKy9RmetrKMLri6xDX4TXARolEDBCAOE1EoSnNRBeT2CNtLUsvHvU1BHeN1cgvEauZ3xrILzxZcqJqg4BhDfurBHeuPO1LLy71xyhAv9vKy5VmetrKG94fZF16IvwOkCjBAJGCCC8RoLwtAbC6wmskbaWhbeFkvC+jfAauZ0RroHwRhgqR6oyBBDeuKNGeOPO17Lw/rrmbSrw/77iMpW5vobyhtcXWYe+CK8DNEogYIQAwmskCE9rILyewBppa1l4d605XIXSP1dcrjLX11CE1xdZh74IrwM0SiBghADCayQIT2sgvJ7AGmlrWXh3URLedxBeI7czwjUQ3ghD5UhVhgDCG3fUCG/c+VoW3p1r3qoC/90VHVTm+hrKG15fZB36IrwO0CiBgBECCK+RIDytgfB6AmukrWXh3anmMBVK763oqDLX11CE1xdZh74IrwM0SiBghADCayQIT2sgvJ7AGmlrWXiblugI7wc/IrxGrmd8ayC88WXKiaoOAYQ37qwR3rjztSy8O5TcogL/3z9eoTLX11De8Poi69AX4XWARgkEjBBAeI0E4WkNhNcTWCNtLQvvdiVDVSh99GMnlbm+hiK8vsg69EV4HaBRAgEjBBBeI0F4WgPh9QTWSFvLwrutkvB+jPAauZ0RroHwRhgqR6oyBBDeuKNGeOPO17Lwbl0yRAX+f37srDLX11De8Poi69AX4XWARgkEjBBAeI0E4WkNhNcTWCNtLQvvVkrC+ynCa+R2RrgGwhthqBypyhBAeOOOGuGNO1/LwrtFyWAV+J//2EVlrq+hvOH1RdahL8LrAI0SCBghgPAaCcLTGgivJ7BG2loW3s1rDVKh9MXyripzfQ1FeH2RdeiL8DpAowQCRgggvEaC8LQGwusJrJG2loV3MyXhnY/wGrmdEa6B8EYYKkeqMgQQ3rijRnjjztey8DapNVAF/oLl3VLNXbV6tVSvVl2qV6+Wqk9WxbzhzYpkBn0Q3gwg0gICSgQQXiXwBRqL8BYItNIYy8LbuNbNKlS+Wn6l89wflq+Q0y/qI+3OPE6OPXxf5z5ZFiK8WdJM2QvhTQmQcggoEkB4FeEXYDTCWwDIiiMsC+8mSsK7yFF4B426X+6874l8mgN6XYTwKt5rs6MRXrPRsBgENkgA4d0goqAfQHiDjm+Dy1sW3ka1Bmxwfx8PLF7e3antkqXfyvIVK6TNJddL53anIbxOFCMvQngjD5jjRU0A4Y06XkF4487XsvA2rNVfBf6S5T1SzW3Vuptcft7JCG8qipEWI7yRBsuxqgQBhDfumBHeuPO1LLwbKQnvN+WEd/qsF+TLhV9XeBGa77yd7N+yRZm/Q3jj/s9MqtMhvKnwUQwBVQIIryp+78MRXu+IVQdYFl5VML8Yfs+DT8mnXyyscJ29dttJDj9wb4TXSljW90B4rSfEfhBYPwGEN+7bgfDGnS/Cm32+vOHNnmk0HRHeaKLkIFWQAMIbd+gIb9z5IrzZ5Zv7/t21a9bKsWdfJe3PPl6O/cO+UqNGcXYDHDvxtWSO4HyUIbw+qNITAoUhgPAWhrPWFIRXi3xh5iK82XHu3GekzJz7SpmGj028SbbfZvPshjh0QngdoP2yZM2atfL1km/y/+2lQf26qbohvKnwUQwBVQIIryp+78MRXu+IVQcgvKr4CzIc4U2B+cW//F069B4u3/+wPN+l5W92ka4Xny4tmm1fYdfZ816XDr1vXefvXp81Rkpq1hCEN0UYlEJAmQDCqxyA5/EIr2fAyu0RXuUACjAe4U0B+aXX/yELv1oiB+67hyxfvkKuG3qX5N743t6/U4Vdn5r3mlx14xiZNqZvmb/fZssmUq1aNYQ3RRaUQkCbAMKrnYDf+QivX77a3RFe7QT8z0d4M2Sc+466HjeOljdnj5PioqJ1OueEt+/gCTLv4eEVTuUNb4Zh0AoCBSaA8BYYeIHHIbwFBl7gcQhvgYErjEN4M4Sek933P/xsnTe4P4/ICW/H3sPlhFb7S0lJTdl7j2bS6uCWpXKM8GYYBq0gUGACCG+BgRd4HMJbYOAFHofwFhi4wjiENyPoP7/dHTuom+y7968r7PrWOx/mf3Ix98Ntn89fJFMefVranHSY9Op4Vv55hDejMGgDAQUCCK8C9AKORHgLCFthFMKrAL3AIxHeDIA//+rb0q7bILm28zly2vGHJO744IxnpffN49f7EYjEjXgQAhCAAAQgAAEIQGC9BBDelJcj98Y2951zN3Q/X0466v8q1W3ey29J++6D5bWZo6VWSc1K1fIwBCAAAQhAAAIQgEAyAghvMk4VPvXIzOel501jpMdlbeTQA/YqfaZRg3pSp3YtmTDlScl9Fdmk4T3zfzf5odnSrOnW0nzn7WTpsm+l23WjpEZxkYwf2j3FFpRCAAIQgAAEIAABCPwvAghvivtx3dCJcv8jc9bp8PPb3oEj75Mp0+fKq0+Myj8z5I4pMu7eGaXP7968qQzs3V622nzTFFuEWbrs2+8l9+sHGzWoH+YB2Pp/Eli7dq2sXrOmwm8rAV3YBFauWi1fLVoiGzfaKP/94fyJi0Du/y5/9fXS/K+GbdK4kRQVVY/rgJymyhJAeAsc/fIfV8jCRUukft060rBBvQJP1x+X+yUd3W+4Q+Y8/0Z+mZz0D7+hgzTeuIH+cmyQGYHcD3EOHTNV5kwdmllPGukTGHPPY3LLmGmli+S+Zebazm2lwUbpfsuk/snYIEcg9wIn9yLn5z+bbdpIbr2hw3p/mRLUIBASAYQ3pLQi2HXs5Mdl6vS5Mml4L6ldq6Zc3GNo/vdrX3/leRGcjiN88tl8ubDrIPn0i4WS+3+WCG9cd2LqY3Nl6y2ayB7Nd5T/fL5Azu88QM5vfYy0Pf3IuA5aRU+T+y+quRcxv929Wf5/gevad6SsWrWaj91V0fsQ27ER3tgSNX6eUy68Nv/dwxeecWx+059/6O/tp+/M/7Y5/oRN4Of/OXTOc2/I2MmPIbxhx7nB7XPfMvPZFwsRog2SCvOBrtfdnv/toUP6XBLmAdgaAr8ggPByHQpKoOVR7fPfaJGT3tyff7z7kZzaro+8MH1E/vuJ+RMHgSfmvCwDb78P4Y0jzgpPkfssb6vWXeWYw/aVLu1Pi/ikVe9oj856XnL/pfXdf/9HhvS5VHbZcZuqB4ETR0cA4Y0uUrsHyv0gU4tDzpWRN3WSg/bdI7/oBx99Jse37SVP3T9YNt9sE7vLs1mlCCC8lcIV5MPXDrpTZsx+WR6f1F+aNG4Y5BlYumICuc9pv/a3d2XBV4vl+ivPl3323AVUEAieAMIbfIRhHSD3hrdfjwvkiIP25g1vWNFValuEt1K4gnt45ISHZcSEh+W+UdfKbrtsH9z+LJyMwB2TpsvdD8ySeQ8PT1bAUxAwTADhNRxOjKvlPsN75CH7yAVtjskfj8/wxpiyCMIbZ665z3MOHnV//usW7xrWI/+d4vyJl8CsZ/4ina69jd8GGm/EVepkCG+Vilv/sLmvNZr22DP5b2moU7tE2ncfwrc06MeS2Qa5j63kfqr7yadfyX8t2czJA6Va9Wp8H29mhHUbXT1gnDz0xDwZNaCL7LDt5qXL5L6Ro7ioSHc5pqcmkHtzv/8+u+V/QdKixd9I7ofWapfU5IcSU5OlgQUCCK+FFKrQDt99vzz/f0SffenN/KlbNNtehvfryGcAI7kD73/4mZxwbq8ypznuiP2kf892kZywah+jVetu+a+cK/9nxt0DZNutNqvacCI4fa/+Y+XhJ58rPcmeLXaS/r3aVclfjhRBnByhHAGElyuhQmDpsu9k5cpV/MIJFfoMhQAEIFAxgRUrVsqCRUukXp3aVfKXI3Ev4iWA8MabLSeDAAQgAAEIQAACEBARhJdrAAEIQAACEIAABCAQNQGEN+p4ORwEIAABCEAAAhCAAMLLHYAABCAAAQhAAAIQiJoAwht1vBwOAhCAAAQgAAEIQADh5Q5AAAIQgAAEIAABCERNAOGNOl4OBwEIQAACEIAABCCA8HIHIAABCEAAAhCAAASiJoDwRh0vh4MABCAAAQhAAAIQQHi5AxCAAAQgAAEIQAACURNAeKOOl8NBAAIQgAAEIAABCCC83AEIQAACEIAABCAAgagJILxRx8vhIAABCEAAAhCAAAQQXu4ABCAAAQhAAAIQgEDUBBDeqOPlcBCAAAQgAAEIQAACCC93AAIQgAAEIAABCEAgagIIb9TxcjgIQAACEIAABCAAAYSXOwABCEAAAhCAAAQgEDUBhDfqeDkcBCAAAQhAAAIQgADCyx2AAAQgAAEIQAACEIiaAMIbdbwcDgIQgAAEIAABCEAA4eUOQAACEIAABCAAAQhETQDhjTpeDgcBCEAAAhCAAAQggPByByAAAQhAAAIQgAAEoiaA8EYdL4eDAAQgAAEIQAACEEB4uQMQgAAEIAABCEAAAlETQHijjpfDQQACEIAABCAAAQggvNwBCEAAAhCAAAQgAIGoCSC8UcfL4SAAAQhAAAIQgAAEEF7uAAQgAAEIQAACEIBA1AQQ3qjj5XAQgAAEIAABCEAAAggvdwACEIAABCAAAQhAIGoCCG/U8XI4CEAAAhCAAAQgAAGElzsAAQhAAAIQgAAEIBA1AYQ36ng5HAQgAAEIQAACEIAAwssdgAAEIAABCEAAAhCImgDCG3W8HA4CEIAABCAAAQhAAOHlDkAAAhCAAAQgAAEIRE0A4Y06Xg4HAQhAAAIQgAAEIIDwcgcgAAEIQAACEIAABKImgPBGHS+HgwAEIAABCEAAAhBAeLkDEIAABCAAAQhAAAJRE0B4o46Xw0EAAhCAAAQgAAEIILzcAQhAAAIQgAAEIACBqAkgvFHHy+EgAAEIQAACEIAABBBe7gAEIAABCEAAAhCAQNQEEN6o4+VwEIAABCAAAQhAAAIIL3cAAhCAAAQgAAEIQCBqAghv1PFyOAhAAAIQgAAEIAABhJc7AAEIQAACEIAABCAQNQGEN+p4ORwEIAABCEAAAhCAAMLLHYAABCAAAQhAAAIQiJoAwht1vBwOAhCAAAQgAAEIQADh5Q5AAAIQgAAEIAABCERNAOGNOl4OBwEIQAACEIAABCCA8HIHIAABCEAAAhCAAASiJoDwRh0vh4MABCAAAQhAAAIQQHi5AxCAAAQgAAEIQAACURNAeKOOl8NBAAIQgAAEIAABCCC83AEIQAACEIAABCAAgagJILxRx8vhIAABCEAAAhCAAAQQXu4ABCAAAQhAAAIQgEDUBBDeqOPlcBCAAAQgAAEIQAACCC93AAIQgAAEIAABCEAgagIIb9TxcjgIQAACEIAABCAAAYSXOwABCEAAAhCAAAQgEDWB/wcAKlm68/HwywAAAABJRU5ErkJggg=="
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "mean = nested_bayes_tree.optimize()\n",
    "px.imshow(mean.vector().reshape((M,N)), zmin=-1, zmax=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2e33aca3",
   "metadata": {},
   "source": [
    "And of course, now ancestral sampling can also be done in parallel (here in the Bayes net):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "7e0f2cc5",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[-0.17 -0.09 -0.27]] [[-0.05 -0.48 -0.49]] [[0.15 0.13 0.17]] [[ 0.2  -0.16  0.41]]\n",
      "[[-0.85  0.16  0.1 ]] [[-0.28 -0.38 -0.46]] [[ 0.17 -0.1   0.1 ]] [[ 0.35 -0.46  0.58]]\n",
      "[[-0.54 -0.44  0.38]] [[-0.05 -0.27  0.16]] [[ 0.68 -0.07  0.21]] [[ 0.38 -0.04  0.05]]\n"
     ]
    }
   ],
   "source": [
    "gbn = graph.eliminateSequential(nested)\n",
    "samples = sample_bayes_net(gbn, 3)\n",
    "for row in range(M):\n",
    "    print(\" \".join([str(np.round(samples[keys[(row,col)]],2)) for col in range(N)]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Discussion\n",
    "\n",
    "Direct (multi-frontal) elimination), as shown above, works like a charm in Gaussian Markov random fields with arbitrary connectivity. However, this will *only* scale up to large graphs when small \"separators\" are to be found. In particular, in image applications, the smallest separator that can be found will always be the smallest image dimension, here $n=\\min(M,N)$, and hence the root clique of the Bayes tree will have that dimension as well. When images are large, this will be prohibitively expensive, as a dense $n \\times n$ matrix will have to be factorized, and we might need to resort to different methods. Still, $256\\times256$ images for example are well within reach for exact inference. \n",
    "\n",
    "In the next section we discuss loopy belief propagation and Gibbs sampling as alternatives."
   ]
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "9f7376ced4243bb13dfcffa8a3ba834e0602aa8334cd3a1d8ba8d285f4628083"
  },
  "kernelspec": {
   "display_name": "Python 3.8.12 ('gtbook')",
   "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.12"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}