notebooks/research/benchmarking_di_approximations_pymc.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "from math import tanh\n",
    "import numpy as np\n",
    "import scipy.stats as st \n",
    "import matplotlib.pyplot as plt \n",
    "import pymc as pm\n",
    "import arviz as az\n",
    "\n",
    "def a(S_r, b, c, d):\n",
    "    return S_r - b * np.tanh(c * d)\n",
    "\n",
    "def b(S_r, F_bar):\n",
    "    return (S_r - max(S_r - F_bar, 0)) / F_bar\n",
    "\n",
    "def B(g, a, b, c, d):\n",
    "    return a + b * np.tanh(-c * (g - d))\n",
    "\n",
    "def B_equiv(g, S_r, b):\n",
    "    _c = 1.0\n",
    "    _d = 1.0\n",
    "    _a = a(S_r, b, _c, _d)\n",
    "    return B(g, _a, b, _c, _d)\n",
    "\n",
    "def B_approx(g, S_r, b):\n",
    "    return S_r - b * (1523/2000 - (3/4) * (1 - g))\n",
    "\n",
    "def F(credit_supply, space_pledged, history_size, min_replication_factor):\n",
    "    return credit_supply / max((space_pledged/min_replication_factor - history_size), 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"media/beta.png\" height=200>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Multiprocess sampling (4 chains in 4 jobs)\n",
      "CompoundStep\n",
      ">NUTS: [g, S_r, C, S, H]\n",
      ">Metropolis: [Rf]\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "\n",
       "    <div>\n",
       "        <style>\n",
       "            /* Turns off some styling */\n",
       "            progress {\n",
       "                /* gets rid of default border in Firefox and Opera. */\n",
       "                border: none;\n",
       "                /* Needs to be in here for Safari polyfill so background images work as expected. */\n",
       "                background-size: auto;\n",
       "            }\n",
       "            .progress-bar-interrupted, .progress-bar-interrupted::-webkit-progress-bar {\n",
       "                background: #F44336;\n",
       "            }\n",
       "        </style>\n",
       "      <progress value='44000' class='' max='44000' style='width:300px; height:20px; vertical-align: middle;'></progress>\n",
       "      100.00% [44000/44000 00:06<00:00 Sampling 4 chains, 0 divergences]\n",
       "    </div>\n",
       "    "
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Sampling 4 chains for 1_000 tune and 10_000 draw iterations (4_000 + 40_000 draws total) took 6 seconds.\n"
     ]
    }
   ],
   "source": [
    "# N_YEARS = 4\n",
    "# BLOCK_TIME_IN_SECONDS = 10\n",
    "# BLOCKS_PER_DAY = 24 * 60 * 60 / BLOCK_TIME_IN_SECONDS\n",
    "# BLOCKS_PER_YEAR = BLOCKS_PER_DAY * 365.25\n",
    "# BLOCKS_PER_N_YEARS = BLOCKS_PER_YEAR * N_YEARS\n",
    "# MAX_DISBURSAL_DURING_N_YEARS = 500_000_000 # 500 Million\n",
    "# MAX_AVG_BLOCK_REWARD_DURING_N_YEARS = MAX_DISBURSAL_DURING_N_YEARS / BLOCKS_PER_N_YEARS\n",
    "# MIN_DISBURSAL_DURING_N_YEARS = MAX_DISBURSAL_DURING_N_YEARS / 10\n",
    "# MIN_AVG_BLOCK_REWARD_DURING_N_YEARS = MIN_DISBURSAL_DURING_N_YEARS / BLOCKS_PER_N_YEARS\n",
    "\n",
    "model = pm.Model()\n",
    "\n",
    "KB = 2 ** (1 * 10)\n",
    "MB = 2 ** (2 * 10)\n",
    "GB = 2 ** (3 * 10)\n",
    "TB = 2 ** (4 * 10)\n",
    "PB = 2 ** (5 * 10)\n",
    "EB = 2 ** (6 * 10)\n",
    "ZB = 2 ** (7 * 10)\n",
    "world_datasphere_size = 15 * ZB\n",
    "max_credit_supply = 1_000_000_000\n",
    "max_block_size = 1000 # HACK\n",
    "with model:\n",
    "\n",
    "    utilization_ratio = pm.Beta('g', alpha=1.1, beta=10)\n",
    "    reference_subsidy = pm.Triangular('S_r', c=15.0, lower=0.0, upper=1_000.0)\n",
    "\n",
    "    credit_supply_in_credits = pm.Uniform('C', 0.0, max_credit_supply)\n",
    "    space_pledged_in_bytes = pm.Uniform('S', 0.0, world_datasphere_size)\n",
    "    history_size_in_bytes = pm.Uniform('H', 0.0, world_datasphere_size)\n",
    "    min_replication_factor = pm.DiscreteUniform('Rf', 0, 1000)\n",
    "\n",
    "    _F = pm.math.switch(space_pledged_in_bytes/min_replication_factor - history_size_in_bytes > 0, space_pledged_in_bytes/min_replication_factor - history_size_in_bytes, 1.0)\n",
    "    F = pm.Deterministic('F', _F)\n",
    "\n",
    "    Fbar = pm.Deterministic('Fbar', F * max_block_size)\n",
    "\n",
    "    _b = pm.math.switch(reference_subsidy - Fbar > 0, 1.0, reference_subsidy / Fbar)\n",
    "\n",
    "    linear_sensitivity = pm.Deterministic('b', _b)\n",
    "\n",
    "\n",
    "    reward_approx = pm.Deterministic('B_approx', B_approx(utilization_ratio, reference_subsidy, linear_sensitivity))\n",
    "    reward_true = pm.Deterministic('B_true', B_equiv(utilization_ratio, reference_subsidy, linear_sensitivity))\n",
    "    error = pm.Deterministic('B_error', reward_true - reward_approx)\n",
    "\n",
    "    idata = pm.sample(draws=10_000)\n",
    "\n",
    "\n",
    "\n",
    "# credit_supply_samples = None\n",
    "# space_pledged_sampels = None\n",
    "# history_size_samples = None\n",
    "# min_replication_factor_samples = None\n",
    "\n",
    "# Fbar_samples = F(credit_supply, space_pledged, history_size, min_replication_factor)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [],
   "source": [
    "ds = idata.posterior[['B_error', 'g', 'S_r', 'b']]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = ds.to_dataframe().reset_index().drop(columns=['chain', 'draw'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [],
   "source": [
    "import seaborn as sns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [
    {
     "ename": "AttributeError",
     "evalue": "QuadMesh.set() got an unexpected keyword argument 'index'",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mAttributeError\u001b[0m                            Traceback (most recent call last)",
      "Cell \u001b[0;32mIn[49], line 1\u001b[0m\n\u001b[0;32m----> 1\u001b[0m sns\u001b[38;5;241m.\u001b[39mheatmap(df, index\u001b[38;5;241m=\u001b[39m\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mg\u001b[39m\u001b[38;5;124m'\u001b[39m, columns\u001b[38;5;241m=\u001b[39m\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mS_r\u001b[39m\u001b[38;5;124m'\u001b[39m, values\u001b[38;5;241m=\u001b[39m\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mB_error\u001b[39m\u001b[38;5;124m'\u001b[39m)\n",
      "File \u001b[0;32m~/anaconda3/lib/python3.11/site-packages/seaborn/matrix.py:459\u001b[0m, in \u001b[0;36mheatmap\u001b[0;34m(data, vmin, vmax, cmap, center, robust, annot, fmt, annot_kws, linewidths, linecolor, cbar, cbar_kws, cbar_ax, square, xticklabels, yticklabels, mask, ax, **kwargs)\u001b[0m\n\u001b[1;32m    457\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m square:\n\u001b[1;32m    458\u001b[0m     ax\u001b[38;5;241m.\u001b[39mset_aspect(\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mequal\u001b[39m\u001b[38;5;124m\"\u001b[39m)\n\u001b[0;32m--> 459\u001b[0m plotter\u001b[38;5;241m.\u001b[39mplot(ax, cbar_ax, kwargs)\n\u001b[1;32m    460\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m ax\n",
      "File \u001b[0;32m~/anaconda3/lib/python3.11/site-packages/seaborn/matrix.py:306\u001b[0m, in \u001b[0;36m_HeatMapper.plot\u001b[0;34m(self, ax, cax, kws)\u001b[0m\n\u001b[1;32m    303\u001b[0m     kws\u001b[38;5;241m.\u001b[39msetdefault(\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mvmax\u001b[39m\u001b[38;5;124m\"\u001b[39m, \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mvmax)\n\u001b[1;32m    305\u001b[0m \u001b[38;5;66;03m# Draw the heatmap\u001b[39;00m\n\u001b[0;32m--> 306\u001b[0m mesh \u001b[38;5;241m=\u001b[39m ax\u001b[38;5;241m.\u001b[39mpcolormesh(\u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mplot_data, cmap\u001b[38;5;241m=\u001b[39m\u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mcmap, \u001b[38;5;241m*\u001b[39m\u001b[38;5;241m*\u001b[39mkws)\n\u001b[1;32m    308\u001b[0m \u001b[38;5;66;03m# Set the axis limits\u001b[39;00m\n\u001b[1;32m    309\u001b[0m ax\u001b[38;5;241m.\u001b[39mset(xlim\u001b[38;5;241m=\u001b[39m(\u001b[38;5;241m0\u001b[39m, \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mdata\u001b[38;5;241m.\u001b[39mshape[\u001b[38;5;241m1\u001b[39m]), ylim\u001b[38;5;241m=\u001b[39m(\u001b[38;5;241m0\u001b[39m, \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mdata\u001b[38;5;241m.\u001b[39mshape[\u001b[38;5;241m0\u001b[39m]))\n",
      "File \u001b[0;32m~/anaconda3/lib/python3.11/site-packages/matplotlib/__init__.py:1446\u001b[0m, in \u001b[0;36m_preprocess_data.<locals>.inner\u001b[0;34m(ax, data, *args, **kwargs)\u001b[0m\n\u001b[1;32m   1443\u001b[0m \u001b[38;5;129m@functools\u001b[39m\u001b[38;5;241m.\u001b[39mwraps(func)\n\u001b[1;32m   1444\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21minner\u001b[39m(ax, \u001b[38;5;241m*\u001b[39margs, data\u001b[38;5;241m=\u001b[39m\u001b[38;5;28;01mNone\u001b[39;00m, \u001b[38;5;241m*\u001b[39m\u001b[38;5;241m*\u001b[39mkwargs):\n\u001b[1;32m   1445\u001b[0m     \u001b[38;5;28;01mif\u001b[39;00m data \u001b[38;5;129;01mis\u001b[39;00m \u001b[38;5;28;01mNone\u001b[39;00m:\n\u001b[0;32m-> 1446\u001b[0m         \u001b[38;5;28;01mreturn\u001b[39;00m func(ax, \u001b[38;5;241m*\u001b[39m\u001b[38;5;28mmap\u001b[39m(sanitize_sequence, args), \u001b[38;5;241m*\u001b[39m\u001b[38;5;241m*\u001b[39mkwargs)\n\u001b[1;32m   1448\u001b[0m     bound \u001b[38;5;241m=\u001b[39m new_sig\u001b[38;5;241m.\u001b[39mbind(ax, \u001b[38;5;241m*\u001b[39margs, \u001b[38;5;241m*\u001b[39m\u001b[38;5;241m*\u001b[39mkwargs)\n\u001b[1;32m   1449\u001b[0m     auto_label \u001b[38;5;241m=\u001b[39m (bound\u001b[38;5;241m.\u001b[39marguments\u001b[38;5;241m.\u001b[39mget(label_namer)\n\u001b[1;32m   1450\u001b[0m                   \u001b[38;5;129;01mor\u001b[39;00m bound\u001b[38;5;241m.\u001b[39mkwargs\u001b[38;5;241m.\u001b[39mget(label_namer))\n",
      "File \u001b[0;32m~/anaconda3/lib/python3.11/site-packages/matplotlib/axes/_axes.py:6227\u001b[0m, in \u001b[0;36mAxes.pcolormesh\u001b[0;34m(self, alpha, norm, cmap, vmin, vmax, shading, antialiased, *args, **kwargs)\u001b[0m\n\u001b[1;32m   6223\u001b[0m     C \u001b[38;5;241m=\u001b[39m C\u001b[38;5;241m.\u001b[39mravel()\n\u001b[1;32m   6225\u001b[0m kwargs\u001b[38;5;241m.\u001b[39msetdefault(\u001b[38;5;124m'\u001b[39m\u001b[38;5;124msnap\u001b[39m\u001b[38;5;124m'\u001b[39m, mpl\u001b[38;5;241m.\u001b[39mrcParams[\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mpcolormesh.snap\u001b[39m\u001b[38;5;124m'\u001b[39m])\n\u001b[0;32m-> 6227\u001b[0m collection \u001b[38;5;241m=\u001b[39m mcoll\u001b[38;5;241m.\u001b[39mQuadMesh(\n\u001b[1;32m   6228\u001b[0m     coords, antialiased\u001b[38;5;241m=\u001b[39mantialiased, shading\u001b[38;5;241m=\u001b[39mshading,\n\u001b[1;32m   6229\u001b[0m     array\u001b[38;5;241m=\u001b[39mC, cmap\u001b[38;5;241m=\u001b[39mcmap, norm\u001b[38;5;241m=\u001b[39mnorm, alpha\u001b[38;5;241m=\u001b[39malpha, \u001b[38;5;241m*\u001b[39m\u001b[38;5;241m*\u001b[39mkwargs)\n\u001b[1;32m   6230\u001b[0m collection\u001b[38;5;241m.\u001b[39m_scale_norm(norm, vmin, vmax)\n\u001b[1;32m   6232\u001b[0m coords \u001b[38;5;241m=\u001b[39m coords\u001b[38;5;241m.\u001b[39mreshape(\u001b[38;5;241m-\u001b[39m\u001b[38;5;241m1\u001b[39m, \u001b[38;5;241m2\u001b[39m)  \u001b[38;5;66;03m# flatten the grid structure; keep x, y\u001b[39;00m\n",
      "File \u001b[0;32m~/anaconda3/lib/python3.11/site-packages/matplotlib/collections.py:1939\u001b[0m, in \u001b[0;36mQuadMesh.__init__\u001b[0;34m(self, coordinates, antialiased, shading, **kwargs)\u001b[0m\n\u001b[1;32m   1936\u001b[0m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_bbox\u001b[38;5;241m.\u001b[39mupdate_from_data_xy(\u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_coordinates\u001b[38;5;241m.\u001b[39mreshape(\u001b[38;5;241m-\u001b[39m\u001b[38;5;241m1\u001b[39m, \u001b[38;5;241m2\u001b[39m))\n\u001b[1;32m   1937\u001b[0m \u001b[38;5;66;03m# super init delayed after own init because array kwarg requires\u001b[39;00m\n\u001b[1;32m   1938\u001b[0m \u001b[38;5;66;03m# self._coordinates and self._shading\u001b[39;00m\n\u001b[0;32m-> 1939\u001b[0m \u001b[38;5;28msuper\u001b[39m()\u001b[38;5;241m.\u001b[39m\u001b[38;5;21m__init__\u001b[39m(\u001b[38;5;241m*\u001b[39m\u001b[38;5;241m*\u001b[39mkwargs)\n\u001b[1;32m   1940\u001b[0m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mset_mouseover(\u001b[38;5;28;01mFalse\u001b[39;00m)\n",
      "File \u001b[0;32m~/anaconda3/lib/python3.11/site-packages/matplotlib/_api/deprecation.py:454\u001b[0m, in \u001b[0;36mmake_keyword_only.<locals>.wrapper\u001b[0;34m(*args, **kwargs)\u001b[0m\n\u001b[1;32m    448\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;28mlen\u001b[39m(args) \u001b[38;5;241m>\u001b[39m name_idx:\n\u001b[1;32m    449\u001b[0m     warn_deprecated(\n\u001b[1;32m    450\u001b[0m         since, message\u001b[38;5;241m=\u001b[39m\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mPassing the \u001b[39m\u001b[38;5;132;01m%(name)s\u001b[39;00m\u001b[38;5;124m \u001b[39m\u001b[38;5;132;01m%(obj_type)s\u001b[39;00m\u001b[38;5;124m \u001b[39m\u001b[38;5;124m\"\u001b[39m\n\u001b[1;32m    451\u001b[0m         \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mpositionally is deprecated since Matplotlib \u001b[39m\u001b[38;5;132;01m%(since)s\u001b[39;00m\u001b[38;5;124m; the \u001b[39m\u001b[38;5;124m\"\u001b[39m\n\u001b[1;32m    452\u001b[0m         \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mparameter will become keyword-only \u001b[39m\u001b[38;5;132;01m%(removal)s\u001b[39;00m\u001b[38;5;124m.\u001b[39m\u001b[38;5;124m\"\u001b[39m,\n\u001b[1;32m    453\u001b[0m         name\u001b[38;5;241m=\u001b[39mname, obj_type\u001b[38;5;241m=\u001b[39m\u001b[38;5;124mf\u001b[39m\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mparameter of \u001b[39m\u001b[38;5;132;01m{\u001b[39;00mfunc\u001b[38;5;241m.\u001b[39m\u001b[38;5;18m__name__\u001b[39m\u001b[38;5;132;01m}\u001b[39;00m\u001b[38;5;124m()\u001b[39m\u001b[38;5;124m\"\u001b[39m)\n\u001b[0;32m--> 454\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m func(\u001b[38;5;241m*\u001b[39margs, \u001b[38;5;241m*\u001b[39m\u001b[38;5;241m*\u001b[39mkwargs)\n",
      "File \u001b[0;32m~/anaconda3/lib/python3.11/site-packages/matplotlib/collections.py:201\u001b[0m, in \u001b[0;36mCollection.__init__\u001b[0;34m(self, edgecolors, facecolors, linewidths, linestyles, capstyle, joinstyle, antialiaseds, offsets, offset_transform, norm, cmap, pickradius, hatch, urls, zorder, **kwargs)\u001b[0m\n\u001b[1;32m    198\u001b[0m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_offset_transform \u001b[38;5;241m=\u001b[39m offset_transform\n\u001b[1;32m    200\u001b[0m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_path_effects \u001b[38;5;241m=\u001b[39m \u001b[38;5;28;01mNone\u001b[39;00m\n\u001b[0;32m--> 201\u001b[0m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_internal_update(kwargs)\n\u001b[1;32m    202\u001b[0m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_paths \u001b[38;5;241m=\u001b[39m \u001b[38;5;28;01mNone\u001b[39;00m\n",
      "File \u001b[0;32m~/anaconda3/lib/python3.11/site-packages/matplotlib/artist.py:1223\u001b[0m, in \u001b[0;36mArtist._internal_update\u001b[0;34m(self, kwargs)\u001b[0m\n\u001b[1;32m   1216\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21m_internal_update\u001b[39m(\u001b[38;5;28mself\u001b[39m, kwargs):\n\u001b[1;32m   1217\u001b[0m \u001b[38;5;250m    \u001b[39m\u001b[38;5;124;03m\"\"\"\u001b[39;00m\n\u001b[1;32m   1218\u001b[0m \u001b[38;5;124;03m    Update artist properties without prenormalizing them, but generating\u001b[39;00m\n\u001b[1;32m   1219\u001b[0m \u001b[38;5;124;03m    errors as if calling `set`.\u001b[39;00m\n\u001b[1;32m   1220\u001b[0m \n\u001b[1;32m   1221\u001b[0m \u001b[38;5;124;03m    The lack of prenormalization is to maintain backcompatibility.\u001b[39;00m\n\u001b[1;32m   1222\u001b[0m \u001b[38;5;124;03m    \"\"\"\u001b[39;00m\n\u001b[0;32m-> 1223\u001b[0m     \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_update_props(\n\u001b[1;32m   1224\u001b[0m         kwargs, \u001b[38;5;124m\"\u001b[39m\u001b[38;5;132;01m{cls.__name__}\u001b[39;00m\u001b[38;5;124m.set() got an unexpected keyword argument \u001b[39m\u001b[38;5;124m\"\u001b[39m\n\u001b[1;32m   1225\u001b[0m         \u001b[38;5;124m\"\u001b[39m\u001b[38;5;132;01m{prop_name!r}\u001b[39;00m\u001b[38;5;124m\"\u001b[39m)\n",
      "File \u001b[0;32m~/anaconda3/lib/python3.11/site-packages/matplotlib/artist.py:1197\u001b[0m, in \u001b[0;36mArtist._update_props\u001b[0;34m(self, props, errfmt)\u001b[0m\n\u001b[1;32m   1195\u001b[0m             func \u001b[38;5;241m=\u001b[39m \u001b[38;5;28mgetattr\u001b[39m(\u001b[38;5;28mself\u001b[39m, \u001b[38;5;124mf\u001b[39m\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mset_\u001b[39m\u001b[38;5;132;01m{\u001b[39;00mk\u001b[38;5;132;01m}\u001b[39;00m\u001b[38;5;124m\"\u001b[39m, \u001b[38;5;28;01mNone\u001b[39;00m)\n\u001b[1;32m   1196\u001b[0m             \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;129;01mnot\u001b[39;00m \u001b[38;5;28mcallable\u001b[39m(func):\n\u001b[0;32m-> 1197\u001b[0m                 \u001b[38;5;28;01mraise\u001b[39;00m \u001b[38;5;167;01mAttributeError\u001b[39;00m(\n\u001b[1;32m   1198\u001b[0m                     errfmt\u001b[38;5;241m.\u001b[39mformat(\u001b[38;5;28mcls\u001b[39m\u001b[38;5;241m=\u001b[39m\u001b[38;5;28mtype\u001b[39m(\u001b[38;5;28mself\u001b[39m), prop_name\u001b[38;5;241m=\u001b[39mk))\n\u001b[1;32m   1199\u001b[0m             ret\u001b[38;5;241m.\u001b[39mappend(func(v))\n\u001b[1;32m   1200\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m ret:\n",
      "\u001b[0;31mAttributeError\u001b[0m: QuadMesh.set() got an unexpected keyword argument 'index'"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAi4AAAGiCAYAAADA0E3hAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAA9hAAAPYQGoP6dpAAAbnElEQVR4nO3dbWzV5f348U9paasu1AhaQTJ+4EBRoo42VkqY0WkNGg0PjBgXEafJmukQOp1UFm+ISSObbt6Bd6AxQUe8nQ86pQ8mViHb6IpxQqIRtKKtpCwq3qwIfP8PDP3vrEU5hRau9vVKzoNeXN+e6+Rad95+v+emIMuyLAAAEjDsUC8AAGB/CRcAIBnCBQBIhnABAJIhXACAZAgXACAZwgUASIZwAQCSIVwAgGQIFwAgGXmHy2uvvRYXX3xxjBkzJgoKCuLFF1/83mPWrFkTFRUVUVpaGhMmTIiHHnqoL2sFAIa4vMPlyy+/jNNPPz0eeOCB/Zq/ZcuWuPDCC2PGjBnR2toat9xyS8ybNy+ee+65vBcLAAxtBQfyJYsFBQXxwgsvxKxZs/Y55+abb46XXnopNm3a1D1WW1sbb775Zqxbt66vdw0ADEH9/hqXdevWRU1NTc7YBRdcEOvXr49vvvmm12O6urri888/z7l1dXX191IBgMNcv4dLR0dHlJeX54yVl5fHrl27orOzs9djGhoaoqysLOfW0NDQ30sFAA5zRQNxJwUFBTk/77069b/je9XX10ddXV3OWElJSf8sDgBIRr+Hy/HHHx8dHR05Y9u2bYuioqIYOXJkr8eUlJQIFQCgh36/VDRt2rRoamrKGVu9enVUVlbG8OHD+/vuAYBBJO9w+eKLL2LDhg2xYcOGiPj27c4bNmyItra2iPj2Ms+cOXO659fW1sYHH3wQdXV1sWnTplixYkUsX748brzxxoPzCACAISPvt0O/+uqrcc455/QYv+qqq+KJJ56IuXPnxvvvvx+vvvpq97+tWbMmFixYEG+//XaMGTMmbr755qitrT3gxQMAQ8sBfY4LAMBA8l1FAEAyhAsAkAzhAgAkQ7gAAMkQLgBAMoQLAJAM4QIAJEO4AADJEC4AQDKECwCQDOECACRDuAAAyRAuAEAyhAsAkAzhAgAkQ7gAAMkQLgBAMoQLAJAM4QIAJEO4AADJEC4AQDKECwCQDOECACRDuAAAyRAuAEAyhAsAkAzhAgAkQ7gAAMkQLgBAMoQLAJAM4QIAJEO4AADJEC4AQDKECwCQDOECACRDuAAAyRAuAEAyhAsAkAzhAgAkQ7gAAMkQLgBAMoQLAJAM4QIAJEO4AADJEC4AQDKECwCQDOECACRDuAAAyRAuAEAyhAsAkAzhAgAkQ7gAAMkQLgBAMoQLAJAM4QIAJEO4AADJEC4AQDKECwCQDOECACRDuAAAyRAuAEAyhAsAkAzhAgAko0/hsnTp0hg/fnyUlpZGRUVFNDc3f+f8lStXxumnnx5HHnlkjB49Oq6++urYvn17nxYMAAxdeYfLqlWrYv78+bFo0aJobW2NGTNmxMyZM6Otra3X+a+//nrMmTMnrrnmmnj77bfjmWeeiX/84x9x7bXXHvDiAYChpSDLsiyfA6qqqmLq1KmxbNmy7rHJkyfHrFmzoqGhocf83//+97Fs2bJ47733usfuv//+WLJkSXz44YcHsHQAYKjJ64zLzp07o6WlJWpqanLGa2pqYu3atb0eU11dHVu3bo3GxsbIsiw++eSTePbZZ+Oiiy7a5/10dXXF559/nnPr6urKZ6kAwCCUV7h0dnbG7t27o7y8PGe8vLw8Ojo6ej2muro6Vq5cGbNnz47i4uI4/vjj4+ijj477779/n/fT0NAQZWVlObfezuYAAENLn16cW1BQkPNzlmU9xvbauHFjzJs3L2699dZoaWmJl19+ObZs2RK1tbX7/P319fXx2Wef5dzq6+v7slQAYBApymfyqFGjorCwsMfZlW3btvU4C7NXQ0NDTJ8+PW666aaIiDjttNPiqKOOihkzZsSdd94Zo0eP7nFMSUlJlJSU5LM0AGAIyOuMS3FxcVRUVERTU1POeFNTU1RXV/d6zFdffRXDhuXeTWFhYUR8e6YGAGB/5X2pqK6uLh577LFYsWJFbNq0KRYsWBBtbW3dl37q6+tjzpw53fMvvvjieP7552PZsmWxefPmeOONN2LevHlx5plnxpgxYw7eIwEABr28LhVFRMyePTu2b98eixcvjvb29pgyZUo0NjbGuHHjIiKivb095zNd5s6dGzt27IgHHnggfv3rX8fRRx8d5557btx1110H71EAAENC3p/jAgBwqPiuIgAgGcIFAEiGcAEAkiFcAIBkCBcAIBnCBQBIhnABAJIhXACAZAgXACAZwgUASIZwAQCSIVwAgGQIFwAgGcIFAEiGcAEAkiFcAIBkCBcAIBnCBQBIhnABAJIhXACAZAgXACAZwgUASIZwAQCSIVwAgGQIFwAgGcIFAEiGcAEAkiFcAIBkCBcAIBnCBQBIhnABAJIhXACAZAgXACAZwgUASIZwAQCSIVwAgGQIFwAgGcIFAEiGcAEAkiFcAIBkCBcAIBnCBQBIhnABAJIhXACAZAgXACAZwgUASIZwAQCSIVwAgGQIFwAgGcIFAEiGcAEAkiFcAIBkCBcAIBnCBQBIhnABAJIhXACAZAgXACAZwgUASIZwAQCSIVwAgGQIFwAgGcIFAEiGcAEAktGncFm6dGmMHz8+SktLo6KiIpqbm79zfldXVyxatCjGjRsXJSUlceKJJ8aKFSv6tGAAYOgqyveAVatWxfz582Pp0qUxffr0ePjhh2PmzJmxcePG+OEPf9jrMZdddll88sknsXz58vjRj34U27Zti127dh3w4gGAoaUgy7IsnwOqqqpi6tSpsWzZsu6xyZMnx6xZs6KhoaHH/Jdffjkuv/zy2Lx5cxxzzDEHvmIAYMjK61LRzp07o6WlJWpqanLGa2pqYu3atb0e89JLL0VlZWUsWbIkTjjhhJg0aVLceOON8fXXX+/zfrq6uuLzzz/PuXV1deWzVABgEMorXDo7O2P37t1RXl6eM15eXh4dHR29HrN58+Z4/fXX41//+le88MIL8cc//jGeffbZuO666/Z5Pw0NDVFWVpZz6+1sDgAwtOT9GpeIiIKCgpyfsyzrMbbXnj17oqCgIFauXBllZWUREXHPPffEpZdeGg8++GAcccQRPY6pr6+Purq6nLGSkpK+LBUAGETyCpdRo0ZFYWFhj7Mr27Zt63EWZq/Ro0fHCSec0B0tEd++JibLsti6dWtMnDixxzElJSVCBQDoIa9LRcXFxVFRURFNTU05401NTVFdXd3rMdOnT4+PP/44vvjii+6xd955J4YNGxZjx47tw5IBgKEq789xqauri8ceeyxWrFgRmzZtigULFkRbW1vU1tZGxLeXeebMmdM9/4orroiRI0fG1VdfHRs3bozXXnstbrrppvj5z3/e62UiAIB9yfs1LrNnz47t27fH4sWLo729PaZMmRKNjY0xbty4iIhob2+Ptra27vk/+MEPoqmpKX71q19FZWVljBw5Mi677LK48847D96jAACGhLw/xwUA4FDxXUUAQDKECwCQDOECACRDuAAAyRAuAEAyhAsAkAzhAgAkQ7gAAMkQLgBAMoQLAJAM4QIAJEO4AADJEC4AQDKECwCQDOECACRDuAAAyRAuAEAyhAsAkAzhAgAkQ7gAAMkQLgBAMoQLAJAM4QIAJEO4AADJEC4AQDKECwCQDOECACRDuAAAyRAuAEAyhAsAkAzhAgAkQ7gAAMkQLgBAMoQLAJAM4QIAJEO4AADJEC4AQDKECwCQDOECACRDuAAAyRAuAEAyhAsAkAzhAgAkQ7gAAMkQLgBAMoQLAJAM4QIAJEO4AADJEC4AQDKECwCQDOECACRDuAAAyRAuAEAyhAsAkAzhAgAkQ7gAAMkQLgBAMoQLAJAM4QIAJEO4AADJEC4AQDKECwCQDOECACSjT+GydOnSGD9+fJSWlkZFRUU0Nzfv13FvvPFGFBUVxRlnnNGXuwUAhri8w2XVqlUxf/78WLRoUbS2tsaMGTNi5syZ0dbW9p3HffbZZzFnzpz46U9/2ufFAgBDW0GWZVk+B1RVVcXUqVNj2bJl3WOTJ0+OWbNmRUNDwz6Pu/zyy2PixIlRWFgYL774YmzYsGGfc7u6uqKrqytnrKSkJEpKSvJZKgAwyOR1xmXnzp3R0tISNTU1OeM1NTWxdu3afR73+OOPx3vvvRe33Xbbft1PQ0NDlJWV5dy+K4oAgKGhKJ/JnZ2dsXv37igvL88ZLy8vj46Ojl6Peffdd2PhwoXR3NwcRUX7d3f19fVRV1eXM+ZsCwCQV7jsVVBQkPNzlmU9xiIidu/eHVdccUXccccdMWnSpP3+/S4LAQC9yStcRo0aFYWFhT3Ormzbtq3HWZiIiB07dsT69eujtbU1rr/++oiI2LNnT2RZFkVFRbF69eo499xzD2D5AMBQktdrXIqLi6OioiKamppyxpuamqK6urrH/BEjRsRbb70VGzZs6L7V1tbGSSedFBs2bIiqqqoDWz0AMKTkfamorq4urrzyyqisrIxp06bFI488Em1tbVFbWxsR374+5aOPPoonn3wyhg0bFlOmTMk5/rjjjovS0tIe4wAA3yfvcJk9e3Zs3749Fi9eHO3t7TFlypRobGyMcePGRUREe3v7936mCwBAX+T9OS4AAIeK7yoCAJIhXACAZAgXACAZwgUASIZwAQCSIVwAgGQIFwAgGcIFAEiGcAEAkiFcAIBkCBcAIBnCBQBIhnABAJIhXACAZAgXACAZwgUASIZwAQCSIVwAgGQIFwAgGcIFAEiGcAEAkiFcAIBkCBcAIBnCBQBIhnABAJIhXACAZAgXACAZwgUASIZwAQCSIVwAgGQIFwAgGcIFAEiGcAEAkiFcAIBkCBcAIBnCBQBIhnABAJIhXACAZAgXACAZwgUASIZwAQCSIVwAgGQIFwAgGcIFAEiGcAEAkiFcAIBkCBcAIBnCBQBIhnABAJIhXACAZAgXACAZwgUASIZwAQCSIVwAgGQIFwAgGcIFAEiGcAEAkiFcAIBkCBcAIBnCBQBIhnABAJIhXACAZAgXACAZfQqXpUuXxvjx46O0tDQqKiqiubl5n3Off/75OP/88+PYY4+NESNGxLRp0+KVV17p84IBgKEr73BZtWpVzJ8/PxYtWhStra0xY8aMmDlzZrS1tfU6/7XXXovzzz8/Ghsbo6WlJc4555y4+OKLo7W19YAXDwAMLQVZlmX5HFBVVRVTp06NZcuWdY9Nnjw5Zs2aFQ0NDfv1O0499dSYPXt23Hrrrb3+e1dXV3R1deWMlZSURElJST5LBQAGmbzOuOzcuTNaWlqipqYmZ7ympibWrl27X79jz549sWPHjjjmmGP2OaehoSHKyspybvsbRQDA4FWUz+TOzs7YvXt3lJeX54yXl5dHR0fHfv2Ou+++O7788su47LLL9jmnvr4+6urqcsacbQEA8gqXvQoKCnJ+zrKsx1hvnn766bj99tvjz3/+cxx33HH7nOeyEADQm7zCZdSoUVFYWNjj7Mq2bdt6nIX5X6tWrYprrrkmnnnmmTjvvPPyXykAMOTl9RqX4uLiqKioiKamppzxpqamqK6u3udxTz/9dMydOzeeeuqpuOiii/q2UgBgyMv7UlFdXV1ceeWVUVlZGdOmTYtHHnkk2traora2NiK+fX3KRx99FE8++WREfBstc+bMiXvvvTfOOuus7rM1RxxxRJSVlR3EhwIADHZ5h8vs2bNj+/btsXjx4mhvb48pU6ZEY2NjjBs3LiIi2tvbcz7T5eGHH45du3bFddddF9ddd133+FVXXRVPPPHEgT8CAGDIyPtzXAAADhXfVQQAJEO4AADJEC4AQDKECwCQDOECACRDuAAAyRAuAEAyhAsAkAzhAgAkQ7gAAMkQLgBAMoQLAJAM4QIAJEO4AADJEC4AQDKECwCQDOECACRDuAAAyRAuAEAyhAsAkAzhAgAkQ7gAAMkQLgBAMoQLAJAM4QIAJEO4AADJEC4AQDKECwCQDOECACRDuAAAyRAuAEAyhAsAkAzhAgAkQ7gAAMkQLgBAMoQLAJAM4QIAJEO4AADJEC4AQDKECwCQDOECACRDuAAAyRAuAEAyhAsAkAzhAgAkQ7gAAMkQLgBAMoQLAJAM4QIAJEO4AADJEC4AQDKECwCQDOECACRDuAAAyRAuAEAyhAsAkAzhAgAkQ7gAAMkQLgBAMoQLAJAM4QIAJEO4AADJEC4AQDL6FC5Lly6N8ePHR2lpaVRUVERzc/N3zl+zZk1UVFREaWlpTJgwIR566KE+LRYAGNryDpdVq1bF/PnzY9GiRdHa2hozZsyImTNnRltbW6/zt2zZEhdeeGHMmDEjWltb45Zbbol58+bFc889d8CLBwCGloIsy7J8DqiqqoqpU6fGsmXLuscmT54cs2bNioaGhh7zb7755njppZdi06ZN3WO1tbXx5ptvxrp163q9j66urujq6soZKykpiZKSknyWCgAMMnmdcdm5c2e0tLRETU1NznhNTU2sXbu212PWrVvXY/4FF1wQ69evj2+++abXYxoaGqKsrCzndsEFF/SIGQZeV1dX3H777fbiMGAvDh/24vBiPw4f/bEXeYVLZ2dn7N69O8rLy3PGy8vLo6Ojo9djOjo6ep2/a9eu6Ozs7PWY+vr6+Oyzz7pvH374YaxZs8b/CA8DXV1dcccdd9iLw4C9OHzYi8OL/Th89MdeFPXloIKCgpyfsyzrMfZ983sb38tlIQCgN3mdcRk1alQUFhb2OLuybdu2HmdV9jr++ON7nV9UVBQjR47Mc7kAwFCWV7gUFxdHRUVFNDU15Yw3NTVFdXV1r8dMmzatx/zVq1dHZWVlDB8+PM/lAgBDWd5vh66rq4vHHnssVqxYEZs2bYoFCxZEW1tb1NbWRsS3r0+ZM2dO9/za2tr44IMPoq6uLjZt2hQrVqyI5cuXx4033rjf91lSUhK33Xaby0eHAXtx+LAXhw97cXixH4eP/tiLvN8OHfHtB9AtWbIk2tvbY8qUKfGHP/whfvKTn0RExNy5c+P999+PV199tXv+mjVrYsGCBfH222/HmDFj4uabb+4OHQCA/dWncAEAOBR8VxEAkAzhAgAkQ7gAAMkQLgBAMg6bcFm6dGmMHz8+SktLo6KiIpqbm79z/po1a6KioiJKS0tjwoQJ8dBDDw3QSge/fPbi+eefj/PPPz+OPfbYGDFiREybNi1eeeWVAVzt4Jbv38Veb7zxRhQVFcUZZ5zRvwscQvLdi66urli0aFGMGzcuSkpK4sQTT4wVK1YM0GoHt3z3YuXKlXH66afHkUceGaNHj46rr746tm/fPkCrHbxee+21uPjii2PMmDFRUFAQL7744vcec1Ceu7PDwJ/+9Kds+PDh2aOPPppt3Lgxu+GGG7Kjjjoq++CDD3qdv3nz5uzII4/Mbrjhhmzjxo3Zo48+mg0fPjx79tlnB3jlg0++e3HDDTdkd911V/b3v/89e+edd7L6+vps+PDh2T//+c8BXvngk+9e7PXpp59mEyZMyGpqarLTTz99YBY7yPVlLy655JKsqqoqa2pqyrZs2ZL97W9/y954440BXPXglO9eNDc3Z8OGDcvuvffebPPmzVlzc3N26qmnZrNmzRrglQ8+jY2N2aJFi7Lnnnsui4jshRde+M75B+u5+7AIlzPPPDOrra3NGTv55JOzhQsX9jr/N7/5TXbyySfnjP3iF7/IzjrrrH5b41CR71705pRTTsnuuOOOg720IaevezF79uzst7/9bXbbbbcJl4Mk3734y1/+kpWVlWXbt28fiOUNKfnuxe9+97tswoQJOWP33XdfNnbs2H5b41C0P+FysJ67D/mlop07d0ZLS0vU1NTkjNfU1MTatWt7PWbdunU95l9wwQWxfv36+Oabb/ptrYNdX/bif+3Zsyd27NgRxxxzTH8sccjo6148/vjj8d5778Vtt93W30scMvqyFy+99FJUVlbGkiVL4oQTTohJkybFjTfeGF9//fVALHnQ6steVFdXx9atW6OxsTGyLItPPvkknn322bjooosGYsn8l4P13N2nb4c+mDo7O2P37t09vqSxvLy8x5cz7tXR0dHr/F27dkVnZ2eMHj2639Y7mPVlL/7X3XffHV9++WVcdtll/bHEIaMve/Huu+/GwoULo7m5OYqKDvmf9qDRl73YvHlzvP7661FaWhovvPBCdHZ2xi9/+cv497//7XUuB6Ave1FdXR0rV66M2bNnx3/+85/YtWtXXHLJJXH//fcPxJL5LwfrufuQn3HZq6CgIOfnLMt6jH3f/N7GyV++e7HX008/HbfffnusWrUqjjvuuP5a3pCyv3uxe/fuuOKKK+KOO+6ISZMmDdTyhpR8/i727NkTBQUFsXLlyjjzzDPjwgsvjHvuuSeeeOIJZ10Ognz2YuPGjTFv3ry49dZbo6WlJV5++eXYsmWLr505RA7Gc/ch/8+yUaNGRWFhYY9a3rZtW48y2+v444/vdX5RUVGMHDmy39Y62PVlL/ZatWpVXHPNNfHMM8/Eeeed15/LHBLy3YsdO3bE+vXro7W1Na6//vqI+PbJM8uyKCoqitWrV8e55547IGsfbPrydzF69Og44YQToqysrHts8uTJkWVZbN26NSZOnNivax6s+rIXDQ0NMX369LjpppsiIuK0006Lo446KmbMmBF33nmnM/QD6GA9dx/yMy7FxcVRUVERTU1NOeNNTU1RXV3d6zHTpk3rMX/16tVRWVkZw4cP77e1DnZ92YuIb8+0zJ07N5566inXjQ+SfPdixIgR8dZbb8WGDRu6b7W1tXHSSSfFhg0boqqqaqCWPuj05e9i+vTp8fHHH8cXX3zRPfbOO+/EsGHDYuzYsf263sGsL3vx1VdfxbBhuU91hYWFEfH//2ufgXHQnrvzeilvP9n79rbly5dnGzduzObPn58dddRR2fvvv59lWZYtXLgwu/LKK7vn731L1YIFC7KNGzdmy5cv93bogyTfvXjqqaeyoqKi7MEHH8za29u7b59++umhegiDRr578b+8q+jgyXcvduzYkY0dOza79NJLs7fffjtbs2ZNNnHixOzaa689VA9h0Mh3Lx5//PGsqKgoW7p0afbee+9lr7/+elZZWZmdeeaZh+ohDBo7duzIWltbs9bW1iwisnvuuSdrbW3tfmt6fz13HxbhkmVZ9uCDD2bjxo3LiouLs6lTp2Zr1qzp/rerrroqO/vss3Pmv/rqq9mPf/zjrLi4OPu///u/bNmyZQO84sErn704++yzs4jocbvqqqsGfuGDUL5/F/9NuBxc+e7Fpk2bsvPOOy874ogjsrFjx2Z1dXXZV199NcCrHpzy3Yv77rsvO+WUU7IjjjgiGz16dPazn/0s27p16wCvevD561//+p3//99fz90FWeZcGQCQhkP+GhcAgP0lXACAZAgXACAZwgUASIZwAQCSIVwAgGQIFwAgGcIFAEiGcAEAkiFcAIBkCBcAIBn/D5DFpefk3GIvAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "sns.heatmap(df, index='g', columns='S_r', values='B_error')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}