{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# [Damped Oscillator - CDE](https://github.com/Ziaeemehr/vbi_paper/blob/main/docs/examples/do_cpp_cde.ipynb)\n",
"\n",
"- Inference without torch dependency:\n",
" - DMN\n",
" - MAF\n",
"\n",
"![\"Open](\"https://colab.research.google.com/assets/colab-badge.svg\")
"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"import pickle\n",
"import corner\n",
"import numpy as np\n",
"from vbi import BoxUniform\n",
"import autograd.numpy as anp\n",
"import matplotlib.pyplot as plt\n",
"from multiprocessing import Pool\n",
"from vbi.cde import MDNEstimator, MAFEstimator\n",
"from sklearn.preprocessing import StandardScaler\n",
"from vbi.models.numba.damp_oscillator import DO\n",
"from vbi.utils import posterior_shrinkage_numpy, posterior_zscore_numpy\n",
"# switch to C++ implementation\n",
"# from vbi.models.cpp.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)"
]
},
{
"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_cut\": 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": [
"Damped Oscillator Model (Numba)\n",
"-------------------------------\n",
"a = 0.1\n",
"b = 0.05\n",
"dt = 0.01\n",
"t_start = 0.0\n",
"t_end = 100.0\n",
"t_cut = 20.0\n",
"method = rk4\n",
"output = output\n",
"initial_state = [0.5 1. ]\n",
"\n"
]
}
],
"source": [
"ode = DO(parameters)\n",
"print(ode())"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYYAAAEiCAYAAAD9DXUdAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAA9hAAAPYQGoP6dpAABQYklEQVR4nO2dd3hUVfrHv9Myk0mZ9EY6NRBqAAGRpoKAuioq6iq66m+xi7iui25Rdl10dV3bihVxrayroigWbCAghJJA6C2k9zIlZTKZub8/ztwpyZR770wyQN7P88xzhjs3JyeX5H7vW877yjiO40AQBEEQduShXgBBEARxZkHCQBAEQbhBwkAQBEG4QcJAEARBuEHCQBAEQbhBwkAQBEG4QcJAEARBuEHCQBAEQbihDPUC+hubzYbq6mpERUVBJpOFejkEQRD9AsdxMBqNSEtLg1zu2yYYcMJQXV2NjIyMUC+DIAgiJFRUVCA9Pd3nOQNOGKKiogCwixMdHR3i1RAEQfQPBoMBGRkZjnugLwacMPDuo+joaBIGgiAGHEJc6BR8JgiCINwgYSAIgiDcIGEgCIIg3BhwMQaCIAYOVqsVFosl1MvoF1QqFRQKRVDmImEgCOKcg+M41NbWorW1NdRL6VdiYmKQkpIS8B4tEgaCIM45eFFISkqCVqs95zezchyH9vZ21NfXAwBSU1MDmo+EoQ/5355KbDnWgIfnj8CgmPBQL4cgBgRWq9UhCvHx8aFeTr8RHs7uMfX19UhKSgrIrUTC0EeUNrbhof/tA8cBnRYrXlsyMdRLIogBAR9T0Gq1IV5J/8P/zBaLJSBhoKykPmLToVpwHHv/07EGdFqsoV0QQQwwznX3kSeC9TOTMPQRe8paHO+7um3YV9EausUQBEGIgIShjyip1AMAojTMW3eszhjK5RAEQQiGhKEP6LRYUWPoBADMz08BABwlYSAIQiDPP/88cnJyoNVqccUVV0Cv1/fr9ydh6AMqW9rBcUCkWonJOSwrorSxLcSrIgjibOCRRx7BSy+9hLfffhtbt25FUVERHn/88X5dAwlDH1DW1A4AyIzTOtJUq1o6QrkkgiDOAnbt2oWnnnoK69atw4wZMzBhwgQsXboUX3zxRb+ug4ShDzhtF4aseC3SY5kwVOs7YbNxoVwWQRBnOM888wzmzJmDCRMmOI4lJiaisbGxX9dB+xj6gJpWZh2kx4YjRaeBXMYykxrbzEiK0oR4dQQx8OA4Dh0hSBkPVykEp5CazWZs2LABzzzzjNvxjo4O6HS6vlieV0gY+oB6oxkAkBytgUohR3K0BjX6TlS1dJAwEEQI6LBYMfLP3/T79z20ch60YcJus3v37kVHRwcefPBB/P73v3cct1gsmD17dl8t0SMkDH1Ag10YEqPUAIAUHROGOoM5lMsiCOIM5tixY9BoNCgpKXE7fvnll+P888/v17WQMPQBDSa7MEQyYUiwj40mEgaCCAXhKgUOrZwXku8rFIPBgKSkJAwZMsRxrLy8HEeOHMGiRYv6YnleIWHoA+rtexiSopkg8JYDb0kQBNG/yGQywS6dUJGQkACDwQCO4xxxiSeeeAILFizAyJEj+3UtZ/aVOgvptFhh6OwGACRGsngCWQwEQfhjzpw56OzsxJNPPonrr78e77//Pj7//HMUFhb2+1ooXTXI8FZBmFKO6HCmu4mRYQBIGAiC8E5ycjLWrl2L1atXY+TIkdi+fTu2bt2KjIyMfl8LWQxBxjW+wJuDvMVAriSCIHyxePFiLF68ONTLIIsh2DT2yEhyfd9o6grJmgiCIMRAwhBkWttZk5C4iDDHMYoxEARxNkHCEGRa2plVEBOuchyLt8cY2rusaO/qDsm6CIIghELCEGRaO5jFoNM6hSFSrYRSzuINvEVBEARxpkLCEGT4G3+s1ulKkslkiLH/m7coCIIgzlRIGIJMK+9KcrEYACDW/m89WQwEQZzhkDAEGd5i0IW7CwMvFC0kDARBnOGQMAQZPsbg6koC4HAltXaQK4kgiDMbEoYg482VxGcpUfCZIIgzHRKGIMPf+GPC3S2GWPu+htYgBZ9b27vw/HfH8fWB2qDMRxAEwUMlMYJIp8Xq6BKl62Ex8DGHYMUY7v2gCD8fZ+3+Pvi/KZg6OD4o8xIEQZDFEET09viCQi5DtMZdc/mYQzBcSScbTA5RAIA120oDnpMgiDODe+65B9OnT/f4WXZ2Np544ok+XwNZDEHENSOpZ59XPuYQDFfSNweZ+ygpSo16oxlbjzfC3G2FWim8KQhBEGcehw4dwurVq7FlyxaPn+fl5aG4uLjP10HCEERaPZTD4HEIQ0fgFsPeshYAwG9n5OKVzafQaDKjqLwVU3LJnUQQHuE4wNLe/99XpQV6PCT64umnn8akSZO8tvKMi4tDRUVFsFbnFRKGIMLf9HtmJAGurqTALAaO41BU3goAGJ8Zi0nZsfjqQC1KKvUkDAThDUs78Pe0/v++j1QDYRGCTu3u7sbHH3+MP/3pT45jS5cuxeTJk3HbbbcBAIxGIyIihM0XCBRjCCJ6L5vbAFdXkgUcx0n+HlWtHWhq64JKIcOotGjkD9IBAA5U6yXPSRBE6Dl58iSMRiNGjx4NALDZbPjoo48QGRnpOGf//v3Iy8vr87WQxRBEjGZWOTVK40EY7Omr3TYOJnO3x3OEcLKhDQCQHR8BjUqBkWnRAIADVSQMBOEVlZY9vYfi+wqktbUVABxC8M0336ClpQVhYezeUVhYiLKyMlxxxRXBXmUvSBiCiLGTWQyRmt6XVaOSI0whR5fVBmOndGE41WACAOQmMnNyREoUAOB0UzssVhtUCjICCaIXMplgl06oyMrKgkwmwwcffICIiAg8+OCDWLBgAT777DNkZ2dj6dKlmDNnDmbMmNHnawnpXWTLli247LLLkJaWBplMhvXr1/s8v6amBjfccAOGDx8OuVyOZcuW9cs6hWLs5C2G3sIgk8kcPaD1AQSgT9kthtxE9lSRHKVBuEoBq41DZUuH5HkJgggtKSkpeOKJJ/Duu+9i/vz5WL58OVatWoXNmzdj+vTpGD58OD766KN+WUtIhaGtrQ1jx47FSy+9JOh8s9mMxMREPProoxg7dmwfr048JrswRHuxBvjjhkCEodFuMSSwpx+5XIZs+/tS+2cEQZydrFixAnq9HrW1tbj99tsxZswYlJaWoq2tDR9++CHi4uL6ZR0hdSXNnz8f8+fPF3x+dnY2nn/+eQDAmjVr+mpZkjGa7a4ktefLGmUPShs6pXdx62kxAEwkDtcYcKqhDXNGSJ6aIAgCwACIMZjNZpjNzl7LBoOhz76XL1cS4MxWkmoxmLutqNF3AgCy4p1BrRyHxdAmaV6CIAhXzvlI5apVq6DT6RyvjIyMPvteTmHw5kpigmHolCYMdXomcGqlHPERziJ9vCupvDkEG3gIgjjnOOeFgffZ8a++3DXoyEry4kqKtlsMUoPP1XoWXE7VadxKbgyKCQcAVFHwmSCIIHDOu5LUajXUanW/fC9/riRn8FlajKHGIQzhbscdwtDaAY7jetVpIgiCEMM5bzH0J0Z/WUnhgbmSqltZfCEtxl0YUnQayGSAuduGprbAi/S1mbtx/4dFuOS5Lfj6QE3A8xFEKAikwsDZSrB+5pAKg8lkQnFxsaNaYGlpKYqLi1FeXg6AuYGWLFni9jX8+SaTCQ0NDSguLsahQ4f6e+m96LbaHL0YPG1wAwIPPle3MoshLUbjdjxMKUdSFLOKguFOWvXVYXxWXI0jtUbc+0ERBbWJswqViv2dtbcPvJgb/zPz10AqIXUl7d69G7Nnz3b8e/ny5QCAm2++GWvXrkVNTY1DJHjGjx/veL9nzx68//77yMrKwunTp/tlzd4wmZ3uIb+uJIkWA5+R1NOVBDB3Up3BjOrWDozNiJE0PwA0msz4oJDFYaLUShjN3Xh180k8uWiM5DkJoj9RKBSIiYlBfX09AECr1Z7z7lWO49De3o76+nrExMRAoQisBH9IhWHWrFk+TZ+1a9f2Onammoe8G0mjknstS+EMPkuLMfAWQ2oPiwFg7qW95a2oag3MYthYUgOrjcPYdB1WLMjDda/twJclNXjs8lHQqKjfA3F2kJKSAgAOcRgoxMTEOH72QDjng8/9BS8MkWrvJpwjXTVAV9KgGA8WQ6wzAB0IfBOgy8amYXJ2HFJ1GtToO7H9ZCPmjEgOaG6C6C9kMhlSU1ORlJQEiyU47XTPdFQqVcCWAg8JQ5DgU1V7tvR0JTpcuiupvavbsWM6RdfbYghGympXtw177E2AZg1PhFwuw6zhifigsALbTjSRMBBnHQqFImg3y4EEZSUFCX+pqoAz+Gwyd8NmE+cSazCyzW0alRxRHvZJpNnjDnwcQgolVXp0WmyIiwjDYHvJjWmDEwAA2040+vpSgiDOIUgYggQffPaWkQQ4RYPjnL0bhNJoYsKQGKX2GEjjrYhag3RhKCxtBgBMyo51fI9pg1lXuCO1RkcjIoIgzm1IGIIE70qK8hFjUCsV0KjYJRcbZ+AthoRIz5v1kqOZMDSazLBYbaLm5impagUATMiMdRyLj1QjM47VZdpv/5wgiHMbEoYgYRDgSgKcKatiy2I0mNjGNW/CEB8RBqVcBo5ziohYDtcYAcDRFY5nTDprH7q/Mjhd4jq6rPjmYC11nSOIMxQShiDhyEryJwwSA9CNRqcryRNyucyxya1Ogjupvasbp5vYRra8VHdhGGffF7GvolX0vD1pMJqx8IWfsfSdPbj0xa144fvjAc9JEERwIWEIEiZ7LwZ/LTudu5/FxRgaTL5dSQCQbI8zSBGGI7VGcBwTnp7fY0x6DIDgWAx/+fwATjW2IUzJfvWe3XQMO081BTwvQRDBg4QhSDjrJPlzJUmrl+TPYgCAFHucoVZCZtLhGtanoqe1AACj0qIhk7HAdpNJmpsKAEoq9dhYUgu5DPjs7vNx/WRWAv3Jr49InpMgiOBDwhAkhKSrAi6uJJExBkdWUmSY13P4AHStQfzN2ykMUb0+i1ArHQHoo7VG0XPzvF/IyptcOiYNeanReODiYQhTyFFU3oq95S2S5yUIIriQMAQJk4Cdz4D0vs9CXEkpAbiSTtaz+MLQpN7CAADDk9nxIxKFoaPLis+LqwAA10/OBAAkRWlw2dg0AMBHu/uuTwZBEOIgYQgSvGvIv8XAu5JE7mMwsqykvnIl8RVUcxMjPH4+IoUXBmmtUbedaERblxWDYsIxJdfZ0HzRhEEAgI0ltejqlpZm60pHlxXv7CjDU18fcezLIAhCHFQSI0gIdSVJKb3dZu52lPT2GXyOlmYxtJm7HRvjchO8CIM99iDVlfT9EVbM7MK8JLcNeuflxiMpSo16oxlbjjXgopHSy240msz49es7cbSOrXH1Tydxz+wheHDusHO+uiZBBBOyGIIEv/NZ6D4GMcFnfl9CuEqBCC9tQwH33c9iqtDyaaqxWhVitJ5jGMPtFsOxOpPoch4cx+FHuzDMGZHk9plCLsOC0akAgE2H6kTN2/N73PdBEY7WGZEQqcZC+5wv/XgCH+2plDwvQQxESBiCAMdxLsLgJ8YgIV3VtRyGL5Kj2eftXVZRJTdON7LmHtlerAUAyI6PgFopR4fFivJmcQ1QjtQaUWvoRLhKgSm58b0+58Xix6P1ksuqf7SnEttPNkGjkuPD307Bv389AQ9cNAwAsHLDIUlxF2+cqaXfCSJYkDAEgfYuK6z2p+i+2PnsLIfhPSMJALRhSsf3rxdxIyxtNAEAcnwIg0Iuw9BkVlhPbJxh12nm65+YHeuxp8N5uXHQhilQbzTjYLX4GIbFanNslFt20TAMSWLrvGfOEIzPjIHJ3I1/bTomet6eFJW34Ndv7MDwP32NCX/dhL98dgCt7YG3UiWIMw0ShiDAxxcUchnC/TSzkdL3WajFALgGoIWnrJ7iA88+hAEARqSwOIPYzCRncb44j5+rlQpMH8KquH5/WHxjlfVFVahs6UBCZBhunprtOK6Qy/DHhXkAgP/ursCJepPouXk+K67Cta/+gm0nmtDVbUNzWxfe/qUMC1/YilMN0ucliDMREoYgwO96jlQr/QY5pQSf/dVJckVKldXTdmHISYj0eR6fmSQmAM1xnJvF4A3enfTDUfHC8J9fygAAt03PRXiYuzAXZMXhorxk2Djg1c0nRc8NANtPNOLB/+6DxcphwegUbHpgBt6+dTKy4rWoau3A9a/vCLhBUmt7F17+6QSuf20HFjz/M5a+sxvri6okF0QkiEAgYQgCQgvoAU5XUluXVfAfvb/Kqq5IyUziU1WzE7Q+z5NiMVS2dKDOYIZSLsP4DO/CMNsuDPsrWx0WkhAOVOlRUqWHSiHD4kkZHs+5a/ZgAMD64irU6MXdwPUdFty/rhjdNg6Xj03DS9dPwNDkKMwcloiP75yGYcmRqDOYsfSd3ei0Z46J5duDtZj1zE/4x9dH8cupJhyqMeCbg3VYtq4Y85//GbtPB552W6PvwOf7qvHK5pN4fcspfHeozlERmCB6QumqQcDZ1lOAMIQ7g9OGDgviBdzspbmShAlDS1sXWux9FrLj/biS7LuiTze1ob2rG9ow/z8vby3kD9L1epp3JTlag1Fp0ThYbcBPRxtwdUG6oPX/174xbu6oFMRFeI7BTMiMxXk5cdhZ2ow3fy7FHy8dKWhuAHjmm6NoMJqRmxCBf1w9BnK50yJMiFRjzS2TcPlL23CgyoBHPi3BP68ZKyo19rPiKiz/7z5YbRyGJUfi5mnZSNOFo7iiFe/sKMOJehMWv7YDv583HL+dkSs67bawtBkv/XgCW4419PpMpZDhorxk3H5BLgqyvIu2L8zdVmw70YiDVQaUNbejq9uGCLUS6bHhGJ8Rg/GZsT7/34XS3NYFU2c3bBwHXbgKMVoVpSD3ISQMQcDkqJPkOyMJYH7vKI0Sxs5utIoUBkEWg0hXUqk9VTUlWuMzFZb//gmRYWg0deF4nQlj7VVXfbHrNCt1MTnHc3zBlTkjknCw2oAfj9YLEoZOixXri9hu6sUTPVsLPHfOGoydpc14v7Acd88eglgvIuJKUXkL3t3J3FR/uzLfY+A8PVaLl24Yj5veLMQne6swPiMGN7nEOXzxYWE5VnxaAo4DFk1Ix1OLRkOpYEb87BFJuPX8HPz58wP4rLgaq746gr3lLXj6mrGCfs/qDZ34+8bDWF9cDQCQyYDRg3QYkhSJrm4bDlUbcKqxDV8dqMVXB2oxKTsWd84ajNnDk/zecG02DjtLm/FZcRU2ltT43KyplMswITMW04cmYPrQBIwZpHP8jJ7gOA51BjMOVOlxoFqPA1UGHKzW9+pMGKaQIzNeixEpUchLjcbw5CgMSYpEemy41/nbu7pR2dKBypZ2VDSzsbKlA7WGTpgtNpi7rVAp5IjWqBClUSJGG4b4yDDERYQhTmsfI8MQH8HeR6qVsHFs3jazFcZOC/Qd7GXotEDfboG+o5u977Cgq9sGjUqOcJUCmjAFIsKU0IYpEKlWIkKtdIz8MYVchk6LFeZuGzotVrR1WWHq7IbJbMGQpCjJYi4EEoYgYBS465lHF66CsbNbcGZSg4ACejwpIl1JzviCb2uBZ0RKNLaeaMSRWoNAYbDHFwT8Es8ekYQXf2BPtxarDSofNxAA+PpALQyd3RgUE+4IXntj5rBEjEyNxqEaA97aVorlc4f7PL/basMjnx4AxwFXTRjkaHHqiWmDE/DwJcPx941HsPKLQxiZpvP7R/vWtlI8vuEQAODGKZlYeXm+mzUCADqtCs8tHodJ2XFYueEQvjlYh2N12/DKjQWOfSU9sdo4/OeX0/jnt8dgMndDJgOum5SJO2bmIquHRXi4xoA1W0uxvrgKu063YNfa3RiREoWlM3Mxb1SKm0XIcRz2V+qxsaQGG/ZVo9rlRp0SrcG0IfHITYiARqWAobMbpY1t2H26GTX6ThSebkbh6WY8u+kYojRKTM2Nx8TsWCRGqSGXyWDo7EZZYxuO1ZtwqFqPRlPvTC+ZDNDahbmty4ouqw0n6k04UW/CF/trHOepFDJkxGmhC1chIkyJrm4bTOZu1Bs7Pc4bCAq5zJGN2N/cPj2HhOFMR2gvBp4YrQqVLR2ChIHjOJcCesF3JTnjC0KFIQpbTzQ6mvr4ormty5EJNNFLRpIrY9NjEBcRhua2Luwpa/G458GVj/eyjWtXF6T3uqn2RCaT4Z45Q3DXe3uxdvtp3D4j1+eT91vbTuNwjQExWhUeXZDnd+3/d0Eu9lXo8WVJDe56bw++uPcCr0L+8k8n8I+vj9q/LgePLMjz+pQuk8lw45Qs5A/S4a5396C0sQ1X/Hsbll88DDdNzXJYMRzH4aejDfjnpqM4UMVSfsdmxOCvvxrlKJvek7zUaDx9zVg8OHc41mwrxXs7ynCk1ogH1u1DmKIEealRiI9Uw2TuxrE6I1pdWrtGaZRYkJ+KK8YPwnk5cR6vP8dxKG9ux9YTjdh6vBHbTzZB32HBt4fq8K2PzYwKuQxDkyIxKk2H/EHRyB+kQ15qtMNVa+62ot5gxskGE47UGnGkxoAjtUaUNrbB3G3DqYY2r3NHa5RIj9UiIy4c6bFapMeGI1UXjvAwBcIUclisTEgMHRY0t3eh2dSF5rYu9r6tC032f3dYrG6ioJDLEK1RQheuQnS4io0a9j46nB0PU8gdT//tXezVZu5mL7vV4freauOgUckdnR/D7enoUWqlIyW7ryBhCAJGgbueefjMJCE9lNu6rOi0sCB1QpR/90eyjt2MGk1mdFttPs12wJmqmuMn8MwzXERm0p4y5kYakhTp1f/vikIuw8xhifi0qAo/Hqn3KQw1+g5sPdEIgLlhhHDJqBQMSYrEiXoT3vmlDHfPHuLxvMqWdvzrO7bvYcX8EYLcfTKZDE9dPQZH64w4UW/Cb9/ZjbW/mez4vwaYC2bVV4fx+s+lAID7LxyKZRcNFeQrH5cRgy/uuwD3fVCErSca8cTGw3jhh+OYmBWLMKUc+yud7pZojRJ/mJ+H6yZl+BVMgGWyPbIgD3fPGoJ3dpzGB4UVqGrtwL4e/TfCVQrMyUvCpaNTMXtEkkfXWs9rkhUfgaz4CPz6vCxYbRwOVOmx9UQjjtYa0dzWBRvHIUKtREasFkOSIpGXylxDvuZWKxXIiNMiI06LWcOdO+ltNg7V+g6UN7fD1NmN9i4rwpRyaMMUSIhUIyNWC53WvxtOCB1dVug7LI751Ur5ORXzIGEIArwryV9lVR6HMAiwGHg3UkSYQlCwNz5CDaVchm4bh3qjGWkx4T7PF5qqysP3azhSawDHcT7/GPhsmkk+0lR7MntEEj4tqsIPR+qxwseT+vqianAcMDk7DpnxwkRNLpfhrlmDsfy/+/DGz6dw43lZvW4UHMfhj+sPoL3LisnZcbimwHfswpVItRKv3FiAq17ehqLyVlzzynY8ceVoTMyKxbE6E/725SH8fJyJ2SMLRuC3MwYLnhsA4iLC8Patk/Hxnko8990xVOs78eNRZ1A5IkyBX0/Jwm9n5AqKR/VEp1XhnjlDcffsIShtbMPxehNa2roQHqZAbkIkhqVEQq2UHkhWyGUYmxEjyAUpBblcZrcChP0+BEJ4mCIoQfUzFRKGIGASka4KALpw9vTcKsBicASeBcQXAPbHl6LToLKlAzX6Dp/CwHEcyppYeQuhFsOQpEjIZUBLuwX1RrMjPdYTzviCfzcSz8yhiVDIZTheb0JFczsy4nqvi+M4hxtpUcEgwXMDwOVj0/DyTydxot6EZ749ir9eke/2+Sd7q/DT0QaEKeRYtWi0oCduV4YkRWLd0qm46c1CHKsz4ZpXfoFKIYPFytwOaqUcTy4ajSvHC7NyeqKQy3DtpAxcXZCOoooWHK01wWqzISs+ApNz4vw+xQtBJpMhNzESuYl9664gzlxoH0MQENq9jUeMxdAoYg8DDy8GVa2+4wyNpi5HgFLoU5ZGpXAEqn3tZ+i0WFFSxVwR3nY8e0KnVaEgk1kY3x/27IfeX6nHiXoTNCq5owCfUJQKOVb+ahQA4N2dZfjJZUPdgSo9Hl1fAgC478IhGCzxxpiXGo2vl12AayemQ6OSw2LloJDLcMmoFGy8/wLJouCKXC5DQVYcbjgvEzdNzcaMYYlBEQWCAMhiCAp8AT0xwWdAoCvJJKxOkiuD7MJQ7Wc3bpk9VTVNFy7qpjIiNRonG9pwpMaAmcMSPZ6zv1IPi5VDYpQaGXG+3Vk9mTsqGYWnm/FpcTVuOT+n1+fv7mAppJeMSvFbtNAT0wYn4PrJmfigsBx3v7cXjy4cCZVChr99eRidFhsuGJqAO2d5jj8IJSFSjX9cPRYrf5WPRpMZCZFqunETZw1kMQQBR7qq6BiD//Q5aRYDc+/4E4ZSkamqPHkp/ru57S5zxhfEBuWuGD8ISrkM+ypacaLe/Xs0mcz4bB/LzRe6X8ATj10+EtOHJKCty4pHPi3BQ//bD32HBeMzY/DSDROgEOlC8oZGpUB6rJZEgTirIGEIAkaxFoOY4LPJf+e2nqQJtBj4PgxZAoO3PHwAmncVeWK3fWNbgYj4Ak9CpNqRbfL+TveWn//5pQxd3TaMHqTDhMwY0XPzqJUKvPWbSfj9JcMxIiUKI1KisPziYfjg/6a4ZRIRxECEXElBQGj3Nh7+xiMq+NwHMYbTjsCzOIthnD2r5ES9Cfp2S6/MHpuNc6SqislIcuWmqVn47nAd3i8sw52zBiMxSo0mkxlvbmWpnlLKQ/REpZDjrllDcFeAbiOCONcgiyEIOLKSBLqSosUEn6UIg06gxdDIWwzihCE+Uo1su5VRVNHS6/NDNQboOyyICFM4rAuxzBiagHEZMei02PD4hoPottqw4pMSmMzdyB8U7ejQRhBE8CFhCBCL1eboxyzUYhATfHYW0BMefOZjDPoOC9q8dHKTkqrqygR75tDe8tZen/Ebz6bkxvsta+ENmUyGxy4fBYVchi/212DSE9/h20N1UClkWHXlGNFppARBCIeEIUBcb7xCYwy8K4nfHu+LRqPwXgw8UfYiYAC8lpl2TVX1tFfAHxPsdVr2lvW2GH4+zjZdXTDUd/0if4zLiMFTi8ZApZChpd0CjUqO5xaPx+h0XUDzEgThG4oxBAgfX9Co5IKfjvnKiVYbB32HxWvGSpu522GNiN3JOigmHEdqjahq7cSQpN4F1067pKpK2c3KN93ZU9aCTovV8TN0dFkdFVWnD/WcyiqGqwvSMX1IAg7XGDA6XSdpRy9BEOIgiyFAnL0YhGeyyGQyQQFo3o0UrlL4LYndE3+ZSVJTVXmGJ0chJVqDDosVO0udjWR+PFqPrm4b0mPDMThR2tw9SdFpMHtEEokCQfQTJAwBwm9uE7rrmUdIyqqjc5uI+AIPv8mtornd4+cn7X2KcyXevGUymaPrmusO5S/tJZAXjkk9p4qKEcRAgoQhQBwF9EQKg5DMJCkZSTx8GW3eZdST43VMGIYme67rL4SL8pgwfLm/BuZuK5pMZnxnF4lLR6dJnpcgiNBCMYYAcZTDEOnqcbqSvO9+5je3SREGPtPIW2364/YdxUMDqOs+c1giUqI1qDV04rOiapxuYvXwx6SzOvoEQZydkMUQIAaRm9t4hKSsSimHwcOX0S5raoetR5cpvsUhEJgwKBVy/Ob8bADAI5+WYPXmkwCAu2YNITcSQZzFhFQYtmzZgssuuwxpaWmQyWRYv36936/ZvHkzCgoKoNFokJubi1deeaXvF+oDk4TgMyCswqqzc5v4GEN6bDgUchk6LFbUGd13QJ+sbwPHAfERYYKa0Pji1uk5GJ8Zg24bB44DrhiXhnmjkgOakyCI0BJSYWhra8PYsWPx0ksvCTq/tLQUCxYswAUXXICioiI88sgjuO+++/Dxxx/38Uq9I7bfM0+Mlt3sm9u8u5LE9mJwRaWQI9O+P4HPQOLh3UjBaA+oUsjx/u1T8I9FY/DvGybgn9eOI2uBIM5yQhpjmD9/PubPny/4/FdeeQWZmZl47rnnAAB5eXnYvXs3nnnmGSxatKiPVukbk8i2ngCA9makaNjX+RYGewE9iU/12fFalDa2obSxza2ZPd+Wc2hycBqxhIcpcO0k4Z3OCII4szmrYgy//PIL5s6d63Zs3rx52L17NywW/+Ul+gKxBfRQewB4fhyu/elCjJKdRlMfWQwAHI1mjvUoj73f3st39CDaQUwQRG/OKmGora1FcrK7/zo5ORnd3d1obGz0+DVmsxkGg8HtFUxEb3D75SXArIfS2o47lJ+jyX7z90QgwWcAGGXPDDpY7fyZbfaG7AAwelCMpHkJgji3OauEAUAv/zXHcR6P86xatQo6nc7xysgIrsvDZBYRY7DZgGPfOP45Q74fLSbPO5ONnRa0dbFyGMnR0oQhP41ZBIdqDI7MpNNNbTCau6FWyoPmSiII4tzirBKGlJQU1NbWuh2rr6+HUqlEfHy8x69ZsWIF9Hq941VRUeHxPKk4LAYhwtB6GuhoBuRKcEotdLJ2JHeWottq63VqnYFlEkVplNCGSQsF5SZGQqOSo73LilL7RrciezXUkWnRkiufEgRxbnNW3RmmTp2KTZs2uR379ttvMXHiRKhUnl05arUa0dHRbq9gIqokRt0hNiaNBNLGAQBGyMrR7GGTW62euZFSojWS16aQyzDKbjXwVVC3nXSWxCYIgvBESIXBZDKhuLgYxcXFAFg6anFxMcrLywGwp/0lS5Y4zr/jjjtQVlaG5cuX4/Dhw1izZg3efPNN/O53vwvF8gGIjDHUH2Zj0kjIUvIBACPk5WgyeRAGu8WQHIAwAMBUuwD8fLwRHMdhm71XwvmDAyuJTRDEuUtIhWH37t0YP348xo8fDwBYvnw5xo8fjz//+c8AgJqaGodIAEBOTg42btyIn376CePGjcNf//pXvPDCCyFLVQVcurcJsRgajrAxKY9ZDQCGyyo9CkNdkIRh5nBW+vrn4w3YWdqMOoMZ2jCFo2w2QRBET0K6j2HWrFmO4LEn1q5d2+vYzJkzsXfv3j5clXA6LVZ02eMDgmIMent8IzYLCI8DAGTI6nGgrXdmEi8MKbrAdiaPz4hBYpQaDUYzrnttBwDgkvwUrz0gCIIgzqoYw5mGybV7m5AAsb6KjboMIDYbAJAua0BTj5IVAFCrtwtDgBaDUiHHLdOyHf+WyYBbz88JaE6CIM5tqLpqADjjC0r/PYit3YCxmr3XpQPaBFihgFrWDXNLJYDBbqcHy5UEAP93QS5O1JtQWNqMO2cNRj5tbCMIwgckDAEgKr5gqgU4GyBXARFJgFwOoyYVMZ2VkLWU9To9WMFnAAhTyvGvxeMCnocgiIEBuZICwNGkR0gvBn0lG6PTADm77J2R6QAAhd5dGKw2ztG9LUUXuDAQBEGIQbQwfPfdd14/e/XVVwNazNmGwS4MfDc2n/DCoEt3Hotm75Vt7pv2ag2dsHGASiGjPscEQfQ7ooVh4cKFePDBB9HV5UyxbGhowGWXXYYVK1YEdXFnOnwvBZ0QYTDab/5RqY5DYbGs/WW4ud7tVL5Pc1oM66lAEATRn4gWhi1btmDDhg2YNGkSDh48iC+//BL5+fkwmUzYt29fX6zxjMXQIWLXc1sDGyMSHYe0CaxuU5y1yS3DiReGjFhtkFZKEAQhHNHCcN5556GoqAhjxoxBQUEBrrzySjz44IP44Ycfgl6g7kxHlMXQZq/+GuHccayJHQQASJK1oFbvLKZXYW+7mREXHqSVEgRBCEdS8Pno0aPYtWsX0tPToVQqceTIEbS3twd7bWc8ooShvbcw8G6lFFkLavTOvQyVdoshnSwGgiBCgGhhePLJJzF16lRcfPHFOHDgAHbt2uWwIH755Ze+WOMZi6jgswdXEi8MiWhFbauz/WZFi92VFEfCQBBE/yNaGJ5//nmsX78eL774IjQaDUaNGoXCwkJcddVVmDVrVh8s8cyFtxgkC0NkEmyQQyHj0Fxf7Thc1sSEIZOEgSCIECB6g1tJSQkSEtwrc6pUKjz99NO49NJLg7awswFxMYYmNmpdyl3LFegIi0dEVwP0dWUApkHfbkG9fQ/D4MSIIK+YIAjCP6Ithp6i4MrMmTMDWszZhmBh6GoHLHZXkavFAMAawVqVtjexOkrH61l/5lSdBlEage1CCYIgggjtfA4AZ7qqnxs4H3hWqAF1lNtHyugUAIDFWAebjcPxehMAYGiy+3mC2fkq8MxwYMMy1kqUIAhCJCQMEuE4DgbeYtD6EQZHfCGBlTd1QRPDLAadtRVVrR04UmMAAAxNktCPueEY8NXDrC7TnreAg5+In4MgiAEPCYNEOi02Ry8Gv64kD3sYeOSRzLWUIDNgf6Uee+09mcdmxIhf1J61AFz6W+x6Q/wcBEEMeEgYJMLHFxRyGSLC/DS94YVB6yE+Y485xMsM2HKsAYfsFsPELAkd1k7Y61hd8iQby3c4g94EQRACIWGQiGMPg0YJmcxPPaOOFjaGe7jZRyYBAOKhx7rdFbDaOGTFa5EWI3LXs6EaaDwKQAaMvQ5IHg2AA079KG4egiAGPCQMEhGVqtqpZ2N4TO/P7O6lRLnRcejSMam9z/NH5S42poxmApR9vvtxgiAIgZAwSETfLkYYWtmo8dA5ze5KylCzdNY0nQa3Tc8Vv6C6g2xMHcPG9ElsrNwtfi6CIAY01MFNIqLKYfAWgyam92d2YdBaWrD1oRlI1GmhVvqJWXii9gAbk/PZmD6RjTX7gO4uQBkmfk6CIAYkZDFIpEWMxdDRykZPFgO/E5qzIV1jliYKAFDXQxhisgC1DrBZgKYT0uYkCGJAQsIgkZY21qgoPkLAk7ivGINCBYTHsff8fgexmE1Aq709aPIoNspkQNII9r7+kLR5CYIYkJAwSKS5nQlDrCBhaGWjJ1cS4CyTIVUYWkrZGB4HaOOcx5Py2Fh/WNq8BEEMSEgYJNJsYsIQJ8Zi8ORKAoIgDKfZGJfjfjzRLgwNR6TNSxDEgISEQSIOi0ErQBj4GIMnVxLg3BHNb4QTvRi7xRCb7X6ctxj4jKVg0GkAjHXBm48giDMOEgaJCI4xdJuBbnvbzr62GGJ7WgzD2dhaxjKTAuXkD8CzecA/hwHf/inw+QiCOCMhYZBIi9AYA+9GgoxlCXmCFwZTvcTF2C2Gnq6kyGRAFQFwNmdwWiqdBuDTO4AuVv0V218Ajn0T2JwEQZyRkDBIwGbjHOmqfmMMvDCoowG5l8vdV64kmQyIs2+Waz4lbW6efR8Cpjo23+Sl7NiWZwKbkyCIMxISBgkYOi2w2lgV0xh/Jbcd8QUv1gLgtBjaJQiDtRvQV7D3PV1JABBvF4amk+LndqX4PTZOXgpc8CAgVwGVhUBtSWDzEgRxxkHCIIFme3whUq30vyHNX0YSEFiMwVgN2LoBRRgQ5aHGUjAsBn0lUFMMyOTA6GuAqGRg2Dz22aHPpc9LEMQZCQmDBHhhEJaq2spGb3sYABdhkGAxGGrYGJXq2VUVDGE4aa/QOqgAiLDv1M67jI1HvpQ+L0EQZyQkDBJoNJkBiBUGXxaDPcZgNgCWTnGLMVazMTrN8+dxg9nYHIAriS/dnTvbeWzYPGZB1B9kFgVBEOcMJAxisVlRZ2DCkBKt8X++vz0MABMNuT1WITbOYLALgyc3EuC0GFrLpaeslm1nY+5M57HwWCB1HHt/equ0eXvS1gjsWwdUUKlwggglJAxCqTsI/C0F+Nco1BvZU31ytNr/1/mqrMojk0nPTDL4sRiiUgCV1p6yWi5uboB1gDPa3VWpY90/y7mAjaU/i5+3J9VFwIsFwKe/Bd68CPjqDwDH+f86giCCDgmDUFThbKOa2eiwGJKEWAxCYgyAdGEwusQYPOGassrvdxBDvX3XdGw2oI5y/yzbLgynt4if15WudmDdEnat+J9j52qg+P3A5iUIQhIkDEJRR7Oxy4R6PWuqkyxIGARkJQHSM5P44HO0j65v/P4GKQFovpwGX87blcwpgEzBLJFA4gw7/g3oy4HodODuQmCOfVf1pj+zyrEEQfQrJAxCcXlaNhrYzV6QK0lIjAGQLgx88DnKiysJcMlMkmAx8H0ekkb2/kwd5SzzLbWFaHcXUPgGe3/hnwFNNHD+/WzN7Y1A4avS5iUIQjIkDEJRqtleAQDtxlYAQbYYtLwrSYQwcJyLxeBLGOwb3yRZDPZeDrwA9CTQFqJHNgCmWla+Y9SV7JhCBcz4PXu/6022iS9QWsqYBbL+LmD/fwGbLfA5CeIchYRBDHarwdppACBUGFrZ2Bcxho4WwMriHV5jDID0GIPN6uzl4MmVBAQuDPwGufE3urcfzb+KdbczVAHHvpY2N8+pn4DV04Btz7Md3J/8H/D+NYClI7B5CeIcJeTC8PLLLyMnJwcajQYFBQX4+WffGS7//ve/kZeXh/DwcAwfPhz/+c9/+mmlAMIiAQBR6IAuXCWsrafYGIOYdFVDFRu1Cb57OvOlMlpOs5u9UJpLWcBdGd67QB8PLww1xeLTYa0W5+a5YfPdP1OqmVgAwN63xc3rSuMJ4MMbWfG/jPOAafexn+fEd8DHt1PmE0F4IKTCsG7dOixbtgyPPvooioqKcMEFF2D+/PkoL/ecVrl69WqsWLECjz32GA4ePIjHH38cd999NzZs2NA/C7YHoCNlHciM0/o/32bz3dbTFSkxBiGBZwDQpbN9EtYuZ3qrEPiMpKQRgNxL6Y/4wcwa6u50xiOEUlEImPXMMhg0offn4+zCcOJ7abvCOQ74YhnQZQQypwE3bwDm/hW46RNAoQaOfAHsflP8vD3RVwLfPQa8cyXw3yVA8QdM9AjiLCWkwvDss8/itttuw+233468vDw899xzyMjIwOrVqz2e/84772Dp0qVYvHgxcnNzcd111+G2227DU0891T8LtruSItGBzHgBwtBlYvsHABFZSSJugEICzwC7qcdmsfdi4gyOjCQv8QWApcNKdScd/5aNQy7yLDyJw9jeCc4KHPxU3NwAcPAT4PTPzEK48hVmhQBA1jTg4sfZ+2//BOirxM/NU/I/4KXJwNZ/sX4Vhz4D1t8BvDYbaDwufV6CCCEhE4auri7s2bMHc+fOdTs+d+5cbN++3ePXmM1maDTufv3w8HAUFhbCYumHJzReGIRaDHx8QaFm+yB8EeESfBbq3hBqMQDS4gy8MCT5EAbARRhEZiYd38TGoXO9nzP6WjaWfCRubo4Dfv4Xez99mVMYeSYvBTKmAJZ2YJPEpkMnvmfuKEsbc1Nd9gIw6xHWe7uuBHjz4sB2cXMcc3n9dwnwzzzg74OAlyYBXz1MokP0KcpQfePGxkZYrVYkJye7HU9OTkZtba3Hr5k3bx7eeOMNXHHFFZgwYQL27NmDNWvWwGKxoLGxEampvW+QZrMZZrPZ8W+DwSB90XZhiEIHRqRE+TkZwuMLgFMYujuZpdFzM5knhFoMgDPOEGyLAQDSJ7KxslD43PpK5qqSyYHBc7yfl78I+PaPQMVOFiPp2XPCGye/ZzdnVQQw+be9P5fLgQX/AF6dCRz4GJh4K5A9Xfj6W8qAj28DwDGX1+UvOK2eib8BPrgOqNoDvLsIuOULIHWM8LkB5vL7/F4mDK40HmOvna+yGMzFKwFtnLi5ASY65b+wIohVe9lGSaWauR3TJwEjFrKEA5lM/NyudLSyhx1LO4vRRaUAYRGBzUn0OSETBh5Zj188juN6HeP505/+hNraWkyZMgUcxyE5ORm33HIL/vGPf0Ch8OwDX7VqFR5//PGgrNUWFgk5mCtpVJqAm73QPQwA+2NRadkfUFujMGGQYjEI3ctgNjmtC0HCIGM3blM9EJnkf37ejZQ+yfeNLTqV1Wg69ROw/yNg5kMCFg+WgQQABTd7nz91LLuJ714DbPw9sHQLoBDwJ2HpANbdyLLC0iYAlz7r7gqLTGLxjHcXsZvvO1cCt34NJAwVtvaKQja/qY6lSE+8laXyRiSyLLHi94CjG4Gid9g47+/AmMXCbuIcxzrvbXkaqPLg+ms4wsTop1VA4ghg7HVsbl/p0G7XphMo386sqZM/APWHep8TkwWkjAbSxrPYUtp4VnvLH9ZuVgyy/hBQfwRoOsH+D8wGQK5kVnlEIhCTyb5HbBZ7kIge5D1G5ja/BWhvYkLW1sD+Dtsa2DG5kv19hkWwmJpGx/6uNTrnv1VeshQ5jv3OmI32l97lvYk9HClU7P9aafcuqMLZ91Nq2Mj/W8jvZxAImTAkJCRAoVD0sg7q6+t7WRE84eHhWLNmDV599VXU1dUhNTUVr732GqKiopCQkODxa1asWIHly5c7/m0wGJCRkSFpzVsrzJgBIF5lRk6CgKceMRYDwKyG1nL2C+ktC8gVo4A9DDyOvQwChYFPU41Mdloz3tDogKQ89gdbUQjkXep/focb6WL/5465zi4MHwIzfuf/Bli1Fyjdwv6Yp9zl+9w5f2Lxi/qDTCDO82BduMJxwBfLgdr9LBts8TvO2IUrYRHADeuAtZeyc9++DLjlSxas98Xed4Avl7NEgaSRwLXvAAlDnJ/HD2bXt3wH8MUD7Jp/uhTY9wFw6b+cDwA9sVlZsH3L087mSkoNE5zcWexmarUwa+TUT+z/p+EIC6p/v5KdM+7XwPAFQJjWfd66g6yQ4snvgdPbnD3OedQ6dmPrMrFXaxl7HfnCeU5sjl0kJrC1yBWsnayplrnNakvYeqwSCkHKVUBMBhOL8BjHfiR0GtjfKC8GHc3i53ZFqXGKBGdlpV662uyxRhHZgL6Qq5hATLoVuOix4MzpgZAJQ1hYGAoKCrBp0yZceeWVjuObNm3Cr371K59fq1KpkJ6eDgD48MMPcemll0LupW2mWq2GWi1gh7IfKprbsavaghkqYFgsoJALeDoTuoeBJyLRLgwCM5P4dFUhriTXGAPH+b+51gt0I/FkTLYLw07/wtBtZjcfwHd8gSfvMnazbDrB3DO868ob255jY/7V7IbgC20cE4cvlwM//o3tn/AlhLvfBPa9z57yrl7DXC/e0OiAmz5l4tBwGHj7cuA3X3p2h1m7gW8ece70zrscuGI1oI70PHfmFGbhbH8B2PwPdj1fngpMvI25mBJHsP9jYw1zF+18hV0/wO5eux2Yek9v6y53JjD5/9gN89BnrF5V+S/s6f/kD+znjs1hP1t3J/t95fuA80SlMvfg4DmsVDvfwwNgRRnrDwI1+4HqvUzEW0qdrwMfe7+eAHNHJY5gmXIJw9nfjDrKeSM21bI1tZxm7r7WcsBmYS5UIW5UmZwJfkQi+z2ISGS/Izar8yZvNjBvQKfe+QLHroeplr08T84yG9VRLq9I9vdotTDR6+5kL0sH8x5YOtkIe9zRZmEWhy0Imz59EFJX0vLly3HTTTdh4sSJmDp1Kl577TWUl5fjjjvuAMCe9quqqhx7FY4dO4bCwkKcd955aGlpwbPPPosDBw7g7bcDyHMXiL7Dgri4eMAIjEkUYJYCEiwGESmrlg5mRgPCXEkxmeyXvsvE5vfn7hEaX+DJOA/Ys5ZZDP4o28Z+2SNTgBQBvnd1JBOH/evYk7EvYWg66dw0d/59gpaOglvY2mv3A98/Dlz+oufzKgpZ1VeAPa25liH3RkQCcPPnwNqF7Gn89QuBRa+7x1UaTwCf38NuwAALYM94yHuPcB6FirVZHXkFsx5KN7O6Uzv+zTKx5EqWqsuj0bF4y5S7/MclNDpgwhL2aj7Fen7v+4DdaHv29giLZEKVOwsYfCGzHr09eETEAzkz2IunvZntg6nay6rstjWyG586kv2OxGYDKfks5hGT5f+6uGKzMnHkhaLLxB5M+J9RE+0iBInMpSVmfoClpXcZXcSi1cX1FMmsR3UUG6XEbDiOrdnSbheMDmGu5gAIqTAsXrwYTU1NWLlyJWpqapCfn4+NGzciK4tlkNTU1LjtabBarfjnP/+Jo0ePQqVSYfbs2di+fTuys7P7fK35g3TInzMG+AyI4NqFfZGYGAMgriwG70ZShguzSJRqVqROX87cSUKFwV9GEk/GeWysLmIb3XxtuHN1Iwn9Qxl7HROGAx8D81Z5n3/7CwA4ZokIFTW5AljwNLBmHrD3P+xpvaeLS1/J/P42CzDyV2yjnFAik4AlnwPvXcMC4u9cCWRNB9LGsaf445vYE29YJEur5bvjCSV+MLDkMxYbKHydudEc7hwZE9L8q4Hxv5Z2Q4nLBWY/AsxaweIeDUfZTUqhAnQZrBlUIL5vbZzTwgg2cgWz6nTp4pILRH0PuV1kBD4AikUmY/ELbzGMPiDkwee77roLd93l2Q+8du1at3/n5eWhqKioH1blBf6Pymz0fR6PlBgDIGwvg2vgWejNNS7HLgyngMzzvJ/HceIthrhctlGtvYk9eft6qucDz0LcSDw5M5mLwlgDHPuK3Zx7Yqhxluqe/oDwuQH2xDvpdmDXG6xkxq3fAInD2WemeuC9a9lNMWkU8Kt/i3/yi04Fbt/E9k3sXgOUbWUvnqHzgIXPMMtOCjIZE7OhFzO3hL6CPS3r0v2nSov5HlEp7EWc04RcGM4qRAtDKxvFxBgAYWUxHH0YBGaLAEwYSjf738tgqGZrlymcN0d/yGTMaji6kcUZvAlD00n2lCxXMteDUOQKYNwNwM//BH552bMw7Pg389NmTGGb2MQy7+8shlFdBKy5BJj5MLO0tjwDGCqBiCTghg+lm/GqcHbzn3Yvq/+kr2BzDrlQuAALQaHyHoQmCAGEvFbSWQXfk8EscC9EX8YYHJ3bBMQXeBwpq36CcLy1kDDUc8aNN/iNbuU7vJ/Du5EypzL/rhgm/5ZlZVTs6B3LMFSzSqwA87tLQakGbvyEpU92NANfP8xKahgqmZ/7N19Jf6J3JTYLOG8pMPdvLA4STFEgiCBAwiAGexE9wRaD2BiDGFeSv85tnnBscvNjMYjNSOLJnMrGsm3ey1of+4qNw+aJmxtgLoyxi9n7TX9x3yH+3WPM751xnrAUWG9o44BbvwUueZJ1qMs6n2Ut3bndPW2UIM5hyJUkBldXkpCUT8kxBjEWgxhXkkiLwVNzHl+kTwTColicoaa4d2G8Tj3LcwdYPrwUZv4BOPAJ20S181Vgyh3AvnUsMA0ZcMmqwHfrKsOAKXeyF0EMQMhiEAMvDLZulmvsD6kxhrZG/41kJFkM2WzsaHZaM57w1c7TFwqVM4XzxPe9Pz/xPcvqSRjmf6OXN2IygDl/ZO+//gOwZj4rWgcAFywHBhVIm5cgCAckDGLgXUmAsF7EYi0GPl2VszpFxRtCOrf1RB3JdjID3gPQ3WaWbw8AySItBoAFUoHeNX4AttEKAIZdIn5eV6bcBUy5GwDHLAfOBky4GZj9aGDzEgQBgFxJ4pDLmauky8gC0JGJ3s/t7rLvWITwGIMyjIlIp55ZDd42IXGci8UgMnUwNoelXTaXsiBrTxqOMItIE8Ny1MUyxO7fryxk7i5euDoNTmEYdYX4eV2RyYBL/g6MXsQ2RaWOAzImBTYnQRAOyGIQi9CUVd5agIzVihGKkMyk9ibmkgHYzlAx8C4cb2Wba+3NdlJGS/PVx2SwIDRncy+VffBTtukqYTirhxMMBhWw8g0kCgQRVEgYxKIWmJnEu4LU0eK22AvZ/cwHniMSfe8w9kRSHhv5zKOe8AXWUkaLm9eVsdezcc/brAaQzQbseJkdG3dD4MFhgiD6FBIGsQi1GBypqiK3yQvJTJISeObhA8q8ZdCTYAhD/lWsWU3zSdavuegd5qJS61iZa4IgzmhIGMQi1pUkNCOJR0iLTympqjy8MDSfYtUiXeG44AiDOoqVxwZY1dIN9rpCMx/qu3oyBEEEDRIGsTiEwc/uZ0eqqsgbIV/czlTn/ZxALIbIRHtmEufsucDTWsZK+spVLBYQCOfdybqv8eQv8t8bgSCIMwLKShKLoyyGP1eSvSS20IwkHv5mz9/8PRGIxQCwHc2mOqDugHtNI77MROpY8bGLnsjlwKI3Wc1/mYxlDlFsgSDOCshiEAtvMfRsTtITsXsYeKIHsZFvwuMJqamqPHwPhKo97scrdrKRL6EdKDKZs3UjiQJBnDWQMIhFcIyhlY1iYwx8UTyDL4tBQmVVVzKnsLHsF/fjDmGYLG1egiDOCUgYxCK0kJ7YAno8vMXQ3ujsNNUTQyUbdYPEzc3DWwRNx1mvAYB10eJLYZAwEMSAhoRBLHypaMcGNi9IzUoKj2VNxQHPcQZHj1lI25kMsB3VfIG8MntRu5M/sE1pSSOlxy4IgjgnIGEQC3+j91WEDpDuSpLJnAFoPsjsSmsFG8PjvDeKFwLfRvHwBjYetZfDFtNVjSCIcxISBrHwwWR/FoNUVxLgEoD2IAx6uzAE2jBm1FVsPPo1a/DOC0Te5YHNSxDEWQ8Jg1j4G72/6qdSs5IApyvHo8VQzsYYiW4knkETWPlrSxvwYgFgNbNNbT17KBAEMeAgYRBLX7uSAJfMJF/CkCV+XldkMuDCv7D31i42XvQYpZUSBEEb3EQTHstGSxtgtbDmND2x2ViZaSBAV1Jl7894YZAaeHYl71LgqteZG2nUlcCQiwKfkyCIsx4SBrG4uoY6Wj33ZDAbAHC9zxcKbw20nO79WbBiDDxjrmUvgiAIO+RKEotc4SyL4S3OwB9XhgNKtfjv4ejNfNq94T3gFItgCQNBEEQPSBikwMcNvGUmBZKRBACxWQBkrFOca5XVtiZnDSZePAiCIIIMCYMU+B4L3gLQgWQkAczK0KWz9669mZvsXdd0GUCYVtrcBEEQfiBhkILDYmj1/HkgGUk8cTlsbD7lPMa344wfIn1egiAIP5AwSIG3BHi3Tk/am9mojZP+PWI9CANvMSQMlT4vQRCEH0gYpOBvk1t7ExsDEQb+5u/aTKfhGBvJYiAIog8hYZCCv01ujiY9AQgD31qzdr/zWM0+988IgiD6ANrHIAXBFkO89O/BN9NpOc2C2d1mwFgNQEbCQBBEn0LCIAV/6arBEAZtHMs+0lcANfsBSwc7Hj/E2SyIIAiiDyBhkAJfFsObKykYwWeAFbTTV7CeCXxjoGC13SQIgvACxRik4C/GEAyLAQAGX8jG498CJ763H5sd2JwEQRB+IItBCrwlwAtATxwWQ4DCMPRiQCYHqvawfyvUzgY7BEEQfQRZDFKISGBje2PvWkZWC2C2xx4CyUoCWF8GvqEOAIy7PnD3FEEQhB/IYpCC1i4M1i7m++f7QAMum95k0mslubLgaVb+QhHG+iUQBEH0MSQMUgjTAqoI1pOhrcFdGHj3Ungsq8QaKNo44PIXA5+HIAhCIORKkorDndQjzhCsjCSCIIgQQcIgFV4Y2hrcj7fby2QHGngmCIIIESEXhpdffhk5OTnQaDQoKCjAzz//7PP89957D2PHjoVWq0Vqaip+85vfoKnJS3ZQXxJh79zm2i8BAEz1bIxM6t/1EARBBImQCsO6deuwbNkyPProoygqKsIFF1yA+fPno7y83OP5W7duxZIlS3Dbbbfh4MGD+Oijj7Br1y7cfvvt/bxyOAPQ7T2EwVjLxsiU/l0PQRBEkAipMDz77LO47bbbcPvttyMvLw/PPfccMjIysHr1ao/n79ixA9nZ2bjvvvuQk5OD6dOnY+nSpdi9e3c/rxwurqSeFkMdGyOT+3c9BEEQQSJkwtDV1YU9e/Zg7ty5bsfnzp2L7du3e/yaadOmobKyEhs3bgTHcairq8P//vc/LFy4sD+W7I5XYSBXEkEQZzchE4bGxkZYrVYkJ7s/WScnJ6O2ttbj10ybNg3vvfceFi9ejLCwMKSkpCAmJgYvvug9ndNsNsNgMLi9goIjxtAj+Gyyrz2KXEkEQZydhDz4LJPJ3P7NcVyvYzyHDh3Cfffdhz//+c/Ys2cPvv76a5SWluKOO+7wOv+qVaug0+kcr4yMjOAsnLcIjD1EjCwGgiDOckImDAkJCVAoFL2sg/r6+l5WBM+qVatw/vnn46GHHsKYMWMwb948vPzyy1izZg1qamo8fs2KFSug1+sdr4qKiuD8ANHpbDRUO4/ZbC7CQDEGgiDOTkImDGFhYSgoKMCmTZvcjm/atAnTpk3z+DXt7e2Qy92XrFCw3cVcz5pFdtRqNaKjo91eQSE6jY1mvbMkdnsTwFkByJyuJoIgiLOMkLqSli9fjjfeeANr1qzB4cOH8cADD6C8vNzhGlqxYgWWLFniOP+yyy7DJ598gtWrV+PUqVPYtm0b7rvvPkyePBlpaWn9u3h1JKDRsff6Kvtot0aiUgCFqn/XQxAEESRCWitp8eLFaGpqwsqVK1FTU4P8/Hxs3LgRWVlZAICamhq3PQ233HILjEYjXnrpJTz44IOIiYnBnDlz8NRTT4XmB4hOZ13cDFVA0gig1b5WXZDiGARBECFAxnnzwZyjGAwG6HQ66PX6wN1K714NnNjEitxNWAJsewHY9Ccg/2rg6jeDs2CCIIggIObeF/KspLMa3SA28q4k3mKIyQzNegiCIIIACUMgxDCXF1pK2cjHGEgYCII4iyFhCISEYWxsOMrGxuNsjMsJzXoIgiCCAAlDICQOZ2PjccDS4bQcEvNCtyaCIIgAIWEIhNhsQK5kndxO/ghwNta5jXY9EwRxFkPCEAgKFZA4gr3f9QYbk/MBLyU9CIIgzgZIGAIl63w2nvyejZlTQ7cWgiCIIEDCECi5s9z/PXhOSJZBEAQRLEgYAmXoXCBuMHufOhbInBLa9RAEQQRISEtinBMolMBvvgJKNwNDL6b4AkEQZz0kDMEgKhkYc22oV0EQBBEUyJVEEARBuEHCQBAEQbhBwkAQBEG4QcJAEARBuEHCQBAEQbhBwkAQBEG4QcJAEARBuDHg9jHwnUwNBkOIV0IQBNF/8Pc8Id2cB5wwGI1GAEBGRkaIV0IQBNH/GI1G6HQ6n+fIOCHycQ5hs9lQXV2NqKgoyESWrzAYDMjIyEBFRYXfZtoDEbo+/qFr5Bu6Pv6Reo04joPRaERaWhrkct9RhAFnMcjlcqSnpwc0R3R0NP3S+oCuj3/oGvmGro9/pFwjf5YCDwWfCYIgCDdIGAiCIAg3SBhEoFar8Ze//AVqtTrUSzkjoevjH7pGvqHr45/+uEYDLvhMEARB+IYsBoIgCMINEgaCIAjCDRIGgiAIwg0Shh6sWrUKkyZNQlRUFJKSknDFFVfg6NGjbudwHIfHHnsMaWlpCA8Px6xZs3Dw4MEQrTi0rFq1CjKZDMuWLXMco+sDVFVV4cYbb0R8fDy0Wi3GjRuHPXv2OD4fyNeou7sbf/zjH5GTk4Pw8HDk5uZi5cqVsNlsjnMG2vXZsmULLrvsMqSlpUEmk2H9+vVunwu5HmazGffeey8SEhIQERGByy+/HJWVldIWxBFuzJs3j3vrrbe4AwcOcMXFxdzChQu5zMxMzmQyOc558sknuaioKO7jjz/mSkpKuMWLF3OpqamcwWAI4cr7n8LCQi47O5sbM2YMd//99zuOD/Tr09zczGVlZXG33HILt3PnTq60tJT77rvvuBMnTjjOGcjX6G9/+xsXHx/PffHFF1xpaSn30UcfcZGRkdxzzz3nOGegXZ+NGzdyjz76KPfxxx9zALhPP/3U7XMh1+OOO+7gBg0axG3atInbu3cvN3v2bG7s2LFcd3e36PWQMPihvr6eA8Bt3ryZ4ziOs9lsXEpKCvfkk086zuns7OR0Oh33yiuvhGqZ/Y7RaOSGDh3Kbdq0iZs5c6ZDGOj6cNzDDz/MTZ8+3evnA/0aLVy4kLv11lvdjl111VXcjTfeyHEcXZ+ewiDkerS2tnIqlYr78MMPHedUVVVxcrmc+/rrr0WvgVxJftDr9QCAuLg4AEBpaSlqa2sxd+5cxzlqtRozZ87E9u3bQ7LGUHD33Xdj4cKFuOiii9yO0/UBPv/8c0ycOBHXXHMNkpKSMH78eLz++uuOzwf6NZo+fTq+//57HDt2DACwb98+bN26FQsWLABA16cnQq7Hnj17YLFY3M5JS0tDfn6+pGs24GoliYHjOCxfvhzTp09Hfn4+AKC2thYAkJyc7HZucnIyysrK+n2NoeDDDz/E3r17sWvXrl6f0fUBTp06hdWrV2P58uV45JFHUFhYiPvuuw9qtRpLliwZ8Nfo4Ycfhl6vx4gRI6BQKGC1WvHEE0/g+uuvB0C/Qz0Rcj1qa2sRFhaG2NjYXufwXy8GEgYf3HPPPdi/fz+2bt3a67OelVk5jhNdrfVspKKiAvfffz++/fZbaDQar+cN1OsDsAq+EydOxN///ncAwPjx43Hw4EGsXr0aS5YscZw3UK/RunXr8O677+L999/HqFGjUFxcjGXLliEtLQ0333yz47yBen28IeV6SL1m5Erywr333ovPP/8cP/74o1s11pSUFADopcL19fW9FP1cZM+ePaivr0dBQQGUSiWUSiU2b96MF154AUql0nENBur1AYDU1FSMHDnS7VheXh7Ky8sB0O/QQw89hD/84Q+47rrrMHr0aNx000144IEHsGrVKgB0fXoi5HqkpKSgq6sLLS0tXs8RAwlDDziOwz333INPPvkEP/zwA3Jyctw+z8nJQUpKCjZt2uQ41tXVhc2bN2PatGn9vdx+58ILL0RJSQmKi4sdr4kTJ+LXv/41iouLkZubO6CvDwCcf/75vVKcjx07hqysLAD0O9Te3t6rH4BCoXCkqw7069MTIdejoKAAKpXK7ZyamhocOHBA2jUTHzM/t7nzzjs5nU7H/fTTT1xNTY3j1d7e7jjnySef5HQ6HffJJ59wJSUl3PXXX39Op9L5wzUriePo+hQWFnJKpZJ74oknuOPHj3Pvvfcep9VquXfffddxzkC+RjfffDM3aNAgR7rqJ598wiUkJHC///3vHecMtOtjNBq5oqIirqioiAPAPfvss1xRURFXVlbGcZyw63HHHXdw6enp3Hfffcft3buXmzNnDqWrBgsAHl9vvfWW4xybzcb95S9/4VJSUji1Ws3NmDGDKykpCd2iQ0xPYaDrw3EbNmzg8vPzObVazY0YMYJ77bXX3D4fyNfIYDBw999/P5eZmclpNBouNzeXe/TRRzmz2ew4Z6Bdnx9//NHjfefmm2/mOE7Y9ejo6ODuueceLi4ujgsPD+cuvfRSrry8XNJ6qLoqQRAE4QbFGAiCIAg3SBgIgiAIN0gYCIIgCDdIGAiCIAg3SBgIgiAIN0gYCIIgCDdIGAiCIAg3SBgIgiAIN0gYCIIgCDdIGAiin5g1a5Zbb2yCOFMhYSAIgiDcoFpJBNEP3HLLLXj77bfdjpWWliI7Ozs0CyIIH5AwEEQ/oNfrMX/+fOTn52PlypUAgMTERCgUihCvjCB6Q609CaIf0Ol0CAsLg1ardXTkIogzFYoxEARBEG6QMBAEQRBukDAQRD8RFhYGq9Ua6mUQhF9IGAiin8jOzsbOnTtx+vRpNDY2wmazhXpJBOEREgaC6Cd+97vfQaFQYOTIkUhMTER5eXmol0QQHqF0VYIgCMINshgIgiAIN0gYCIIgCDdIGAiCIAg3SBgIgiAIN0gYCIIgCDdIGAiCIAg3SBgIgiAIN0gYCIIgCDdIGAiCIAg3SBgIgiAIN0gYCIIgCDdIGAiCIAg3/h9MJ8EKx2IrMQAAAABJRU5ErkJggg==",
"text/plain": [
""
]
},
"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()"
]
},
{
"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(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.02651486 0.02314122 1.0525634 0.88261193]\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",
"prior = BoxUniform(low=prior_min, high=prior_max, seed=seed)\n"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [],
"source": [
"theta = prior.sample(num_sim)\n",
"control_list = [{\"a\": theta[i, 0], \"b\": theta[i, 1]} for i in range(num_sim)]"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [],
"source": [
"stat_vec = batch_run(parameters, control_list, cfg, n_workers=4)"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [],
"source": [
"scaler = StandardScaler()\n",
"stat_vec = scaler.fit_transform(np.array(stat_vec))"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"((2000, 2), (2000, 4))"
]
},
"execution_count": 15,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# stat_vec = np.array(stat_vec)\n",
"theta.shape, stat_vec.shape"
]
},
{
"cell_type": "code",
"execution_count": 31,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Inferred dimensions: param_dim=2, feature_dim=4\n"
]
},
{
"data": {
"application/vnd.jupyter.widget-view+json": {
"model_id": "dceb7b9ee30648e49c44274c0c9ecfcd",
"version_major": 2,
"version_minor": 0
},
"text/plain": [
"Training: 0%| | 0/500 [00:00"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"print(\"Plotting posterior marginals for MDN ...\")\n",
"\n",
"param_labels = [\"a\", \"b\"]\n",
"\n",
"fig = corner.corner(\n",
" samples,\n",
" labels=param_labels,\n",
" truths=theta_true_np,\n",
" quantiles=[0.16, 0.5, 0.84],\n",
" show_titles=True,\n",
" title_kwargs={\"fontsize\": 12},\n",
" truth_color=\"red\",\n",
" bins=50, # Increase number of bins for better resolution\n",
" range=[(0.08, 0.12), (0.08, 0.12)], # Adjust ranges\n",
" fig=plt.figure(figsize=(6, 6)),\n",
")\n",
"\n",
"plt.tight_layout()\n",
"plt.show()"
]
},
{
"cell_type": "code",
"execution_count": 36,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAaAAAAHNCAYAAABYXK/dAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAV3ElEQVR4nO3dfWxd9Zng8edcm7zQxJuQtzpLgHQ266XJDgthAKMGkQ0TCVSkVDsrtJUof6zURdpWTbP8AaIrJLRStCqjohbB4lYritC2aBRAXYGqibpkQ5fSaWhQO6REIZOpTWpqnIQ4b8b29d0/QlwyvMQ2tp97jj8fKWrvxc59chT4+ve7555TNBqNRgDADKtlDwDA7CRAAKQQIABSCBAAKQQIgBQCBEAKAQIghQABkKI1ewAg4i9r/z57BPhUdo7+zYS/xwoIgBQCBEAKAQIghQABkEKAAEghQACkSDsNu9FoxOnh0xERcfFFF0dRFFmjAJAgbQV0evh0LNi+IBZsXzAWIgBmD1twAKQQIABSCBAAKQQIgBQCBEAKAQIghQABkEKAAEghQACkECAAUggQACkECIAUAgRACgECIIUAAZBCgABIIUAApBAgAFIIEAApBAiAFGkBajQaWS8NQBNIC9AfTwxmvTQATcAWHAApBAiAFAIEQAoBAiBFWoBGnQQHMKvlnQV3/EzWSwPQBGzBAZBCgABIIUAApBAgAFIIEAApBAiAFAIEQAoBAiCFAAGQQoAASCFAAKQQIABSCBAAKQQIgBRpATozPJr10gA0gbQAvfH2QNZLA9AEbMEBkEKAAEghQACkECAAUggQACkECIAUAgRACgECIIUAAZAiLUCnh+pZLw1AE0gL0PHTw1kvDUATsAUHQAoBAiCFAAGQQoAASNEUARoYdEICwGzTFAGq1xvZIwAww5oiQADMPgIEQAoBAiCFAAGQQoAASCFAAKQQIABSCBAAKQQIgBQCBEAKAQIgRd4tuYfdkhtgNksL0MCZP10Bu+/kYNYYACRpii24+qirYQPMNk0RIABmHwECIIUAAZBCgABIIUAApBAgAFIIEAApBAiAFAIEQAoBAiCFAAGQQoAASCFAAKRIC9BI3RWwAWaztAD1n3ov66UBaAJpARptWAEBzGZpAaoVRdZLA9AEmuIkhNPv1bNHAGCGNUWABgaHs0cAYIY1RYAAmH0ECIAUAgRACgECIIUAAZBCgABIIUAApBAgAFKkBWhoZDTrpQFoAmkBOnLS1bABZrO8+wFZAAHMamkBaq25GjbAbOYkBABSCBAAKQQIgBQCBEAKAQIgRVMEqJE9AAAzrikC1Hd8MHsEAGZYUwQIgNlHgABIIUAApEgL0HDdxeAAZrOUANVHGzEwOJLx0gA0iZQAjTaceA0w23kPCIAUAgRAiqYI0JmRevYIAMywpgjQu6eHs0cAYIY1RYAAmH0ECIAUKQHyIVQAUgLUf2Io42UBaCK24ABIIUAApBAgAFIIEAApBAiAFAIEQAoBAiCFAAGQQoAASNEUARoYdDVsgNkmJUBD9fPv/3P6PfcDAphtcq4Fd9K14ABmu6bYgnvnxHvZIwAww5oiQCOjjewRAJhhKQGqCw7ArJcSoL4TgxkvC0ATaYotOODTeffS/vi7//h/4t1L+7NHgXETIKiAnusOxsClR6PnuoPZo8C4CRCU3ND896Lv829FRETf59+KofnOKqUcBAhKrGXx4vhjZ180amdP7Gm0RPz9/UX8w3/vTJ4MLqw1ewBg/AYXnomhz/zpJJ6Wtnp0//m+D3xFIwZGfh7zW/4sBj57bOzZOafmxbwT82dwUriwpgjQaMNp2TAef//vfhnvXv5PTjRoRETxp4cjjaNxePA7cfg//em5Rf+4NK794c0zMCGMX1NswTktG8Zn5a9XR224djY65xQf++URjYjacC1W7l093aPBhDXFCqhetwKC8Vj5m8uj7Q+L4zd3vBynLzn5yT9CjkZcfHRh/PnTnbGgv23GZoTxaooV0PHBkewRoDQW9LfF9V23xIrXL/3Er1vx+qq4/vFN4kPTSgnQKwePnvf49JDbMcBEtAy3xuLuZedvxX1AERGLf780WkaaYpMDPlJKgJ7e03Pe49baJ21iAx9loP1YFKPv/7tzLkTn/rdexMDKdxOmgvFrii04YOKOX3o0Gi2NKOpF1EZqcdkv1kRtpBZFvYhGSyOOX3oke0T4RCkB6lixMONloTLqLfU4tfRERETMP7Ygruu6Jf7l314V13XdEvOPfSYiIk4tPRH1FtvbNK+UDeK/WL04ipbB6Dl24a8FPmz0onos6GuLhb2L41+98G/G3us5d4LCG7fujROffTdGW+vRUm9JnhY+mncooYQuGpwT1z9+SxQf8SGgluHWWPuTv4hGND7yn0Oz8B4QlNSF4iI+NDsBAiCFAAGQIuU9oKde6Y7RGIx4/+K8I6ONaDQaURS2DABmi6ZZAe3rHcgeAYAZ1DQBGnZBUoBZJS1Aa5Z/5rzHf/23+5MmASDDpAP0D++cjAf/97741T+evbDo6Ggjeo6evuD3vXXs7Ne01M5/6ZcO9H/UlwNQUZM+CeHf/vX/jYiI//n/Dl3wa//Ddati0cVz4sU3+sZONNj8+RWx8xdn//lll8yPt4424m/29ERRFPHs3rfis23zo/PPlkx2PPhE/3yRv1+QrWg03A8bgJnXNCchADC7CBAAKQQIgBQCBEAKAQIgxaROw240GnHixImpngVm3MKFC12DEJJMKkD9/f2xfPnyqZ4FZlxfX18sW7YsewyYlSYVoDlz5kRERE9PT7S1tU3pQJw1MDAQq1atcoynybnje+7vMjDzJhWgc1sWbW1t/uM4zRzj6WX7DfI4CQGAFAIEQIpJBWju3LnxwAMPxNy5c6d6Ht7nGE8vxxfyuRgpAClswQGQQoAASCFAAKQQIABSCBAAKSYVoEcffTRWr14d8+bNi/Xr18dLL7001XNV3kSO4a5du6Ioig/9euONN2Zw4mrYvXt33H777bFy5cooiiKee+657JFg1ppwgJ5++unYunVr3H///bF3797YsGFD3HrrrdHd3T0d81XSZI/h/v37o7e3d+zXmjVrZmji6jh16lRcddVV8cgjj2SPArPehD8HdP3118c111wTjz322NhzV155ZWzZsiW2b98+5QNW0USP4a5du2Ljxo1x7NixWLRo0QxOWm1FUcSzzz4bW7ZsyR4FZqUJrYCGhobi1Vdfjc2bN5/3/ObNm+Pll1+e0sGq6tMcw6uvvjra29tj06ZN8eKLL07nmADTbkIB6u/vj3q9HitWrDjv+RUrVsTbb789pYNV1WSOYXt7e3R1dcWOHTvimWeeiY6Ojti0aVPs3r17JkYGmBaf6nYM5zQaDZe1n6CJHMOOjo7o6OgYe9zZ2Rk9PT3x0EMPxU033TStcwJMlwmtgJYuXRotLS0f+km9r6/vQz/R89Gm6hjecMMNceDAgakeD2DGTChAc+bMifXr18fOnTvPe37nzp1x4403TulgVTVVx3Dv3r3R3t4+1eMBzJgJb8Ft27Yt7rzzzrj22mujs7Mzurq6oru7O+6+++7pmK+SLnQM77vvvjh8+HA8+eSTERHx8MMPxxVXXBFr166NoaGheOqpp2LHjh2xY8eOzD9GKZ08eTLefPPNsceHDh2K1157LS655JK47LLLEieD2WfCAbrjjjviyJEj8eCDD0Zvb2+sW7cuXnjhhbj88sunY75KutAx7O3tPe8zQUNDQ3HPPffE4cOHY/78+bF27dp4/vnn47bbbsv6I5TWnj17YuPGjWOPt23bFhERd911VzzxxBNJU8Hs5H5AAKRwLTgAUggQACkECIAUAgRACgECIIUAAZBCgABIIUAApBAgAFIIEAApBKjJ/PSnP40vfOELsWjRoliyZEl88YtfjIMHD2aPBTDlBKjJnDp1KrZt2xa/+tWv4mc/+1nUarX40pe+FKOjo9mjAUwpFyNtcu+8804sX748fvvb38a6deuyxwGYMlZATebgwYPx5S9/OT73uc9FW1tbrF69OiLivNszAFTBhO8HxPS6/fbbY9WqVfH9738/Vq5cGaOjo7Fu3boYGhrKHg1gSglQEzly5Ej87ne/i8cffzw2bNgQERE///nPk6cCmB4C1EQWL14cS5Ysia6urmhvb4/u7u649957s8eiyTUajTg9fDoiIi6+6OIoiiJ5Ihgf7wE1kVqtFj/+8Y/j1VdfjXXr1sU3v/nN+Pa3v509Fk3u9PDpWLB9QSzYvmAsRFAGVkBN5pZbbol9+/ad95wTFYEqsgICIIUAAZBCgABIIUAApBAgAFIIEAApBAgq5PU/HM8eAcZNgKBC/up//CJ7BBg3AQIghQABkEKAAEghQACkECAAUggQACkECIAUAgRACgECIIUAAZBCgABIIUAApBAgAFIIEAApBAiAFAIEQAoBAiCFAAGQQoAASCFAAKQQIABSCBAAKQQIgBQCBEAKAYKKaTQa2SPAuAgQVIz+UBYCBEAKAYKKsQCiLAQIKsZ7QJSFAAGQQoCgYqx/KAsBgoqxA0dZCBAAKQQIKqZhE46SECComJODI9kjwLgIEFTMeyOj2SPAuAgQVIwNOMpCgABIIUBQMa6EQFkIEAApBAgqxgKIshAgAFIIEAApBAgqxhYcZSFAAKQQIKgY14KjLAQIgBQCBBXjPSDKQoAASCFAUDEWQJSFAAGQQoCgYlyMlLIQIABSCBBUjPUPZSFAAKQQIKgYbwFRFgIEQAoBgsqxBKIcBAiAFAIEFeM9IMpCgKBi9IeyECAAUggQVIwtOMpCgABIIUBQMW7JTVkIEAApBAgqxntAlIUAAZBCgKBirIAoCwECIIUAQcU4C46yECAAUggQVMyxU8PZI8C4CBBUzKizECgJAQIghQBBxVj/UBYCBEAKAYKKaXgPiJIQIKgY+aEsBAiAFAIEVWMJREkIEAApBAgqxrXgKAsBgopxEhxlIUAApBAgqBgrIMpCgKBi9IeyECAAUggQVIxL8VAWAgRACgGCirH+oSwECCrGDhxlIUAApBAgqBxLIMpBgABIIUBQMd4DoiwECCpGfygLAQIghQBBxdiCoywECCrGDekoCwECIIUAQcXYgqMsBAiAFAIEFWMBRFkIEFTMmaGR7BFgXAQIKub0UD17BBgXAYKKGbUHR0kIEFRIS61wS25KQ4CgQmpFxKgAURICBBVSKwpbcJSGAEGFtNQKKyBKQ4CgQlqKwpUQKA0Bggqp1YoYtQdHSQgQVMjZkxCyp4DxESCokEYj4viZ4ewxYFwECCpk4bzWGKq7EgLlIEBQIfPntGSPAOMmQFAhhbPgKBEBggqpFW7HQHkIEFTI2WvBZU8B4yNAUCFFRFgDURYCBFXiPSBKRICgQooIAaI0BAgqpCiyJ4DxEyComIb3gCgJAYIKsQVHmQgQVEhRFNY/lIYAQYVYAVEmAgQVUhTeA6I8BAgqpAjX4qE8BAiAFAIEFWMBRFkIEFTI2ZMQJIhyECCokMJbQJSIAEGFFOFipJSHAEGFWAFRJgIEFVIU3gOiPAQIKqSIIk4MjmSPAeMiQFAhF7XWYmhkNHsMGBcBggpZMLfVPYEoDQECIIUAAZBCgABI0Zo9ADB1/jgwGL/+/ZnsMWBcrICgQvb8/lj2CDBuAgRACgECIIUAQYV88y/XxOKLL8oeA8ZFgKBCakURoy4FR0kIEFRIrYgYdTFSSkKAoEJqhfsBUR4CBBVSFIUVEKUhQFAhtuAoEwGCCikiYrguQJSDAEGFtNRqUR9txOBwPXsUuCABggpZtnBuRESces9dUWl+AgRN5MAfT8Sh/lOT/v6LWs/eje4P7w5O1UgwbVwNG5rA87/pjf/8v3499rjzc0viv37x8/F3h47E6mULYsO/WBq12kff6vTYqaGx/79i4byIiOg+ejr+9aX/bHqHhk+paDScMgPZrrj3+Ul/72gMRs/8v4qIiHf+y/G49r+99LFfe8ln5kz6dcri6PtBXji3NS5qrX3o+WY+Bkc/8MNExNTPOh3HoFYUsedbt0zqewUIgBTeAwIghQABkEKAAEghQACkECAAUvgcECRrNBpx4sSJ7DHgU1m4cGEUxUd/Vu3jCBAk6+/vj+XLl2ePAZ9KX19fLFu2bELfI0CQbM6csx8K7Onpiba2tuRpqmdgYCBWrVrl+E6Tc8f33N/jiRAgSHZu26Ktrc1/IKeR4zu9Jrr9FuEkBACSCBAAKQQIks2dOzceeOCBmDt3bvYoleT4Tq9Pc3xdjBSAFFZAAKQQIABSCBAAKQQIgBQCBIkeffTRWL16dcybNy/Wr18fL7308bfT5uNN5Dju2rUriqL40K833nhjBicuv927d8ftt98eK1eujKIo4rnnnpvw7yFAkOTpp5+OrVu3xv333x979+6NDRs2xK233hrd3d3Zo5XKZI/j/v37o7e3d+zXmjVrZmjiajh16lRcddVV8cgjj0z693AaNiS5/vrr45prronHHnts7Lkrr7wytmzZEtu3b0+crFwmehx37doVGzdujGPHjsWiRYtmcNLqKooinn322diyZcuEvs8KCBIMDQ3Fq6++Gps3bz7v+c2bN8fLL7+cNFX5fJrjePXVV0d7e3ts2rQpXnzxxekck48hQJCgv78/6vV6rFix4rznV6xYEW+//XbSVOUzmePY3t4eXV1dsWPHjnjmmWeio6MjNm3aFLt3756JkfkAV8OGRP/0CsKNRmNSVxWe7SZyHDs6OqKjo2PscWdnZ/T09MRDDz0UN91007TOyfmsgCDB0qVLo6Wl5UM/pff19X3op3k+3lQdxxtuuCEOHDgw1eNxAQIECebMmRPr16+PnTt3nvf8zp0748Ybb0yaqnym6jju3bs32tvbp3o8LsAWHCTZtm1b3HnnnXHttddGZ2dndHV1RXd3d9x9993Zo5XKhY7jfffdF4cPH44nn3wyIiIefvjhuOKKK2Lt2rUxNDQUTz31VOzYsSN27NiR+cconZMnT8abb7459vjQoUPx2muvxSWXXBKXXXbZuH4PAYIkd9xxRxw5ciQefPDB6O3tjXXr1sULL7wQl19+efZopXKh49jb23veZ4KGhobinnvuicOHD8f8+fNj7dq18fzzz8dtt92W9UcopT179sTGjRvHHm/bti0iIu6666544oknxvV7+BwQACm8BwRACgECIIUAAZBCgABIIUAApBAgAFIIEAApBAiAFAIEVNLNN98cW7duzR6DTyBAAKQQIABSCBBQWSMjI/G1r30tFi1aFEuWLIlvfetb4fKXzUOAgMr64Q9/GK2trfHLX/4yvvvd78Z3vvOd+MEPfpA9Fu9zNWygkm6++ebo6+uL119/fez23Pfee2/85Cc/iX379iVPR4QVEFBhN9xww1h8IiI6OzvjwIEDUa/XE6fiHAECIIUAAZX1yiuvfOjxmjVroqWlJWkiPkiAgMrq6emJbdu2xf79++NHP/pRfO9734tvfOMb2WPxvtbsAQCmy1e+8pU4c+ZMXHfdddHS0hJf//rX46tf/Wr2WLzPWXAApLAFB0AKAQIghQABkEKAAEghQACkECAAUggQACkECIAUAgRACgECIIUAAZDi/wNRW2lWKc5FuAAAAABJRU5ErkJggg==",
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"from vbi.plot import pairplot_numpy\n",
"\n",
"limits = [(0.0, 1.0), (0.0, 1.0)]\n",
"\n",
"fig, ax = pairplot_numpy(\n",
" samples,\n",
" points=theta_true_np.reshape(1, -1),\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"
]
},
{
"cell_type": "code",
"execution_count": 39,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Inferred dimensions: param_dim=2, feature_dim=4\n"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"Training: 0%| | 0/500 [00:00"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"print(\"Plotting posterior marginals for MAF ...\")\n",
"\n",
"param_labels = [\"a\", \"b\"]\n",
"\n",
"fig = corner.corner(\n",
" samples,\n",
" labels=param_labels,\n",
" truths=theta_true_np,\n",
" quantiles=[0.16, 0.5, 0.84],\n",
" show_titles=True,\n",
" title_kwargs={\"fontsize\": 12},\n",
" truth_color=\"red\",\n",
" bins=50, \n",
" fig=plt.figure(figsize=(6, 6)),\n",
")\n",
"plt.tight_layout()\n",
"plt.show()"
]
},
{
"cell_type": "code",
"execution_count": 42,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAaAAAAHNCAYAAABYXK/dAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAa7UlEQVR4nO3dbWyd5Zkn8Ov4HDsvhEwgIRlnCJCOslmazFIIA4QtDFHYSDDNKNVqhVSJ8mGlDtK2Ko34AEtXSEirrFSkohaBmvYDRUgtWgUYJFBXUZcoVKidkgnbFwoDWWYSgtmQFJo3yLGPz35IjmMf2yE29rmf+zy/nxSBj19y5VGcv+/ruZ77rjSbzWYAQIf1pC4AgHISQAAkIYAASEIAAZCEAAIgCQEEQBICCIAkBBAASdRSFwBE/Iee/5S6BPhMdg7/zyl/jhUQAEkIIACSEEAAJCGAAEhCAAGQRLIpuGazGScHT0ZExPze+VGpVFKVAkACyVZAJwdPxoJtC2LBtgUjQQRAeWjBAZCEAAIgCQEEQBICCIAkBBAASQggAJIoRAAd/WQwdQkAdFghAuiF//Ne6hIA6LBCBNDRU0OpSwCgwwoRQKeGhlOXAECHFSKAhhoCCKBsChFAg1ZAAKVTiACqWwEBlE4hAmhwWAABlE0xAsgQHEDpFCKAtOAAyqcQAWQIAaB8ihFAVkAApSOAAEiiEAFU14IDKJ1CBJAxbIDyKUYANZqpSwCgw5IFUGP4bOhowQGUT7IAGj14YDNSgPJJFkCjHz71ICpA+RRiBaQFB1A+yQJoaNTggSEEgPJJ14IbteoZGhZAAGVTjHtAWnAApVOMe0CN4Wg2rYIAyiRdALWterThAMqlECugCG04gLJJdw+obf83O2IDlEshpuAiPIwKUDYJW3Bj7/lowQGUSyEeRI3wMCpA2RRmCME9IIBySRZAp9rvAWnBAZRKYVZAhhAAyiXhPaC2FpwVEECpFGIvuAhDCABlU5wx7EYjUSUApFCcB1GHrIAAyqQwQwjGsAHKpTC7YRvDBigXKyAAkrAbNgBJJGzBjX27fWcEALpbgVpwpuAAyqQwY9hacADlUpgVkCk4gHIpTABZAQGUS2FacHbDBiiXdCugYS04gDIrzGakWnAA5VKYFpwxbIByKcyBdO4BAZRLYQ6kcw8IoFwK1IITQABlYggBgCTS3QMabjuSWwsOoFQK04Krm4IDKJUkAdRsNscNIbSfkApAd0sSQO3ttwhj2ABlkySAJrrfYwgBoFySBNBEYaMFB1AuaVZAEwSQFhxAuRSmBWcMG6BcErXgxg8h2IwUoFwKcw9ICw6gXArTgmsMN6MxwXg2AN2pMEMIEUaxAcokTQtukoEDbTiA8ijWCsgkHEBpFGIIobdaOfO6e0AAZZFoCGFs0PRWT5fhHhBAeRSiBdcKoFNacAClUYghhFpPqwUngADKohD3gObUtOAAyqYgLbjTKyD7wQGURyF2QmjdA/IcEEB5FGIz0t6RFpwxbICyKNYKSAsOoDQKMYTQ5zkggNIRQAAkkSSA2h84rZ2ZgvMgKkB5FGIFZCsegPIpRAD1tabgrIAASqMQU3B9ngMCKJ1CPAdUcxwDQOkUYiuePs8BAZROIe4B2YoHoHwKcQ+ot3amBWcFBFAahVgBeRAVoHwS3QOa+EhuLTiA8ihIC641hGAKDqAstOAASKIQAWQrHoDyKUYLzpHcAKVTiBVQay84QwgA5VGQFZAWHEDZFGIrHkdyA5RPITYjHdkJwWakAKXR8QBqDDejMTw2aIxhA5RPxwNoopDp7dGCAyibjgfQRJNupuAAyqfzK6AJVjlnD6QTQABlkWwF1Hr4NMKBdABllGAFdHoAoTV6HXG2BWcKDqA8Eq6Azv7WjmMAKJ/OB9DQ+BZcbVQLrtm0CgIog2Rj2H2jW3A9Z/9/aFgAAZRBsgAaM4RQ6xn3fgC6W7IW3OjQGR1GJuEAyiHZEEJtVAuu2lOJSmXs+wHobglacOPHsCuVyqj94NwDAiiDdC246tjfeiSAtOAASqEQU3AREb32gwMolYT3gCpjXrcdD0C5FGIKLmL0oXQCCKAMEj4H1BZAVkAApVKYe0Cm4ADKJd1ecLW2e0AjO2JbAQGUQYIhhDPPAfVM3II7pQUHUAqFuQd0tgUngADKoEBTcAIIoEwKsRt2RETfmbdNwQGUQ3FacFZAAKXS8QA6NVkLbuRYbmPYAGVQiN2wR7+tBQdQDp0PoCEtOAASbkY6+U4IAgigDNINIbTthNBrCg6gVNJtxTNJC855QADlUJgWXK8WHECpJGvB1ca14FpHchvDBiiDBFNwpwOmfQU0RwsOoFQK14ITQADlUJghhLMtOAEEUAbp7gG1b0aqBQdQKslacHPG7QV3OpBMwQGUQ3G24jEFB1AqhdmMtNWCO2UFBFAKHQ2gZrM50oLrneQ4BkMIAOXQ0QAaHHXWT2+PnRBgpnx06eH4x//8v+OjSw+nLgXOW4cD6Gy4jNsN2xQcTNuB6/bF0Uv/GAeu25e6FDhv6QKoNtkQggCCqajPOxWHPv9uREQc+vy7UZ93KnFFcH46GkCth1B7KhHVnra94M7sDedIbpiagS/8azQrp79vKtVm/N2zr8T/eu+1tEXBeah18jcbGUCojs+9s0dyNzpZEmTlkws/jvoFn4x57d1r90Wc+XmuGRGP/e5P8Tf98+Lon3848jF9J+bG3GPzOlgpfLqOBlBrCKG9/RYx+kRUKyCYzO/+46/io8vbBg2aMSaA/u+xofjrZ96N+Pt3Rz5k0b8siWt/fEuHqoTzk6QF1z6AEHE2lEzBweSW/9PK6BnsOZ00LZVJPzyiGdEz2BPL966c7dJgyqa9Amo2z3+lMjTcjE8GG3H4+Ombo+dqwQ0NN6Mx3Iyec31TwQyoVPL7S7b8N5fHwvcuit/c8UqcvPj4uX+EHI6Y/8cL4989vT4WHF7YsRrhfE07gFbe/+K0f9Pe2vhv/NFtub/8r9P/2nC+/uV//G3qEqZlweGFcf32W+P1za/G//urdyf9uGW/XxGff35dVIc62mmH89bxrXjm9vbEli/8xbjXL+irxl/9xZ91uhzIUnWwFhftv2RsK260ZsRF/7pE+FBo0/7buefbt573x1Z7KjG3txpzaj0jbY8T9RNjPqZSqcQ//Jd/Hx99PDjy2lTafFA2R/s/jMpwJZrV5tlBhDP/rQxX4ujyjyL2pK0RzmXaAbR4wZyZrCMiInp6KnHxBX0z/nWhG/3p0j9Gs9qMSqMSleFKXPrqX8a71+6LZk8zmtVm/OnSI6lLhHOyPocMNaqNOLHkWEREzPtwwcigwfJ/Wnl6QGHJ8Tix5Fg0qo2oNqqJq4WJCSDI0HBvIxYcWhgXDlwU//bFL4zc62kNKLxx29449ucfxXBNAFFcAggy1PtJX1z/g1ujMsFDQNXBWqx5/q+jGc0J3w9F0fEpOGBmfFq4CB+KTgABkESyFtz83vlx/P7jI/8PQLkkC6BKpRIX9F2Q6rcHIDEtOACSEEAAJCGAAEhCAAGQhAACIIlK05bTACRgBQRAEgIIgCQEEABJCCAAkhBAACQxrb3gms1mHDt2bKZrgY678MILo1JxbAGkMK0AOnz4cCxdunSma4GOO3ToUFxyySWpy4BSmlYA9fX1RUTEgQMHYuHChTNaEKcdPXo0VqxY4RrPktb1bf1dBjpvWgHUalksXLjQP46zzDWeXdpvkI4hBACSEEAAJDGtAJozZ048+OCDMWfOnJmuhzNc49nl+kJ6NiMFIAktOACSEEAAJCGAAEhCAAGQhAACIIlpBdBjjz0WK1eujLlz58a6devi5Zdfnum6ut5UruGuXbuiUqmM+/XGG290sOLusHv37ti8eXMsX748KpVKPPfcc6lLgtKacgA9/fTTcc8998QDDzwQe/fujZtuuiluu+222L9//2zU15Wmew3ffPPNGBgYGPm1atWqDlXcPU6cOBFXXXVVPProo6lLgdKb8nNA119/fVxzzTXx+OOPj7x25ZVXxpYtW2Lbtm0zXmA3muo13LVrV2zYsCE+/PDDWLRoUQcr7W6VSiWeffbZ2LJlS+pSoJSmtAKq1+uxZ8+e2LRp05jXN23aFK+88sqMFtatPss1vPrqq6O/vz82btwYL7300myWCTDrphRAhw8fjkajEcuWLRvz+rJly+L999+f0cK61XSuYX9/f2zfvj127NgRzzzzTKxevTo2btwYu3fv7kTJALPiMx3H0NJsNm1rP0VTuYarV6+O1atXj7y9fv36OHDgQDz88MNx8803z2qdALNlSiugJUuWRLVaHfeT+qFDh8b9RM/EZuoa3nDDDfHWW2/NdHkAHTOlAOrr64t169bFzp07x7y+c+fOuPHGG2e0sG41U9dw79690d/fP9PlAXTMlFtwW7dujTvvvDOuvfbaWL9+fWzfvj32798fd99992zU15U+7Rref//9cfDgwXjyyScjIuKRRx6JK664ItasWRP1ej2eeuqp2LFjR+zYsSPlHyNLx48fj7fffnvk7XfeeSdee+21uPjii+Oyyy5LWBmUz5QD6I477ogjR47EQw89FAMDA7F27dp48cUX4/LLL5+N+rrSp13DgYGBMc8E1ev1uPfee+PgwYMxb968WLNmTbzwwgtx++23p/ojZOvVV1+NDRs2jLy9devWiIi466674oknnkhUFZST84AASMJecAAkIYAASEIAAZCEAAIgCQEEQBICCIAkBBAASQggAJIQQAAkIYAASEIAFczPfvaz+OIXvxiLFi2KxYsXx5e+9KXYt29f6rIAZpwAKpgTJ07E1q1b49e//nX8/Oc/j56envjyl78cw8PDqUsDmFE2Iy24Dz74IJYuXRq//e1vY+3atanLAZgxVkAFs2/fvvjKV74Sn/vc52LhwoWxcuXKiIgxxzMAdIMpnwfE7Nq8eXOsWLEifvjDH8by5ctjeHg41q5dG/V6PXVpADNKABXIkSNH4g9/+EP84Ac/iJtuuikiIn7xi18krgpgdgigArnoooti8eLFsX379ujv74/9+/fHfffdl7osCq7ZbMbJwZMRETG/d35UKpXEFcH5cQ+oQHp6euKnP/1p7NmzJ9auXRvf+ta34jvf+U7qsii4k4MnY8G2BbFg24KRIIIcWAEVzK233hqvv/76mNcMKgLdyAoIgCQEEABJCCAAkhBAACQhgKCLGFghJwIIush/+4ffpS4BzpsAgi7ym3c/Sl0CnDcBBF2kPuTYDvIhgKCLDDbcAyIfAgi6SL1hBUQ+BBB0ESsgciKAoIsMWgGREQEEXUQAkRMBBF1ksNGM4WFtOPIggKDLDA5bBZEHAQRdxiACuRBA0GU8jEouBBB0GYMI5EIAQZexAiIXAgi6jN0QyIUAgi6jBUcuBBB0mcEhU3DkQQBBl6k3GqlLgPMigKDL1K2AyIQAgi7jHhC5EEDQZQQQuRBA0GU8B0QuBBB0Gc8BkQsBBF3GZqTkQgBBl9GCIxcCCLqMIQRyIYCgywggciGAoMuc0oIjEwIIuowVELkQQNBlBBC5EECQuWZz7Ni1KThyIYAgc+3P/XgOiFwIIMhce8vNTgjkQgBB5tpbblpw5EIAQeYGh8cGjiEEciGAIHPtgSOAyIUAgsxpwZErAQSZa596q5uCIxMCCDI3rgVnBUQmBBBkrr3l5h4QuRBAkDnPAZErAQSZax/DNoRALgQQZE4LjlwJIMjc0LgpOAFEHgQQZK597HpwyBg2eRBAkDktOHIlgCBz46bgDCGQCQEEmRsyhk2mBBBkrj1wtODIhQCCzLUH0HBz/KoIikgAQeYmmnpzLDc5EECQucFGY9xr7gORAwEEmZtotWMSjhwIIMjc6LDpq57+ljaIQA4EEGSuPmoz0t5qJSIEEHkQQJC50UMIvWdWQFpw5EAAQeZGj1y3VkCGEMiBAILM1ccEUOsekDFsik8AQebqo1pwc2qGEMiHAILMjQ6bmntAZEQAQeYGR4VNb809IPIhgCBzg8MTPAdkBUQGBBBkbnQLrrfnTAvOCogMCCDI3OgjufsMIZARAQSZGxw6uxnpyBj2BDtkQ9EIIMjc6Gd+WgF0ygqIDAggyNyYe0A1QwjkQwBB5uyGTa4EEGRucHhUC652+r8eRCUHAggyN9EYthUQORBAkLmJWnB1m5GSAQEEGWsMN2NUB25kCEELjhwIIMhYe6utZgiBjAggyNiptpVOnyO5yYgAgoy1B02fFhwZEUCQsfYA6q3ajJR8CCDIWPtKp2YMm4wIIMjYZC24QWPYZEAAQcbqbbte9zmSm4wIIMhY+70e94DIiQCCjI0bQqgZwyYfAggy1n7sQu+Z54C04MiBAIKMtR885zgGciKAIGPtK6CzAWQKjuITQJCx9qDpNQVHRgQQZKzeaIx5u2YKjowIIMjYYPtzQDX3gMiHAIKMta90PIhKTgQQZKw9aDwHRE4EEGRsst2wBxvNaDZNwlFsAggyNlkARRhEoPgEEGSsvQXXNyqAPAtE0QkgyFh9kueAIsY/pApFI4AgY+0tuGpPJao9Z/aD04Kj4AQQZGyiaTcbkpILAQQZmyhkbEhKLgQQZGyiNptjucmFAIKMTRQyNiQlFwIIMlYfaox7rbUCMoRA0QkgyNi5VkDuAVF0AggyNvEUnBYceRBAkLFTE03BOZKBTAggyNhEIdNXtSM2eRBAkLFzteAmWh1BkQggyNiED6J6DohMCCDImCk4ciaAIGPn2orHFBxFJ4AgY+feikcAUWwCCDJ2zt2wBRAFJ4AgYxMdOudBVHIhgCBjWnDkTABBpprN5oRTcGfPAzKGTbEJIMjUZAGjBUcuBBBkarIhg5EA0oKj4AQQZGqiAYSIUfeArIAoOAEEmWoNGVR7KmNe77UZKZkQQJCp1majrcBpcSIquRBAkKnWCqc19dZydiseU3AUmwCCTLWm4NpXQDYjJRcCCDJVH2nBjf027q0ZwyYPAggy1brH07rn09JnBUQmBBBkqhUw7SugvpopOPIggCBTrRZbrb0F50huMiGAIFODWnBkTgBBps6OYbdNwdVsRkoeBBBkauRB1J7JngOyAqLYBBBkauQ5oPYWnPOAyIQAgkxNNgVnN2xyIYAgU5NtxdPaGUELjqITQJCpVsCMm4LTgiMTAggy1Wqx1dp3wz6zIhpuRjSGTcJRXAIIMjU4dO7NSCO04Sg2AQSZqjcaETHBcQyjWnIGESgyAQSZao1htwdQbdQJqe4DUWQCCDI12XEMlUrFdjxkQQBBplrttfYHUSNGHcvtHhAFJoAgU4OTrIBOv+ZIBopPAEGmzu6EUBn3vpHdEIaMYVNcAggy1WrBzZlgBTTSgrMCosAEEGSqtbppP5AuwplA5EEAQaYm24x09GuDhhAoMAEEmZpsL7jRr52yAqLABBBk6twroDNTcFZAFJgAgkyNBFBt/PtGWnCO5abABBBkqnUkd1+1Ou59Z6fgGh2tCaZigp+dgE777y+8PuXPef/oJxER0dsz/jmg1hTcs3vfi9ffO/rZioNP8cDffn5anyeAoAB++PI70/7cP5vfN8FrvRERsfufP4jd//zBtL82nA8BBBn7+7/53LQ+b+XiC2LV0gXjXv/Wrf8mLl00zxQchVZpNpvuUkLGTtRPxIJtp0Po+P3H44K+CxJXBOfHEAIASQggAJIQQAAkIYAASEIAAZCEAAIgCWPYkLlmsxknB09GRMT83vlRqYzfGQGKSAABkIQWHABJCCAAkhBAACQhgABIQgABkIQAAiAJAQRAEgIIgCQEEABJCCAAkhBAACRRS10AlF2z2Yxjx46lLgM+kwsvvHDKG+EKIEjs8OHDsXTp0tRlwGdy6NChuOSSS6b0OQIIEuvr64uIiAMHDsTChQsTV9N9jh49GitWrHB9Z0nr+rb+Hk+FAILEWm2LhQsX+gdyFrm+s2s651AZQgAgCQEEQBICCBKbM2dOPPjggzFnzpzUpXQl13d2fZbr60huAJKwAgIgCQEEQBICCIAkBBAASQggSOixxx6LlStXxty5c2PdunXx8ssvpy4pS1O5jrt27YpKpTLu1xtvvNHBivO3e/fu2Lx5cyxfvjwqlUo899xzU/4aAggSefrpp+Oee+6JBx54IPbu3Rs33XRT3HbbbbF///7UpWVlutfxzTffjIGBgZFfq1at6lDF3eHEiRNx1VVXxaOPPjrtr2EMGxK5/vrr45prronHH3985LUrr7wytmzZEtu2bUtYWV6meh137doVGzZsiA8//DAWLVrUwUq7V6VSiWeffTa2bNkypc+zAoIE6vV67NmzJzZt2jTm9U2bNsUrr7ySqKr8fJbrePXVV0d/f39s3LgxXnrppdksk0kIIEjg8OHD0Wg0YtmyZWNeX7ZsWbz//vuJqsrPdK5jf39/bN++PXbs2BHPPPNMrF69OjZu3Bi7d+/uRMmMYjdsSKh9B+FmszmtXYXLbirXcfXq1bF69eqRt9evXx8HDhyIhx9+OG6++eZZrZOxrIAggSVLlkS1Wh33U/qhQ4fG/TTP5GbqOt5www3x1ltvzXR5fAoBBAn09fXFunXrYufOnWNe37lzZ9x4442JqsrPTF3HvXv3Rn9//0yXx6fQgoNEtm7dGnfeeWdce+21sX79+ti+fXvs378/7r777tSlZeXTruP9998fBw8ejCeffDIiIh555JG44oorYs2aNVGv1+Opp56KHTt2xI4dO1L+MbJz/PjxePvtt0fefuedd+K1116Liy++OC677LLz+hoCCBK544474siRI/HQQw/FwMBArF27Nl588cW4/PLLU5eWlU+7jgMDA2OeCarX63HvvffGwYMHY968ebFmzZp44YUX4vbbb0/1R8jSq6++Ghs2bBh5e+vWrRERcdddd8UTTzxxXl/Dc0AAJOEeEABJCCAAkhBAACQhgABIQgABkIQAAiAJAQRAEgIIgCQEENCVbrnllrjnnntSl8E5CCAAkhBAACQhgICuNTQ0FF//+tdj0aJFsXjx4vj2t78dtr8sDgEEdK0f//jHUavV4le/+lV873vfi+9+97vxox/9KHVZnGE3bKAr3XLLLXHo0KH4/e9/P3I893333RfPP/98vP7664mrI8IKCOhiN9xww0j4RESsX78+3nrrrWg0GgmrokUAAZCEAAK61i9/+ctxb69atSqq1WqiihhNAAFd68CBA7F169Z488034yc/+Ul8//vfj29+85upy+KMWuoCAGbLV7/61fj444/juuuui2q1Gt/4xjfia1/7WuqyOMMUHABJaMEBkIQAAiAJAQRAEgIIgCQEEABJCCAAkhBAACQhgABIQgABkIQAAiAJAQRAEv8fiONpJfA+PtYAAAAASUVORK5CYII=",
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"from vbi.plot import pairplot_numpy\n",
"\n",
"limits = [(0.0, 1.0), (0.0, 1.0)]\n",
"\n",
"fig, ax = pairplot_numpy(\n",
" samples,\n",
" points=theta_true_np.reshape(1, -1),\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"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"language_info": {
"name": "python"
}
},
"nbformat": 4,
"nbformat_minor": 2
}