{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# [Damped Oscillator - Numba/C++](https://github.com/Ziaeemehr/vbi_paper/blob/main/docs/examples/damp_oscillator.ipynb)\n",
    "\n",
    "<a href=\"https://colab.research.google.com/github/Ziaeemehr/vbi_paper/blob/main/docs/examples/damp_oscillator.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.sbi_inference import Inference\n",
    "from sklearn.preprocessing import StandardScaler\n",
    "from vbi.models.cpp.damp_oscillator import DO\n",
    "\n",
    "# change to numba version is compatible with current pipeline\n",
    "# from vbi.models.numba.damp_oscillator import DO "
   ]
  },
  {
   "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": 6,
   "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_cut\": 20,\n",
    "    \"output\": \"output\",\n",
    "    \"initial_state\": [0.5, 1.0],\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "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_cut': 20, 'output': 'output', 'initial_state': [0.5, 1.0]}\n"
     ]
    }
   ],
   "source": [
    "ode = DO(parameters)\n",
    "print(ode())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYYAAAEiCAYAAAD9DXUdAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAA9hAAAPYQGoP6dpAABZ0ElEQVR4nO2dd3wUdf7/X7M9PaQnpBA6iEAAUYqCqPgNiuX8nv2w3Z2eZ8XzTvTOwqno6deven7BO09Fz8bPAjZORVEQGy1Beg0khIT0bMv2z++Pz8zszu7sZjdtSXg/H484m9mZyWcSnNe+u8AYYyAIgiAIEU28F0AQBEGcWJAwEARBEApIGAiCIAgFJAwEQRCEAhIGgiAIQgEJA0EQBKGAhIEgCIJQQMJAEARBKNDFewF9jc/nw7Fjx5CSkgJBEOK9HIIgiD6BMQaLxYKCggJoNJFtgpNOGI4dO4aioqJ4L4MgCCIu1NTUoLCwMOIxJ50wpKSkAOC/nNTU1DivhiAIom8wm80oKiqSn4GROOmEQXIfpaamkjAQBHHSEY0LnYLPBEEQhAISBoIgCEIBCQNBEASh4KSLMRAEcfLg9XrhdrvjvYw+Qa/XQ6vV9si1SBgIghhwMMZQX1+Ptra2eC+lT0lPT0deXl63a7RIGAiCGHBIopCTk4PExMQBX8zKGIPdbkdDQwMAID8/v1vXI2HoJb7YWY/PdtbjzxeMRUaSId7LIYiTBq/XK4tCZmZmvJfTZyQkJAAAGhoakJOT0y23EglDL8AYw81vbAFjQElGEu48d0S8l0QQJw1STCExMTHOK+l7pHt2u93dEgbKSuoFats6wBh//fPRtriuhSBOVga6+0iNnrpnEoZe4ECDVX7dZHXGcSUEQRCxQ8LQCzRZXfLrY+2OOK6EIAgidkgYeoEWmzPgtQtM8isRBEFEwXPPPYfS0lIkJibikksuQXt7e5/+fBKGXqDF5i+o8foYbC5vHFdDEER/4v7778cLL7yA1157DRs2bEBFRQUeeeSRPl0DCUMvEGgxAECb3RXmSIIgCD+bNm3Ck08+iRUrVuCss87CpEmTcPPNN+OTTz7p03WQMPQCLTalELR3nBwl+QRBdI+nn34ac+bMwaRJk+R92dnZaGpq6tN1UB1DL0DCQBAnFowxdLj73qWboNdGnULqdDrx8ccf4+mnn1bs7+joQFpaWm8sLyxxFYb169fjqaeewpYtW1BXV4eVK1fikksuiXiO0+nE4sWL8cYbb6C+vh6FhYV44IEHcOONN/bNoqPA4vAovjeTMBBEXOlwezH2wc/7/OfuWnw+Eg3RPWa3bt2Kjo4O3HPPPfjjH/8o73e73Tj77LN7a4mqxFUYbDYbJkyYgBtuuAGXXXZZVOdcfvnlOH78OF5++WUMHz4cDQ0N8Hg8nZ/Yh9icfD3piXq02d2wOSn4TBBEZPbt2weTyYTt27cr9l900UWYMWNGn64lrsJQXl6O8vLyqI//7LPPsG7dOhw6dAgZGRkAgCFDhvTS6rqOlIWUk2JEm90NexxMWIIg/CTotdi1+Py4/NxoMZvNyMnJwfDhw+V91dXV2LNnT9QfnHuKfhVj+OijjzBlyhT87W9/w7///W8kJSXhoosuwl//+le5gVS8YYzJFkN2ihH7jlthd55YFg1BnGwIghC1SydeZGVlwWw2gzEmxyUee+wxzJs3D2PHju3TtZzYv6kgDh06hA0bNsBkMmHlypVoamrCrbfeipaWFrzyyiuq5zidTjid/vRRs9ncq2t0enzw+HhBW3ayEQBgpzoGgiA6Yc6cOXA4HHjiiSdw1VVX4a233sJHH32EjRs39vla+lW6qs/ngyAIePPNNzF16lTMmzcPzzzzDJYvX46Ojg7Vc5YsWYK0tDT5q6ioqFfXaAuwDrJkYSCLgSCIyOTm5mL58uVYtmwZxo4di++//x4bNmzo9WeWGv1KGPLz8zF48GBF6taYMWPAGMPRo0dVz1m0aBHa29vlr5qaml5doxRoTtBrkWziBhlZDARBRMMVV1yB6upq2O12fPLJJxg2bFhc1tGvhGHGjBk4duwYrFZ/99J9+/ZBo9GgsLBQ9Ryj0YjU1FTFV29iFS2GJKMOSQYSBoIg+h9xFQar1YrKykpUVlYCAKqqqlBZWYnq6moA/NP+ggUL5OOvvvpqZGZm4oYbbsCuXbuwfv163HvvvbjxxhtPmOCzTXQbJRu1SDDwjARyJREE0Z+IqzBs3rwZZWVlKCsrAwAsXLgQZWVlePDBBwEAdXV1skgAQHJyMtasWYO2tjZMmTIF11xzDebPn4/nn38+LutXQ2ExGCVhIIuBIIj+Q1yzkmbPnh2xJfXy5ctD9o0ePRpr1qzpxVV1D1uAMCTodYp9BEEQ/YF+FWPoD0gikEwWA0EQ/RQShh5GykpKNGjlqsd4NO8iCILoKiQMPUygxWAShcHp9sVzSQRBEDFBwtDDSH2SEg06GHX81+v0kMVAEET/gYShh3GIbiOTXiNbDA6yGAiC6EeQMPQwknWQoNfKFoPD442YfUUQBHEiQcLQw0jWgUmvhVG0GBgD3F4SBoIg+gckDD1Mh8vvSpIsBoBbDQRBEJ1x2223YebMmarvDRkyBI899livr6Fftd3uD0gCYAxwJQFiZpIpXqsiCKI/sGvXLixbtgzr169XfX/MmDFyC6HehIShh5GCz9IQcKNOA6fHJ+8nCCIOMAa47X3/c/WJgDh0JxqeeuopnHbaaWFHeWZkZPR6h2iAhKHHCYwxSFunxwenhzKTCCJuuO3A4wV9/3PvPwYYkqI61OPx4P3338df/vIXed/NN9+MqVOn4qabbgIAWCwWJCVFd73uQDGGHiYwXTVwSxYDQRCROHjwICwWC0499VQAfDDZu+++i+TkZPmYn3/+GWPGjOn1tZDF0MP4hYFbDEadWP1MFgNBxA99Iv/0Ho+fGyVtbW0AIAvB559/jtbWVhgMBgDAxo0bceTIEVxyySU9vcoQSBh6GNmVpJNcSWL1M1kMBBE/BCFql068KCkpgSAIePvtt5GUlIR77rkH8+bNw4cffoghQ4bg5ptvxpw5c3DWWWf1+lrIldTDSFlJCQb+qyWLgSCIaMjLy8Njjz2GN954A+Xl5Vi4cCGWLFmCdevWYebMmRg1ahTefffdPlkLWQw9jFTHYAyyGCjGQBBEZyxatAiLFi1S7KuqqurzdZDF0IMwxmTLgGIMBEH0V0gYepDAhz9lJREE0V8hYehBAh/+ZDEQBNFfIWHoQaSMJJ1GgF4rBp/JYiAIop9BwtCDdATVMABkMRAE0f8gYehBgqueA1+TxUAQfcvJOAOlp+45rsKwfv16zJ8/HwUFBRAEAatWrYr63O+++w46nQ4TJ07stfXFSnDVM+C3GGiKG0H0DXq9HgBgt8ehaV6cke5Z+h10lbjWMdhsNkyYMAE33HADLrvssqjPa29vx4IFC3DOOefg+PHjvbjC2AhuoMdfc+11ecliIIi+QKvVIj09HQ0NDQCAxMRECDF0OO2PMMZgt9vR0NCA9PR0aLXazk+KQFyFoby8HOXl5TGfd/PNN+Pqq6+GVquNycrobdRcSXKMweUF3B2APiGmazLG8Pu3tmJ3nQXLbzgNJZkndlk/QZwI5OXlAYAsDicL6enp8r13h35X+fzqq6/i4MGDeOONN/Doo4/GezkKZGHQBbqSuEhcVfMw8MT3wI3/AQZPjvqaFTVtWL29HgCw/PvDeGj+KT23YIIYoAiCgPz8fOTk5MDtdsd7OX2CXq/vtqUg0a+EYf/+/bjvvvvw7bffQqeLbulOpxNOp1P+3mw299byAvokBQiDXoM0WDHJ8jXfUfFmTMKw9Uir/Hrz4dYIRxIEEYxWq+2xh+XJRL/JSvJ6vbj66qvxyCOPYOTIkVGft2TJEqSlpclfRUVFvbZGKcZg1CmDz6OFgIlLTftiuuauY34hO9RoPSkzLQiC6Fv6jTBYLBZs3rwZt912G3Q6HXQ6HRYvXoxt27ZBp9Nh7dq1quctWrQI7e3t8ldvjsWTGugpYwwa5AoBn/Tbj8Z0zaNtHfJrm8uLJqure4skCILohH7jSkpNTcX27dsV+5YuXYq1a9fivffeQ2lpqep5RqMRRqOxL5You5KU6aoa5AQKg7kW8PkATXSaXNvaofi+ptWO7JS+uR+CIE5O4ioMVqsVBw4ckL+vqqpCZWUlMjIyUFxcjEWLFqG2thavv/46NBoNxo0bpzg/JycHJpMpZH+88KerBlgMeq3SYvC6gI4WICmr0+t5fQz1ZgcAIC/VhHqzA40WZydnEQRBdI+4upI2b96MsrIylJWVAQAWLlyIsrIyPPjggwCAuro6VFdXx3OJMSFNaUsIshgyhaCAt70lqus1WBzw+hh0GgHjBqcCAAkDQRC9TlwthtmzZ0cMpi5fvjzi+Q8//DAefvjhnl1UN1CvfNYgBUEVmPbmqK5X386thdxUE3JTTQCABhIGgiB6mX4TfO4PhGuilyIo4wTRCkOrnQeaM5IMclyBLAaCIHobEoYexJ+uGhhjCLAYBFEwOqJzJbXYeGHOoCQDspJJGAiC6BtIGHoQyZUUWOBm0GqQDNFiGDSEb6O1GGyixZCoR0aSAQDQ3kHpqgRB9C4kDD2IQ5r3rFNWPqcI3GJgMQpDi+hKGpRkQFoC75bYZj85yvsJgogfJAw9iMOlEmMIsBh8qYV8p9MS1fX8FoNfGNo7SBgIguhdSBh6EH+BW0CMAS4YBL7fk5wvHhhdv6YWm4rFQMJAEEQvQ8LQgzjU6hi8NgCAjwlwJYrtcKO1GAKyktITuTC4PD6aBkcQRK9CwtCDyFlJAcIgiCJghQkufQrfGaUwyBZDogHJRh20Gj5shOIMBEH0JiQMPUiHyqAeuLgI2JAAlyaR73NG50pqFQUgI8kAQRAC3EmUmUQQRO9BwtCDqFU+w82rlzuYAQ5tMt8XhcXg8zG0SVlJohtJDkCTxUAQRC/Sb7qr9gecoispMMYADxcGJ/SARhzLGYXFYHN54BO7haQmKIWBAtAEQfQmJAw9hNfH4PJK3VUDhYFXKjuhh0+QXEkWgDEgwoByi8MDANBrBbmSmlJWCYLoC8iV1EMEZgopYgweXsPghAF2KcbAfIDLFvF6VicXhhSTHoIoICkmruM28T2CIIjegIShh1AIg07FYmB6dDCjv19SJ3EGi4NbBclGv1EnCYNkTRAEQfQGJAw9hNQOw6DTQKMJcBHJMQYDnF4fYJRSViPHGaSHvyQGgF8krGQxEATRi5Aw9BByRpIu6FcqZiU5oOfBaSMfuNO5xcAf/oEWQ7JRr3ivq+ukAjmCICJBwtBDqKaqAkqLwRO9xRAYY5BINnXPYjjYaMUZS77CzCfXoqbF3vkJBEGclJAw9BDhhcEfY3B6vIBBTFl1RX4wSzGG1ABXUorkSnJ0LStp2TcH0WZ3o8nqwivfVXXpGgRBDHxIGHoIqR2GIiMJUNQxOD2+AGGInJUku5ICYwzdsBgYY/h2f6P8/dd7GmK+BkEQJwckDD1EZ64kBww8xiALgzXi9SIFn7sSY6huseO42T/97XCznSqoCYJQhYShh5AtBl2YGAPTw+X1xm4xGHsmxrC7jge7Tx2chuIMXk+x41h7zNchCGLgQ8LQQ0gWgzHEleSvfFZaDJ0VuPFP8ylqMYYuCMPBRm6hDMtOwqi8FMU+giCIQKglRg/hH9ITZDG4eeWzAwZ4PT4gURQGd3QWg0IYxAwlq8MDxphcER0Nhxr5zxuanSwHtg83UWYSQRChxNViWL9+PebPn4+CggIIgoBVq1ZFPP6DDz7Aeeedh+zsbKSmpmLatGn4/PPP+2axneAPPofJSoKUlSR2WI3SlZSiEnz2+BgPZMfAoSZuHQzNTkJJJhenw82R10AQxMlJXIXBZrNhwoQJeOGFF6I6fv369TjvvPOwevVqbNmyBWeffTbmz5+PioqKXl5p54QtcAuIMcSSlaRWx5Co18p998wxpqweaebWQWlWEoaQMBAEEYG4upLKy8tRXl4e9fHPPvus4vvHH38cH374IT7++GOUlZX18Opiw9lZHQMM8MYQY1DrlaTRCEg26GBxemB1eJCTEt3a3F6fPA0uL9WEVFFsjrZ2xOySIghi4NOvYww+nw8WiwUZGRlhj3E6nXA6/WmaZnN009NiReqVFFrHIMUY9PApXEmxp6sC3J1kcXpiCkA3W7koaDUCHxNq4mt1eXxotbuRkWSI+loEQQx8+nVW0v/8z//AZrPh8ssvD3vMkiVLkJaWJn8VFRX1ylo6rXyWWmJEYTG4PD45hpASkK4KcAtiiFCH9I3/C5iPRbW2RgtfQ1ayARqNAKNOi6xkLgZ17R1RXYMgiJOHfisMb7/9Nh5++GGsWLECOTk5YY9btGgR2tvb5a+amppeWU+nvZJiiDEEWgPJwRaDUYtl+udQ/POzwMd3RrW2RitfQ1ayUd6Xm2oCANS3O6K6BkEQJw/90pW0YsUK3HTTTXj33Xdx7rnnRjzWaDTCaDRGPKYnkLKSjGG7qxrAonQlSfGFJIMWWo3S/z9M24Axmmr+zf4vuMBIYhOGJgt3JWWn+H8P+Wkm7DxmRr2ZhIEgCCX9zmJ4++23cf311+Ott97CBRdcEO/lyHTeXTW4wC18DYFanySJkTii3HGsstO1NVq5Kyk7wGLIS+uexXCgwYLfv7UV//f1Afik4dQEQQwI4moxWK1WHDhwQP6+qqoKlZWVyMjIQHFxMRYtWoTa2lq8/vrrALgoLFiwAM899xzOOOMM1NfXAwASEhKQlpYWl3uQ8Aefw9cxsChdSf7Asz7kvSG+auWOpr3AkBkR1ybFGAIthrxuuJK8Pobfvr4Fh5ps+PTnOmSnGHH5lN6J3RAE0ffE1WLYvHkzysrK5FTThQsXoqysDA8++CAAoK6uDtXV/gfhP/7xD3g8Hvz+979Hfn6+/HXnndH52nsTv8UQqY7BC+hFYfB0AD71gTlqqaoSeZ6ggHPzwU7XJlkMWQqLIQEAuuRK+mZvAw41+YVt+XeHY74GQRAnLnG1GGbPng3Gwrshli9frvj+m2++6d0FdQOpjiEh0GLweQEff8g7YAACLQaAWw2m1JBr+YvbQv886d4mAMCRxHEose8AWg51ujY1iyFHfC29FwtfiS27L55YgE9+rsOuOjNqWuwoEpvzEQTRv+l3MYYTFdWWGB7/p3E5xqAzAoJ4TBh3kuRKSlVxJaW6mwEAB0yn8B1tnWdZNakIQ3YXhYExhvX7+FyHiycWYGJROgDguwNNMV2HIIgTFxKGHsLfRC/gV+rxP3R5HYMXDOi0X1IkV1Kiiz+A92lHiAd3Xsvgr2PwC4P0usXugscbfd+lerMDR1s7oNUIOL00EzOGZwEAfqpqifoaBEGc2JAw9BBy2+3AeQxiZ1Wm0cMHDXyMN8DrbFiPJZwryWWDwcPnKuzAcL7P3qwQILV1SdcLtBgykgzQCABjkNtlRMPOWl45Pjw7GUlGHSYW8aD/9lqa7UAQAwUShh4ioitJ538gR1PkFjZd1cKzsOzMiCp3JqAVr2upC7suyVow6DSK+dFajYBM0WpoiMGdtKuOC8MpBTw2Mm4wF4aDjVbYujAngiCIEw8Shh5CNStJ+iSvM8m7nO6AKW5u9VoGa7h0VRv37TeyNFhdXiA1n+83RxCGgBqG4GZ5kjupyRq9MOwUp76NFYUhJ8WEvFQTGPOLRqy4PD78fLQNHS71LC2CIPoWEoYegDH/fAQ1i0HQmWDQ8l81txgiVz9LMYaU4BiDgz94zUjkmUspBeIJ4eMMUuA5KyW0+rsrAWjp4T82359NJVkPe7ogDO12Ny78+7e46IXvcO4z61DbRr2bCCLekDD0AIFDc1RdSXqT3CojGldS2HRVJ3/wWlgiLA43WEqueEJD2LWpVT1LSI30GqO0GBxuL2pa+INbGg8KAMNzudDtb4h9VOgTn+3BvuP8vNq2Djz04Y6Yr0EQRM9CwtADSG4kIGhQjxxjMMmzoPkUt+hiDCGuJAd341iQCLeXwZvAM4JgC58qqlbDICHtk3opdYY07CfFpFO06h4hDoY4EKMwNFudeH/rUQDAkl+cCq1GwJe7G7CnvndaoxMEER0kDD2AFHjWaQTotGoxBqOcraTslxTOlRQm+CxZDOCFZE7DIL5fjD2o0SRbDKEzFyQrIlqLQZr4NiQzSRGvGJ7TNYvhw8pjcHl8OHVwGq48rQjnjuFdcv/fpqMxXYcgiJ6FhKEHCNtAT0xXhS4hJleSHGMIFgYxxtCh4ec79KIw2JvDri0ai6HREl1bjCOSMGQpu7lKwtBocaLdHv3I0a/2HAcAXFI2GIIg4JeTeb+l1dvrIlbER8Ll8eG42dHl8wmCIGHoEVSL2wCFxWDQqbiSnKGfsBlj/hhDcPBZtBhcWv4gtknC0FVXkpyVFJ0rqaqJu5KGZCpbXyQbdcgXu7UeaLREdS2Lw42fDvGiuHNGc0th5ogsJBq0qDc7sPNY7O6kDfubMG3JVzj98a9w6dLv0UAtxQmiS5Aw9ABSmqWiuA0IijGouJLcoRaD3eWF1MU6NMbAH5ZuPRcGsyad74/gSlJroCcRa1bSkQBXUjCS1RBtnOG7A03w+BiGZifJFohJr8VMsZJ67Z7wAfVwa/vtvzejWSzWq6xpw42vbYI7hqpugiA4JAw9gL+4LbzFoHQliRk9KhaDFF/QaoTQ64kWg0fPzzdrxJRRu7rFwBhTHdIjIYlFe4ebWzKdIAWfh2SFNsuTAtD7j0cnDJsOtwIApg/LVOw/R4wzfBWjMCz+eBfsLi9OL83AlwvPQnqiHjtqzXhlQ1VM1wnE52MxtQshiIFCzMLw5Zdfhn3vH//4R7cW01/xu5KCLQYxxqAPjDFEzkqyOv3xheCCNMliYEb+EG4VxBkUjnbAG+rbt7m86BDjH2oWQ1qCHmO0R3GH9gO0HYvcvtvh9uKYOB+6JJLF0BidMGw5woVhSkmGYv+skVwYfj7ahtYoW3XsPNaOr/Y0QCPw7KbhOSm4v3wMAODFdQdjrsh2e314+vO9mLj4C4z6y2f49WubqL6COKmIWRguuOAC3HPPPXC5/P/TNjY2Yv78+Vi0aFGPLq6/4Aw7vU0lK8njA4zhm+iZHeFbbksWAzNwS6HZm+jv1KoSgJZcREkGLZJUGvJp3Da8rf8rFurfQ9oHVwO+8J+Oa1rsYIzHPTKTQjOc5MykKCwGh9srV1BPLhmkeC8vzYTReSlgDNgQZcfWl7/lVsGF4wswNJuv4xeTBqM0Kwmtdjf+/eORSKcr8PoY7nqnEi98fQBmhwdeH8OXuxtw8Qvf4XBT+OFKkWCMob7dgQYKihP9hJiFYf369fj4449x2mmnYefOnfj0008xbtw4WK1WbNu2rTfWeMIT3pWkUscQ2BLDFRqolVNVjaEttyWLQRBnOFicPiBR/MStEmeQUlXVqp4BAHs+QTr4Gkxt+4Ej36kfB+Cw6EYqyUoMtWQAjBCFobatA3ZX5E/o22vb4fYyZKcYUTgoIeT9s0ZmAwDW7QsfO5GwOT34zw7eQ+q66UPk/TqtBr8/mzcafHlDlaLWJBJv/XQEn26vg0GrwXNXTsTqO87E6LwUNFmduP7VjTA7os+6AoAPK2tx9tPf4IwlX2Hq41/h3GfW4cPKWhII4oQmZmE4/fTTUVFRgfHjx2Py5Mm49NJLcc8992Dt2rUoKjo5xzvK6arBwWd3gDCIriSXNyDGoOZKisJi0Cak82OdHiCJP0TVMpPkjCQVNxIA4NA65ff7v1A/DpA/Lau5kQBgUJJBtiQONUb+ZC25kSYXD1IVmVkBwtDZA/TznfXocHsxJDMRk4rTFe9dNKEA+WkmNFqcWFlRG/E6ABfSpz7fCwC4f95oXDxxMMYWpOL1m6ZicHoCDjfb8dCHOzu9DsAtj0Uf/Iw736nE4WY7tBoBGgE42GjDne9U4va3K2LqDdVocWLJf3bjnP/5BiMeWI1TH/4c1/zrR7y9sTpq0Qukw+XFtpo2bNjfhM2HW2Lql0UMfLo0wW3v3r3YtGkTCgsLcezYMezZswd2ux1JSeoPjYFO2DqGQItBrcBNNfgcpk8SY4CTf7rXJqYBaOMikigGbyMJQziL4dhWAMAn3jNwofZH4OimsPcoFbeVhhEGABiWk4zmqhbsb7DIXVfVkIUhyI0kMWXIICTotWi0OLG7ziI37FNDeuBLtRCBGHQa3DSzFI9+uhv/XH8Il08pglYTKkQST/5nD8wOD04pSMWvpg2R9+ekmPD8VRPxyxd/wMqKWswelY2LJw4Oex3GGP68ajve3lgDjQDccc4I/ObMofAxhle/O4znv9qPT36uw5FmO15aMAV5aaaw1/L5GFZsrsGS1btlNyMAuL0efHegGd8daMYza/bh5rOG4prTS5Bg0Ia9ls3pwertdVhZUYsfDzXL2W8Sw7KTcNbIbMwZnYOppRmhWXYBuDw+7Kk3Y1tNG7YdbcfuOjNcHh8MOg0K0hNQNCgRw3KSMDw7GcNykpGZZFD9EODzMTRanTjSbEd1C/+yOjzwMQaDToMM8QNHdooROSkm5KQakZFogCbo78gYQ5vdjSarE802FywODywON2xODzQaAXqtBkadBokGHRINWiSK7lX+mm+NOg08Pgab0wOr0wOb0wu7ywOnx8e/3HyeilGngVGnhVGv8b/WaWDQaWDQ8q3Hx+Dy+ODy+vjW44PH54PXx+DxMb71MjAwaAQBAgCN+OEB4FtBELcQIAiAIAAaQYBGEDAoSY+clPD/brpLzMLwxBNP4KGHHsJvf/tbPPXUUzh48CCuvfZajB8/Hm+88QamTZvWG+s8oXGIvZKMUWclRQo+h5/FAMYFyJCUDqBNtBjEthgqMYamCKmqcFqBRv7p+E3vOVwY6rbxcaSa0AeClJFUkhl+fOeInGRsrGqJmLLKGMNWURgmhREGo06LacMysXZPA9bvbwwrDA1mhzw57tIy9Qf1VVOL8fe1B1DVZMMXO+tRfmq+6nFbjrTg3S284nrxxeNCBGRySQZunzMCz321Hw+s3IGJRemq1hNjDI9+uhtvb6yBIADPXlmGiyYUyO/fcc4ITC3NwO/e2ILtte246IUNeGnBFEwQJ+EFsv+4Bfev3C5ncI0bnIrfzRqO8YVpsDo9WL+vEa99fxjH2h149NPdeHHdQfz2rKG4+vQSecgTYww7j5mxYlMNVlbUyv++AN4rKyvZCIvDg2PtHTjYaMPBRhte/e4wkgxanDkiG2cMzRDngzM0WJzYW2/Bjtp27K6zcOtXBbUalPREPYZlJ6NwUAK8PgaLw4OaVjtqWzsUvcaiQacRkJVsRKJRC7fXB6fbhxabi8866QYaASFieaJyw4wheGj+Kb12/ZiF4bnnnsOqVatQXl4OADjllFOwceNG3H///Zg9ezaczpPPJI3KYgjslSRmFcFt4wFfjV9QzJ20w4CghSmR+/MtDg+QJQlDjBZD8wEADA7DIPzkGAOHYILJbQea9gE5Y0IOr2pSr3oOJJoAdHWLHc02FwxaDcYNDm8JzBqZjbV7GrBubyNumTVM9ZgPK4/Bx7jlEc7FlWTUYcG0Evx97QG8uO4g/mtcXsgnV6+P4S+ruIvo8imFYS2Z2+cMx/cHm7DpcCtuf7sC790yXS5clHj2y/14WUyRffKy8QpRkDhjaCY+um0mbnptE/Ydt+Lyf/yAe88fhWvPKIFJr0WD2YGXvj2EV787DI+PIdGgxT1zR+G6aSWKlitj8lNxw4xSrKw4ihe+PoCalg48vnoPnvp8L0bnpcKo0+Bws13hJhqSmYj/nlyIiyYMRnGAyLd3uPH9gSZ8s7cRa/c2oNHixGc76/HZznrV3wXAH/bjC9MxoTAN4wanIcWkg8PtRW1rB6qa7DjYaMXBRitq2zrQZndjy5FW2VoMRKsRUJBuQnFGIooGJWJQkgFaQYDD7UWLzYVGqxONFv7VLApAfZjixbQEPTKTDEhJ0CPVxC0BxnimmdPjg93FrQCbk2fs2ZweWZgCRcGg0yDFqINJr4VJz60Cg04DQeDWErcivHC4/RaB0+MNERaDTgOjaEXotAJ0Gg20GgE6jQCNhlsKDICPMYDxrY8BDEzOBfExBhbwHsBUpzv2JDFfffv27ciSHkYier0eTz31FC688MIeW1h/Qg4+hytw05uUWUmGgIeY2+YXCgS2w1AvboMpFcnie1an2+9KipCVpCoMrfzh5UgZAp9ZgyPaEozy7OVWRJAwOD3+VFW14jYJuZlehJRV6cEwbnBqRFeFFGfYfKQFNqdHNasq0I0UieumD8E/1x/CtqPt+OFQM6YPU/77fePHI9hVZ0aqSYc//dfosNfRaTV47soyzHv+W/x8tB33vLsN/3v5BOi0GjDG8NxX+/HcV/sBAA/PH4vLp4SPuRVlJOL9303Hne9UYu2eBjz66W48/cVeZCYZUdfeIT9gzh2Ti0cuPgWD00OD9AB/8FxxWjF+MakQKytq8eI3B3GoyaaYqJeg1+Ls0dm45vQSTB+WqerSSUvQo/zUfJSfmg+fj2HHsXas3dOA3XVmHDc7odUIGJSox4jcFIzJT8XEwnQUZSSoXiuYDpcXh5qsONhoQ317B/RaDZKMOhSmJ6BwUCLy003Qa6MLd7q9PjRZnWgwO+Fwe6EX3TcZSQZkJhsi/psKh8frg93tRYfLC4O4tmDBj+VaLq8PWo0Ag1YT1e/nRCRmYQgWhUBmzZrVrcX0V1SH9ABBMQYpK8nHB/cIWu4acimFwSpnJYWxGIypsjVhcQS4ktRiDBFabqPlEADAlz4EqAUO+AowCnuBpv0hh9a0dIAxnvaapdKMT0KyGI4022V/czCdxRckhmQloTgjEdUtdvxwsBnnjs1VvL+33oJddWbotQIuDOMekshKNuLyKUX4949HsOybgwphqG93yAHne88fJU+1C0dBegKevWIifv3aZny87RiOttoxb1w+vt7bgO8PcnH+w9yRuH5GacTrAFz8/7VgCt7dUoPnvzqA2rYOuV5icskg3DZnOM4eldPpdQBAr9Xg8ilF+OXkQhxptmN/gxVOjxf5aQmdinAwGo2A8YXpGF+YHvU5kUgwaHFKQRpOKQgfd4oWvVaD/LQE5KepC2VX0Gk1SNVqkBr8YayL19JFKXInMnG9g/Xr12P+/PkoKCiAIAhYtWpVp+esW7cOkydPhslkwtChQ/Hiiy/2/kI7wRm2wE0txuDlUSRpWE9QAFpKV00N00APplT5PavT409XVYsxRBjSgxZuMRhyeErnLpf4AGoOFYbA5nmRPgHlphqRbNTB62NysDqYrdVtAIBJxZGFAVBmJwUjWQtnj8rBIJW6imB+c+ZQaATg2/1N2FbD1+DzMdy/cjusTg8mFqXj6tNLOr0OAMwelYOl10xCkkGLiuo2PLZ6N74/2AyDVoPHLz0Vt80ZEdV1AP4QvuK0Ymz409n4cuEsvP+7adj4wDl4/3fToxaFQARBwJCsJJw3NhcXji/A5JJBXfoUTZzcxFUYbDYbJkyYgBdeeCGq46uqqjBv3jyceeaZqKiowP3334877rgD77//fi+vNDJh6xgCu6vqA1xJQNjW21JwMDTGILoGjGlyjQPPSlIPPjPG/BaDqivpMAAgIWcYNAJwwCd+6laxGOT4QgQ3EsAfSpF6JlmdHuwVZy2ECzwHIgnDN/saFGmrXh/DR5VcGMIFnYMpzkzEfNHff8c7FdhW04YHP9qBtXsaYNBq5HkQ0TL3lDysWTgLt509HPNOzcPvZg/DV/fMwtWnF0d9jUCk393kkoxezTYhiGjo3QhGJ5SXl8tB7Gh48cUXUVxcjGeffRYAMGbMGGzevBlPP/00Lrvssl5aZeeEDz6HyUoCePWzBSGZSf501UgxBtGV5PQASWJwM8iV1N7hhtvLH6aq7p9WXg2szRyK7BQLDlnE6zQf4KmxAZZBNBlJEsNzklFZ06YqDNtq2uBjwOD0BOSmdv7wmzYsE3qtgJqWDhxutqNUDHyv3dOAY+0OpCXocfbo6D9VPzT/FGw50oojzXZc/H/+Yr4lvzgVY/LDB8LDUZCegD+cPyrm8wjiRKdfOcN++OEHzJ07V7Hv/PPPx+bNm+F2q1ekOp1OmM1mxVeXaD0MvPHfwP9bEPJW2AI3tRiD1KwujMVg6WSsJ4wpcvzB5fHBaUjn+ztaFC0tpEyUVJMu1JXAGGAVs01S8pCbasIRlgsGgf+cIJE5HKGrqv9eXcDK3+HhA7/EVdqvsFtl/rOUploWVIgWjiSjDqcN4a6yDyv9BWqvfsfdYFeeVhQqxhHISDLg/908DbNHZUOnETA0Kwn/+NVkXDa5MOprEMTJQL8Shvr6euTmKoOQubm58Hg8aGpS76uzZMkSpKWlyV9drs72uIADa0KrheF3JUVVx+CWXElSvyT1GEOIK8kREHwOCEzbtOn8BfMBjjZ5f0OkjKSOVsAr9rpKzkVuqgku6GE3ilXUbdWKw/1dVSMIw6aXgG1vIdl5HI/qXoH5yM8hh2ypFusXoogvSFw5lbtmXv/hCGxODzbsb8L3B5uh1Qj41bToYgKBFKQnYPkNU7H/sXKs/cNsnH9KXszXIIiBTr8SBgAhwU/J9xwuKLpo0SK0t7fLXzU1NV37wXrR9eEJzZ+2ixZDoiHoYa7orhocYwgXfOaWT0iGhFj1DFMqtBoBiWKFq9UtAEYx2yMgzhAxVdUqtrQ2pQN6E3JT+TGtBjHO0HZYPtTl8eFoq/qAHhnGgJ/8SQBageHCjg8Vg3LcXh82VvHBPFNLM0IuEY554/JQnJGIFpsLt7yxBX98j/fj+tUZJSgc1LlrKxz9NY2QIPqCfiUMeXl5qK9XFtw0NDRAp9MhMzNT9Ryj0YjU1FTFV5fQielxHkdIF1KHSxKGKLOSANXqZ7fXJ1sfkdJVAb+ryexwA0mhbTGkqWzZaoFMyY2UzK2vPNHff1wjWmMBFsPRVjt8jOfCh22tUb+dn6NLAK55DwBQrt2IbUf866msaYPd5cWgRD3GxuDP12k1eOKyU+WMomPtDgzNSsLd542M+hoEQcRGvxKGadOmYc2aNYp9X3zxBaZMmQK9vvs5yBHRBzxgg6wGu5u7f6KrfA4IPgMKYbAG9MIJ60oSO6tKwmF2BBa5+R/EDeIcZ9XAs4XPWkYKFwIpEFzNRFdSq79NdWDgOeyn7APi32TYHGDYHFi1aUgXbGja+Y18iNS6YvrwrJA+N50xfVgW3rjpdJw3NhdXTS3GOzefgbSEXv57E8RJTFyzkqxWKw4cOCB/X1VVhcrKSmRkZKC4uBiLFi1CbW0tXn/9dQDALbfcghdeeAELFy7Eb37zG/zwww94+eWX8fbbb/f+YnUBBTUeB2DwuzE6XPxhr7AYvB7AJz7og5voAQExBn/rbSlV1aTXhFaCBlkM6YkGADa0292qKauNZm6tqKY+yhYD969LVbX7naKLp80vDFKqammk+MLRzXw7ZCag0aI5fxaSj34E4fB6ANcCAL7czcXozOHhCyQjMX14FqZ38VyCIGIjrhbD5s2bUVZWhrKyMgDAwoULUVZWhgcffBAAUFdXh+pqv1ujtLQUq1evxjfffIOJEyfir3/9K55//vm+SVXV6gCNqKNu5TSvDnH+QEKgxRBoVahmJYVaDOZw7TAAPqUNAEw8njAokR/T1qHuSpJqGHIixRhEi6Eog4vcdls63x/gSjrUxGMgYQPPjPm7shZOAQBknnI2AKDUvg117R2oarJhR60ZWo2AudEGexmLODiIIIjeI64Ww+zZsyP221++fHnIvlmzZmHr1q29uKoI6BL4J/yAhz5jTB6fqbAYPM6A80ww6vg5IQVuAcFneRaDWoMsVYsBaLW7AlxJLfLhDZLFkKoiDBZljKEgPQE6jYAqbxb/F9FWLTf321fP1zcyNzn0OgC3LmyNgEYP5I3nlx05C/gcmCgcwKubD6HFyUVxxvAsZERRpYy6n4H3bgDaaoCZdwGzFynqKgiC6F36VYwh7khxhgCLwenxyQ3PFL3wJfHQGgCNJjTGoFLHYIk0pCcoxpAu+tjbFK6k0BiDelaSGGMQXUlajYDCQQmoYxlggpanslrrwRjDvgbu6pIa5IVQJ07tyxvn//1kDIXDmAWj4MH6rz/Hv77lfZluCJiwFhaXDXjnGl5o53UC654EKt/s/DyCIHoMEoZYCMxMEgmcwqXqStLxh6XkSuIDOnz+xnmBwuDkrqSQwDNjIRaD1B+oLdBiEF1JLo8PrXZ+LdUYg2QxpPhrQoozk+CFFjaT6Oppq0aj1Yk2uxsawd8gLwRxpgOyAzqyCgIMJVMBAGPYQfgYMGN4JmaPyla/RiBbXwfaq4G0YuD0W/i+NQ+FuO8Igug9SBhiQcVikNxIhuCuirIw8E/sgdXH4Yb1+F1JQTEGd4c/kC1aDFJWTqvdHTKsR6p61msFORahQIoxJPv9/cUZXPSa9eK+1iPyXIXijMTwFcayMCjTRzWFkwAAF2cfx2/OLMXSayZ3XjvAGPDTP/jrMxcCcx/lAmFvAn5eEflcgiB6DBKGWNCFFrnZXZ233AagaEHNhUEt+NzZkB6NfN4gMcbAs5KUMxkaAmY9hzyM3R3+hnzJ/j5DJRlcqGoh7ms7gj31ohspN4wbCQCaRGHICuoZVMATCsZrDuGBC8ZGl17acojPidDogVN/CWj1wOk38/c2/avz89XY+x/gX+cCL54JbFnOxYcgiIiQMMSCXnQlBVgMjnBVz26lMGg1AvRa/pB2eryqlc9hx3pKVc/GFDkIm54oWQyuUGEQK46z1RrVSW4knUnOcAKAUXn84b/HIaasth7Gz0fbAACnhpvf7PMBTWK6cXaQMORzYUDLQaCjTf38YA6u5dviM/x1HhOu4tlg9duBhj3RXUdi5yrg7St51lT9z8DHdwLfLIntGgRxEkLCEAv60BiDZDGEDGEPshgAKGsZVILP7R1h2mE4lPEFIFAYAlxJbjvgsisshhBkN1KuItNHmqu8TUpZbT2CSnFugdo8Yr7gat72Q2sA0oP6FiVl+vfVVaqfH4wkDMPOVl5n+Ln89fZ3o7sOAFgbgY9u568nXgPMuo+/Xvek/+fESkcrt0AOrvULP0EMQEgYYkEXajFIMYaECC23JRStt+XK51BhCHG7yLMY/MIgu5I6XGD6JP5wBgB7kywMqqmqQe0wJLKSjchNNaLGxwPE3tbDctXzxHCTvKTZDZnDeZ1HMKI7Cccq1M8PxOsGqtbz18PmKN879Zd8u/3d6F1Ba//KXXD5E4GL/g6cvQg47Tf8vY/ujP3BvuU14H/HcQvk35cCz40H9n4W2zUIop9AwhALKo305OK2cBaD3l8xLQmDw+31WwxeF+/cCsAcThiCUlUBv8Xg9jLY3D5F9bPUQE+1uC2oHUYgpxSkoYbxGIPGXAs9PBiRk4w0tQA24A88Z4WZWDaYB6BRG0XdydFNXCQTMoC8Ccr3Rs0D9Em8ZiKaa1mOA5Vv8df/9QSgEf825z0CpBRwSyeg6V+n/Pgi8PEdfH2DSnnQ3nqci8SW16K/jsSR74FVvwdeOgd4/RJgw7OKGhSCiDckDLEQwWIIbaCnzEoCgESxcM3uCogxAICbB6AliyE9+EHsDHUlJei1ckC71RaYstqM+na+vjy1GENQO4xApg/LRCPS4BSMEMBQIDTJU9RUCRd4lpAshmge5oFuJE3QP0tDIjB6Hn+9473Or7VlOeBzA4VTgZJpAddJAs75C3/97f+ozskO4fAG4PNF/PWZ9wC3bwXu3AZMug4A43GL3R93fh2AWykf3Ay8Wg5UvgHUbgYOfQ18+RDw3ITYg+MdrcC6p4BXyoG/TwZemgN8eg9w4KuuVY1b6oGajfxvcXQLvz5xUhLXyud+h4rF4M9K6jzGkCSKh83p4Rk3WiMv4nJagYRB4V1JKhaDIAjISDSg3uxAs82FoiR/ALq2jZ8/eJDKwHRreIth9qhsPPqpgCPeLIzU1KJIaMS5Y0OPk2ncx7fBgWeJgjIAAmA+yh86KRHaYcjCMEf9/XGXcVfSjg94GqsmTPqsxwVsfpm/ljKaAhl/JfDjMh6MXvckMO+p8GuytwDv/4bPuph4DTDnLzwuozEB85/jWWJbXgXe/zWw4COg+PTw17I18cK9mh8BQQuUXcNjJ9YGYPOrQMNOLjI73gcuegEYFGHWhLUB+OH/gE0vK3ptAQBqt/AMrrRiYNICoOxaIDVf/To+H3BsK4+b7PsMOL4j9Jj0YqDwNKDodL7NO5X/2w2EMT7Iqq4SOL4LaNgFNB/kH2g8Tv7vNjGLuxyzhgOZI/j9pZcACen+63ic/PdkPc7/vVjrufVnqePV9YKGu0yNKdwVmpwjbnN5IoXexP+f8rl5tp/TytfgNPP/h5xmnsjhtAAQ+Ic2fQL/f9SQyK1SfQL/AKFP5NfzefkHQY+Db70u/trjFLdSIauRX0dn4FuNjv+eNDqeZafV8fUzBoAFbBHwPYLeY0HvgV9D0PDfX05A7VAPQ8IQC2oWQxQttyWSRIvBJrqfYEgCOpxyymr4GEOoxQAAualG1JsdPAtJtBiYrRG1rTyLqCBdRRgkV1Jy6AN/WHYyppQMQvWxHIxELU5Lbcfp4WYnMBZgMYRpgW1M4f94G3bxRntjLlQ/zt7ityqGnq1+zLA5/H9+az13xZSeqX7crg/5gyUlHxh7cej7Gg0XltcvAja/Akz9bXhX2Cd3A5Zj/EE27yllWw5BAOY9zR9a+z4D3rocuOE/QO7Y0Os07gPe+iV/eBrTgCteB4bO9r8/5Ubu2vrqrzzOsmw6cN5iYPINSuuprQb44QVuWUgPpNxxXAAHlfJ6j6r1XFzaq4GvH+VZWKPKgZH/5b/Plip+3MGv/B8UAP7ASSsEDCl86JO5lrdHaavm1wT4/wMFZUByNhdMy3GgaZ9iSFQI9iaeinx0Y+h7OhN/2DNfyNAqIgKn/w4of6LXLk/CEAtqBW6uMMFn6ZiArqxSSqtkZcCYzEdyuqzw+VhMMQZAmrXQzoPNYozBaWmCzVUKAChIU7MYwruSBEHAXy8Zhx3/KgC8FbhyJAtflGZrEl0NQvgHKwAMnsyFoTaCMFStA8CA7NFA2mD1Y3RGYMxFQMW/+UMqnDD8tIxvp9wU+slWYugs/qDc9xnw5cPAlSotN3a8D+xaxT/xXfaSPyYUiFYH/PcrPBhd8xPfLlil/CR34Evg3Rt5AkF6CXDNu6EWlkYLTPs9X9OHvweqfwA+XcitgtHz+KfXYxX8Wkx0EQ2eApz1B35O4N/olEuB8x8Hdn3ErZnqH4A9n/AvNQwpwPBzuHiMmAskBnwQ6GjjVkDNJn5/RzdxAaj+XuV3YQByxgL543kVfPZI/mFFa+QfbCz1QPN+nrDQfIC3drc3KT91A9yaSs7lFm1Kvvg6D0jK5vfpcfKGktYGLmrS1mnxX0tr4K5aQyIXYmMK/3/HmMI/XBmTAQj+490dckYf39rEfR38b6wz8S/JutAZxS+TP+nD6+JrkywKr5tbGz63+Noj/u0E8e8lAAIg/idgX+AWyn2BVkR6FydRRgkJQyzILTFUspLCBp8DXEnGAFcSEFDLYIHV5ZF7LqXGYDEAYt2CmLLa0co/AWYkGULXBIR0Vg1mTH4qRp87A/j8U+R66lSPAeC3FtKLFQH2EAqn8Ie51JpbjYNf8204N5LEuMv4tXZ9yD/BBz/4j27mrhStAZh8feRrnbcY2L+GPzAPfaP8BN9ey331AHDmH/yxEjUMScDVK4BXL+DuoJfOAabfDuSewl00294GwLgr5oo3+SftcGQOA65fza2Hrx/nNSDf/115TOksXhVeOit8Y0F9AjDhCv7VsJtXjdds4hYA8wFpRTwxYNgcoGQGd3+okZDOfy/S78bn4w/3YxX+2prkXP5vIGds+OuEw2XjHzCkqv7EDD5VkBomxh0ShliQLYbQGEP44HOgMIiuJKdoMUgFZk4zr2AGz1wKiVfILbeVwiD1QWqwOIFB/JOe28If/AXpKoFnn5f7agFVi0FCGDSEvwiYyxCC3AojTHxBYjBvxY1jFfznB8cGGIteGIacyT852hr5w3zEecr3pUyjcf8d+QEsrXvKjXxW9Qc3Azev459MnRZgxTXcGsqfwD+Vd0bCIOC6j3lH2Kp1wLogE3/yDUD5kwq3Ylg0GmDarcCkX/Gg9rFKHofKGMbvN1a/cs4Y4NyHYzsn0tqyR3X+N48WQ5K6JUbEHRKGWJBbYkTjSgoffLZLMQbJAnCYw8cXgACLQVmBLFkMx80O2ZXEbLz6WdWNZGvknxgFjb8oTg1JGFojCEOTGHgOF1+QyBnDg3ouKxeTYB9880HuD9cagJLpka+l1XFXycZ/AtveUQqDuQ7YuZK/Vgs6q3HeI9zX3rQXeOV87n76eQUPwiZkAJe/Ht4dFUxSJvCrlXwN29/jv+vsUdxyKZoa3TUCMaYAE6/mXwTRx1C6aizILTEC6hhkV1KQxkriESAMUoxBan0hWwCO9vCpqoDfbA+2GCRXksUpB5+1Di4MqhlJUjuMpOzwWT2Av2K5o8Uf3wgmWotBo/W7YqSBPoEc/Ipvi8+I7tPjhKv4dudKpXB9+z/cJVE8HSiY2Pl1AP7zrnqbu0JaDwNr/sJFITETuPY9v0BGi0YLnPrfwNXvAL/5CrhkaddEgSDiDAlDLERoohe28lklxiAHnwNdSZEsBpWWGIDflXTc7JQtAKOrDQAwVG3imjyHIUIKKsCDc1JdRDh3UrQWA+CvJZAqmwPZL82LPqfz6wDcNz70bIB5uR8e4CmSW5bz12ffH911JDKHAb9dx1tmjLoAmHk3cOuPPGhOECcp5EqKhYhN9MJlJanFGIJdSe1RupKUXU6lWc1NVic6dOlIAJDks0ADH4Zmq8xPkOcwRKgnkBhUypvyNR/gueuK9Vh4IBOIThiGzQHWP8WLucTJcAB4Fsjhb/nrked3fh2JOX/mMYaf3+GxhD2refbHyPLw2UqRSMzgLTMIggBAFkNsqFoM/CEfXYFbUB2DZDEExBhCMpLE9/nxSothUJJBFpLDHdytpAHDIFgwNLsbFgPAU0cBfxFbIJK1kJStTG8MR+FpPAPL3swLyyQOf8t/T2lF/p8XDYVTgFl/5K+//zvP3kkr4kVnBEF0GxKGWFCxGDrcPK88bFZSQCpnolz5LLmSxAd9JFeSx8mzUoAQVxIADBFdRodbnPAa0wEA+Xo7clUnt4npp9FYDFLsoFGl1bUkFuFaYQSj1QOlZ/HX+wIaz+1cxbcj5saeojh7ETD/eZ62OWkBcNMXYVNwCYKIDRKGWIhgMYS6kkIthmS5V1KoK6nNHq64TUxVhaAqDKWZiQCAqmYbOvTpAIDxGW5oNCoP2ghVzyHIFsPe0PekGobsKNxIElIV8rZ3eIqqywbs/ojvG3959NeREARg8nXAdR/x7qmpBbFfgyAIVUgYYkHFYpDiBVL8QEYtKymkjiGdbx3taLFxqyAzKahIKDDwHNxcDkBpFo8lHDhuRSu4cIwf5FZfvzWGGIP00G/eD3g9yveO7xSPiSGnfsx87k5qreJN3ra+7u9WWhShxxBBEH1O3IVh6dKlKC0thclkwuTJk/Htt99GPP7NN9/EhAkTkJiYiPz8fNxwww1obm7um8WqDOqRUk+TQ4RBJStJciW5QtNVW2y89XZGUlARlFzcpj5FbXwR319R04ajHv56VGKYnjOyxRCFMKQV80pvr4uncgZSt41v8yeEnBYWQ5LYkRTAh7cCax/lr2fcSZWuBHGCEVdhWLFiBe666y488MADqKiowJlnnony8nJUV1erHr9hwwYsWLAAN910E3bu3Il3330XmzZtwq9//eu+WXBgEz3G4PH64BBjDCEWg+xK8scYpGPsKpXPzbIwBFsMbeKxoW4kAJhUNAgAUNVkw04rtx6GmSyhBzIW0Fk1CmHQaPxWQ8NO/35LvXgdAcgb1/l1AjnrD0DqYH6+y8rbYpddG9s1CILodeIqDM888wxuuukm/PrXv8aYMWPw7LPPoqioCMuWLVM9/scff8SQIUNwxx13oLS0FDNnzsTNN9+MzZsj9OHpSQyJ4gsGeBywi6mqQHTzGCRhcHl9fO5zQOVzi1V0JSUHCYOUqhrGYkhL1GOiOHqzjvEMoVRXQ+iB9hae0glEF2MAAuYpbPHvqxOzirJGxt7OIDGDdyCd+ltgxl28oVy0lcUEQfQZcRMGl8uFLVu2YO7cuYr9c+fOxfffq3RvBDB9+nQcPXoUq1evBmMMx48fx3vvvYcLLrigL5bMu1xKuGxyfEGnEeTpbDIqWUkpRp3sNWnvcPutAOaFq4O7f0IthtCxnsH85syhAIDjojDAfCz0ICm+kJARfbOzwtP4NrAB3jGxPXb++OiuEcygEt4A77xHlL34CYI4YYhbgVtTUxO8Xi9yc5WfXnNzc1FfX696zvTp0/Hmm2/iiiuugMPhgMfjwUUXXYS///3vqscDgNPphNPplL83m8O0eIgGjZa7hjwdgMsKm5tbA0lGXWh7arnAzW8xaDQCUk16tHe4Ye5wIyc5mbd19nmQAjs6BJM8y1mmkxgDAFwwPh85qdNgqNUCa/7uLz4LxBxDqqqE1ACvdisPQGt1fKIZABRPC38eQRD9mrgHn4MfqIyFnwGwa9cu3HHHHXjwwQexZcsWfPbZZ6iqqsItt9wS9vpLlixBWlqa/FVU1M0+5pL7xGWHVYwVhASepT7sgCLGAPjTUdvsbh50FS2BVMGO9AQ9tMFppo7IriSJ04ZkYMIpYoM6S13oiMj2GnEBMdx/1kj+cz0d3J3kdvj7HUl1CQRBDDjiJgxZWVnQarUh1kFDQ0OIFSGxZMkSzJgxA/feey/Gjx+P888/H0uXLsUrr7yCujr12QGLFi1Ce3u7/FVTU9O9hcvCYINdTlUNE18AFFlJgL9JnlTQJj3wU2APdSMBYVtuq5KcB0DgmUT2oEwtWRgKO7+OhEbDx08CwN5PeRsKj4MPUMkcHv11CILoV8RNGAwGAyZPnow1a9Yo9q9ZswbTp6u3X7bb7dAE5fJrtfyhzMIMUTcajUhNTVV8dW/hojC4bXKqamJwZ9WA7quBdQyA32LwC4NkMdiQGZyqCnQafFagM/A5uECoO6lNFIZYJz+NFuM3O1byITkAMPYSSjEliAFMXF1JCxcuxL/+9S+88sor2L17N+6++25UV1fLrqFFixZhwYIF8vHz58/HBx98gGXLluHQoUP47rvvcMcdd2Dq1KkoKOijytcAi0GqRwitYRCFQaMPaW+dGiwMkiupU4shCmEA+Kd5gE8hC6T9KN/G4koCeGO6pGw+M0EaD1l2TWzXIAiiXxHX7qpXXHEFmpubsXjxYtTV1WHcuHFYvXo1Skr4PIC6ujpFTcP1118Pi8WCF154Affccw/S09MxZ84cPPnkk323aCkzyWWTYwxhXUkqIy/TA2MMgNyELl2wwpcSQRgiZCUpf0ARn9PbFlQL0pUYA8BTdM/7K7BKjONMuTG02ypBEAOKuLfdvvXWW3Hrrbeqvrd8+fKQfbfffjtuv/32Xl5VBKQ5zYoYQ7ArKTQjSSLElZTAhSFDsECXqtL4Lsrgs0zGML5tOejf5/X4U1i7MkR84lW8qV5HK5+FQBDEgCbuwtDvCHQlddYOQxdqMUjCYJaEQRyIkw4rklSFIYbgM8AHzwB8ZKaE+SgfbKM1RF/cFszgSV07jyCIfgcJQ6xI1c9uf7pqSPBZaqCnD33QS1lJbUHCkCFYMCgtkjCkR7e+DF7shpZD/n1NB/zvRRrpSRAEgROgjqHfIbuSrAEWQ7iW26GupKxkvq/RIloVojAMggV5wRaD1wO4pHnPMbqS2qoBryg+8hjOEdFdgyCIkxoShliRg892OSsptOV2aAM9Cf+cZn6M08Af+IMEK3KDLYaOVvGFEL3FkJIH6JO460jqitq8n2+jGcNJEMRJDwlDrKjEGMILQ6jFkJvK9zVZnfD6GJp9fI5zpmBBSvB1pCI1UxpvRxENggDkiEN2pDGa0sS1TLIYCILoHBKGWJGFwSoP3AkJPrvtymMDyEw2QiMAPgY025w46uRWxSDBEtoKpKOFb0V3U9TkT+TbY5WAz+cXiNxTYrsOQRAnJSQMsSJXPtthCWcxuERhCOzGKqLVCHKcocHsxCE7f22Cy3+ehGQxxCoMBRP5tq4SaD7Aq6d1CUDO2NiuQxDESQkJQ6wEuJKklNOQOc1um3hsqDAAQI7oTjpuduBAK+BiYvBashAk7JLFkBHbGgdP5tujm3l/I4CLRbTuKIIgTmpIGGIlUBgcXBhSTeEsBvVBNgVp3H1U02LHkdYOtILHGUIa33XVYsgZC6QWcpfWf+7l+4bNie0aBEGctJAwxIr4sGcufxO91BCLQYoxqFsMw3J4yuvBRhsONFjRyiRhCLYYJGGI0WIQBGDsRcp9Yy+O7RoEQZy0kDDEikESBqs88iAlxGIQXUkqMQYAGJbNhaGypg1VTTY0M7Gq2daoPFBKV02IURgAYMadQGIWfz3hat7SgiAIIgrI6RwrAYN6AMCo08CoCy5wCx98BoDhosWwvZZXNVv0WYAPfMBOIF11JQG8nuG2Tby4rXBq7OcTBHHSQhZDrIjCILhtAFioGwnwxxjCuJLG5qciyeAXE8MgsWW4JWikqa2Jb2N1JUkkZgDFZ/CBOwRBEFFCT4xYEdtfC8yHJDhCA8+APyspTPDZoNNg9ugc+fvCYrG/UbDFYD3Ot8kxzGkmCILoJuRKihV9AqDRAT4PUmDvksUAAH++YAxcHh9G56VgRIELqIDSYmDM/30KCQNBEH0HCUOsCAK3GjpakCJ0INWkIgydxBgAID8tAS8tmMK/qRaDzIEWg70F8IlN8LraKpsgCKILkCupK4izEcJbDFKBm7orKQTJIrDUQ051sorWQmImn+VMEATRR5AwdAVpTrNgDxNj6NxiUCDNafY4AEcbfy1ZDxRfIAiijyFh6AribITOYwxRWgw6oz8ltf0o31J8gSCIOEHC0BVEi0E1xsBYQFZSlBYDAAwq5Vtp8pokEKn53VgoQRBE7JAwdIWAGENIAz2PE2A+/jpCVlIIwSM5W6r4VhIMgiCIPoKEoStIriTBjoykoMCwO6B1dpg6BlUyxZGczQf5tlUUhgwSBoIg+pa4C8PSpUtRWloKk8mEyZMn49tvv414vNPpxAMPPICSkhIYjUYMGzYMr7zySh+tVsTotxhChMFl5VutMbY217LFIAqCNJaTLAaCIPqYuNYxrFixAnfddReWLl2KGTNm4B//+AfKy8uxa9cuFBcXq55z+eWX4/jx43j55ZcxfPhwNDQ0wOPx9O3CTf6spIykIFeS08K3xpTYrpkhWgyNe4CONn/VM1kMBEH0MXEVhmeeeQY33XQTfv3rXwMAnn32WXz++edYtmwZlixZEnL8Z599hnXr1uHQoUPIyOD9g4YMGdKXSwYAeAwp0EGyGILmOjvMfCuKR9TkjuUV1fYmYO9qvi+tGEgY1O31EgRBxELcXEkulwtbtmzB3LlzFfvnzp2L77//XvWcjz76CFOmTMHf/vY3DB48GCNHjsQf/vAHdHR09MWSZWzgQeVUoSM0+NxVi0GfAOSO4683vcy3+eO7sUqCIIiuETeLoampCV6vF7m5ynYPubm5qK+vVz3n0KFD2LBhA0wmE1auXImmpibceuutaGlpCRtncDqdcDqd8vdms7nbazezRKQBSNd0QKsRgn6geH1jjBYDABRN5XOaazf7vycIguhj4h58FgTlg5UxFrJPwufzQRAEvPnmm5g6dSrmzZuHZ555BsuXLw9rNSxZsgRpaWnyV1FRUbfX3Oo1AQDSBHvom90RhlMuVX4/al7s1yAIgugmcROGrKwsaLXaEOugoaEhxIqQyM/Px+DBg5GWlibvGzNmDBhjOHr0qOo5ixYtQnt7u/xVU1PT7bU3+/ignVRYQ9/saowBAIqnASPL+euJ1wJZI7q4QoIgiK4TN2EwGAyYPHky1qxZo9i/Zs0aTJ8+XfWcGTNm4NixY7Ba/Q/kffv2QaPRoLCwUPUco9GI1NRUxVd3Oe7hMYYE1gG4Hco3uxpjAHjn1qveBu7eCVz8QjdXSRAE0TXi6kpauHAh/vWvf+GVV17B7t27cffdd6O6uhq33HILAP5pf8GCBfLxV199NTIzM3HDDTdg165dWL9+Pe69917ceOONSEhI6LN119j18DDxV9fRonxTFoYuCpAgAGmFfEsQBBEH4pquesUVV6C5uRmLFy9GXV0dxo0bh9WrV6OkpAQAUFdXh+rqavn45ORkrFmzBrfffjumTJmCzMxMXH755Xj00Uf7dN3HLS60IgXZaOdzmVML/G/KMYYuWAwEQRAnAHEf1HPrrbfi1ltvVX1v+fLlIftGjx4d4n7qa46bHWhhKcgWRGEIRLIYuhJjIAiCOAGIe1ZSf+S42YFWiBZBsDA42vm2q64kgiCIOEPC0AXq27nFAICP4AxEGrRjSgNBEER/hIQhRjpcXpgdHrSyMBaDXZzfnJjRtwsjCILoIUgYYqSmlRe12XSiRRBsMUhCkUDCQBBE/4SEIUYON/HpbJokcRRnoMXgsgMesQJbGtVJEATRzyBhiJHDzVwY9CnZfIetwf+mVNOg0VG6KkEQ/RYShhipEi2GhMzBfIcloKWH5FZKzKQCNYIg+i0kDDGyo5YXsOUUiBPXzMf8b1J8gSCIAQAJQww43F7sqefCMHz4SL7TZfU3zusIsBgIgiD6KSQMUeL1MZz1t6/h9jLkpZowOCcTMIqZSZY6vrWK8YakrPgskiAIogcgYYiS/+yoQ4OFD/y5uKyAz4yQeiSZa8Wt6FYK7J1EEATRz4h7r6T+wgWn5sP+314cbrLh9jninITUfKBxN2AWLQbJckjJj88iCYIgegAShigRBAGXTwma/iZZBu3ikCCyGAiCGACQK6k7ZAzj2+YDfEvCQBDEAICEoTtIozeb9gE+rz/WkDo4fmsiCILoJiQM3SFTFIbmA0DrYcDrAnQmPoGNIAiin0LC0B0ySgFBw2sZqtaL+4YBGm1810UQBNENSBi6g84I5Izlr7e+xreSe4kgCKKfQsLQXYqn8e2xCr4tmhq/tRAEQfQAJAzdZehs5fdDzozLMgiCIHoKEobuMvK/gAyxoV7hVCDv1PiuhyAIoptQgVt30eqA61cD+z4DTrmU2m0TBNHvIWHoCVLzgSk3xHsVBEEQPULcXUlLly5FaWkpTCYTJk+ejG+//Taq87777jvodDpMnDixdxdIEARxkhFXYVixYgXuuusuPPDAA6ioqMCZZ56J8vJyVFdXRzyvvb0dCxYswDnnnNNHKyUIgjh5EBhjLF4//PTTT8ekSZOwbNkyed+YMWNwySWXYMmSJWHPu/LKKzFixAhotVqsWrUKlZWVUf9Ms9mMtLQ0tLe3IzU1tTvLJwiC6DfE8uyLm8XgcrmwZcsWzJ07V7F/7ty5+P7778Oe9+qrr+LgwYN46KGHovo5TqcTZrNZ8UUQBEGEJ27C0NTUBK/Xi9zcXMX+3Nxc1NfXq56zf/9+3HfffXjzzTeh00UXN1+yZAnS0tLkr6Kios5PIgiCOImJe/BZCErvZIyF7AMAr9eLq6++Go888ghGjhwZ9fUXLVqE9vZ2+aumpqbbayYIghjIxC1dNSsrC1qtNsQ6aGhoCLEiAMBisWDz5s2oqKjAbbfdBgDw+XxgjEGn0+GLL77AnDlzQs4zGo0wGo29cxMEQRADkLgJg8FgwOTJk7FmzRpceuml8v41a9bg4osvDjk+NTUV27dvV+xbunQp1q5di/feew+lpaVR/Vwp1k6xBoIgTiakZ140+UZxLXBbuHAhfvWrX2HKlCmYNm0a/vnPf6K6uhq33HILAO4Gqq2txeuvvw6NRoNx48Ypzs/JyYHJZArZHwmLxQIAFGsgCOKkxGKxIC0tLeIxcRWGK664As3NzVi8eDHq6uowbtw4rF69GiUlJQCAurq6TmsaYqWgoAA1NTVISUlRjWVEwmw2o6ioCDU1NQM61ZXuc2BB9zmw6Op9MsZgsVhQUND56OG41jH0N06WGgi6z4EF3efAoi/uM+5ZSQRBEMSJBQkDQRAEoYCEIQaMRiMeeuihAZ/+Svc5sKD7HFj0xX1SjIEgCIJQQBYDQRAEoYCEgSAIglBAwkAQBEEoIGGIkq5OmjtRWbJkCU477TSkpKQgJycHl1xyCfbu3as4hjGGhx9+GAUFBUhISMDs2bOxc+fOOK24Z1iyZAkEQcBdd90l7xso91lbW4trr70WmZmZSExMxMSJE7Flyxb5/YFwnx6PB3/+859RWlqKhIQEDB06FIsXL4bP55OP6Y/3uX79esyfPx8FBQUQBAGrVq1SvB/NPTmdTtx+++3IyspCUlISLrroIhw9erRrC2JEp7zzzjtMr9ezl156ie3atYvdeeedLCkpiR05ciTeS+sy559/Pnv11VfZjh07WGVlJbvgggtYcXExs1qt8jFPPPEES0lJYe+//z7bvn07u+KKK1h+fj4zm81xXHnX2bhxIxsyZAgbP348u/POO+X9A+E+W1paWElJCbv++uvZTz/9xKqqqtiXX37JDhw4IB8zEO7z0UcfZZmZmeyTTz5hVVVV7N1332XJycns2WeflY/pj/e5evVq9sADD7D333+fAWArV65UvB/NPd1yyy1s8ODBbM2aNWzr1q3s7LPPZhMmTGAejyfm9ZAwRMHUqVPZLbfcotg3evRodt9998VpRT1PQ0MDA8DWrVvHGGPM5/OxvLw89sQTT8jHOBwOlpaWxl588cV4LbPLWCwWNmLECLZmzRo2a9YsWRgGyn3+6U9/YjNnzgz7/kC5zwsuuIDdeOONin2/+MUv2LXXXssYGxj3GSwM0dxTW1sb0+v17J133pGPqa2tZRqNhn322Wcxr4FcSZ3Q1Ulz/Y329nYAQEZGBgCgqqoK9fX1ivs2Go2YNWtWv7zv3//+97jgggtw7rnnKvYPlPv86KOPMGXKFPzyl79ETk4OysrK8NJLL8nvD5T7nDlzJr766ivs27cPALBt2zZs2LAB8+bNAzBw7jOQaO5py5YtcLvdimMKCgowbty4Lt13XJvo9Qe6Mmmuv8EYw8KFCzFz5ky5U610b2r3feTIkT5fY3d45513sHXrVmzatCnkvYFyn4cOHcKyZcuwcOFC3H///di4cSPuuOMOGI1GLFiwYMDc55/+9Ce0t7dj9OjR0Gq18Hq9eOyxx3DVVVcBGDh/z0Ciuaf6+noYDAYMGjQo5JiuPKdIGKIk2klz/ZHbbrsNP//8MzZs2BDyXn+/75qaGtx555344osvYDKZwh7X3+/T5/NhypQpePzxxwEAZWVl2LlzJ5YtW4YFCxbIx/X3+1yxYgXeeOMNvPXWWzjllFNQWVmJu+66CwUFBbjuuuvk4/r7farRlXvq6n2TK6kTYp0019+4/fbb8dFHH+Hrr79GYWGhvD8vLw8A+v19b9myBQ0NDZg8eTJ0Oh10Oh3WrVuH559/HjqdTr6X/n6f+fn5GDt2rGLfmDFj5Lb1A+Xvee+99+K+++7DlVdeiVNPPRW/+tWvcPfdd2PJkiUABs59BhLNPeXl5cHlcqG1tTXsMbFAwtAJgZPmAlmzZg2mT58ep1V1H8YYbrvtNnzwwQdYu3ZtyAS80tJS5OXlKe7b5XJh3bp1/eq+zznnHGzfvh2VlZXy15QpU3DNNdegsrISQ4cOHRD3OWPGjJB043379smzTQbK39Nut0OjUT62tFqtnK46UO4zkGjuafLkydDr9Ypj6urqsGPHjq7dd8zh6pMQKV315ZdfZrt27WJ33XUXS0pKYocPH4730rrM7373O5aWlsa++eYbVldXJ3/Z7Xb5mCeeeIKlpaWxDz74gG3fvp1dddVVJ3zaXzQEZiUxNjDuc+PGjUyn07HHHnuM7d+/n7355pssMTGRvfHGG/IxA+E+r7vuOjZ48GA5XfWDDz5gWVlZ7I9//KN8TH+8T4vFwioqKlhFRQUDwJ555hlWUVEhp8RHc0+33HILKywsZF9++SXbunUrmzNnDqWr9jb/93//x0pKSpjBYGCTJk2S0zr7KwBUv1599VX5GJ/Pxx566CGWl5fHjEYjO+uss9j27dvjt+geIlgYBsp9fvzxx2zcuHHMaDSy0aNHs3/+85+K9wfCfZrNZnbnnXey4uJiZjKZ2NChQ9kDDzzAnE6nfEx/vM+vv/5a9f/H6667jjEW3T11dHSw2267jWVkZLCEhAR24YUXsurq6i6th7qrEgRBEAooxkAQBEEoIGEgCIIgFJAwEARBEApIGAiCIAgFJAwEQRCEAhIGgiAIQgEJA0EQBKGAhIEgCIJQQMJAEARBKCBhIIg+Yvbs2YpZ0wRxokLCQBAEQSigXkkE0Qdcf/31eO211xT7qqqqMGTIkPgsiCAiQMJAEH1Ae3s7ysvLMW7cOCxevBgAkJ2dDa1WG+eVEUQoNNqTIPqAtLQ0GAwGJCYmyhO5COJEhWIMBEEQhAISBoIgCEIBCQNB9BEGgwFerzfeyyCITiFhIIg+YsiQIfjpp59w+PBhNDU1yQPsCeJEg4SBIPqIP/zhD9BqtRg7diyys7NRXV0d7yURhCqUrkoQBEEoIIuBIAiCUEDCQBAEQSggYSAIgiAUkDAQBEEQCkgYCIIgCAUkDARBEIQCEgaCIAhCAQkDQRAEoYCEgSAIglBAwkAQBEEoIGEgCIIgFJAwEARBEAr+P4X9r6FEIk8WAAAAAElFTkSuQmCC",
      "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": 9,
   "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": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "def wrapper(par, control, cfg, verbose=False):\n",
    "    ode = DO(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": 11,
   "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": 12,
   "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": 13,
   "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": 14,
   "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": 15,
   "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": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "stat_vec = batch_run(parameters, control_list, cfg, n_workers=4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "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": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(torch.Size([2000, 2]), torch.Size([2000, 4]))"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "theta.shape, stat_vec_st.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " Neural network successfully converged after 315 epochs.train Done in 0 hours 0 minutes 39.261273 seconds\n"
     ]
    }
   ],
   "source": [
    "posterior = obj.train(theta, stat_vec_st, prior, method=\"SNPE\", density_estimator=\"maf\", num_threads=4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "with open(\"output/posterior.pkl\", \"wb\") as f:\n",
    "    pickle.dump(posterior, f)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "# with open(\"output/posterior.pkl\", \"rb\") as f:\n",
    "#     posterior = pickle.load(f)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "xo = wrapper(parameters, theta_true, cfg)\n",
    "xo_st = scaler.transform(xo.reshape(1, -1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "e911a18aac1a4f168b9de025b70011b6",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "  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": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAbgAAAHcCAYAAACkr7//AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAc6ElEQVR4nO3de4yd9Xkn8O+ZMxcb22MbG1g7AQMNSxNbmyZQLmlIQaRkyyYruqsqadg2aJsm+aNVqZVdkUvrqq3ESo02VRq1G4dVQsSmyVbQqlLSbFAaQhAkIQRSmtKUuGAb4tTYGHzB9sycc/YPe47nzMX2zLwz7znvfD6ShX0unoejYb48v2ut1Wq1AgAV01d2AQCwEAQcAJUk4ACoJAEHQCUJOAAqScABUEkCDoBKEnAAVFJ/mV/8yMiRrLxzZZLk8IcOZ8XgijLLgSXnF/p+uewS4Izub/7lnN6ngwOgkgQcAJUk4ACoJAEHQCUJOAAqScABUEldE3Bfe+pfyy4BgArpmoB7ZMe+sksAoEK6JuBGms2ySwCgQron4MZaZZcAQIV0UcDp4AAojoADoJK6JuCOCzgACtQ1ATcq4AAoUNcEnFWUABSpewJOBwdAgQQcAJXURQHXKLsEACqkiwJOBwdAcbon4BoCDoDidFHAOaoLgOJ0T8AZogSgQF0TcMfHmmm1dHEAFKNrAi5JRg1TAlCQrgo4C00AKEp3BZx5OAAKIuAAqKSuCrjjTjMBoCClBlyj2bmoRAcHQFFKDbjRSYtKXHoKQFG6KuCsogSgKCUHnCFKABZGqQE31jRECcDC6K4hSgEHQEEMUQJQSd3VwTXsgwOgGN0VcDo4AApS7iKTSUOUFpkAUBQdHACV1FUBp4MDoCjlBpyzKAFYIOUG3JijugBYGN21yGRUwAFQjFIDbnLHZh8cAEXpqg7OHBwARemqVZQCDoCilNzBWWQCwMIodw6uaZEJAAtDBwdAJXXVHJyTTAAoSsnbBKyiBGBhdNcQpYADoCBdFXDHx2z0BqAYJc/BTRqitMgEgIKUvE3AECUAC6Pk2wQsMgFgYZQacA33wQGwQLrqNgH74AAoSlddeDrWbKU5qasDgLnoqpNMEispAShGV83BJYYpAShGuXNw04SZhSYAFKHcIcoJ++AG+0+U4jQTAIrQNSeZDNVPlKKDA6AIXbPIZKBeS2KRCQDF6JrDlof6dXAAFKdrhigHBRwABeqaIcpTi0wEHADz1zUd3IAODoAClTsHN2Gj92BdBwdAcbpmiHKoXk9iFSUAxeiagDNECUCRSgu4Vqs17SpKJ5kAUITSAm5s0kHLtgkAUKTSAm7yVTmO6gKgSOUF3NikDm7g5FFdAg6AApQXcM3OIBvsO9nBWUUJQAG6ZojSSSYAFKl7higtMgGgQKUF3OShyIGTG711cAAUocRtApNWUQ7o4AAoTvcMUdYtMgGgOF00RHlykcmok0wAmL+uW0WpgwOgCOXNwTU6hyiHrKIEoEBd08EJOACK1D1zcP21aR8HgLnomg5ucHwf3KiAA2D+umYOzjYBAIrUNUOUjuoCoEjdM0TpsGUAClTiSSYzdXA2egMwfyWeRTn9HJwODoAidN8cXKOZVqs13VsA4Kx13WHLrdbU7g4AZqtrFpkM9dfbv7eSEoD5Ki/gmtMPUSYCDoD565ohynpfLfW+E8d1WWgCwHx1zRBlMuE0EwEHwDx1V8C1V1LaCwfA/JQYcFNXSg45zQSAgnRnByfgAJinrgw4HRwA89VdAWeRCQAF6co5OAEHwHx1VQc3fpqJS08BmK+uCjiLTAAoSom3CUwdojy1yMQ+OADmp7z74CwyAWABdeUQpW0CAMxXd66itMgEgHkqbw5umi7NIhMAilLeHFzTECUAC6erhih1cAAUpcQLT6fZ6G0VJQAFKXEf3DQBN3DyJBMBB8A8lTgHN80QZd0qSgCKUUrANZqtNKYLOCeZAFCQUgJuuk3eiUUmABSnuwKubpsAAMUoJeDGptkikCRDAzo4AIrRX8YXHe/g+mqdj1tkAl2kr37q903z4vSeUjq48QAbqHd++fYik1EBB8D8lDQHd2KIcnCGgNPBATBfJc3BnQiw+qQxyiGrKAEoSClzcKeGKGvJhKH9oX4nmUC3qA+vbP++8dLLJVYCc1PqEOVMc3CGKAGYr1L3wY0H2rj2PrhRK7YAmJ9Stwn0T5qD08FBF+kv5ccDFKbUIcr+GYYoRxutNKc5qxIAzlY5AXdyEcnApCHKlUOn/o/x4LHRRa0JgGopdw6u3jlEuWygnvNWDSVJdr34yqLXBUB1lBNwJ4cf+2tTv/ymc89JkuzcL+BgodUGBtu/Jmu+9HL7F/SiUoco+/unfvmL1p0IOB0cAPPRVUOUSbLp3BVJkp37jyxqTQBUS8nbBKYZolxniBIWS2t0ZObnxsYWsRIoXrknmRiiBGCBlDpEOTDNEOVFJxeZ/OTgsRxzogkAc1TqEOV0AbduxWBWDNZzZKSR5w68ktecv2qxy4Mlqf/Sizv+PPYvz5ZSBxSlpAtPT24TmHyld5JarZaL1o0vNDFMCcDclHof3OTbBMbZCwfAfJW8TWCGgLPQBIB5KmkObvr74MZd1N4qYC8cLJbWASeWUC0lzcGd6ODq0ywySSZs9tbBATBH3TkHd7KDe+7Fo2m4NgeAOSh1iHKmObgNq5elv6+WkUYzPzl4LK9as3wxy4MlqXHgQMef/9+Pn8jDPzmaDz6yL407r82a59aXVBnMTalDlNPtg0tOXIT66rUnQm2XlZRQmk/+w8v59t7j2X3VjrJLgVkr9TaB+jRnUY4b3wu360ULTaAM+442cu+/HE6S7H3dcxlZfrzkimB2ypmDa55+iDKxFw7K9rl/PpjxKfBWrZU9r99ZbkEwS+Ue1dU//RBlMuFWASspYcEdW3U0IyuOdTz2e/ePpbkySe3Er+d+dkfWPntex2sGjyzLskPmyOlOpQTcyNjpV1Empw5dNgcHC+8f/vO389KmfZ0PtnIi3HLin0fXHsl33v+1jpeseXZ9rrz7+kWoEGav1JNMprsPbtymdS4+hcWy8XuXpG+070SojZs8wDLxz62kb7QvGx+/ZBGqg7kppYMbn4MbOM1XH+/gDh4by0uvjGTNOYOLURosSRv/flOGf7w2f//Oh/PKuYdP/7++zeScF1fl333x2qzcN7xoNcJslbNN4CyGKJcP1nP+qqEkFprAYli5bzhXb39rLvjBq0/7ugt+cGGu/tSNwo2uV+6Fp6cZokwsNIHFVh/tz9pd53UOVU7UStbuXJ/6WCmDPzArJQXc+GHLM6+iTJKLT87D/Z9v7czREbd7w2I4uOFAas2T/22OB93Jf9aatRzc+FIZZcGslXoWZf9phiiT5D1vujgrBuv59jMv5tfvflTIwSJ4+dUvplVvpdaopW+sLxc9cln6xvpSa9TSqrfy8qv3l10inJVSb/Q+3Rxckmx51ep87tevyorBeh7esV/IwQJr1Bs5sv5QkmT5gZW5avtb82+/+vpctf2tWX7gxIjKkfWH0qj775DuV+4c3BkCLkmu2HRuPvfrV2XlUH8e3rE///GTD+V/fWNHdpuXg8I1BxpZuXc4Gx6/uGMhyfgClA2Pb8rKfx1Os1/A0f1qrVZr0e+jed3vfSWvjDTyt7f/bF73qQuSJIc/dDgrBlfM+J7Hdh7IbZ/5Tg4dG2s/9toNw/mp81Zk45rl2bh6WS5evyI//W+Gc8HwUGq108/vAckv9P3ylMdaaaU2ZRPc2T8PRbu/+Zdzet+8l0Jt/eITs37PsdET//d3Nh3cuCs2rc03/tsN+dt/2JMvP7knj+zYn6f2HMxTew5Oee3q5QO5+OQKzLFmK41mK4P9fVk2UM85g/Us66+n3ldLX18t9VrSJwwr73++82fKLqFnnCm8hBu9Yt4Bd9/jz8/pfQP1WlYOze7Ln7tiMLdevSm3Xr0p+w8fz6PPvpjnXzqWH790NM8fOJofvXA4z+w7kpePjub7z708p7qoJgEHS8+8A+7DN//0nN635VWrs2rZwJy/7rqVQ/n3WzZMefzYaCM7XjicH790LPW+E1fy1Gu1jDQaOTrSzCsjYzk21kyz2cpYs5Vms5XWjJt+AOhV8w64973lp+b83iMjxZ8zuWygns0bV2fzxtWF/90A9I5SVlECwEITcABUkoADoJIEHACVJOAAqKRSTjIBgIWmgwOgkgQcAJUk4ACoJAEHQCUJOAAqaV5nUbZarRw6dKioWmBBrVq1yj2BsITMK+D27duX888/v6haYEHt3bs35513XtllAItkXgE3ODiYJNm9e3eGh4cLKQiKdvDgwVx44YXt71dgaZhXwI0P9wwPDws4up7hSVhaLDIBoJIEHACVNK+AGxoayrZt2zI0NFRUPVA436ewNDlsGYBKMkQJQCUJOAAqScABUEkCDoBKEnAAVNKcAu7RRx/NzTffnLVr12bFihW56qqr8vnPf77o2mDO7rnnnrz//e/PlVdemaGhodRqtXz2s58tuyxgEc36qK4HHnggb3vb2zI4OJh3vetdWb16de67777ceuutefbZZ/PhD394IeqEWfnoRz+anTt3Zv369dmwYUN27txZdknAIptVBzc2Npb3vve9qdVqefDBB/PpT386H/vYx/L9738/mzdvzrZt2/L0008vVK1w1u666648++yzeeGFF/KBD3yg7HKAEswq4P7u7/4uO3bsyLvf/e684Q1vaD++atWq/O7v/m7Gxsbymc98pvAiYbbe+ta3ZtOmTWWXAZRoVgH3wAMPJEluuummKc+NP/aNb3xj/lUBwDzNKuDGhx8vu+yyKc+tXbs269evN0QJQFeYVcC9/PLLSZLVq1dP+/zw8HD7NQBQJvvgAKikWQXceOc2U5d28ODBGbs7AFhMswq48bm36ebZDhw4kH379k07PwcAi21WAffzP//zSZKvfvWrU54bf2z8NQBQplkF3I033phLL700n//85/PEE0+0Hz906FD+8A//MP39/bntttsKLhEAZm/WN3p//etfz9ve9rYMDQ3lV37lVzI8PJz77rsvzzzzTP7oj/4oH/nIRxaqVjhrd911Vx566KEkyZNPPpnvfe97+bmf+7m85jWvSZLccsstueWWW0qsEFhosw64JPnOd76Tbdu25ZFHHsnIyEg2b96c22+/PbfeeutC1Aizdtttt+Xuu++e8flt27bl93//9xevIGDRzSngAKDb2QcHQCUJOAAqScABUEkCDoBKEnAAVJKAA6CSBBwAlSTgAKgkAQdAJQk4ACpJwAFQSQKui3zlK1/Jm9/85qxZsybr1q3L29/+9uzYsaPssgB6koDrIkeOHMnWrVvz6KOP5mtf+1r6+vryS7/0S2k2m2WXBtBz3CbQxV544YWcf/75efLJJ7Nly5ayywHoKTq4LrJjx468+93vzqWXXprh4eFccsklSZJdu3aVXBlA7+kvuwBOecc73pELL7wwn/70p7Nx48Y0m81s2bIlIyMjZZcG0HMEXJfYv39/nnrqqXzqU5/KddddlyR56KGHSq4KoHcJuC6xdu3arFu3Ltu3b8+GDRuya9eu3HHHHWWXBdCzzMF1ib6+vnzhC1/IY489li1btuR3fud38sd//MdllwXQs6yiBKCSdHAAVJKAA6CSBBwAlSTgAKgkAQdAJQk4ACpJwAFQSQIOgEoScABUkoADoJIEHACVJOAAqCQBB0AlCTgAKknAAVBJAg6AShJwAFSSgAOgkgQcAJXUX3YBQHdqtVp5ZfSVJMk5A+ekVquVXBHMjg4OmNYro69k5Z0rs/LOle2gg14i4ACoJAEHQCUJOAAqScABUEkCDoBKsk0AmOLI8bF865l9ZZcB8yLggCk+fv8/Z/tDTyXLy64E5s4QJTDFj18+WnYJMG8CDphiZKxVdgkwbwIOmGKk0Sy7BJg3AQdMMTLWKLsEmDcBB0wx2jBESe8TcMAUI2OGKOl9Ag6YQsBRBQIOmMIiE6pAwAFT6OCoAgEHTKGDowoEHDCFDo4qEHDAFKM6OCpAwAFT6OCoAgEHdGg2Wxlr2uhN7xNwQAcLTKgKAQd0EHBUhYADOph/oyoEHNBBwFEVAg7oYIsAVSHggA46OKpCwAEdjgs4KkLAAR2soqQqBBzQYVQHR0UIOKCDDo6qEHBAB4tMqAoBB3SwTYCqEHBAB6soqQoBB3QwRElVCDigg0UmVIWAAzrYJkBVCDiggw6OqhBwQAdzcFSFgAM6jDRaZZcAhRBwQAcdHFUh4IAOAo6qEHBAh5FGo+wSoBACDugwOmYOjmoQcEAH2wSoCgEHdDAHR1UIOKCDDo6qEHBABx0cVSHggA7jAbd8oF5yJTA/Ag7oMD5EuWJIwNHbBBzQYfxG75VD/SVXAvMj4IAO7SHKQR0cvU3AAR3GA27loA6O3ibggA6n5uAEHL1NwAEdxju4cwQcPU7AAR2soqQqBBzQ4dQcnICjtwk4oMP4NgFDlPQ6AQe0NZutjDZOXJezwipKepyAA9omHrRsFSW9TsABbaMCjgoRcEDbxJsEzhn044He5jsYaBsfouzvq2Wo3ypKepuAA9rGO7jB/r4M1v14oLf5DgbaxufgBvv7MlCvlVwNzI+AA9qOn+zgBup96dfB0eN8BwNt7SHKel+G+v14oLf5Dgbaxjd5D/X3ZUAHR4/zHQy0TVxkMqCDo8f5DgbaRhqNJCfm4CZ2cK1Wq6ySYM4EHNDWsU2g79SPh0ZTwNF7BBzQNnJyDm6w3peB/tqEx5szvQW6loAD2sY7uIFJG71Hx3Rw9B4BB7RN3CZQ75vQwTV1cPQeAQe0jZ9kMtTfl1qtNuVx6CUCDmibuMhkotExAUfvEXBA2/hiksnnUOrg6EUCDmg7PkMHd1wHRw8ScEBb+zaBeuddcGM6OHqQgAPaTm0TmDxEaZsAvUfAAW3jATc06aBlG73pRQIOaJt44el0j0MvEXBA24zbBAxR0oMEHNB2vHHqRu+JRqyipAcJOKBt5g5OwNF7BBzQdmqbgICj9wk4oM0cHFUi4IC2ibcJdDyug6MHCTigzTYBqkTAAW0znUXpNgF6kYAD2kZm2iagg6MHCTigzTYBqkTAAW0zbhMYs4qS3iPggLaZOriRRqOMcmBeBBzQNtM2Afvg6EUCDmgbD7LJHdyYgKMHCTggSdJqtdqrJScH3HGLTOhBAg5I0rkVYPI2Afvg6EUCDkjSeSXOkG0CVICAA5J0LiSZ0sEJOHqQgAOSnOrg6n211PtqHc9ZRUkvEnBAkpm3CCSO6qI3CTggyanN3JNXUCad83PQKwQckCQZGZt+D1ySjOng6EECDkhyahhyuiFKi0zoRQIOSDLzOZQTn4NeIuCAJDPfJJAkI1ZR0oMEHJDk9B3cWFPA0XsEHJAkOT42fpt3bcpzhijpRQIOSJIZD1pOLDKhNwk4IMmpA5UH++tTnhtrttI0TEmPEXBAkonbBKYOUU58HnqFgAOSnH6RSWKYkt4j4IAkp98mcOJ5Q5T0FgEHJDm1inKmDs5KSnqNgAOSnAqwyXfBjTNESa8RcECS028TmPg89AoBBySZuE1AB0c1CDggyelvE0jMwdF7BByQ5PQ3eic6OHqPgAOSnMUc3JhtAvQWAQckOfNGb4tM6DUCDkhyFtsEzMHRYwQckOTMQ5Tm4Og1Ag5IcirAhgxRUhECDkhy5iFK2wToNQIOSHI22wSsoqS3CDggSTJyMsBm3ibQWMxyYN4EHJDkVIDNvMhEB0dvEXBAklOLSGacg7PIhB4j4IAkp+bgZlxFaZEJPUbAAUlODUHaB0dVCDggiQtPqR4BByQ5i7MoDVHSYwQckFardeb74KyipMcIOKBjC4A5OKpCwAEdWwDc6E1VCDigI7x0cFSFgANy8OhokqTeV0u9rzbtawQcvUbAAfm/392dJPmZC9fM+JrjhijpMQIOlrjDx8dyz7d2Jkl+47pLZ3ydDo5eI+Bgifvio7tz8NhYLl2/Ir/wugtmfJ3Dluk1Ag6WsNFGM//7m/+SJHnvdZfOOP+WWEVJ7xFwsIR96e/35McvH8v6lYP5T2981Wlfa4iSXiPgYAn71IMnurfb3nRxlg3UT/ta1+XQawQcLGFP7TmYcwbr+S/XbDrjaw1R0mv6yy4AKNe7fvairDln8Iyve+7A0Vx8x5cWoSLo9Oz/+A9zep8ODpaw37jukvzXN198xtedv+rMAQjdptZqtaz9BaY4MnIkK+9cmSQ58N8PptEQcpRj3cqhOb3PECVwRgP1vqxZPrcfMlAWQ5QAVJKAA6CSBBwAlSTgAKgkAQdAJdkmAEyr1WrlldFXkiTnDJyTWm3mg5ihGwk4ACrJECUAlSTgAKgkAQdAJQk4ACpJwAFQSQIOgEoScABUkoADoJIEHACVJOAAqCQBB0Al9ZddAFCOVquVQ4cOlV0GnJVVq1bN+sBvAQdL1L59+3L++eeXXQaclb179+a8886b1XsEHCxRg4ODSZLdu3dneHi45Gqq4eDBg7nwwgt9pgUa/0zHv19nQ8DBEjU+3DM8POyHccF8psWby32EFpkAUEkCDoBKEnCwRA0NDWXbtm0ZGhoqu5TK8JkWbz6faa3VarUWoCYAKJUODoBKEnAAVJKAA6CSBBwAlSTgYIn6sz/7s1xyySVZtmxZrrjiinzzm98su6SeMJvP7YEHHkitVpvy65/+6Z8WseLe9OCDD+Yd73hHNm7cmFqtlr/+67+e9d8h4GAJ+uIXv5jbb789H/nIR/L444/nuuuuyy/+4i9m165dZZfW1eb6uf3whz/Mnj172r8uu+yyRaq4dx05ciSvf/3r88lPfnLOf4dtArAEXX311XnjG9+YP//zP28/9trXvja33HJL7rzzzhIr626z/dweeOCB3HDDDTlw4EDWrFmziJVWS61Wy1/91V/llltumdX7dHCwxIyMjOSxxx7LTTfd1PH4TTfdlIcffrikqrrffD63N7zhDdmwYUNuvPHGfP3rX1/IMplAwMESs2/fvjQajVxwwQUdj19wwQX5yU9+UlJV3W8un9uGDRuyffv23Hvvvbnvvvty+eWX58Ybb8yDDz64GCUveW4TgCVq8unsrVZrTie2LzWz+dwuv/zyXH755e0/X3vttdm9e3c+9rGP5S1vecuC1okODpac9evXp16vT+k69u7dO6U74ZSiPrdrrrkmTz/9dNHlMQ0BB0vM4OBgrrjiitx///0dj99///1505veVFJV3a+oz+3xxx/Phg0bii6PaRiihCVo69at+dVf/dVceeWVufbaa7N9+/bs2rUrH/jAB8ouraud6XP70Ic+lOeffz6f+9znkiR/8id/kosvvjibN2/OyMhI7rnnntx777259957y/zX6AmHDx/Oj370o/afn3nmmTzxxBM599xzc9FFF53V3yHgYAl65zvfmf379+cP/uAPsmfPnmzZsiVf/vKXs2nTprJL62pn+tz27NnTsSduZGQkH/zgB/P8889n+fLl2bx5c770pS/l5ptvLutfoWd897vfzQ033ND+89atW5Mk73nPe/LZz372rP4O++AAqCRzcABUkoADoJIEHACVJOAAqCQBB0AlCTgAKknAAVBJAg6AShJwAPN0/fXX5/bbby+7DCYRcABUkoADoJIEHEABxsbG8pu/+ZtZs2ZN1q1bl49+9KNx1G+5BBxAAe6+++709/fn29/+dj7xiU/k4x//eO66666yy1rS3CYAME/XX3999u7dmx/84Aep1WpJkjvuuCN/8zd/k3/8x38subqlSwcHUIBrrrmmHW5Jcu211+bpp59Oo9EosaqlTcABUEkCDqAA3/rWt6b8+bLLLku9Xi+pIgQcQAF2796drVu35oc//GH+4i/+In/6p3+a3/7t3y67rCWtv+wCAKrg137t13L06NFcddVVqdfr+a3f+q28733vK7usJc0qSgAqyRAlAJUk4ACoJAEHQCUJOAAqScABUEkCDoBKEnAAVJKAA6CSBBwAlSTgAKgkAQdAJf1/eEUfK1bommMAAAAASUVORK5CYII=",
      "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)"
   ]
  }
 ],
 "metadata": {
  "language_info": {
   "name": "python"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}