{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# [Damped Oscillator - C++](https://github.com/Ziaeemehr/vbi_paper/blob/main/docs/examples/do_cpp.ipynb)\n",
    "\n",
    "<a href=\"https://colab.research.google.com/github/Ziaeemehr/vbi_paper/blob/main/docs/examples/do_cpp.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import torch\n",
    "import pickle\n",
    "import numpy as np\n",
    "import sbi.utils as utils\n",
    "import matplotlib.pyplot as plt\n",
    "from multiprocessing import Pool\n",
    "from sbi.analysis import pairplot\n",
    "from vbi.inference import Inference\n",
    "from sklearn.preprocessing import StandardScaler\n",
    "from vbi.models.cpp.damp_oscillator import DO_cpp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "from vbi import report_cfg\n",
    "from vbi import extract_features\n",
    "from vbi import get_features_by_domain, get_features_by_given_names"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "seed = 2\n",
    "np.random.seed(seed)\n",
    "torch.manual_seed(seed);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "parameters = {\n",
    "    \"a\": 0.1,\n",
    "    \"b\": 0.05,\n",
    "    \"dt\": 0.01,\n",
    "    \"t_start\": 0,\n",
    "    \"method\": \"rk4\",\n",
    "    \"t_end\": 100.0,\n",
    "    \"t_transition\": 20,\n",
    "    \"output\": \"output\",\n",
    "    \"initial_state\": [0.5, 1.0],\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Damp Oscillator model\n",
      "{'a': 0.1, 'b': 0.05, 'dt': 0.01, 't_start': 0, 'method': 'rk4', 't_end': 100.0, 't_transition': 20, 'output': 'output', 'initial_state': [0.5, 1.0]}\n"
     ]
    }
   ],
   "source": [
    "ode = DO_cpp(parameters)\n",
    "print(ode())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYYAAAEiCAYAAAD9DXUdAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjAsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvlHJYcgAAAAlwSFlzAAAPYQAAD2EBqD+naQAAWlhJREFUeJztnXecFPX9/1+z9W6vcv0O7jh6kXaAIKARFYNIsOUbG1GCxkQDCUo0EQtYopifhmgSI9FYY8ESwQLREBAQRfohHY52x3G9bbvtn98fn5nZnd3Zvd1ry8H7+Xisszc7M3zmwHntuwuMMQaCIAiCENHEewEEQRDE2QUJA0EQBKGAhIEgCIJQQMJAEARBKCBhIAiCIBSQMBAEQRAKSBgIgiAIBSQMBEEQhAJdvBfQ3fh8Ppw5cwYpKSkQBCHeyyEIgugWGGOwWCwoKCiARhPZJjjvhOHMmTMoLCyM9zIIgiDiQkVFBfr06RPxmPNOGFJSUgDwX05qamqcV0MQBNE9mM1mFBYWys/ASJx3wiC5j1JTU0kYCII474jGhU7BZ4IgCEIBCQNBEAShgISBIAiCUHDexRgIgjh/8Hq9cLvd8V5Gt6DX66HVajvlWiQMBEGcczDGUF1djebm5ngvpVtJT09HXl5eh2u0SBgIgjjnkEQhJycHJpPpnC9mZYzBbrejtrYWAJCfn9+h65EwdBH/3V+NL/ZX45GZw5GRZIj3cgjivMHr9cqikJmZGe/ldBuJiYkAgNraWuTk5HTIrUTC0AUwxvDLt3eCMaBvRhIWTBsU7yURxHmDFFMwmUxxXkn3I92z2+3ukDBQVlIXUNncCsb4++9PN8d1LQRxvnKuu4/U6Kx7JmHoAspqrfL7eqszjishCIKIHRKGLqDe6pLfn2lxxHElBEEQsUPC0AU02pwB711gkl+JIAiiDRhjWLZsGfr16weTyYTrrrsOLS0t3boGEoYuoNHmL6jx+hhsLm8cV0MQRE/igQcewEsvvYQ333wTX3/9NXbu3InHHnusW9dAwtAFBFoMANBsd4U5kiAIws/WrVuxbNkyvP/++/jBD36AcePG4a677sKaNWu6dR0kDF1Ao00pBC2t50dJPkEQHeO5557DFVdcgbFjx8r7cnNzUV9f363roDqGLoCEgSDOLhhjaHV3v0s3Ua+NOoXU6XRi9erVeO655xT7HQ4H0tLSumJ5YYmrMGzatAnPPvssdu7ciaqqKqxcuRLXXXddxHOcTieeeOIJvP3226iurkZ+fj4WL16MO+64o3sWHQUWh0fxs5mEgSDiSqvbi+GLv+z2P/fAE9NhMkT3mN21axdaW1vx29/+Fr/73e/k/W63G5dddllXLVGVuAqDzWbD6NGjcccdd+CGG26I6pwbb7wRNTU1ePXVVzFw4EBUVVXB5/N18Upjw+bkwpBu0qPZ7obNScFngiAic+TIESQlJaG0tFSxf+bMmZgyZUq3riWuwjBjxgzMmDEj6uO/+OILbNy4EcePH0dGRgYAoLi4uItW136kLKScFCOa7W7Y42DCEgThJ1GvxYEnpsflz40Ws9mMrKwsDBw4UN536tQpHD16FD/+8Y+7Ynlh6VExhk8//RTjx4/H//t//w//+te/kJSUhGuuuQZPPvmk3EAq3jDGZIshO8WIIzVW2J2eNs4iCKIrEQQhapdOvMjKykJLSwsYY3Jc4qmnnsLVV1+N4cOHd+tazu7fVBDHjx/H5s2bkZCQgJUrV6K+vh6/+tWv0NDQgNdff131HKfTCafTnz5qNpu7dI1Ojw8eHy9oy042AgDsVMdAEEQbXH755XA4HHjmmWdw880345133sFnn32Gbdu2dftaelS6qs/ngyAIeOeddzBhwgRcffXVWLZsGd588020traqnrN06VKkpaXJr8LCwi5doy3AOsiShYEsBoIgIpObm4s33ngDL730Ei644AJ899132Lx5c5c/s9ToUcKQn5+P3r17K1K3hg0bBsYYTp8+rXrOokWL0NLSIr8qKiq6dI1SoDlRr0VyAjfIyGIgCCIabrrpJpSXl8Nut+Pzzz/HgAED4rKOHiUMU6ZMwZkzZ2C1+ruXHjlyBBqNBn369FE9x2g0IjU1VfHqSqyixZBk1CHJQMJAEETPI67CYLVaUVpaKqdnnThxAqWlpSgvLwfAv+3ffvvt8vG33norMjMzMXfuXBw4cACbNm3CAw88gDvuuOOsCT7bRLdRslGLRAPPSCBXEkEQPYm4CsOOHTtQUlKCkpISAMDChQtRUlKCxYsXAwCqqqpkkQCA5ORkrF27Fs3NzRg/fjxmz56NWbNm4S9/+Utc1q+GwmIwSsJAFgNBED2HuGYlTZ06NWJL6jfeeCNk39ChQ7F27douXFXHsAUIQ6Jep9hHEATRE+hRMYaegCQCyWQxEATRQyFh6GSkrCSTQStXPcajeRdBEER7IWHoZAIthgRRGJzus6uXE0EQRCRIGDoZqU+SyaCDUcd/vU4PWQwEQfQcSBg6GYfoNkrQa2SLwUEWA0EQPQgShk5Gsg4S9VrZYnB4vBGzrwiCIM4mSBg6Gck6SNBrYRQtBsYAt5eEgSCIngEJQyfT6vK7kiSLAeBWA0EQRFvcc889uPjii1U/69OnD5555pkuX0OParvdE5AEwBjgSgLEzKSEeK2KIIiewP79+/Hyyy/j66+/Vv182LBhIRPeugIShk5GCj5LQ8CNOg2cHp+8nyCIOMAY4LZ3/5+rNwHi0J1oePbZZ3HhhRdi8uTJqp9nZGSgurq6s1YXFhKGTiYwxiBtnR4fnB7KTCKIuOG2A08XdP+f+9AZwJAU1aEejwcff/wxHn30UXnfL3/5S0yYMAF33nknAMBisXRLw1CKMXQygemqgVuyGAiCiMSxY8dgsVgwcuRIAHww2YcffoiUlBT5mO+//75bxnySxdDJ+IWBWwxGnVj9TBYDQcQPvYl/e4/Hnxslzc3NAHgXaQD48ssv0dTUhIQEHpz87rvvUFlZieuvv77TlxkMCUMnI7uSdJIrSax+JouBIOKHIETt0okXffv2hSAIeO+995CUlIT7778fM2fOxCeffILCwkLcfffdmDZtWtiMpc6EXEmdjJSVlGjgv1qyGAiCiIa8vDw89dRTePvttzFjxgz89re/xVNPPYV169bhkksuwbBhw/DBBx90y1rIYuhkpDoGY5DFQDEGgiDaYtGiRVi0aJFi38mTJ7t9HWQxdCKMMdkyoBgDQRA9FRKGTiTw4U9ZSQRB9FRIGDqRwIc/WQwEQfRUSBg6ESkjSacRoNeKwWeyGAiC6GGQMHQirUE1DABZDARB9DxIGDqR4KrnwPdkMRBE93I+zkDprHuOqzBs2rQJs2bNQkFBAQRBwKpVq6I+95tvvoFOp8OYMWO6bH2xElz1DPgtBpriRhDdg16vBwDY7XFomhdnpHuWfgftJa51DDabDaNHj8Ydd9yBG264Ierzmpubcfvtt+OKK65ATU1NF64wNoIb6PH3XHtdXrIYCKI70Gq1SE9PR21tLQDAZDJBiKHDaU+EMQa73Y7a2lqkp6dDq9W2fVIE4ioMM2bMwIwZM2I+7+6778att94KrVYbk5XR1ai5kuQYg8sLuFsBfWydERljmPfuLhyssuCNuReib+bZXdZPEGcDeXl5ACCLw/lCenq6fO8docdVPr/++us4fvw43n77bfzhD3+I93IUyMKgC3QlcZG4peIx4JlvgTv+A/QeF/U1d1c0Y81e3n/9jW9PYsmsCzpvwQRxjiIIAvLz85GTkwO32x3v5XQLer2+w5aCRI8ShqNHj+LBBx/E119/DZ0uuqU7nU44nU75Z7PZ3FXLC+iTFCAMeg3SYMVYy1d8x+53YhKGXaea5Pc7TjZFOJIgiGC0Wm2nPSzPJ3pMVpLX68Wtt96Kxx9/HIMHD476vKVLlyItLU1+FRYWdtkapRiDUacMPg8VKvwH1R+J6ZoHzviF7Hid9bzMtCAIonvpMcJgsViwY8cOzJ8/HzqdDjqdDk888QT27NkDnU6H9evXq563aNEitLS0yK+KigrV4zoDqYGeMsagQa4Q8E2/5XRM1zzd3Cq/t7m8qLe6OrZIgiCINugxrqTU1FTs3btXse/vf/871q9fj48++gj9+vVTPc9oNMJoNHbHEmVXkjJdVYOcQGEwVwI+H6CJTpMrm1oVP1c02ZGd0j33QxDE+UlchcFqtaKsrEz++cSJEygtLUVGRgaKioqwaNEiVFZW4q233oJGo8GIESMU5+fk5CAhISFkf7zwp6sGWAx6rdJi8LqA1kYgKavN63l9DNVmBwAgLzUB1WYH6izONs4iCILoGHF1Je3YsQMlJSUoKSkBACxcuBAlJSVYvHgxAKCqqgrl5eXxXGJMSFPaEoMshkwhKOBtb4zqerUWB7w+Bp1GwIjeqQBAwkAQRJcTV4th6tSpEYOpb7zxRsTzH3vsMTz22GOdu6gOoF75rEEKgiow7Q1RXa+6hVsLuakJyE3lc19rSRgIguhiekzwuScQroleiqCME0QrDE12HmjOSDLIcQWyGAiC6GpIGDoRf7pqYIwhwGIQRMFojc6V1GjjhTm9kgzISiZhIAiieyBh6EQkV1JggZtBq0EyRIuhVzHfRmsx2ESLwaRHRpIBANDSSumqBEF0LSQMnYhDmvesU1Y+pwjcYmAxCkOj6ErqlWRAWiLvlthsPz/K+wmCiB8kDJ2Iw6USYwiwGHypffhOpyWq6/ktBr8wtLSSMBAE0bWQMHQi/gK3gBgDXDAIfL8nOV88MLp+TY02FYuBhIEgiC6GhKETcajVMXhtAAAfE+Ayie1wo7UYArKS0k1cGFweH02DIwiiSyFh6ETkrKQAYRBEEbAiAS59Ct8ZpTDIFoPJgGSjDloNHzZCcQaCILoSEoZOpFVlUA9cXARsSIRLY+L7nNG5kppEAchIMkAQhAB3EmUmEQTRdZAwdCJqlc9w8+rlVmaAQ5vM90VhMfh8DM1SVpLoRpID0GQxEATRhfSY7qo9AafoSgqMMcDDhcEJPaARx3JGYTHYXB74xG4hqYlKYaAANEEQXQkJQyfh9TG4vFJ31UBh4JXKTujhEyRXkgVgDIgwoNzi8AAA9FpBrqSmlFWCILoDciV1EoGZQooYg4fXMDhhgF2KMTAf4LJFvJ7VyYUhJUEPQRSQlASu4zbxM4IgiK6AhKGTUAiDTsViYHq0MqO/X1IbcQaLg1sFyUa/UScJg2RNEARBdAUkDJ2E1A7DoNNAowlwEckxBgOcXh9glFJWI8cZpIe/JAaAXySsZDEQBNGFkDB0EnJGki7oVypmJTmg58FpIx+407bFwB/+gRZDslGv+Ky966QCOYIgIkHC0EmopqoCSovBE73FEBhjkEhO6JjFcKzOiouWrsPFf1yPikZ72ycQBHFeQsLQSYQXBn+MwenxAgYxZdUV+cEsxRhSA1xJKZIrydG+rKSXNhxDs92NeqsLr31zol3XIAji3IeEoZOQ2mEoMpIARR2D0+MLEIbIWUmyKykwxtABi4Exhq+P1sk/f3WoNuZrEARxfkDC0Em05UpywMBjDLIwWCNeL1LwuT0xhvJGO2rM/ulvJxvsVEFNEIQqJAydhGwx6MLEGJgeLq83dovB2DkxhoNVPNg9sncaijJ4PcW+My0xX4cgiHMfEoZOQrIYjCGuJH/ls9JiaKvAjX+bT1GLMbRDGI7VcQtlQHYShuSlKPYRBEEEQi0xOgn/kJ4gi8HNK58dMMDr8QEmURjc0VkMCmEQM5SsDg8YY3JFdDQcr+N/Xv/sZDmwfbKeMpMIggglrhbDpk2bMGvWLBQUFEAQBKxatSri8R9//DGuvPJKZGdnIzU1FZMmTcKXX37ZPYttA3/wOUxWEqSsJLHDapSupBSV4LPHx3ggOwaO13ProH92EvpmcnE62RB5DQRBnJ/EVRhsNhtGjx6NF198MarjN23ahCuvvBJr1qzBzp07cdlll2HWrFnYvXt3F6+0bcIWuAXEGGLJSlKrYzDptXLfPXOMKaunGrh10C8rCcUkDARBRCCurqQZM2ZgxowZUR///PPPK35++umn8cknn+Czzz5DSUlJJ68uNpxt1THAAG8MMQa1XkkajYBkgw4WpwdWhwc5KdGtze31ydPg8lITkCqKzemm1phdUgRBnPv06BiDz+eDxWJBRkZG2GOcTiecTn+aptkc3fS0WJF6JYXWMUgxBj18CldS7OmqAHcnWZyemALQDVYuClqNwMeEJvC1ujw+NNndyEgyRH0tgiDOfXp0VtJzzz0Hq9WKG2+8MewxS5cuRVpamvwqLCzskrW0WfkstcSIwmJweXxyDCElIF0V4BZEsVCF9G1/BsxnolpbnYWvISvZAI1GgFGnRVYyF4OqltaorkEQxPlDjxWGd999F48//jg++OAD5OTkhD1u0aJFaGlpkV8VFRVdsp42eyXFEGMItAaSgy0GoxYv6V9A0ffPA58tiGptdVa+hqxko7wvNzUBAFDd4ojqGgRBnD/0SFfSihUr8POf/xwffvghpk2bFvFYo9EIo9EY8ZjOQMpKMobtrmoAi9KVJMUXkgxaaDVK//8AbS2Gacr5D0f/ywVGEpsw1Fu4Kyk7xf97yE9LwP4zZlSbSRgIglDS4yyG9957D3PnzsV7772HmTNnxns5Mm13Vw0ucAtfQ6DWJ0liME4pd5wpbXNtdVbuSsoOsBjy0jpmMZTVWjDv3V148asy+KTh1ARBnBPE1WKwWq0oKyuTfz5x4gRKS0uRkZGBoqIiLFq0CJWVlXjrrbcAcPfRnDlz8MILL2DixImorq4GACQmJiItLS0u9yDhDz6Hr2NgUbqS/IFnfchnxb5y5Y76w0DxlIhrk2IMgRZDXgdcSV4fwy/e2onj9Tas/r4K2SlG3Di+a2I3BEF0P3G1GHbs2IGSkhI51XThwoUoKSnB4sWLAQBVVVUoL/c/CF9++WV4PB7MmzcP+fn58mvBguh87V2J32KIVMfgBfSiMHhaAZ/6wBy1VFWJPE9QwLnhWJtrkyyGLIXFkAgA7XIlbThci+P1fmF745uTMV+DIIizl7haDFOnTgVj4d0Qb7zxhuLnDRs2dO2COoBUx5AYaDH4vICPP+QdMACBFgPArYaE1JBr+YvbQv960r31AIBTphHoa98HNB5vc21qFkOO+F76LBbWiS27rx1TgM+/r8KBKjMqGu0oFJvzEQTRs+lxMYazFdWWGB7/t3E5xqAzAoJ4TBh3kuRKSlVxJaW6GwAAZQkX8B3NbWdZ1asIQ3Y7hYExhk1H+FyHa8cUYExhOgDgm7L6mK5DEMTZCwlDJ+FvohfwK/X4H7q8jsELBrTZLymSK8nk4g/gI9pB4sFt1zL46xj8wiC9b7S74PFG33ep2uzA6aZWaDUCJvbLxJSBWQCArScao74GQRBnNyQMnYTcdjtwHoPYWZVp9PBBAx/jDfDaGtZjCedKctlg8PC5CvswkO+zNygESG1d0vUCLYaMJAM0AsAY5HYZ0bC/kleOD8xORpJRhzGFPOi/t5JmOxDEuQIJQycR0ZWk8z+QoylyC5uuauFZWHZmxAl3JqAVr2upCrsuyVow6DSK+dFajYBM0WqojcGddKCKC8MFBTw2MqI3F4ZjdVbY2jEngiCIsw8Shk5CNStJ+iavS5B3Od0BU9zc6rUM1nDpqjbu269jabC6vEBqPt9vjiAMATUMwc3yJHdSvTV6YdgvTn0bLgpDTkoC8lITwJhfNGLF5fHh+9PNaHWpZ2kRBNG9kDB0Aoz55yOoWQyCLgEGLf9Vc4shcvWzFGNICY4xOPiD1wwTz1xKKRBPCB9nkALPWSmh1d/tCUBLD//h+f5sKsl6ONQOYWixu/Gjv36Na/72DaYt24jKZurdRBDxhoShEwgcmqPqStInyK0yonElhU1XdfIHr4WZYHG4wVJyxRNqw65NrepZQmqkVxelxeBwe1HRyB/c0nhQABiYy4XuaG3so0Kf+eIQjtTw8yqbW7Hkk30xX4MgiM6FhKETkNxIQNCgHjnGkCDPguZT3KKLMYS4khzcjWOBCW4vgzeRZwTBFj5VVK2GQULaJ/VSagtp2E9Kgk7RqnuQOBiiLEZhaLA68e9dpwEAS28YCa1GwP8O1uJQdde0RicIIjpIGDoBKfCs0wjQadViDEY5W0nZLymcKylM8FmyGMALyZyGXny/GHtQo162GEJnLkhWRLQWgzTxrTgzSRGvGJjTPovhk9IzcHl8GNk7DTdfWIhpw3iX3A+2n47pOgRBdC4kDJ1A2AZ6YroqdIkxuZLkGEOwMIgxhlYNP9+hF4XB3hB2bdFYDHWW6NpinJKEIUvZzVUShjqLEy326EeOrjtUAwC4rqQ3BEHAT8bxfktr9lZFrIiPhMvjQ43Z0e7zCYIgYegUVIvbAIXFYNCpuJKcod+wGWP+GENw8Fm0GFxa/iC2ScLQXleSnJUUnSvpRD13JRVnKltfJBt1yBe7tZbVWaK6lsXhxtbjvCjuiqHcUrh4UBZMBi2qzQ7sPxO7O2nz0XpMWroOE59eh+v//i1qqaU4QbQLEoZOQEqzVBS3AUExBhVXkjvUYrC7vJC6WIfGGPjD0q3nwmDWpPP9EVxJag30JGLNSjoV4EoKRrIaoo0zfFNWD4+PoX92kmyBJOi1uFispF5/KHxAPdzafvGvHWgQi/VKK5pxx5vb4Y6hqpsgCA4JQyfgL24LbzEoXUliRo+KxSDFF7QaIfR6osXg0fPzzRoxZdSubjEwxlSH9EhIYtHS6uaWTBtIwefirNBmeVIA+mhNdMKw/WQTAGDygEzF/ivEOMO6GIXhic8OwO7yYmK/DPxv4Q+QbtJjX6UZr20+EdN1AvH5WEztQgjiXCFmYfjqq6/CfvaPf/yjQ4vpqfhdScEWgxhj0AfGGCJnJVmd/vhCcEGaZDEwI38INwniDApHC+AN9e3bXF60ivEPNYshLVGPYdrT+I32YzSfidy+2+H24ow4H7pvJIuhLjph2HmKC8P4vhmK/ZcO5sLw/elmNEXZqmP/mRasO1QLjcCzmwbmpOChGcMAAMs3Hou5Itvt9eG5Lw9jzBP/xZBHv8DP39xO9RXEeUXMwnDVVVfhgQcegNvtfxDV19dj1qxZePDBBzt1cT0FZ9jpbSpZSR4fYAzfRM/sCN9yW7IYmIFbCg1ek79Tq0oAWnIRJRm0SFJpyKdx2/Ce/kks1H+EtI9vBXzhvx1XNNrBGI97ZCaFZjjJmUlRWAwOt1euoB7Xt5fis7y0BAzNSwFjwOYoO7a++jW3Cn40qgD9s/k6bhjbG/2yktBkd+Nf352KdLoCr4/h3hWl+NtXZTA7PPD6GP53sBbX/u0bnKwPP1wpEowxVLc4UEtBcaKH0C6LYeXKlbjwwgtx4MABrF69GiNGjIDZbEZpaWkXLPHsJ7wrSaWOIbAlhis0UCunqhpDW25LFoMgznCwOH2ASfzGrRJnkFJV1aqeAQCHPkc6+BoSmo8Cp75RPw7ASdGN1DfLFGrJABgkCkNlcyvsrsjf0PdWtsDtZchOMaJPr8SQz38wOBsAsPFI+NiJhM3pwX/28R5ScyYXy/t1Wg3mXcYbDb66+YSi1iQS7249hdV7q2DQavDCzWOw5jeXYGheCuqtTvzs9W0wO6LPugKAT0orcdlzG3DR0nWY8PQ6TFu2EZ+UVpJAEGc1MQvD5MmTUVpaihEjRmDs2LG4/vrrcd9992HDhg3o27dvV6zxrEdOVw0OPrsDhEF0Jbm8ATEGNVdSFBaDNjGdH+v0AEn8IaqWmSRnJKm4kQAAxzcqfz76X/XjAPnbspobCQB6JRlkS+J4XeRv1pIbaVxRL1WRuTRAGNp6gH65vxqtbi+KM00YW5Su+Oya0QXIT0tAncWJlbsrI14H4EL67JeHAQAPXT0U147pjeEFqXjrzgnonZ6Ikw12LPlkf5vXAbjlsejj77FgRSlONtih1QjQCMCxOhsWrCjFr9/bHVNvqDqLE0v/cxBX/GkDBj28BiMf+xKz//kd3ttWHrXoBdLq8mJPRTM2H63HjpONMfXLIs592jXB7ciRI9ixYwf69OmDM2fO4PDhw7Db7UhKUn9onOuErWMItBjUCtxUg89h+iQxBjj5t3utKQ1AMxcRkxi8jSQM4SyGM7sAAJ97L8KPtN8Bp7eHvUepuK1fGGEAgAE5yWg40YijtRa566oasjAEuZEkxhf3QqJeizqLEwerLHLDPjWkB75UCxGIQafBnRf3wx9WH8TLm47jxvGF0GpChUjij/85BLPDgwsKUnHbpGJ5f05KAv5yyxj8ZPkWrNxdialDsnHtmN5hr8MYwyOr9uK9bRXQCMBvrhiEuy7pDx9jeP2bk/jLuqP4/PsqnGqw45XbxyMvLSHstXw+hvd3VGDpmoOymxEA3F4PvilrwDdlDVi29gh++YP+mD2xLxIN2rDXsjk9WLO3Cit3V+K74w1y9pvEgOwk/GBwNi4fmoMJ/TJCs+wCcHl8OFRtxp6KZuw53YKDVWa4PD4YdBoUpCeisJcJA3KSMDA7GQNykpGZZFD9EuDzMdRZnTjVYEd5I39ZHR74GINBp0GG+IUjO8WInJQE5KQakWEyQBP098gYQ7PdjXqrEw02FywODywON2xODzQaAXqtBkadBiaDDiaDFibRvcrf861Rp4HHx2BzemB1emBzemF3eeD0+PjLzeepGHUaGHVaGPUa/3udBgadBgYt33p8DC6PDy6vj289Pnh8Pnh9DB4f41svAwODRhAgANCIXx4AvhUEcQsBggAIAqARBGgEAb2S9MhJCf/vpqPELAzPPPMMlixZgl/84hd49tlnUVZWhttuuw2jRo3C22+/jUmTJnXFOs9qHGKvJGPUWUmRgs/hZzGAcQEyJKUDaBYtBrEthkqMoT5CqiqcVqCOfzt+x3sFF4aqPXwcqSb0gSBlJPXNDD++c1BOMradaIyYssoYwy5RGMaGEQajTotJAzKx/lAtNh2tCysMtWaHPDnu+hL1B/UtE4rw1/VlOFFvw3/3V2PGyHzV43aeasSHO3nF9RPXjggRkHF9M/DrywfhhXVH8fDKfRhTmK5qPTHG8IfVB/HetgoIAvD8zSW4ZnSB/PlvrhiECf0ycM/bO7G3sgXX/G0zXrl9PEaLk/ACOVpjwUMr98oZXCN6p+KeSwdiVJ80WJ0ebDpShze/PYkzLQ78YfVBLN94DL/4QX/cOrGvPOSJMYb9Z8x4f3sFVu6ulP99AbxXVlayERaHB2daWnGszoZjdTa8/s1JJBm0uGRQNi7qnyHOB2eotThxuNqCfZUtOFhl4davCmo1KOkmPQZkJ6NPr0R4fQwWhwcVTXZUNrUqeo1Fg04jICvZCJNRC7fXB6fbh0abi8866QAaASFiebYyd0oxlsy6oMuuH7MwvPDCC1i1ahVmzJgBABgxYgS2bduGhx56CFOnToXTef6ZpFFZDIG9ksSsIrhtPOCr8QuKuY12GBC0SDBxf77F4QGyJGGI0WJoKAPA4DD0wlbHMDiEBCS47UD9ESBnWMjhJ+rVq54DiSYAXd5oR4PNBYNWgxG9w1sClw7OxvpDtdh4uA53XzpA9ZhPSs/Ax7jlEc7FlWTU4fZJffHX9WVYvvEYrhqRF/LN1etjeHQVdxHdOL5PWEvm15cPxLfH6rH9ZBN+/d5ufHT3ZLlwUeL5/x3Fq2KK7B9/PEohChIX9c/Ep/Mvxp1vbseRGitu/McWPDB9CH56UV8k6LWoNTvwytfH8fo3J+HxMZgMWvz2h0MwZ1JfRcuVYfmpmDulH1buPo2/fVWGisZWPL3mEJ798jCG5qXCqNPgZINd4SYqzjTh/8b1wTWje6MoQORbWt34tqweGw7XYf3hWtRZnPhifzW+2F+t+rsA+MN+VJ90jO6ThhG905CSoIPD7UVlUytO1NtxrM6KY3VWVDa3otnuxs5TTbK1GIhWI6AgPQFFGSYU9jKhV5IBWkGAw+1Fo82FOqsTdRb+ahAFoDpM8WJaoh6ZSQakJOqRmsAtAcZ4ppnT44Pdxa0Am5Nn7NmcHlmYAkXBoNMgxahDgl6LBD23Cgw6DQSBW0vcivDC4fZbBE6PN0RYDDoNjKIVodMK0Gk00GoE6DQCNBpuKTAAPsYAxrc+BjAwORfExxhYwGcAU53u2JnEfPW9e/ciS3oYiej1ejz77LP40Y9+1GkL60nIwedwBW76BGVWkiHgIea2+YUCge0w1IvbkJCKZPEzq9PtdyVFyEpSFYYm/vBypBTDZ9bglLYvhngOcysiSBicHn+qqlpxm4TcTC9Cyqr0YBjROzWiq0KKM+w41Qib06OaVRXoRorEnMnFeHnTcew53YItxxsweYDy3+/b353CgSozUhN0+P1VQ8NeR6fV4IWbS3D1X77G96db8NsP9+DPN46GTqsBYwwvrDuKF9YdBQA8Nms4bhxfGPZahRkm/PueyViwohTrD9XiD6sP4rn/HkZmkhFVLa3yA2basFw8fu0F6J0eGqQH+IPnpguLcMPYPli5uxLLNxzD8XqbYqJeol6Ly4ZmY/bEvpg8IFPVpZOWqMeMkfmYMTIfPh/DvjMtWH+oFgerzKgxO6HVCOhl0mNQbgqG5adiTJ90FGYkql4rmFaXF8frrThWZ0N1Syv0Wg2SjDr0SU9En14m5KcnQK+NLtzp9vpQb3Wi1uyEw+2FXnTfZCQZkJlsiPhvKhwerw92txetLi8M4tqCBT+Wa7m8Pmg1AgxaTVS/n7ORmIUhWBQCufTSSzu0mJ6K6pAeICjGIGUl+fjgHkHLXUMupTBY5aykMBaDMVW2JiyOAFeSWowhQsttNB4HAPjSi4FKoMxXgCE4DNQfDTm0orEVjPG01yyVZnwSksVwqsEu+5uDaSu+IFGclYSiDBPKG+3YcqwB04bnKj4/XG3BgSoz9FoBPwrjHpLISjbixvGF+Nd3p/DShmMKYahuccgB5wemD5Gn2oWjID0Rz980Bj9/cwc+23MGp5vsuHpEPr46XItvj3Fxvv+Hg/GzKf0iXgfg4v/P28fjw50V+Mu6MlQ2t8r1EuP69sL8ywfisiE5bV4HAPRaDW4cX4ifjOuDUw12HK21wunxIj8tsU0RDkajETCqTzpG9UmP+pxIJBq0uKAgDRcUhI87RYteq0F+WiLy09SFsj3otBqkajVIDf4y1s5r6aIUubOZuN7Bpk2bMGvWLBQUFEAQBKxatarNczZs2ICxY8fCaDRi4MCBeOONN7p8nW3hDFvgphZj8PIokjSsJygALaWrpoZpoIeEVPkzq9PjT1dVizFEGNKDRm4xGHJ4SucBl/gAaggVhsDmeZG+AeWmGpFs1MHrY3KwOphd5c0AgLFFkYUBUGYnBSNZC5cNyUEvlbqKYO66pD80AvD10XrsqeBr8PkYHlq5F1anB2MK03HrxOiy6qYOycHfZ49FkkGL3eXNeGrNQXx7rAEGrQZPXz8S8y8fFNV1AP4QvunCImz+/WX438JL8e97JmHbw1fg3/dMjloUAhEEAcVZSbhyeC5+NKoA4/r2ate3aOL8Jq7CYLPZMHr0aLz44otRHX/ixAnMnDkTl112GUpLS3Hvvffi5z//Ob788ssuXmlkwtYxBHZX1Qe4koCwrbel4GBojEF0DRjT5BoHnpWkHnxmjPktBlVX0kkAQGLOAGgEoMwnfutWsRjk+EIENxLAH0qReiZZnR4cFmcthAs8ByIJw4YjtYq0Va+P4dNSLgzhgs7BFGWaMEv09/9mxW7sqWjG4k/3Yf2hWhi0GnkeRLT88II8rF14KeZfNhBXj8zDPVMHYN1vL8WtE4uivkYg0u9uXN+MLs02IYho6NoIRhvMmDFDDmJHw/Lly9GvXz/86U9/AgAMGzYMmzdvxp///GdMnz69q5bZJuGDz2GykgBe/WxBSGaSP101UoxBdCU5PUCSGNwMciW1tLrh9vKHqar7p4lXA2sz+yM7xYLjFvE6DWU8NTbAMogmI0liYE4ySiuaVYVhT0UzfAzonZ6I3NS2H36TBmRCrxVQ0diKkw129BMD3+sP1eJMiwNpiXpcNjT6b9VLZl2AnaeacKrBjmtf9BfzLb1hJIblhw+Eh6MgPRH3Tx8S83kEcbbTo5xhW7ZswbRp0xT7pk+fji1btoQ9x+l0wmw2K17toukk8Pb/AR/cHvJR2AI3tRiD1KwujMVgaWOsJ4wpcvzB5fHBaUjn+1sbFS0tpEyU1ARdqCuBMcAqZpuk5CE3NQGnWC4YBP7nBInMyQhdVf336gJW3oPHyn6CW7TrcFBl/rOUploSVIgWjiSjDhcWc1fZJ6X+ArXXv+FusJsvLAwV4whkJBnwwS8nYeqQbOg0AvpnJeEft43Dj8f1ifoaBHE+0KOEobq6Grm5yiBkbm4uzGYzWlvVm5wtXboUaWlp8quwMHymSEQ8LqBsbWi1MPyupKjqGNySK0nql6QeYwhxJTkCgs8BgWmbNp2/YT7A0Szvr42UkdTaBHjFBnXJuchNTYALetiNYhV1c7nicH9X1QjCsP0VYM+7SHbW4A+612A+9X3IITvLxfqFKOILEjdP4K6Zt7acgs3pweaj9fj2WAO0GgG3TYq90r4gPRFvzJ2Ao0/NwPr7p2L6BXkxX4MgznV6lDC0h0WLFqGlpUV+VVRUtO9CetH14QnNn7aLFoPJEPQwV3RXDY4xhAs+c1dSSIaEWPWMhFRoNQJMYoWr1S0ARjHbIyDOEDFV1Sq2tE5IB/QJyE3lxzQZxDhD80n5UJfHh9NN6gN6ZBgDti6Xf9QKDD9q/UQxKMft9WHbCT6YZ0K/jJBLhOPqEXkoyjCh0ebC3W/vxO8+2gMAuO2ivujTq23XVjh6ahohQXQHPUoY8vLyUFNTo9hXU1OD1NRUJCaqp68ZjUakpqYqXu1CJ17f4wjpQupwScIQZVYSoFr97Pb6ZOsjUroq4Hc1mR1uICm0LYY0lS1bLZApuZGSufWVJ/r7azSiNRZgMZxussPHeC582NYa1Xv5ObpEYPZHAIAZ2m3Yc8q/ntKKZthdXvQy6TE8Bn++TqvBMz8eKWcUnWlxoH9WEu67cnDU1yAIIjZ6lDBMmjQJ69atU+xbu3Zt97Th0Ac8YIOsBrubu3+iq3wOCD4DCmGwBvTCCetKEjurSsJhdgQWufkfxLXiHGfVwLNFFNcULgRSILicia6kJn+b6sDAc9hv2WVr+XbA5cCAy2HVpiFdsKF+/wb5EKl1xeSBWSF9btpi8oAsvH3nRFw5PBe3TCjCil9ehLTEjuecEwShTlyzkqxWK8rKyuSfT5w4gdLSUmRkZKCoqAiLFi1CZWUl3nrrLQDA3Xffjb/97W/43e9+hzvuuAPr16/HBx98gNWrV3f9YnUBFonHARj8boxWF3/YKywGrwfwiQ/64CZ6QECMwd96W0pVTdBrQitBgyyGdJMBgA0tdrdqymqdmVsrqqmPssXA/etSVe1Rp+jiafYLg5Sq2i9SfOH0Dr4tvhjQaNGQfymST38K4eQmAD8FAPzvIBejSwaGL5CMxOSBWZjcznMJgoiNuFoMO3bsQElJCUpKSgAACxcuRElJCRYvXgwAqKqqQnm5363Rr18/rF69GmvXrsXo0aPxpz/9Cf/85z+7J1VVqwM0oo66lYHuVnH+QGKgxRBoVahmJYVaDOZw7TAAPqUNABJ4PKGXiR/T3KruSpJqGHIixRhEi6Ewg4vcXls63x/gSjpez2MgYQPPjPm7svYZDwDIvOAyAEA/+x5UtbTiRL0N+yrN0GoE/DDaYC9jEQcHEQTRdcTVYpg6dWrEfvtqVc1Tp07F7t27u3BVEdAl8m/4AQ99xpg8PlNhMUjxBUAUBn5OSIFbQPBZnsWg1iBL1WIAmuyuAFdSo3x4rWQxpKoIg0UZYyhIT4ROI+CEN4v/i2gul5v7Hanm6xucmxx6HYBbF7Y6QKMH8kbxyw6+FPgSGCOU4fUdx9Ho5KI4ZWAWMqKoUkbV98BHc4HmCuDie4GpixR1FQRBdC09KsYQd6Q4Q4DF4PT45IZnil74knhoDYBGExpjUKljsEQa0hMUY0gXfezNCldSaIxBPStJjDGIriStRkCfXomoYhlggpanslqrwRjDkVru6pIa5IVQxbOEkDfC//vJ6A+HMQtGwYNNX32Jf37N+zLNDZiwFhaXDVgxmxfaeZ3Axj8Cpe+0fR5BEJ0GCUMsBGYmiQRO4VJ1Jen4w1JyJfEBHT5/47xAYXByV1JI4JmxEItB6g/UHGgxiK4kl8eHJju/lmqMQbIYUvw1IUWZSfBCC1uC6OppLked1Ylmuxsawd8gLwRxpgOyAzqyCgIMfScAAIaxY/AxYMrATEwdkq1+jUB2vQW0lANpRcDEu/m+tUtC3HcEQXQdJAyxoGIxSG4kQ3BXRVkY+Df2wOrjcMN6/K6koBiDu9UfyBYtBikrp8nuDhnWI1U967WCHItQIMUYkv3+/qIMLnoNenFf0yl5rkJRhil8hbEsDMr0UU2fsQCAa7NrcNcl/fD32eParh1gDNj6D/7+koXAD//ABcJeD3z/fuRzCYLoNEgYYkEXWuRmd7XdchuAogU1Fwa14HNbQ3o08nm9xBgDz0pSzmSoDZj1HPIwdrf6G/Il+/sM9c3gQlUJcV/zKRyqFt1IuWHcSABQLwpDVlDPoAKeUDBKcxwPzxweXXpp43E+J0KjB0b+BNDqgYm/5J9t/2fb56tx+D/AP6cByy8Bdr7BxYcgiIiQMMSCXnQlBVgMjnBVz26lMGg1AvRa/pB2eryqlc9hx3pKVc/GFDkIm26SLAZXqDCIFcfZao3qJDeSLkHOcAKAIXn84X/IIaasNp3E96ebAQAjw81v9vmAejHdODtIGPK5MKDxGNDarH5+MMfW823RRf46j9G38Gyw6r1A7aHoriOxfxXw3s08a6r6e+CzBcCGpbFdgyDOQ0gYYkEfGmOQLIaQIexBFgMAZS2DSvC5pTVMOwyHMr4ABApDgCvJbQdcdoXFEILsRspVZPpIc5X3SCmrTadQKs4tUJtHzBdcztt+aA1AelDfoqRM/76qUvXzg5GEYcBlyusMFBsn7v0wuusAgLUO+PTX/P2Y2cClD/L3G//o/3NipbWJWyDH1vuFnyDOQUgYYkEXajFIMYbECC23JRStt+XK51BhCHG7yLMY/MIgu5JaXWD6JP5wBgB7vSwMqqmqQe0wJLKSjchNNaLCxwPE3qaTctXzmHCTvKTZDZkDeZ1HMKI7CWeiSC/2uoETm/j7AZcrPxv5E77d+2H0rqD1T3IXXP4Y4Jq/ApctAi68i3/26YLYH+w73wT+PIJbIP+6HnhhFHD4i9iuQRA9BBKGWFBppCcXt4WzGPT+imlJGBxur99i8Lp451YA5nDCEJSqCvgtBreXweb2KaqfpQZ6qsVtQe0wArmgIA0VjMcYNOZK6OHBoJxkpKkFsAF/4DkrzMSy3jwAjcpd6p8Hcno7F8nEDCBvtPKzIVcD+iReMxHNtSw1QOm7/P1VzwAa8e/myseBlAJu6QQ0/WuT75YDn/2Gr69XPx60t9Zwkdj5ZvTXkTj1LbBqHvDKFcBb1wGbn1fUoBBEvCFhiIUIFkNoAz1lVhIAmMTCNbsrIMYAAG4egJYshvTgB7Ez1JWUqNfKAe0mW2DKagOqW/j68tRiDEHtMAKZPCATdUiDUzBCAEOBUC9PUVMlXOBZQrIYonmYB7qRNEH/LA0mYOjV/P2+j9q+1s43AJ8b6DMB6BvQR8uQBFzxKH//9Z9U52SHcHIz8OUi/v6S3wK/3gUs2AOMnQOA8bjFwc/avg7ArZSPfwm8PgMofRuo3AEc/wr43xLghdGxB8dbm4CNzwKvzQD+Og545XJg9W+BsnXtqxq3VAMV2/jfxemd/PrEeUlcK597HCoWgz8rqe0YQ5IoHjanh2fcaI28iMtpBRJ7hXclqVgMgiAgw2RAtdmBBpsLhUn+AHRlMz+/dy+VjrPW8BbD1CHZ+MNqAae8WRisqUShUIdpw0OPk6k7wrfBgWeJghIAAmA+zR86KRHaYcjCcLn65yN+zF1J+z7maayaMOmzHhew41X+XspoCmTUzcB3L/Fg9MY/Alc/G35N9kbg33fxWRdjZgOXP8rjMpoEYNYLPEts5+vAv38O3P4pUDQx/LVs9bxwr+I7QNACJbN57MRaC+x4Hajdz0Vm37+Ba/4G9Iowa8JaC2x5Edj+qqLXFgCgcifP4EorAsbeDpT8FEjNV7+Ozwec2cXjJke+AGr2hR6TXgT0uRAonMi3eSP5v91AGOODrKpKgZoDQO0BoOEY/0LjcfJ/t6Ys7nLMGghkDuL3l94XSEz3X8fj5L8naw3/92Kt5tafpYpX1wsa7jI1pnBXaHKOuM3liRT6BP7/lM/Ns/2cVr4Gp5n/P+Q080QOpwWAwL+06RP5/6MGE7dK9Yn8C4TexK/n8/Ivgh4H33pd/L3HKW6lQlYjv47OwLcaHf89aXQ8y06r4+tnDAAL2CLgZwR9xoI+A7+GoOG/v5yA2qFOhoQhFtQshihabkskiRaDTXQ/wZAEtDrllNXwMYZQiwEAclONqDY7eBaSaDEwWx0qm3gWUUG6ijBIrqTk0Af+gOxkjO/bC+VncjAYlbgwtQUTw81OYCzAYgjTAtuYwv/x1h7gjfaG/Uj9OHuj36rof5n6MQMu5//zW6u5K6bfJerHHfiEP1hS8oHh14Z+rtFwYXnrGmDHa8CEX4R3hX1+H2A5wx9kVz+rbMshCMDVz/GH1pEvgHdvBOb+B8gdHnqduiPAuz/hD09jGnDTW0D/qf7Px9/BXVvrnuRxlpcmA1c+AYybq7SemiuALX/jloX0QModwQWwVz9e73FiExeXlnLgqz/wLKwhM4DBV/nvs/EEP+7YOv8XBYA/cNL6AIYUPvTJXMnbozSX82sC/P+BghIgOZsLpqUGqD+iGBIVgr2epyKf3hb6mS6BP+yZL2RoFRGBifcAM57pssuTMMSCWoGbK0zwWTomoCurlNIqWRkwJvORnC4rfD4WU4wBkGYttPBgsxhjcFrqYXP1AwAUpKlZDOFdSYIg4MnrRmDfPwsA727cPJiFL0qz1YuuBiH8gxUAeo/jwlAZQRhObATAgOyhQFpv9WN0RmDYNcDuf/GHVDhh2PoS346/M/SbrUT/S/mD8sgXwP8eA25Wabmx79/AgVX8G9+PX/HHhALR6oD/e40Hoyu28u3tq5Tf5Mr+B3x4B08gSO8LzP4w1MLSaIFJ8/iaPpkHlG8BVi/kVsHQq/m31zO7+bWY6CLqPR74wf38nMC/owuuB6Y/DRz4lFsz5VuAQ5/zlxqGFGDgFVw8Bv0QMAV8EWht5lZAxXZ+f6e3cwEo/1bld2EAcoYD+aN4FXz2YP5lRWvkX2ws1UDDUZ6w0FDGW7vb65XfugFuTSXncos2JV98nwckZfP79Dh5Q0lrLRc1aeu0+K+lNXBXrcHEhdiYwv/fMabwL1fGZACC/3h3q5zRx7c2cV8r/zvWJfCXZF3ojOIrwZ/04XXxtUkWhdfNrQ2fW3zvEf/uBPHvSwAEQPxPwL7ALZT7Aq2I9HZOoowSEoZYkFtiqGQlhQ0+B7iSjAGuJCCglsECq8sj91xKjcFiAMS6BTFltbWJfwPMSDKErgkI6awazLD8VAydNgX4cjVyPVWqxwDwWwvpRYoAewh9xvOHudSaW41jX/FtODeSxIgf82sd+IR/gw9+8J/ewV0pWgMw7meRr3XlE8DRtfyBeXyD8ht8SyX31QPAJff7YyVqGJKAW98HXp/J3UGvXAFM/jWQewF30ex5DwDjrpib3uHftMOROQD42RpuPXz1NK8B+favymP6XcqrwvtdGr6xoD4RGH0Tf9Ue5FXjFdu5BcB8QFohTwwYcDnQdwp3f6iRmM5/L9LvxufjD/czu/21Ncm5/N9AzvDw1wmHy8a/YEhV/aYMPlWQGibGHRKGWJAthtAYQ/jgc6AwiK4kp2gxSAVmTjOvYAbPXAqJV8gtt5XCIPVBqrU4gV78m57bwh/8BekqgWefl/tqAVWLQULoVczfBMxlCEFuhREmviDRm7fixpnd/M8Pjg0wFr0wFF/Cvzna6vjDfNCVys+lTKMR/xf5ASyte/wdfFb1x78EfrmRfzN1WoD3Z3NrKH80/1beFom9gDmf8Y6wJzYCG4NM/HFzgRl/VLgVw6LRAJN+BYy9jQe1z5TyOFTGAH6/sfqVc4YB0x6L7ZxIa8se0vbfebQYktQtMSLukDDEgtwSIxpXUvjgs12KMUgWgMMcPr4ABFgMygpkyWKoMTtkVxKz8epnVTeSrY5/YxQ0/qI4NSRhaIogDPVi4DlcfEEiZxgP6rmsXEyCffANx7g/XGsA+k6OfC2tjrtKtr0M7FmhFAZzFbB/JX+vFnRW48rHua+9/jDw2nTufvr+fR6ETcwAbnwrvDsqmKRM4LaVfA17P+K/6+wh3HIpnBDdNQIxpgBjbuUvguhmKF01FuSWGAF1DLIrKUhjJfEIEAYpxiC1vpAtAEdL+FRVwG+2B1sMkivJ4pSDz1oHFwbVjCSpHUZSdvisHsBfsdza6I9vBBOtxaDR+l0x0kCfQI6Jo1qLLoru2+PoW/h2/0qlcH39J+6SKJoMFIxp+zoA//NueY+7QppOAmsf5aJgygR++pFfIKNFowVG/h9w6wrgrnXAdX9vnygQRJwhYYiFCE30wlY+q8QY5OBzoCspksWg0hID8LuSasxO2QIwupoBAP3VJq7JcxgipKACPDgn1UWEcydFazEA/loCqbI5kKPSvOgr2r4OwH3j/S8DmJf74QGeIrnzDf7+soeiu45E5gDgFxt5y4whM4GL7wN+9R0PmhPEeQq5kmIhYhO9cFlJajGGYFdSS5SuJGWXU2lWc73ViVZdOhIBJPks0MCH/tkq8xPkOQwR6gkkevXjTfkaynjuumI9Fh7IBKIThgGXA5ue5cVc4mQ4ADwL5OTX/P3gGMazXv4IjzF8v4LHEg6t4dkfg2eEz1aKhCmDt8wgCAIAWQyxoWox8Id8dAVuQXUMksUQEGMIyUgSP+fHKy2GXkkGWUhOtnK3kgYMvWBB/+wOWAwATx0F/EVsgUjWQlK2Mr0xHH0u5BlY9gZeWCZx8mv+e0or9P950dBnPHDp7/j7b//Ks3fSCnnRGUEQHYaEIRZULIZWN88rD5uVFJDKaZIrnyVXkvigj+RK8jh5VgoQ4koCgGLRZXSy0QmvMR0AkK+3I1d1cpuYfhqNxSDFDupUWl1LYhGuFUYwWj3Q7wf8/ZGAxnP7V/HtoB/GnqI4dREw6y88bXPs7cCd/w2bgksQRGyQMMRCBIsh1JUUajEky72SQl1JzfZwxW1iqioEVWHol2kCAJxosKFVnw4AGJXhhkaj8qCNUPUcgmwxHA79TKphyI7CjSQhVSHvWcFTVF024OCnfN+oG6O/joQgAOPmAHM+5d1TUwtivwZBEKqQMMSCisUgxQuk+IGMWlZSSB1DOt86WtBo41ZBZlJQkVBg4Dm4uRyAflk8llBWY0UTuHCM6uVWX781hhiD9NBvOAp4PcrPavaLx8SQUz9sFncnNZ3gTd52veXvVloYoccQQRDdzlkhDC+++CKKi4uRkJCAiRMnYts2lZ4qATz//PMYMmQIEhMTUVhYiPvuuw8ORzcMTlEZ1COlniaHCINKVpLkSnKFpqs22njr7YykoCIoubhNfYraqEK+f3dFM057+PshpjA9Z2SLIQphSCvild5eF0/lDKRqD9/mjw45LSyGJLEjKYBPfgWs/wN/P2UBVboSxFlG3IXh/fffx8KFC7FkyRLs2rULo0ePxvTp01FbW6t6/LvvvosHH3wQS5YswcGDB/Hqq6/i/fffx0MPxZim2B4Cm+gxBo/XB4cYYwixGGRXkj/GIB1jV6l8bpCFIdhiaBaPDXUjAcDYwl4AgBP1Nuy3cuthQIIl9EDGAjqrRiEMGo3faqjd799vqRavIwB5I9q+TiA/uB9I7c3Pd1l5W+ySn8Z2DYIgupy4C8OyZctw1113Ye7cuRg+fDiWL18Ok8mE1157TfX4b7/9FlOmTMGtt96K4uJi/PCHP8Qtt9zSppXRKRhM4hsGeBywi6mqQHTzGCRhcHl9fO5zQOVzo1V0JSUHCYOUqhrGYkgz6TFGHL1ZxXiGUKpLRVTtjTylE4guxgAEzFPY6d9XJWYVZQ2OvZ2BKYN3IJ3wC2DKvbyhXLSVxQRBdBtxFQaXy4WdO3di2rRp8j6NRoNp06Zhy5YtqudMnjwZO3fulIXg+PHjWLNmDa6++uquX7De5H/vssnxBZ1GkKezyahkJaUYdbLXpKXV7bcCmBeuVu7+CbUYQsd6BnPXJf0BADWiMMB8JvQgKb6QmBF9s7M+F/JtYAO8M2J77PxR0V0jmF59eQO8Kx9X9uInCOKsIa4FbvX19fB6vcjNVX6Dzc3NxaFDKmmSAG699VbU19fj4osvBmMMHo8Hd999d1hXktPphNPplH82m8O0eIgGjZa7hjytgMsKm5tbA0lGXWh7arnAzW8xaDQCUhP0aGl1w9zqRk5yMm/r7PMgBXa0CgnyLGeZNmIMADBzVD5yUifBUKkF1v7VX3wWiDmGVFUJqQFe5S4egNbq+EQzACiaFP48giB6NHF3JcXKhg0b8PTTT+Pvf/87du3ahY8//hirV6/Gk08+qXr80qVLkZaWJr8KCzvYx1xyn7jssIqxgpDAs9SHHVDEGAB/Omqz3c2DrqIlkCrYkZ6ohzY4zdQR2ZUkcWFxBkZfIDaos1SFjohsqRAXEMP9Zw3mf66nlbuT3A5/vyOpLoEgiHOOuApDVlYWtFotampqFPtramqQl6f+zfbRRx/Fbbfdhp///OcYOXIkrr/+ejz99NNYunQpfCpzbhctWoSWlhb5VVFR0bFFy8Jgg11OVQ0TXwAUWUmAv0meVNAmPfBTYA91IwFhW26rkpwHQOCZRPYG5WeyMPRp+zoSGg0fPwkAh1fzNhQeBx+gkjkw+usQBNGjiKswGAwGjBs3DuvWrZP3+Xw+rFu3DpMmqbsq7HY7NEH5/FotfzAzlUHqRqMRqampilfHFi0Kg9smp6qagjurBnRfDaxjAPwWg18YJIvBhszgVFWgzeCzAp2Bz8EFQt1JzaIwxDr5aehMvt23kg/JAYDh11GKKUGcw8S9id7ChQsxZ84cjB8/HhMmTMDzzz8Pm82GuXPnAgBuv/129O7dG0uXLgUAzJo1C8uWLUNJSQkmTpyIsrIyPProo5g1a5YsEF1KgMUg1SOE1jCIwqDRh7S3Tg0WBsmVBDsSIloMUQgDwL/NW2v4FLLAOoOW03wbiysJ4I3pkrL5zISWcr6vZHZs1yAIokcRd2G46aabUFdXh8WLF6O6uhpjxozBF198IQeky8vLFRbCI488AkEQ8Mgjj6CyshLZ2dmYNWsWnnrqqe5ZsJSZ5LLJMYawriSVkZfpgTEGQG5Cly5Y4UuJIAwRspKUf0Ahn9PbXK7c354YA8BTdK98Elh1N/95/B2h3VYJgjiniLswAMD8+fMxf/581c82bNig+Fmn02HJkiVYsmRJN6xMBWlOsyLGEOxKCs1IkghxJSVyYcgQLNClqjS+izL4LJMxgG8bj/n3eT3+FNb2DBEfcwtvqtfaxGchEARxTnNWCEOPItCV1FY7DF2oxSAJg1kSBnEgTjqsSFIVhhiCzwAfPAPwkZkS5tN8sI3WEH1xWzC9x7bvPIIgehwkDLEiVT+7/emqIcFnqYGePvRBL2UlNQcJQ4ZgQa+0SMKQHt36MnixGxqP+/fVl/k/izTSkyAIAj2wjiHuyK4ka4DFEK7ldqgrKSuZ76uziFaFKAy9YEFesMXg9QAuad5zjK6k5nLAK4qPPIZzUHTXIAjivIaEIVbk4LNdzkoKbbkd2kBPwj+nmR/jNPAHfi/Bitxgi6G1SXwjRG8xpOQB+iTuOpK6ojYc5dtoxnASBHHeQ8IQKyoxhvDCEGox5KbyffVWJ7w+hgYfn+OcKViQEnwdqUgtIY23o4gGQQByxCE70hhNaeJaJlkMBEG0DQlDrMjCYJUH7oQEn9125bEBZCYboREAHwMabE6cdnKropdgCe231NrIt6K7KWryx/DtmVLA5/MLRO4FsV2HIIjzEhKGWJErn+2whLMYXKIwBHZjFdFqBDnOUGt24ridv0+Ay3+ehGQxxCoMBWP4tqoUaCjj1dO6RCBneGzXIQjivISEIVYCXElSymnInGa3TTw2VBgAIEd0J9WYHShrAlxMDF5LFoKEXbIYMmJbY+9xfHt6B+9vBHCxiNYdRRDEeQ0JQ6wECoODC0NqQjiLQX2QTUEadx9VNNpxqqkVTeBxhpDGd+21GHKGA6l9uEvrPw/wfQMuj+0aBEGct5AwxIr4sGcufxO91BCLQYoxqFsMA3J4yuuxOhvKaq1oYpIwBFsMkjDEaDEIAjD8GuW+4dfGdg2CIM5bSBhixSAJg1UeeZASYjGIriSVGAMADMjmwlBa0YwT9TY0MLGq2VanPFBKV02MURgAYMoCwJTF34++lbe0IAiCiAJyOsdKwKAeADDqNDDqggvcwgefAWCgaDHsreRVzRZ9FuADH7ATSHtdSQCvZ5i/nRe39ZkQ+/kEQZy3kMUQK6IwCG4bABbqRgL8MYYwrqTh+alIMvjFxNCrgL+xVCsPtNXzbayuJAlTBlB0ER+4QxAEESX0xIgVsf21wHxIgiM08Az4s5LCBJ8NOg2mDs2Rf+5TJPY3CrYYrOJku+QY5jQTBEF0EHIlxYo+EdDoAJ8HKbC3y2IAgEdmDoPL48PQvBQMKnABu6G0GBjz/5xCwkAQRPdBwhArgsCthtZGpAitSE1QEYY2YgwAkJ+WiFduH89/KBeDzIEWg70R8IlN8NrbKpsgCKIdkCupPYizEcJbDFKBm7orKQTJIrBUQ051sorWgimTz3ImCILoJkgY2oM0p1mwh4kxtG0xKEjJ51uPA3A08/eS9UDxBYIguhkShvYgzkZoO8YQpcWgM/pTUltO8y3FFwiCiBMkDO1BtBhUYwyMBWQlRWkxAECvfnwrTV6TBCI1vwMLJQiCiB0ShvYQEGMIaaDncQLMx99HyEoKIXgkZ+MJvpUEgyAIopsgYWgPkitJsCMjKSgw7A5onR2mjkGVTHEkZ8Mxvm0ShSGDhIEgiO7lrBCGF198EcXFxUhISMDEiROxbdu2iMc3Nzdj3rx5yM/Ph9FoxODBg7FmzZpuWi38riSoCIPLyrdaY2xtrmWLQRQEaSwnWQwEQXQzca9jeP/997Fw4UIsX74cEydOxPPPP4/p06fj8OHDyMnJCTne5XLhyiuvRE5ODj766CP07t0bp06dQnp6evctOsGflZSRFORKclr41pgS2zUzRIuh7hDQ2uyveiaLgSCIbibuwrBs2TLcddddmDt3LgBg+fLlWL16NV577TU8+OCDIce/9tpraGxsxLfffgu9nj+Ui4uLu3PJ8BhSoINkMQTNdXaY+VYUj6jJHc4rqu31wGHR+kkrAhJ7dXi9BEEQsRBXV5LL5cLOnTsxbdo0eZ9Go8G0adOwZcsW1XM+/fRTTJo0CfPmzUNubi5GjBiBp59+Gl6vt7uWDRt4UDlVaA0NPrfXYtAnArkj+Pvtr/Jt/qgOrJIgCKJ9xNViqK+vh9frRW6usuVDbm4uDh06pHrO8ePHsX79esyePRtr1qxBWVkZfvWrX8HtdmPJkiUhxzudTjidTvlns9nc4XWbmQlpANI1rdBqhKA/ULy+MUaLAQAKJ/A5zZU7/D8TBEF0M2dF8DkWfD4fcnJy8PLLL2PcuHG46aab8PDDD2P58uWqxy9duhRpaWnyq7CwsMNraPImAADSBHvohx0RhguuV/485OrYr0EQBNFB4ioMWVlZ0Gq1qKmpUeyvqalBXp56xW9+fj4GDx4MrdY/z2DYsGGorq6Gy+UKOX7RokVoaWmRXxUVFR1ed4OPD9pJhTX0w/bGGACgaBIweAZ/P+anQNagdq6QIAii/cRVGAwGA8aNG4d169bJ+3w+H9atW4dJkyapnjNlyhSUlZXB5/PJ+44cOYL8/HwYDKHN5oxGI1JTUxWvjlLj4TGGRNYKuB3KD9sbYwB459Zb3gPu2w9c+7cOrpIgCKJ9xN2VtHDhQrzyyit48803cfDgQdxzzz2w2WxyltLtt9+ORYsWycffc889aGxsxIIFC3DkyBGsXr0aTz/9NObNm9dta66w6+Fh4q+utVH5oSwM7RQgQQDS+vAtQRBEHIh7uupNN92Euro6LF68GNXV1RgzZgy++OILOSBdXl4OTcBoysLCQnz55Ze47777MGrUKPTu3RsLFizA73//+25bc43FhSakIBstfC5zaoH/QznG0A6LgSAI4iwg7sIAAPPnz8f8+fNVP9uwYUPIvkmTJuG7777r4lWFp8bsQCNLQbYgCkMgksXQnhgDQRDEWUDcXUk9kRqzA00QLYJgYXC08G17XUkEQRBxhoShHVS3cIsBAB/BGYg0aEdstEcQBNHTIGGIkVaXF2aHB00sjMVgF+c3mzK6d2EEQRCdBAlDjFQ08aI2m060CIItBkkoEkkYCILomZAwxMjJej6dTZMkjuIMtBhcdsDTyt9LozoJgiB6GCQMMXKygQuDPiWb77DV+j+Uaho0OkpXJQiix0LCECMnRIshMbM332Gp9n8ouZVMmVSgRhBEj4WEIUb2VfICtpwCceKa+Yz/Q4ovEARxDkDCEAMOtxeHqrkwDBw4mO90Wf2N81oDLAaCIIgeCglDlHh9DD/4f1/B7WXIS01A75xMwChmJlmq+NYqxhuSsuKzSIIgiE6AhCFK/rOvCrUWPvDn2pICCILg75FkrhS3olspsHcSQRBED+Os6JXUE5g5Mh/2//PiZL0Nv75cnJOQmg/UHQTMosUgWQ4p+fFZJEEQRCdAwhAlgiDgxvFB098ky6DlNN+SxUAQxDkAuZI6QsYAvm0o41sSBoIgzgFIGDqCNHqz/gjg8/pjDam947cmgiCIDkLC0BEyRWFoKAOaTgJeF6BL4BPYCIIgeigkDB0hox8gaHgtw4lN4r4BgEYb33URBEF0ABKGjqAzAjnD+ftdb/Kt5F4iCILooZAwdJSiSXx7ZjffFk6I31oIgiA6ARKGjtJ/qvLn4kvisgyCIIjOgoShowy+CsgQG+r1mQDkjYzvegiCIDoIFbh1FK0O+Nka4MgXwAXXU7ttgiB6PCQMnUFqPjB+brxXQRAE0SmcFa6kF198EcXFxUhISMDEiROxbdu2qM5bsWIFBEHAdddd17ULJAiCOI+IuzC8//77WLhwIZYsWYJdu3Zh9OjRmD59OmprayOed/LkSdx///245BIK9hIEQXQmcReGZcuW4a677sLcuXMxfPhwLF++HCaTCa+99lrYc7xeL2bPno3HH38c/fv378bVEgRBnPvEVRhcLhd27tyJadOmyfs0Gg2mTZuGLVu2hD3viSeeQE5ODu688842/wyn0wmz2ax4EQRBEOGJqzDU19fD6/UiNzdXsT83NxfV1dWq52zevBmvvvoqXnnllaj+jKVLlyItLU1+FRYWtn0SQRDEeUzcXUmxYLFYcNttt+GVV15BVlZ04zMXLVqElpYW+VVRUdHFqyQIgujZxDVdNSsrC1qtFjU1NYr9NTU1yMvLCzn+2LFjOHnyJGbNmiXv8/l8AACdTofDhw9jwIABinOMRiOMRmMXrJ4gCOLcJK7CYDAYMG7cOKxbt05OOfX5fFi3bh3mz58fcvzQoUOxd+9exb5HHnkEFosFL7zwQlRuIsYYAFCsgSCI8wrpmSc9AyMR9wK3hQsXYs6cORg/fjwmTJiA559/HjabDXPn8oKx22+/Hb1798bSpUuRkJCAESNGKM5PT08HgJD94bBYLABAsQaCIM5LLBYL0tLSIh4Td2G46aabUFdXh8WLF6O6uhpjxozBF198IQeky8vLodF0XiikoKAAFRUVSElJgRBj+wqz2YzCwkJUVFQgNTW109Z0tkH3eW5B93lu0d77ZIzBYrGgoKDt0cMCi8auIADwv5C0tDS0tLSc8//w6D7PHeg+zy264z57VFYSQRAE0fWQMBAEQRAKSBhiwGg0YsmSJed8+ivd57kF3ee5RXfcJ8UYCIIgCAVkMRAEQRAKSBgIgiAIBSQMBEEQhAIShihp75S5s5WlS5fiwgsvREpKCnJycnDdddfh8OHDimMcDgfmzZuHzMxMJCcn48c//nFIX6uexjPPPANBEHDvvffK+86V+6ysrMRPf/pTZGZmIjExESNHjsSOHTvkzxljWLx4MfLz85GYmIhp06bh6NGjcVxx7Hi9Xjz66KPo168fEhMTMWDAADz55JOKNg898T43bdqEWbNmoaCgAIIgYNWqVYrPo7mnxsZGzJ49G6mpqUhPT8edd94Jq9XavgUxok1WrFjBDAYDe+2119j+/fvZXXfdxdLT01lNTU28l9Zupk+fzl5//XW2b98+Vlpayq6++mpWVFTErFarfMzdd9/NCgsL2bp169iOHTvYRRddxCZPnhzHVXeMbdu2seLiYjZq1Ci2YMECef+5cJ+NjY2sb9++7Gc/+xnbunUrO378OPvyyy9ZWVmZfMwzzzzD0tLS2KpVq9iePXvYNddcw/r168daW1vjuPLYeOqpp1hmZib7/PPP2YkTJ9iHH37IkpOT2QsvvCAf0xPvc82aNezhhx9mH3/8MQPAVq5cqfg8mnu66qqr2OjRo9l3333Hvv76azZw4EB2yy23tGs9JAxRMGHCBDZv3jz5Z6/XywoKCtjSpUvjuKrOpba2lgFgGzduZIwx1tzczPR6Pfvwww/lYw4ePMgAsC1btsRrme3GYrGwQYMGsbVr17JLL71UFoZz5T5///vfs4svvjjs5z6fj+Xl5bFnn31W3tfc3MyMRiN77733umOJncLMmTPZHXfcodh3ww03sNmzZzPGzo37DBaGaO7pwIEDDADbvn27fMx//vMfJggCq6ysjHkN5Epqg/ZOmetptLS0AAAyMjIAADt37oTb7Vbc99ChQ1FUVNQj73vevHmYOXOm4n6Ac+c+P/30U4wfPx4/+clPkJOTg5KSEsUwqxMnTqC6ulpxn2lpaZg4cWKPus/Jkydj3bp1OHLkCABgz5492Lx5M2bMmAHg3LnPQKK5py1btiA9PR3jx4+Xj5k2bRo0Gg22bt0a858Z9yZ6ZzuRpswdOnQoTqvqXHw+H+69915MmTJF7lJbXV0Ng8Egd6+ViDRd72xlxYoV2LVrF7Zv3x7y2blyn8ePH8dLL72EhQsX4qGHHsL27dvxm9/8BgaDAXPmzJHvJZZpiWcjDz74IMxmM4YOHQqtVguv14unnnoKs2fPBoBz5j4DieaeqqurkZOTo/hcp9MhIyOjXfdNwkBg3rx52LdvHzZv3hzvpXQ6FRUVWLBgAdauXYuEhIR4L6fL8Pl8GD9+PJ5++mkAQElJCfbt24fly5djzpw5cV5d5/HBBx/gnXfewbvvvosLLrgApaWluPfee1FQUHBO3We8IVdSG8Q6Za6nMX/+fHz++ef46quv0KdPH3l/Xl4eXC4XmpubFcf3tPveuXMnamtrMXbsWOh0Ouh0OmzcuBF/+ctfoNPpkJube07cZ35+PoYPH67YN2zYMJSXlwOAfC89/d/xAw88gAcffBA333wzRo4cidtuuw333Xcfli5dCuDcuc9AormnvLw81NbWKj73eDxobGxs132TMLRB4JQ5CWnK3KRJk+K4so7BGMP8+fOxcuVKrF+/Hv369VN8Pm7cOOj1esV9Hz58GOXl5T3qvq+44grs3bsXpaWl8mv8+PGYPXu2/P5cuM8pU6aEpBsfOXIEffv2BQD069cPeXl5ivs0m83YunVrj7pPu90eMp9Fq9XKI37PlfsMJJp7mjRpEpqbm7Fz5075mPXr18Pn82HixImx/6HtDp2fR6xYsYIZjUb2xhtvsAMHDrBf/OIXLD09nVVXV8d7ae3mnnvuYWlpaWzDhg2sqqpKftntdvmYu+++mxUVFbH169ezHTt2sEmTJrFJkybFcdWdQ2BWEmPnxn1u27aN6XQ69tRTT7GjR4+yd955h5lMJvb222/LxzzzzDMsPT2dffLJJ+z7779n11577VmfxhnMnDlzWO/eveV01Y8//phlZWWx3/3ud/IxPfE+LRYL2717N9u9ezcDwJYtW8Z2797NTp06xRiL7p6uuuoqVlJSwrZu3co2b97MBg0aROmqXc1f//pXVlRUxAwGA5swYQL77rvv4r2kDgFA9fX666/Lx7S2trJf/epXrFevXsxkMrHrr7+eVVVVxW/RnUSwMJwr9/nZZ5+xESNGMKPRyIYOHcpefvllxec+n489+uijLDc3lxmNRnbFFVeww4cPx2m17cNsNrMFCxawoqIilpCQwPr3788efvhh5nQ65WN64n1+9dVXqv8/zpkzhzEW3T01NDSwW265hSUnJ7PU1FQ2d+5cZrFY2rUe6q5KEARBKKAYA0EQBKGAhIEgCIJQQMJAEARBKCBhIAiCIBSQMBAEQRAKSBgIgiAIBSQMBEEQhAISBoIgCEIBCQNBEAShgISBILqJqVOnKmZNE8TZCgkDQRAEoYB6JRFEN/Czn/0Mb775pmLfiRMnUFxcHJ8FEUQESBgIohtoaWnBjBkzMGLECDzxxBMAgOzsbGi12jivjCBCodGeBNENpKWlwWAwwGQy9dhJYsT5A8UYCIIgCAUkDARBEIQCEgaC6CYMBgO8Xm+8l0EQbULCQBDdRHFxMbZu3YqTJ0+ivr5eHmBPEGcbJAwE0U3cf//90Gq1GD58OLKzs1FeXh7vJRGEKpSuShAEQSggi4EgCIJQQMJAEARBKCBhIAiCIBSQMBAEQRAKSBgIgiAIBSQMBEEQhAISBoIgCEIBCQNBEAShgISBIAiCUEDCQBAEQSggYSAIgiAUkDAQBEEQCv4/IVh6yEmWgkcAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 400x300 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "sol = ode.run()\n",
    "t = sol[\"t\"]\n",
    "x = sol[\"x\"]\n",
    "plt.figure(figsize=(4, 3))\n",
    "plt.plot(t, x[:, 0], label=\"$\\\\theta$\")\n",
    "plt.plot(t, x[:, 1], label=\"$\\omega$\")\n",
    "plt.xlabel(\"t\")\n",
    "plt.ylabel(\"x\")\n",
    "plt.legend()\n",
    "plt.tight_layout()\n",
    "plt.savefig(\"output/damp_oscillator_ts.jpeg\", dpi=300)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Selected features:\n",
      "------------------\n",
      "■ Domain: statistical\n",
      " ▢ Function:  calc_std\n",
      "   ▫ description:  Computes the standard deviation of the signal.\n",
      "   ▫ function   :  vbi.feature_extraction.features.calc_std\n",
      "   ▫ parameters :  {'indices': None, 'verbose': False}\n",
      "   ▫ tag        :  all\n",
      "   ▫ use        :  yes\n",
      " ▢ Function:  calc_mean\n",
      "   ▫ description:  Computes the mean of the signal.\n",
      "   ▫ function   :  vbi.feature_extraction.features.calc_mean\n",
      "   ▫ parameters :  {'indices': None, 'verbose': False}\n",
      "   ▫ tag        :  all\n",
      "   ▫ use        :  yes\n"
     ]
    }
   ],
   "source": [
    "cfg = get_features_by_domain(domain=\"statistical\")\n",
    "cfg = get_features_by_given_names(cfg, names=[\"calc_std\", \"calc_mean\"])\n",
    "report_cfg(cfg)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def wrapper(par, control, cfg, verbose=False):\n",
    "    ode = DO_cpp(par)\n",
    "    sol = ode.run(control)\n",
    "\n",
    "    # extract features\n",
    "    fs = 1.0 / par[\"dt\"] * 1000  # [Hz]\n",
    "    stat_vec = extract_features(\n",
    "        ts=[sol[\"x\"].T], cfg=cfg, fs=fs, n_workers=1, verbose=verbose\n",
    "    ).values\n",
    "    return stat_vec[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "def batch_run(par, control_list, cfg, n_workers=1):\n",
    "    stat_vec = []\n",
    "    with Pool(processes=n_workers) as pool:\n",
    "        stat_vec = pool.starmap(\n",
    "            wrapper, [(par, control, cfg) for control in control_list]\n",
    "        )\n",
    "    return stat_vec"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.12421611 0.10675827 1.0513536  0.8769057 ]\n"
     ]
    }
   ],
   "source": [
    "control = {\"a\": 0.11, \"b\": 0.06}\n",
    "x_ = wrapper(parameters, control, cfg)\n",
    "print(x_)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "num_sim = 2000\n",
    "num_workers = 10\n",
    "a_min, a_max = 0.0, 1.0\n",
    "b_min, b_max = 0.0, 1.0\n",
    "prior_min = [a_min, b_min]\n",
    "prior_max = [a_max, b_max]\n",
    "theta_true = {\"a\": 0.1, \"b\": 0.05}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "prior = utils.torchutils.BoxUniform(\n",
    "    low=torch.as_tensor(prior_min), high=torch.as_tensor(prior_max)\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "obj = Inference()\n",
    "theta = obj.sample_prior(prior, num_sim)\n",
    "theta_np = theta.numpy().astype(float)\n",
    "control_list = [{\"a\": theta_np[i, 0], \"b\": theta_np[i, 1]} for i in range(num_sim)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "stat_vec = batch_run(parameters, control_list, cfg, n_workers=4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "scaler = StandardScaler()\n",
    "stat_vec_st = scaler.fit_transform(np.array(stat_vec))\n",
    "stat_vec_st = torch.tensor(stat_vec_st, dtype=torch.float32)\n",
    "torch.save(theta, \"output/theta.pt\")\n",
    "torch.save(stat_vec_st, \"output/stat_vec_st.pt\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(torch.Size([2000, 2]), torch.Size([2000, 4]))"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "theta.shape, stat_vec_st.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " Neural network successfully converged after 297 epochs.train Done in 0 hours 0 minutes 40.652620 seconds\n"
     ]
    }
   ],
   "source": [
    "posterior = obj.train(theta, stat_vec_st, prior, method=\"SNPE\", density_estimator=\"maf\", num_threads=4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "with open(\"output/posterior.pkl\", \"wb\") as f:\n",
    "    pickle.dump(posterior, f)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "# with open(\"output/posterior.pkl\", \"rb\") as f:\n",
    "#     posterior = pickle.load(f)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "xo = wrapper(parameters, theta_true, cfg)\n",
    "xo_st = scaler.transform(xo.reshape(1, -1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "3b235dd44d2a4a6f8e8ecf06ba059edc",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Drawing 10000 posterior samples:   0%|          | 0/10000 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "samples = obj.sample_posterior(xo_st, 10000, posterior)\n",
    "torch.save(samples, \"output/samples.pt\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAbgAAAHcCAYAAACkr7//AAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjAsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvlHJYcgAAAAlwSFlzAAAPYQAAD2EBqD+naQAAHMZJREFUeJzt3XuMneV9J/DvmTMXg+0xNuZiB2NCIG1SUuhCuASyIYIN2mxQU1UVbVeCtsqlu2pE6vzDoiaWWrXsbi7bbpOVAqiCLNs2WXWbRNstTZc0oVkghFBaShtCKdimxRhzsz1gz8w5Z/+w58ycGQ9mPO/M+553Ph/JSnzm2P4xmvHXv+d9fs/T6HQ6nQBAzQyUXQAALAUBB0AtCTgAaknAAVBLAg6AWhJwANSSgAOglgQcALU0WOYfPjY+ljW3rEmSHPgPB7J6eHWZ5cCK868GfqbsEuCY/rz9P4/r1+ngAKglAQdALQk4AGpJwAFQSwIOgFoScADUUmUC7oF/3Ft2CQDUSGUC7rtPvVh2CQDUSGUCbmKyXXYJANRIZQLukIADoECVCbjxtoADoDjVCbiJTtklAFAj1Qm4lg4OgOJUJ+AmW2WXAECNVCbgbDIBoEiVCThjAgAUqTIB5xkcAEUScADUUmUC7tCEgAOgOJUJuAkdHAAFqkzAjdtkAkCBKhNwxgQAKFJlAs4mEwCKVJ2A08EBUKDqBFyrnU7HgcsAFKPUgGu3pwOt00kmWgIOgGKUGnATs+6A8xwOgKKUG3CzOjbP4QAoSqkBNzmrYxNwABSl1IAbn7VEecidcAAUpOQOzhIlAEuj3Gdwk7M7OAEHQDHK7eDsogRgiZT7DG6yd4nSlTkAFKVauyh1cAAUpOQ5OGMCACyNkk8ymbVEaUwAgILo4ACoJQEHQC1Va9DbJhMAClKpw5aNCQBQlGotUergAChIpQLOUV0AFKVSS5Q2mQBQlEqdRWkODoCiVGqJUgcHQFEsUQJQSyXfJmCTCQBLo9SAa82+D07AAVCQai1RmoMDoCCVWqLUwQFQlJKvyzEmAMDSqNZhyzo4AApSqYCzixKAopT7DM6gNwBLpFInmejgAChKpZYojQkAUJRKdXCWKAEoSqUCzpgAAEUpeQ7OmAAAS6PkZ3C9gdbuzH0NAI5HyWMCnaO8JuAAWLxKdXBJcmhCwAGweJXZZNIcaCTRwQFQjMrMwQ03D5diowkARajMUV3Dg4dLMSoAQBEqc+HpcPPwEqXjugAoQmU2mUx1cJYoAShCZQa9BRwARSot4DqdTs8uyqHm1DM4AQfA4pUWcK12J50Zc94jQ80kOjgAilFawE3MOsVkZGpMwBwcAAUoL+DavUFmDg6AIpUXcLOCzBwcAEWqzBKlXZQAFKnEgJuvgxNwACxeZQLOmAAARarMEuXI4JHbBAQcAAWoTAc33DwyB2dMAIAClBZws4PMJhMAilRawE3OWqIcMiYAQIEqs0Q5ooMDoECVWaIUcAAUqTonmTiLEoAClfcMrj3rGdzUHNyEgANg8arzDG5IBwdAccp7BjdnifLwHJyTTAAoQmVOMhk6cpKJgAOgCCU+g7OLEoClU6ElyqmAM+gNwOJVZomye1SXTSYAFKAyuyiNCQBQpBLPouwNslXGBAAoUIlHdc1aopy6LscmEwAKUJklStflAFCkygTc1DO4yXYnrVnHeAHAQlVmF+XUHFyiiwNg8SrTwQ0LOAAKVJmAGxxopHH4tC63egOwaJUJuEaj0T3NxHmUACxWZZ7BJTPOozQLB8AiVaaDS5LhQbNwABSjUgE31cFZogRgscoLuMnXWaIUcAAsUnkB1z7aEqWAA6AYlVqinL4yx5gAAItTqSXKYVfmAFCQSnVwI67MAaAg1XoGZ9AbgIJUa4nSmAAABanUEqVBbwCKUuKN3vMPegs4ABartICbPMpZlObgAChKtZYou5tMzMEBsDilBFy73clk21FdACydUgLuaCMCietyAChOKQF3tOdvyYwxASeZALBI5XRw83Rowzo4AApSSsBNBVij0fv6iDk4AApSUgd3eIlyaKD3j3eSCQBFKekZ3OEAG2r2tnDGBAAoSqnP4AabR+/gLFECsFjlPIM7ctDy7A7OmAAARSm1gxsamLVEqYMDoCDlPIM7Mug9NGiTCQBLY7CMP3R6ibI34BzVBdXRGJz+66EzOVliJXB8KrXJxBwcAEUp9xnc7DGBQWMCABSj3EHv2WMCTUuUABSjUrso1584nCQZG2/l1XFr/gAcv5KXKHv/+HUnDmV01eEH27tefG3Z6wKgPiq1ySRJzjz5xCTJzhdfXdaaAKiXkp/BNeZ8bOuG1UmSHS+MLWtNsCINNKd/zNKZnOz+gH5UqSXKJNmy4XAHt0sHB8AilBpww0cJuK2WKAEoQCknmUwtUTaPskR55pEOboeAg6XXMZJDfVVq0DuZDrhnXnwt7XZnWesCoD5KnoOb+8dvWrcqgwONjLfa2b3v4HKXBkBNlLqLcnhw7h8/2BzIm9afkMRzOFhync70j0aj9wf0uZIuPD0yBzdw9G+iqWVKAQfA8Sr1PrijDXonMwLuBQEHwPEpZ4lycv4lykQHB8DilTQmMP8uymR6Fs6oACyfl66/tOfn6++8v6RKoBjlPIObOovyKLsoE6eZALB45TyDm+c+uClTS5Qvjo1n/8GJZasLgPqo5BLl2lVD2bB6OC+OjWfni6/mxzavW87yYEXa8D++1/Pzu3b9v3zvufF86oF9GfiPl+ekZzaWVBkcn1KXKOfr4BLLlFAFtz82loefn8iui58suxRYsMrdBzdl69SZlEYFoBQvHGznfz91+DShPW9/JuMnHCq5IliYUp/BDc+zRJkYFYCyfeWHr+bIt2o6jU6ePX9HuQXBApX6DG5QwEElHFz7WsZX9579+qW/f236J43kmXc+mfVPn9LznuGxVVm1/4TlKBEWrJSAGz/GLsokOdO9cLBs/vanv5uXt+7tea2xP+ne59FIXls/lgc/ek/Pe056emMuuvPK5SgRFqzkZ3DH7uD+6aXXMtlyZxUspc0PvzkDEwMzEq3n/x4289u1kwxMDGTzX715GaqD41NKBzcVWMPzDHonyemjqzLcHMh4q51nXznY3VUJFG/z32zN6D+vz99cd19e3XDg9f/p205OfHFtfvzLl2XN3tFlqxEWqtTrcl5viXJgoJEzNrg2B5bLmr2jueTWq3PaY2e87vtOe2xLLvniVcKNyiv3qK7XCbhkepny7r/dnU7H7d6w1JoTg1m/85SjrE8e0UnW79iY5mQpiz+wIOXe6P06z+CS5Kf/xeF/Sf73B3bkv/zfJ5a8LiDZt+mlNNpHvjengu7I/zbajezb/HIZZcGCVfIsyinXnr85n/rA25Mk//WeJ/L5bwo5WGqvnPFiOs1OGq1GBiYHcub952ZgciCNViOdZievnPFC2SXCG1LSmMCRDm6e++Bm+qUr3pyJVju3/OkP8plv/DD7D07m3135lpx04vBSlwkrTqvZytjG/UmSE15a091IsvnhNx/egLLxQMY27k+r2Uqz1Sy5Wnh9yx5wnU7nDY0JzPTR97wl45PtfPbPf5gv3vuP+dL9O/IzF52RX7r8zTlr4+qlLBdWlPZQK2v2jGbts+vzo//ngu6ztqkNKD/413+V/ae/nPaggKP6lj3gWu1OpvaLDDXe+Arpx646N2dtXJ3/9q0n8/fP7suX7t+RL92/I2/euDoXn7Uhl5y9IeecuiYnrxnJyauHs2rINx8s1NDB4VzyxavTyNx/fDYnBvNjX39nOukc9eNQNYsOuG1ffmRB72/N2A35RpYoZ7r2/M35wI9vyn1PvpDb/vIf8+0fPp+n9o7lqb1j+fJDu3reu3q4mbWrhrJm1WDWrhrMicPNjAw2s2poIEPNgTSSNBqHv0l9q9bf5667oOwS+saxwku40S8WHXD/66/+6bh+3erhZkYWGHDJ4VC6/JyNufycjXnltYk89PSLefCpF/O9p1/Ms68czAsHxjPeamdsvJWx8Vay77jKo2YEHKw8iw64m9//o8f16y7cuuGYuyiPZd0JQ7nqbaflqred1n2t0+lk38HJvDQ2nv0HJ7P/4ET2H5rMwYnWkR/tTLTa3WXSzrwDPwD0s0UH3Ef+5VuO+9eOjY8t9o+fo9FoZN0JQ1l3wlDhvzcA/aOUOTgAWGoCDoBaEnAA1JKAA6CWBBwAtdTouIcGgBrSwQFQSwIOgFoScADUkoADoJYEHAC1tKizKDudTvbv319ULbCk1q5d270iCai/RQXc3r17c+qppxZVCyypPXv25JRTTim7DGCZLCrghoeHkyS7du3K6OhoIQVB0fbt25ctW7Z0v16BlWFRATe13DM6OirgqDzLk7Cy2GQCQC0JOABqaVEBNzIyku3bt2dkZKSoeqBwvk5hZXLYMgC1ZIkSgFoScADUkoADoJYEHAC1JOAAqKXjCrjvfe97ef/735+TTjopq1evzqWXXpqvfOUrRdcGx+2uu+7KRz/60Vx00UUZGRlJo9HIHXfcUXZZwDJa8FFdf/EXf5Frrrkmq1atys/+7M9m7dq1+aM/+qNcd9112bVrVz7xiU8sRZ2wIL/2a7+WHTt2ZOPGjdm0aVN27NhRdknAMltQBzc5OZkPf/jDGRgYyL333ptbb701n/3sZ/PXf/3Xeetb35qbb77ZXyRUwu23356nn346zz//fH75l3+57HKAEiwo4L75zW/mySefzM///M/nggsu6L6+bt263HzzzRkfH8+dd95ZdI2wYFdffXW2bt1adhlAiRYUcN/61reSJO973/vmfOyaa65Jknz7299efFUAsEgLCrgnnngiSXLuuefO+djpp5+eNWvWdN8DAGVaUMC98sorSQ4vSR7N6Oho9z0AUCZzcADU0oICbqpzm69L27dv37zdHQAspwUF3NSzt6M9Z9u9e3cOHDhw1OdzALDcFhRw73nPe5Ik3/jGN+Z87M/+7M963gMAZVpQwF111VU5++yz8/u///t55JFHuq+/8sor+a3f+q0MDw/n+uuvL7pGAFiwBd/oPd9RXTt27MhnPvMZR3VRCbfffnu+853vJEkeffTRPPzww7n88stzzjnnJEmuuOKKfOhDHyqzRGCJLTjgkuTBBx/M9u3bc99992ViYiLveMc7sm3btlx33XVLUSMs2C/8wi+87qk6N9xwg8OXoeaOK+AAoOrMwQFQSwIOgFoScADUkoADoJYEHAC1JOAAqCUBB0AtCTgAaknAAVBLAg6AWhJwANSSgKuQu+++O1dccUVOOumknHzyyfnABz6QJ598suyyAPqSgKuQsbGxbNu2LQ899FDuueeeDAwM5Kd+6qfSbrfLLg2g77hNoML27t2bU045JY8++mjOO++8sssB6Cs6uAp54okn8nM/93M5++yzMzo6mrPOOitJsnPnznILA+hDg2UXwLRrr702W7duzW233ZbNmzen3W7nvPPOy/j4eNmlAfQdAVcRL7zwQh5//PHcdtttefe7350k+c53vlNyVQD9S8BVxPr163PyySfn1ltvzaZNm7Jz587cdNNNZZcF0Lc8g6uIgYGB/OEf/mG+//3v57zzzsuv/uqv5tOf/nTZZQH0LbsoAaglHRwAtSTgAKglAQdALQk4AGpJwAFQSwIOgFoScADUkoADoJYEHAC1JOAAqCUBB0AtCTgAaknAAVBLAg6AWhJwANSSgAOglgQcALUk4ACoJQEHQC0Nll0AUE2dTievTryaJDlx6MQ0Go2SK4KF0cEBR/XqxKtZc8uarLllTTfooJ8IOABqScABUEsCDoBaEnAA1JKAA6CWjAkAcxw4NJnvPr237DJgUQQcMMd/+tMf5M4HHk9OKLsSOH6WKIE5nn3ltbJLgEUTcMAchybbZZcAiybggDnGBRw1IOCAOSZaAo7+J+CAOcYFHDUg4IA5LFFSBwIOmGOi1Sm7BFg0AQfMoYOjDgQcMIcxAepAwAFz2EVJHQg4YA5LlNSBgAPm0MFRBwIO6NFudzLZtouS/ifggB6GvKkLAQf0EHDUhYADethgQl0IOKCHgKMuBBzQww5K6kLAAT10cNSFgAN62GRCXQg4oIcOjroQcEAPAUddCDigh7vgqAsBB/QYb7XKLgEKIeCAHpYoqQsBB/QYt0RJTQg4oIcOjroQcEAPAUddCDigh6O6qAsBB/TQwVEXAg7o4agu6kLAAT10cNSFgAN66OCoCwEH9JjQwVETAg7ooYOjLgQc0MMzOOpCwAE9dHDUhYADeujgqAsBB/QQcNSFgAN6OKqLuhBwQA/P4KgLAQf0mJg8fB/c8KC/HuhvvoKBHoeOdHCrh5slVwKLI+CAHlObTFaPCDj6m4ADekxtMjlxeLDkSmBxBBzQY6qDWzMi4OhvAg7oMRVwOjj6nYADekwtUerg6HcCDujR7eBsMqHPCTighzEB6kLAAV2dTqe7RLnaEiV9TsABXZPtTjqHDzIRcPQ9AQd0zTxoebVdlPQ5AQd0zbwqx0km9DsBB3RNBdxAIzlhSMDR3wQc0DV1Vc5QcyCDzUbJ1cDiCDiga6qDGx4cyFDTXw/0N1/BQNdUBzcyOJBhAUef8xUMdE1ddjrU1MHR/3wFA13jrVYSS5TUg69goOvQ1DO45kCGB/31QH/zFQx0TbRmLlHaRUl/E3BAl12U1ImvYKBr6qiu4eZAhixR0ud8BQNdMzu44YHpvx5a7U5ZJcFxE3BAV88S5YwObuYhzNAvBBzQNX1UV6Nnk8m4gKMPCTiga7qDa/acZCLg6EcCDugan7HJpNGY7uAmBRx9SMABXRPdDq4x63WbTOg/Ag7omtnB9b7eKqMcWBQBB3R1A27WDNy4Do4+JOCArqlNJrNPMfEMjn4k4ICumXNwPa8LOPqQgAO6JlpH7+AMetOPBBzQNdXBjejgqAEBB3TNt8nEmAD9SMABXVO7JecuUQo4+o+AA7rmm4PzDI5+JOCArvHJwwPddlFSBwIO6Jpaipy9RDm1+QT6iYADuubbRTnZFnD0HwEHdE3Mt4vSEiV9SMABXfMd1WWJkn4k4ICuQ/Me1WVMgP4j4ICu6aO6eu+Dc9gy/UjAAV1T4wBzjuqyREkfEnBAV/c2gWaz53WbTOhHAg7o6i5RDvYuUU4YE6APCTggSdJud7qD3rOP6nKjN/1IwAFJeo/jmjsHJ+DoPwIOSNL7nG3ObQI2mdCHBByQpHen5JwlSptM6EMCDkgy86DlRgYGzMHR/wQckGT+Y7pmfgz6iYADkiTjraPfBZck48YE6EMCDkgyPQpwtA7OoDf9SMABSaY3kszeYJIkE+bg6EMCDkgy/2WniQ6O/iTggCQzbxKwyYR6EHBAkhkHLevgqAkBBySZ/7LTxIWn9CcBBySZ/7LTxKA3/UnAAUlmLlE253zMYcv0IwEHJJnu4I42JjDeaqfTEXL0FwEHJJkxBzc4d4ky0cXRfwQckGTGEuVROrjETkr6j4ADkkx3cEebg0vMwtF/BByQ5PXn4BIdHP1HwAFJjh1wLj2l3wg4IMnr76I8/HGbTOgvAg5I8gY6OM/g6DMCDkgyfRyXXZTUhYADkkx3aEOewVETAg5I8voXniaWKOk/Ag5IkkwYE6BmBByQ5NgdnICj3wg4IIldlNSPgAOSvIGjuszB0WcEHJDkDRzVpYOjzwg4IImjuqgfAQckmd5EMtSc7z44AUd/EXBAkukObcQmE2pCwAFJZszBNZtH/bglSvqNgAOSzNhFOTjPEuWkXZT0FwEHJEkOTRr0pl4EHJBkxn1wdlFSEwIOSDJjTMBhy9SEgAMy2WqnfeQRm8OWqQsBB2RixjFc8x7VpYOjzwg4oCe8dHDUhYADcqjVSpI0GsngwHwnmRgToL8IOKAbXkPNgTQaRw+4Q5Yo6TMCDsg3f7AnSbLhxOF532OJkn4j4GCF27P/YP7z3T9Ikvz7975l3vcJOPqNgIMV7jf/5O+z/+BkfvyMdfm3l2yd9312UdJvBBysYH/5xPP52iP/nIFG8psffEea82wwSXRw9B8BByvYJ7/6t0mS6y87K+84Y93rvnfcLkr6jICDFezpF17NqWtH8on3vfWY7x2fbC1DRVCcwbILAMq1/dofy9pVQ8d835PPj+Wsm/5kGSqCXk//x39zXL9OwMEK9qc3vjs/evraY77v5NVDeWlsGQqCAgk4WMHetmn0Db3vm594b9rt+WfkoIoEHHBMw4MDWT08UnYZsCA2mQBQSwIOgFoScADUkoADoJYEHAC11Oh0Os7fAebodDp5deLVJMmJQyfOe08cVJWAA6CWLFECUEsCDoBaEnAA1JKAA6CWBBwAtSTgAKglAQdALQk4AGpJwAFQSwIOgFoScADU0mDZBQDl6HQ62b9/f9llwBuydu3aBR/4LeBghdq7d29OPfXUssuAN2TPnj055ZRTFvRrBBysUMPDw0mSXbt2ZXR0tORq6mHfvn3ZsmWLz2mBpj6nU1+vCyHgYIWaWu4ZHR31l3HBfE6Ldzz3EdpkAkAtCTgAaknAwQo1MjKS7du3Z2RkpOxSasPntHiL+Zw2Op1OZwlqAoBS6eAAqCUBB0AtCTgAaknAAVBLAg5WqC984Qs566yzsmrVqlxyySV58MEHyy6pLyzk83bHHXek0Wj0/Fi1atUyVtu/7r333lx77bXZvHlzGo1GvvrVry749xBwsAJ9+ctfzrZt27J9+/Y8/PDDOf/883PNNddkz549ZZdWacfzeRsdHc2zzz7b/bFjx45lrLh/jY2N5fzzz88XvvCF4/49jAnACnTJJZfkne98Zz7/+c8nSdrtdrZs2ZKPfexjuemmm0qurroW+nm744478vGPfzwvv/zyMldaL41GI3/8x3+cD37wgwv6dTo4WGHGx8fz/e9/P1dffXX3tYGBgVx99dW5//77S6ys2o7383bgwIFs3bo1W7ZsyU/+5E/mscceW45yiYCDFWfv3r1ptVo57bTTel4/7bTTsnv37pKqqr7j+bz9yI/8SH7v934vX/va13LXXXel3W7nXe96V5555pnlKHnFc5sAwBK57LLLctlll3V//q53vStve9vb8sUvfjG/8Ru/UWJlK4MODlaYjRs3ptls5rnnnut5/bnnnsvpp59eUlXVV8TnbWhoKD/xEz+Rf/iHf1iKEplFwMEKMzw8nAsvvDD33HNP97V2u5177rmnp9ugVxGft1arlUcffTSbNm1aqjKZwRIlrEDbtm3LDTfckIsuuigXX3xxfvu3fztjY2P5xV/8xbJLq7Rjfd6uv/76vOlNb8ott9ySJPn1X//1XHrppTnnnHPy8ssv59Of/nR27NiRD33oQ2X+Z/SFAwcO9HS6Tz31VB555JFs2LAhZ5555hv6PQQcrEDXXXddnn/++XzqU5/K7t27c8EFF+Tuu++es4GCXsf6vO3cuTMDA9MLYy+99FI+/OEPZ/fu3Vm/fn0uvPDC3HfffXn7299e1n9C33jooYfy3ve+t/vzbdu2JUluuOGG3HHHHW/o9zAHB0AteQYHQC0JOABqScABUEsCDoBaEnAA1JKAA6CWBBwAtSTgAKglAQewSFdeeWU+/vGPl10Gswg4AGpJwAFQSwIOoACTk5P5lV/5laxbty4bN27MJz/5yTjqt1wCDqAAd955ZwYHB/Pggw/md37nd/K5z30ut99+e9llrWhuEwBYpCuvvDJ79uzJY489lkajkSS56aab8vWvfz1/93d/V3J1K5cODqAAl156aTfckuSyyy7LE088kVarVWJVK5uAA6CWBBxAAb773e/2/PyBBx7Iueeem2azWVJFCDiAAuzcuTPbtm3L448/nj/4gz/I7/7u7+bGG28su6wVbbDsAgDq4Prrr89rr72Wiy++OM1mMzfeeGM+8pGPlF3WimYXJQC1ZIkSgFoScADUkoADoJYEHAC1JOAAqCUBB0AtCTgAaknAAVBLAg6AWhJwANSSgAOglv4/GvwkomP8QMMAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 500x500 with 4 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "limits = [[i, j] for i, j in zip(prior_min, prior_max)]\n",
    "fig, ax = pairplot(\n",
    "    samples,\n",
    "    points=[list(theta_true.values())],\n",
    "    figsize=(5, 5),\n",
    "    limits=limits,\n",
    "    labels=[\"a\", \"b\"],\n",
    "    upper=\"kde\",\n",
    "    diag=\"kde\",\n",
    "    fig_kwargs=dict(\n",
    "        points_offdiag=dict(marker=\"*\", markersize=10),\n",
    "        points_colors=[\"g\"],\n",
    "    ),\n",
    ")\n",
    "ax[0, 0].tick_params(labelsize=14)\n",
    "ax[0, 0].margins(y=0)\n",
    "plt.tight_layout()\n",
    "plt.savefig(\"output/do_cpp.jpeg\", dpi=100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "language_info": {
   "name": "python"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}