Examples/NV Play.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "5c8736f7-3fcd-49e0-9367-6a3627c14b26",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Populating the interactive namespace from numpy and matplotlib\n"
     ]
    }
   ],
   "source": [
    "%load_ext autoreload\n",
    "%autoreload 2\n",
    "%pylab inline\n",
    "pylab.rcParams['figure.figsize'] = (12, 8)\n",
    "import warnings; warnings.filterwarnings(\"ignore\")\n",
    "\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "from sympy import symbols\n",
    "\n",
    "import sys; sys.path.append('../')\n",
    "from PDE_FIND import *\n",
    "\n",
    "import sys; sys.path.append('../../parametric-discovery/')\n",
    "from best_subset import *\n",
    "from frols import frols\n",
    "from p_linear_regression import PLinearRegression\n",
    "from r_pca import R_pca\n",
    "from pde_diff_new import RPCA\n",
    "from RobustPCA.rpca import RobustPCA\n",
    "\n",
    "import scipy.io as sio\n",
    "from scipy.signal import savgol_filter\n",
    "import itertools\n",
    "\n",
    "from sklearn.preprocessing import normalize\n",
    "from sklearn.linear_model import Ridge, BayesianRidge, LinearRegression as SkLinearRegression\n",
    "from sklearn.base import BaseEstimator\n",
    "from sklearn.pipeline import make_pipeline\n",
    "from sklearn.metrics import mean_squared_error\n",
    "from abess.linear import LinearRegression\n",
    "from bess import PdasLm\n",
    "import pysindy as ps\n",
    "\n",
    "from tqdm import trange\n",
    "from os.path import join as join_path\n",
    "FIGURE_EXPORT_PATH = \"/Users/pongpisit/Documents/figure_export/\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "f428bf63-f8aa-4de3-97c8-dd17f45a8c55",
   "metadata": {},
   "outputs": [],
   "source": [
    "def nonzero(wei):\n",
    "    return tuple(np.nonzero(wei)[0])\n",
    "\n",
    "def check_percent(bic_scores, complexities):\n",
    "    slope = (bic_scores[1:]-bic_scores[:-1])/(complexities[1:]-complexities[:-1])\n",
    "    slope_index = np.argmin(slope)\n",
    "    percent_improve = 100*np.abs(bic_scores[slope_index+1]-bic_scores[slope_index])/np.abs(bic_scores[slope_index])\n",
    "    percent_from_1 = 100*np.abs(bic_scores[slope_index+1]-bic_scores[0])/np.abs(bic_scores[0])\n",
    "    return slope_index+1, percent_improve, percent_from_1\n",
    "\n",
    "def percent_coeff(pred, ground):\n",
    "    return 100*(np.abs(pred-ground)/np.abs(ground))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "0e366b04-cdf3-40e8-9cfc-b0e0c25ea165",
   "metadata": {},
   "outputs": [],
   "source": [
    "U = np.load('../Datasets/course_timestep_cylinder/Cylinder_U.npy')\n",
    "V = np.load('../Datasets/course_timestep_cylinder/Cylinder_V.npy')\n",
    "W = np.load('../Datasets/course_timestep_cylinder/Cylinder_W.npy')\n",
    "\n",
    "dt = 0.2\n",
    "dx = 0.02\n",
    "dy = 0.02"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "38789883-e9ed-4262-b484-0a86bf76991f",
   "metadata": {},
   "outputs": [],
   "source": [
    "xmin = 100\n",
    "xmax = 425\n",
    "ymin = 15\n",
    "ymax = 185\n",
    "\n",
    "x = (np.arange(449)*dx)[xmin:xmax]\n",
    "y = (np.arange(199)*dy)[ymin:ymax]\n",
    "t = np.arange(151)*dt\n",
    "\n",
    "W = W[xmin:xmax,ymin:ymax,:]\n",
    "U = U[xmin:xmax,ymin:ymax,:]\n",
    "V = V[xmin:xmax,ymin:ymax,:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "450112c3-065c-4a0c-94b2-4801932031fc",
   "metadata": {},
   "outputs": [],
   "source": [
    "def col_vec(vec): return vec.reshape(-1, 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "435ad57b-c77f-46da-bd05-7ce6344148a9",
   "metadata": {},
   "outputs": [],
   "source": [
    "fd = ps.SmoothedFiniteDifference(axis=-1)\n",
    "wt = fd._differentiate(W, t=t).reshape(-1, 1)\n",
    "\n",
    "fd = ps.SmoothedFiniteDifference(axis=0)\n",
    "wx = fd._differentiate(W, t=x)\n",
    "wxx = fd._differentiate(wx, t=x)\n",
    "fd = ps.SmoothedFiniteDifference(axis=1)\n",
    "wy = fd._differentiate(W, t=y)\n",
    "wyy = fd._differentiate(wy, t=y)\n",
    "wxy = fd._differentiate(wx, t=y)\n",
    "\n",
    "fd = ps.SmoothedFiniteDifference(axis=-1)\n",
    "wt = fd._differentiate(W, t=t).reshape(-1, 1)\n",
    "\n",
    "w = col_vec(W)\n",
    "u = col_vec(U)\n",
    "v = col_vec(V)\n",
    "wt = col_vec(wt)\n",
    "wx = col_vec(wx)\n",
    "wy = col_vec(wy)\n",
    "wxx = col_vec(wxx)\n",
    "wxy = col_vec(wxy)\n",
    "wyy = col_vec(wyy)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "a0908539-3e15-486d-a973-120c7c4730e7",
   "metadata": {},
   "outputs": [],
   "source": [
    "theta = np.hstack([w, u, v, wx, wy, wxx, wxy, wyy])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "8ddec8ed-1509-4ef7-9f77-e20059e0beb6",
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.preprocessing import PolynomialFeatures\n",
    "poly_transformer = PolynomialFeatures(interaction_only=True, include_bias=False)\n",
    "theta = poly_transformer.fit_transform(theta)\n",
    "basis_names = ['w', 'u', 'v', 'w_{x}', 'w_{y}','w_{xx}','w_{xy}','w_{yy}']\n",
    "description = [''.join([basis_names[int(e)] for e in fn.replace('x', '').split()]) for fn in poly_transformer.get_feature_names_out()]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "e77fa604-381a-4851-96a9-50d06ba8ead7",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Form Theta using up to quadratic polynomials in all variables.\n",
    "# X_data = np.hstack([w,u,v])\n",
    "# X_ders = np.hstack([np.ones((w.shape[0],1)), wx, wy, wxx, wxy, wyy])\n",
    "# X_ders_descr = ['','w_{x}', 'w_{y}','w_{xx}','w_{xy}','w_{yy}']\n",
    "# X, description = build_Theta(X_data, X_ders, X_ders_descr, 2, data_description = ['w','u','v'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "c76fc38b-a2f2-4853-8c2b-5810550b518e",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(3,)\n",
      "-22756481.79869665\n",
      "(3, 22)\n",
      "-23900992.342734784\n",
      "(7, 16, 22)\n",
      "-49397232.076707594\n",
      "(5, 7, 16, 22)\n",
      "-52873990.62114105\n",
      "(2, 5, 7, 16, 22)\n",
      "-54806790.22714589\n",
      "(2, 5, 7, 15, 16, 22)\n",
      "-55178672.09783288\n",
      "(2, 5, 7, 13, 15, 16, 22)\n",
      "-55340928.83784835\n",
      "(5, 7, 15, 16, 22, 24, 26, 28)\n",
      "-55858475.88903728\n",
      "(2, 5, 7, 15, 16, 22, 24, 26, 28)\n",
      "-56784891.6495166\n"
     ]
    }
   ],
   "source": [
    "lrs = []\n",
    "ebics = []\n",
    "eff_indices = []\n",
    "max_complexity = 10\n",
    "for s in range(1, max_complexity):\n",
    "    lr = LinearRegression(support_size=s, ic_type='bic', always_select=None)\n",
    "    lr.fit(theta, wt)\n",
    "    lrs.append(lr)\n",
    "    eff_indices.append(nonz(lr.coef_))\n",
    "    ebics.append(lr.ic_)\n",
    "    print(nonz(lr.coef_))\n",
    "    print(lr.ic_)\n",
    "abess_feature_hierarchy = set(eff_indices)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "5d275fc9-19ab-4a50-84ff-814e77fe749f",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "RFE_PLinearRegression\n",
      "(5, 7, 15, 16, 22, 24, 26, 28)\n",
      "(5, 7, 15, 16, 22, 24, 28)\n",
      "(5, 7, 16, 22, 24, 28)\n",
      "(5, 7, 16, 22, 24)\n",
      "(5, 7, 16, 22)\n",
      "(7, 16, 22)\n",
      "(16, 22)\n",
      "(16,)\n",
      "--------------------------------------------------\n",
      "(3, 22)\n",
      "(3,)\n",
      "--------------------------------------------------\n",
      "(2, 5, 7, 15, 16, 22, 24, 26, 28)\n",
      "(2, 5, 7, 16, 22, 24, 26, 28)\n",
      "(2, 5, 7, 16, 22, 24, 28)\n",
      "(2, 5, 7, 16, 22, 28)\n",
      "(2, 5, 7, 16, 22)\n",
      "(5, 7, 16, 22)\n",
      "--------------------------------------------------\n",
      "(2, 5, 7, 13, 15, 16, 22)\n",
      "(2, 5, 7, 15, 16, 22)\n",
      "(2, 5, 7, 16, 22)\n",
      "--------------------------------------------------\n",
      "(2, 5, 7, 16, 22)\n",
      "--------------------------------------------------\n",
      "(5, 7, 16, 22)\n",
      "--------------------------------------------------\n",
      "(2, 5, 7, 15, 16, 22)\n",
      "--------------------------------------------------\n",
      "(7, 16, 22)\n",
      "--------------------------------------------------\n",
      "SelectKBest\n",
      "(5, 7, 15, 16, 22, 24, 26, 28)\n",
      "(5, 7, 15, 16, 22, 24, 28)\n",
      "(5, 7, 16, 22, 24, 28)\n",
      "(5, 7, 16, 22, 24)\n",
      "(5, 7, 16, 22)\n",
      "(7, 16, 22)\n",
      "(16, 22)\n",
      "(16,)\n",
      "--------------------------------------------------\n",
      "(3, 22)\n",
      "(3,)\n",
      "--------------------------------------------------\n",
      "(2, 5, 7, 15, 16, 22, 24, 26, 28)\n",
      "(2, 5, 7, 16, 22, 24, 26, 28)\n",
      "(2, 5, 7, 16, 22, 24, 28)\n",
      "(2, 5, 7, 16, 22, 28)\n",
      "(2, 5, 7, 16, 22)\n",
      "(5, 7, 16, 22)\n",
      "--------------------------------------------------\n",
      "(2, 5, 7, 13, 15, 16, 22)\n",
      "(2, 5, 7, 15, 16, 22)\n",
      "(2, 5, 7, 16, 22)\n",
      "--------------------------------------------------\n",
      "(2, 5, 7, 16, 22)\n",
      "--------------------------------------------------\n",
      "(5, 7, 16, 22)\n",
      "--------------------------------------------------\n",
      "(2, 5, 7, 15, 16, 22)\n",
      "--------------------------------------------------\n",
      "(7, 16, 22)\n",
      "--------------------------------------------------\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "{1: ((3,), 6701971.368732434), 2: ((16, 22), -8899707.81245153), 3: ((7, 16, 22), -19938776.68029077), 4: ((5, 7, 16, 22), -23415533.17342277), 5: ((2, 5, 7, 16, 22), -25348336.844904445), 6: ((2, 5, 7, 15, 16, 22), -25720216.463609237), 7: ((2, 5, 7, 16, 22, 24, 28), -26705352.991840046), 8: ((2, 5, 7, 16, 22, 24, 26, 28), -26997442.985712484), 9: ((2, 5, 7, 15, 16, 22, 24, 26, 28), -27326435.69842902)}"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "br1 = backward_refinement(abess_feature_hierarchy, (theta, wt), mode='rfe', ic_type='bic', verbose=True)\n",
    "br2 = backward_refinement(abess_feature_hierarchy, (theta, wt), mode='k best', ic_type='bic', verbose=True)\n",
    "br_u=br1+br2; br_u"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "6cf2cc6d-1461-46bc-85cd-70d791b12514",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(3,) 6701971.368732434\n",
      "(16, 22) -8899707.81245153\n",
      "(7, 16, 22) -19938776.68029077\n",
      "(5, 7, 16, 22) -23415533.17342277\n",
      "(2, 5, 7, 16, 22) -25348336.844904445\n",
      "(2, 5, 7, 15, 16, 22) -25720216.463609237\n",
      "(2, 5, 7, 16, 22, 24, 28) -26705352.991840046\n",
      "(2, 5, 7, 16, 22, 24, 26, 28) -26997442.985712484\n",
      "(2, 5, 7, 15, 16, 22, 24, 26, 28) -27326435.69842902\n"
     ]
    }
   ],
   "source": [
    "complexities = np.array(sorted(list(br_u.track.keys())))\n",
    "bic_scores = []\n",
    "eff_feature_names = []\n",
    "for com in complexities:\n",
    "    indices, bic = br_u.track[com]\n",
    "    bic_scores.append(bic)\n",
    "    print(indices, bic)\n",
    "bic_scores = np.array(bic_scores)\n",
    "min_bic_score = np.min(bic_scores)\n",
    "div_scale = 1e3\n",
    "rel_bic_scores = (bic_scores-min_bic_score)/div_scale"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "6f8bb984-1bd0-46f8-b2e4-931e6485a9a4",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAUoAAADyCAYAAAArx4ypAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAA9hAAAPYQGoP6dpAAApbUlEQVR4nO3df3ATdf4/8Oemgig/sg2CP+qPNqGKhzo2adXz9DOcTcGRo/5oSqlKxJYm1PPOc4SmKMiJCqR4pzeno01p0YqWJoE5cPTQBJy70/NHaMRT1F7Jth7ij1Hyoyg6nnS/f/SbvabNtkmbzebH6zGTIcl7N/vaO3nx3vdPhud5HoQQQkQp5A6AEEJSHSVKQggZAyVKQggZAyVKQggZwylyB5AK5s2bB41GM+oxR48eRV5e3rh+n85Nj2vTuelxbanP9fl8OHToUOSXPOEXL16ckGPo3ImfK+e16dz0uLbU50Y7hh69Y1RdXS3LuRMhV8xy3e9Er52O9yxnzFl1z+NOzRlkov+ySimVY5NKtt1ztt0vz6f2PUeLjdoox9DT04O2tjb09fUhPz8fNTU1KCwsTNr15aydySXb7jnb7hdIv3tmeJ5m5pSXl2PPnj0jvt+2bRtWrFgBhmHA87zwZ2trK5YvX578QAkhkouWD6iNUkRPTw9WrFiBgYEBnDx5MuLP2tpaHD58WO4QCSFJQolSRFtbGxiGiVrGMAxaW1uTHBEhRC6UKEX09fVBrFWC53n09fUlNyBCiGwoUWJwEGp5eTk6OjqE7/Lz80etUebn5ycpOkJIMnR0dKC8vBxHjx4dUUadOYjeeNvT04O5c+diYGBgxPEKhQLd3d2YM2dOskIkhCQJdebEobCwEK2trVAoFFAoBv9nYhgGCoUCra2tlCQJySKUKEexfPlydHd3o6GhAZdffjkUCgU8Hg8NDSIky1CiHMOcOXOwadMmvPzyy+B5Hl1dXXKHRAhJMkqUMTrnnHNQVlaG9vZ2uUMhhCQZJco4GI1GvPHGG/D5fHKHQghJIkqUcbjpppswffp0PP/883KHQghJIkqUcTj99NNRWVmJ9vZ20cHohJDMQ4kyTkajEb29vXjzzTflDoUQkiSUKBF9Zo6Ya6+9FhdccAF16hCSYUabmUOJEkBeXh727NkT0xp5CoUCy5Ytg91ux/fff5+E6AghyVBdXY09e/ZE3VOHEuU4LFu2DKFQCC+99JLcoRBCkoAS5ThceOGFuOqqq+jxm5AsQYlynIxGI/bu3YuvvvpK7lAIIRKjRDlOVVVVyMnJwYsvvih3KIQQiSVtczG32w0ACAaD8Hg8qKqqglarBQBwHAen0wm1Wg2O42AymcCyrGRliaBSqbB48WK0t7fj3nvvTdjvEkJSULK2gGRZlu/q6uJ5nuebm5t5tVotlGm1WuG9z+fjDQaDpGXDjXfrzN27d/MA+Pfff39c5xNCUk+0fJC0R2+HwyHUIAFE1PyGUqvVQu1TirJEuv7663HGGWfQlEZCMlzSEqVerxfeOxwOmM1mAIOP5CqVKuJYlUoFr9crSVkiTZ48GdXV1di+fTt++umnhP42ISR1JLUzx+v1wmKxoKysDCaTCcBgm2U0fr9fkrJEMxqN+PLLL7Fv376E/zYhJDUkrTMHALRaLdRqNSwWC5xOJwwGg+ixYslOirLwFMaw6urqmGbpAIBOp8PFF1+M9vZ2LFy4MKZzCCGpo6OjI2L6crQpjElNlMBg22RlZSXKysoQCATAsuyImp7f7wfLspKURROewjgeDMPAaDRiw4YN6O/vx4wZM8b1O4QQeQyvGA2tNIUl5dHb7XYjNzdX+KxWqwEMdroMbbscqri4WJIyKdx+++344Ycf4HQ6Jfl9Qoi8kpIoVSpVRPLyer1gWVZ4FB+K4zgUFxeDZVlJyqRw7rnnorS0lKY0EpKhkvLordVqUVVVBZvNBgBwuVwRm3Q5HA5YLBaUlJTA4/HA4XBIWiYFo9EIo9GIvr4+5OfnS3otQkhyMTxPS3VH2/A8Xt9++y3OOussNDY2Yu3atQmKjBCSbNHyAc31TpBp06ahoqKCtokgJANRokwgo9GInp4evPPOO3KHQghJoDHbKPv6+uBwOOByuRAIBITvVSoVysrKYDAYqE3u/5s/fz7OPfdctLe346qrrpI7HEJIgoyaKBsbG8EwDJYsWYLVq1ePKH/vvffwzDPPgGEYbNq0SbIg00VOTg5uv/12NDc34/HHH8epp54qd0iEkAQQ7czZsmULTCYTlErlmD8SCoWwefPmtE2WOp0OeXl5cc3IEfPxxx/jZz/7GZxOJyoqKhIUISFEauEZOkePHo0YlQNQrzeAxPR6D3XFFVfg7LPPxu7duxP2m4SQ5KBe7yQxGo145ZVX8PXXX8sdCiEkAUYkyt7eXtTX1+PgwYPo6+vDggULUFhYiIMHD8oQXnpaunQpAGDHjh0yR0IISYQRibK5uRl6vR7PPPMMmpubYbVa8eqrr6Kzs1OO+NLSGWecgUWLFtGURkIyxIhEWVZWhoqKCjzzzDPQ6/UoKiqCWq0WXWiCRGc0GnHgwAF89NFHcodCCJmgEYkyGAxizZo1AIDS0lIAwGOPPYZQKJTcyNLcokWLkJubS9tEEJIBRiTKiooKNDY2RnynVCqpRhmnU089FUuXLsX27dtx8uRJucMhhExA1F7v4WMnS0pKcODAAezatQv79+9Hf39/UoJLd0ajEZ999hlef/11uUMhhEyA6MycUCiEuro6uN1uzJw5U1jowe/3g2EYFBcXw+Fw0Ireo7jyyitx4YUXor29nWrkhKQx0XGULS0tsFqt8Pv96OnpweHDh3H48GH4/X4cO3YMTz/9dNrOxEmW8DYRO3fuxLfffit3OISQcRJNlAUFBSgoKBA9Ua1WS7a1QrKFNxcbusFQotx+++04ceIEdu3alfDfJoQkTkdHB8rLy6NuLiaaKDmOG/OHPR7PxCJLEeHNxSY6zzuaCy64APPnz6cxlYSkuOrqauzZswd5eXkjykTbKA0GA4qLi8EwDFQqlfD90N0Npd5eIVMYjUbU1tbiyJEjOO+88+QOhxASJ9FEWVBQgAMHDqC3txdut1vYE5tlWej1+lEfy0mkiooK/PrXv8YLL7wwYugVIST1jblwb0FBAerq6pIRS8aaMWMGbr75ZrS3t8NisYBhGLlDIoTEYUKrB1EHReyMRiM+/vhjHDhwQO5QCCFxGnei7O/vz5jOnGTQ6/U4++yzqVOHkDQkmiiLi4uRk5Mj+mJZFk1NTcmMNa2Ft4no6OjAjz/+KHc4hJA4jDrg/MCBAzh58mTU18DAQNR9dIg4o9GIY8eO4a9//avcoRBC4iCaKIuKitDb2zvqyWVlZQkPKJNdcsklKCoqosdvQtLMqG2Ut9xyy6gnh5dhS3dSzswZzmg04qWXXooYj0oIkd+4ZuZkEyln5gxXXV2NgYEBWjGekBQz2swcSpRJduaZZ+L666+nx29C0khMiXLr1q1xfU9GZzQa8fbbb6O7u1vuUAghMYgpUep0OtTX10d8V19fnzGrByVbeXk5lEolbRNBSJqIKVEWFRWhoaFBSJYrV66EyWTC5ZdfLmVsGWvKlCmoqqrC888/j4GBAbnDIYSMIeY2yoKCAphMJsycORMrV65EUVGRlHFlPKPRiP/85z/4+9//LncohJAxxNWZs2nTJtjtdlrZPAGuvvpqqNVq6tQhJA3ElChDoRCqqqrQ0tKC0tJSNDc3o6qqijYZm4DwNhEOhwMnTpyQOxxCyChiSpSbNm2CzWYTdmdkWRY2m41qlhO0bNkyfPvtt/jLX/4idyiEkFEwfHh7xSxWXl6OPXv2yHLta6+9FlOnTsXevXtluT4hJFK0fEADzpHcKYzDGY1GuFwufP7550m/NiHkf+KewhgKhdDX1xfXRfr7+9O2zTKZUxiHq6ysxKRJk/DCCy8k/dqEkP+JewqjUqmEy+WKeQXznTt3wm63Y8aMGROLNAuxLIubbroJzz33HKgVhJDUJLpnTl1dHd577z0sWbIEGo0GJSUlUKvVYFkWwWAQHMfh3XffRW9vL8xmMyoqKpIZd0YxGo1YtGgRDh48SONTCUlBo24uVlRUBLvdjlAoBLvdjnfffRfBYBAsy0Kj0cBsNtNujAmwYMECzJ49G+3t7ZQoCUlBY+7CCAw+itNOjNI55ZRTcNttt+GFF15AU1MTJk2aJHdIhJAhYkqUieD1euF2uwEAHo8HLS0tYFkWAMBxHJxOJ9RqNTiOg8lkkrQsFRmNRjz++ON47bXXsGjRIrnDIYQMxSeJ1WqNeK/VaoXPQ9/7fD7eYDBIWjbc4sWL47kVSQwMDPCXXnopv2TJErlDISSrRcsHSRlH6fV6I2bxGAwGeL1ecBwHjuMijlWr1ULNU4qyVBWe0rh7924EAgG5wyGEDJGURKnVatHS0iJ8DgaDAACVSgW32w2VShVxvEqlEh7VE12Wym699Vb897//hcPhkDsUQsgQom2Uwwec5+fnAxgcWD6e8ZIGg0F439nZCb1eLww1isbv90tSFk14Zk5YdXW1LIPPzznnHJSVlaG9vR0mkynp1yckG3V0dETMyos2M0c0UbpcLpjNZpjNZlRWVgqJ8tixY3C5XGAYZsxdGqMJBoNwOp3o6uoa87hklYVn5qSCO+64A7feeit8Ph80Go3c4RCS8YZXjIZWmsJEE6VarUZXV9eIcX0FBQUoKChAKBTCrl274k6WFosFLpdL6IFmWXZETc/v94NlWUnKUt2NN96I6dOn4/nnn8fvf/97ucMhhGCUNspQKDTq4GelUhn3lLumpiZYLBao1WoEg0EEg0Ho9fqoxxYXF0tSlupOP/10VFZWor29naY0EpIiRBOlWHveUPH0zjqdTmi1WiFJ2u12sCwLtVodcRzHcSguLpasLB0YjUb09vbizTfflDsUQghGefT2+XxjnhzLMcBgoqqsrIz4jmVZocPC4XDAYrGgpKQEHo8notdXirJUd+211yIvLw933303Lr74YuTn56OmpgaFhYVyh0ZIVhJduLelpQUMw2DFihVRT9y6dSsCgQBWr14taYDJIOfCvdFs27YNNTU1AAbHVyoUCvA8j9bWVixfvlze4AjJcNHywagrnC9YsAAMw2DlypXC4hccx6G5uRkA8Oqrr0oYbvKkUqLs6enB3Llzo25jq1Ao0N3djTlz5sgQGSHZIe4Vzl977TWUlpaitrYWOp0OWq0WBoMBWq02Y5JkqmlrawPDMFHLGIZBa2trkiMihIy5KEZDQwMaGhrQ29sLALSsmsT6+vpEe7t5no975XlCyMSJ1ij3798f8Tk8fnKogwcPShJUNsvPzx+1Rhke+E8ISR7RROlyuXD8+HFhL5xor87OzmTGKhk5NxcbrqamZtQaZW1tbZIjIiQ7jLa5mGhnjkKhEK3ZAIN/aRmGwcmTJxMXqUxSqTMHAJ599lnU1taCYRgMDAwIibOtrQ133nmnzNERktmi5QPRNkqTyQSLxTJiFZ6wY8eOoampKbEREgDA8uXLcc0116C1tVVok9yxYwdmzpwpb2CEZCnRRFlZWTlqx41SqRwxiJwkzpw5c4Q1PHmex5dffol169bhV7/6FRQK2o6dkGQS/RtXWlo65smxHEMmjmEYPProo/jXv/4Fu90udziEZJ1Rqyb79+9HfX09Xn/9dQCDQ1cWLFiAmTNn4q677kpKgGTQ1VdfjRtuuAHr16/HTz/9JHc4hGQV0US5c+dO6PV6eDweGAwGHDx4UNigq7GxEYcPH0Z9fX0yY816jzzyCP7973+jvb1d7lAIySqibZR2ux2BQABKpRLBYFDo3Ak/bq9evRpLlixJWqBkcJ91g8GAhx56CLfddhtOPfVUuUMiJCuI1iiLi4uhVCoBDK70U1ZWNqJNsqSkRNroyAgbNmzAZ599FrEHESFEWqKJMjc3N+Lz8DUeox1DpHfxxRfj9ttvx6OPPooTJ07IHQ4hWUE0Ufp8Phw/flz4HAqFoh5Dkm/9+vX45ptv8OSTT8odCiFZIeaZOeGZOMM/Z8LMHJ1Oh7y8PNl2XxyP+vp62O12cBwnNJEQQsYvvBvj0aNHR2x+KFqjNJlM8Pv9wisQCER89vv9qKurkzz4ZAjvwpguSRIA1q5dixMnTuDxxx+XOxRCMkJ1dTX27NmDvLy8EWWiidJsNkOpVIq+WJaF2WyWNHAiLi8vD3fddRf++Mc/4tixY3KHQ0hGE02Uo+3AGM8xRDqNjY3geR5Wq1XuUAjJaDRpOI3NmjUL9957L5588kl88cUXcodDSMaiRJnm7rvvPkyZMgWPPvqo3KEQkrEoUaY5pVKJhoYG2Gw22iaCEIlQoswAv/nNb6BSqbBhwwa5QyEkI8WVKLds2YKqqioAwL59+9Df3y9JUCQ+U6dOxf3334/nnnsO3d3dcodDSMaJOVE2NjaCZVno9XoAg2tRut1uyQIj8TGbzcjLy8P69evlDoWQjBNzoiwpKUFdXV3UOd/pLpU2FxuvU089FQ8++CA6Ozvx/vvvyx0OIWlntM3FYk6U4X29h05j9Hg8CQhPfuk4MyeaO+64A3PmzMG6devkDoWQtDPazBzR9SiHKyoqQnFxMWbOnAmXywW3200DnVPMpEmThLUq3377bVx11VVyh0RIRoi5RllaWgq73Y6ioiLwPA+bzYbrrrtOytjIOCxduhSXXHIJ1q5dK3cohGSMmGuU9fX1ePrpp7F582Yp4yETpFAo8PDDD+Pmm2/G66+/jl/+8pdyh0RI2ou5RulyubB161bs379fynhIAtx4440oKSnBAw88AJFV9AghcYi5RtnV1QWlUolQKIStW7dCpVJBr9djxowZUsZHxoFhGDzyyCNYuHAhXnnlFSxatEjukAhJazHXKMOLwyqVSvA8j4aGhoxZjzITlZWV4f/+7/+wdu1aDAwMyB0OIWkt5kRZVVWF+vp6qFQqeL1euFwudHZ2ShkbmQCGYfDoo4/i4MGD2Llzp9zhEJLWYk6UXV1d0Ol08Pv9ePrpp1FQUCBlXCQBrrnmGlx//fV48MEHM2LLDkLkEnOitFqtWLFihZSxEAk88sgj+OSTT7B9+3a5QyEkbYluLhaLvr4+5OfnJzAceaTj5mLxqKiogNfrRXd3NyZPnix3OISkpNE2FxNNlLt27Yro1d66dWtEeTAYhMvlwquvvipR2MlTXl6OPXv2yB2GZA4dOoRLL70UTz31FOrr6+UOh5CUFi0fiD56b9y4EQcOHBA+P/PMMwgEAsKL53na1CpNzJs3D7fddhsefvhhfP/993KHQ0jaER1HOTRJAkBLS8uIzcTCS66R1Pf73/8eO3bswFNPPYVVq1bJHQ4haSXmzpzc3FzhfSgUws6dOyO+I6lNo9GgpqYGmzdvpgWXCYlTzIly6CK9SqUSFRUVtHBvmlm3bh2+/fZbPPHEE3KHQkhaGXUKYygUgt1uB8MwcLlcI8q7urpiHjLk9XpRV1c3ojeJ4zg4nU6o1WpwHAeTyQSWZSUry2bnnnsu6uvr8Yc//AF33303VCqV3CERkh74MXAcx5vNZl6n0/EWiyXixXHcWKfzPM/zDoeD7+rq4qNdTqvVCu99Ph9vMBgkLYtm8eLFMd1HJvjqq6/4qVOn8haLRe5QCElJ0fLBmIkyzO12TziA4YnS5/NFJDWe53mWZSUrE5NNiZLnef7+++/nTzvtNP6LL76QOxRCUk60fBDXwr1D7d+/H7t27ZpQbdbtdo94/AvPJZeijAxatWoVJk+ejI0bN8odCiFpIeZl1oDBQegcxwEAeJ7HgQMHcMstt4z74sFgMOr3fr9fkjIx4c3FwjJ1hk5Ybm4uVq9ejQ0bNmDVqlU4//zz5Q6JENmEZ+SERdtcLOZE2djYiGAwCL/fD7VajWAwCLPZnJhIhxFLdlKVhTcXyyb33HMP/vSnP2HDhg0jZl0Rkk2GV4yGVprCYk6UGo0GdXV16O3tBcMwyM/Pn/Bq5yzLjqjp+f1+sCwrSRn5n2nTpmHNmjVYvXo1LBYLCgsL5Q6JkJQVcxulWq3Gp59+ioKCAjidzoRcXGxmT3FxsSRlJFJ9fT3OOussrF+/Xu5QCElpMSfKYDAItVqN/v5+fPPNN1i4cCGam5vjvuDQR2C1Wh1RxnEciouLwbKsJGUk0pQpU7Bu3Trs2LEDH3zwgdzhEJK6xtuF7na7+WAwGNOxLpeLb2ho4AHwDQ0NvMPhEMp8Pp/wXUNDAx8IBCQtiybbhgcN9eOPP/JqtZq/8cYb5Q6FkJQQLR9MaD3KgwcP4vLLL09c1pZJpi+zNpbt27dj2bJleOedd3DFFVfIHQ4hsoqWD0Q7c8bqCQ0EArDb7fB4PImJjsimuroamzZtwtq1a/Haa6/JHQ4hKUe0jXL4+pPDXwBoz+gMkZOTg4cffhgulwt/+9vf5A6HkJQjWqO0Wq0jZuMMR+tRZo6bb74ZOp0ODzzwAP7xj3+AYRi5QyIkZYjWKKMlyS1btqCqqgoAsG/fPmg0GukiI0nFMAweeeQRvPnmm6iqqkJ1dTXWrFmDnp4euUMjRHYxDw9qbGwEy7JCLbK0tDRj1qMMT2EcOo0pG33++ecAAIfDgc7OTmzZsgVz587Fs88+K29ghCRBR0cHysvLo05hjDlRlpSUoK6ubsQ4xUwQnsKYyfO7x9LT04O6ujrhM8/zOHnyJAYGBlBbW4vDhw/LGB0h0quursaePXuQl5c3oizmRNnb2wsAEW1X1OOdOdra2kTbJRmGQWtra5IjIiR1xDzXu6ioCMXFxZg5cyZcLhfcbjesVquUsZEk6uvrEx3FMDAwgL6+vuQGREgKiWs9SofDgaKiIvA8D5vNhuuuu07K2EgS5efni9YoaRgYyXYxJ0oAKCgowObNm7F582YAYw9KJ+mjpqZm1IS4Y8cOVFdXCx0+hGSTURNlf38/du3aFfWxq6ioaNQ1Hkl6KSwsRGtrKxQKBXJyciL+bGtrQ3t7O/bv34+5c+fi8ccfx08//SR3yIQkjWii7O3tRX5+PgwGAzQaDfr6+tDf34+VK1eisLAQOTk51JmTYZYvX47u7m6sXr0aS5YswerVq9Hd3Y0777wTy5YtQ3d3N4xGI1atWgWtVos33nhD7pAJSQrRRTFWrlyJyspKlJaWwuv1wm63Y9++fdDpdGBZVljINxNk+6IY8fJ6vaivr8e7776LO+64A01NTZg9e7bcYRGSEHEtiqHRaITZOVqtFvv27YPNZkNRUZG0UZKUp9Vq8dZbb2Hr1q1Ys2YNdu/ejY0bN8JkMiEnJ0fu8AhJONFH7+E9oFqtdkSSPHjwoCRBJRvNzImfQqGAyWRCd3c3KioqcNddd+Gqq66i5hiStsY1M+fYsWM4fvw4+vv70d/fj1AoJLwPvzJlCiPNzBm/M844A1u3bsU///lP/PTTT7jyyiuxcuXKUXe9JCQVjWtmjtVqBcuyyM3NRW5uLgwGg/A+/D3VHkjYz3/+c3g8HjzxxBPo6OjARRddhG3btmFgYEDu0AiZMNFEaTKZ4Pf7hVcgEIh4f+zYMRQUFCQzVpLiTjnlFPz2t79Fd3c3Fi5ciJqaGlx77bV4//335Q6NkAkRTZRmsxlKpVL0xbKssOQaIUOdddZZ2L59O15//XUEg0FotVr87ne/Q39/v9yhETIuookylt5t6gEno5k/fz7ee+89bN68GVu3bsVFF12Ejo4OmhJJ0k5cUxgJidfkyZOxevVqfPzxx/jFL36BW2+9FaWlpfj444/lDo2QmFGiJElx3nnnwel0Yu/evThy5Aguu+wyNDY24rvvvhOO6enpwZo1a2h1dZJyYl5mjZBEWLhwIT744ANs2bIFGzduxIsvvognnngCwWAQdXV1YBgGPM+DYRg0NTWhtbUVy5cvlztskuWoRkmSbsqUKVi3bh0OHTqEyy67DBUVFaitrcXAwICwqjqtrk5SCSVKIhu1Wo2XXnoJt9xyi+gxtLo6SQWUKEFTGOXEMAwmT54sumgwra5OkiUhm4tlMprCKK/8/HwoFNH/U+R5Hh6PB8899xytf0oklZDNxQiRymirqzMMg9zcXNx5552YPXs2Fi9ejOeff54Gr5OkokRJZDfW6uoejwdHjhzBY489hkAgAKPRiNmzZ+Omm27Ciy++iOPHj8t9CyTDiS7cm01o4d7UcPjwYbS2tqKvrw/5+fmora3FnDlzRhx35MgROJ1O2O12vP3225gyZQpuuOEGLFmyBIsWLcK0adNkiJ5kimj5gBIlKFGms08//RQOhwN2ux0ejwennXYaFi1ahKqqKtxwww04/fTT5Q6RpJlo+YAevUlau+CCC7Bq1Sq8++674DgO69evB8dxqKysxKxZs7B06VLs2rUL33//vdyhkjRGiZJkjIKCAlgsFnR1daGnpwdr167FJ598goqKCsyePRu33XYbdu/ejR9++CHiPJo6ScZCj96gR+9M193dLTyef/DBB5g+fTpuvPFGVFVV4fPPP0d9fX3E1Eme52nqZBajNkoRlCizx0cffQSHw4HOzs5RVzBSKBTo7u6O2pmUKD09PWhraxM6r2pqalBYWCjZ9UhsKFGK0Ol0yMvLQ3V1NQ06zyIrVqxAW1ub6BjOM844A3PnzhW2Phm6NYrYd9OmTROdZTTUtm3bsGLFCqrJppCOjg50dHTg6NGj6Orqiiij1YPwv5k5JLt89913QoKKRqVSoaCgAIFAAL29vQgEAggGgwgEAjhx4kTUc3JycsZMqD/++CPWr18f9bq1tbW45pprJK3JkujCFaXy8vIRZZQoSdbKz88Xrf3l5OTglltuwaZNm6KW//jjj0LSFPsz/P7YsWPw+XwRZWLJmed53HXXXbj//vsxb948zJo1K2H3S8aPHr1BbZTZqqenB3Pnzo26U6SUbZRLly6F3W4XTZYKhUKI6cwzz8Qll1yCSy65BPPmzRP+nDFjRsLjIoOi5QOqUZKsFZ46WVtbG7WtUKrH34KCAigUCpw8eXJEWU5ODu677z7ceeed+PDDD4XX3r178ec//1lIoOeff76QQMOvuXPn4rTTThvz+tSJFD+qUYJqlNku1qmTiTLemuwPP/yATz75JCKBfvjhh/j000+Fc+fMmTMigRYWFuKUUwbrRHJ3IqVDko6aD/gM5fP5eKvVyjscDt5qtfKBQED02MWLFycvMEJ4nt+2bRuvUCj4nJyciD+3bdsW92+FQiH+rbfe4ltaWvh77rmHLy0t5c8880weAA+Anzx5Mn/ZZZfxixcv5hmGEb4f+lIoFHxPT0/ib3SItra2hN2zlKLlg4ytUep0OqGLn+M4WCwWOByOqMdSjZLIQeqa7Ndff41Dhw4JNc+XX34Zn332mejx5557Li677DJMnz497tfUqVNHHRYlV3vw0OvHWpPNmjZKjuMiPqvVarjdbpmiISS6OXPmiPaqJ8KsWbMwf/58zJ8/H8Dg8JfOzk7RTqScnBxMmjQJX3/9NTiOw/HjxyNeo9WpGIbBtGnTRBPpRx99NOr5W7ZswUMPPQSlUokpU6bENBY1VtGaG+LduC4jE6Xb7YZKpYr4TqVSwev1QqvVyhQVIfIKryQv1olUXV0tmrh5nseJEydGJM9YXkeOHEFfX59oohwYGIDNZoPNZgMATJo0CUqlEjNmzIBSqRz1vVjZ9OnToVAo0NPTgxUrVkStycYzZjUjE6XYlgF+vz+5gRCSQmpqatDU1BS1jOd51NbWip7LMAymTp2KqVOn4qyzzor72mvWrMGWLVuiJmmFQoGlS5di6dKlCIVCCIVC6O/vH/G+t7c34nN/f3/UBBg2ffp0ABA9JrxxXSy1+oxMlGLEEmh4c7EwmspIMpFcw6GA0ZM0ADz00ENxX5/neXz33XdRk2v487Zt23Do0CHR8/v6+oSpi2HRNhfLyETJsuyI2qPf7wfLslGPpymMJFssX74c11xzTVKHQwHSJOlwu+i0adOibggGAN988w0++eSTqDVZhmGQn58/omIUbQpjRvZ6hxduHTqxPTc3F729vVGTJfV6E5Ic6TBmNWt6vdVqdcRnjuNQXFwsWqMkhCSH1D39wyWqJpuRiRIAHA4HLBYLSkpK4PF4RMdQEkIyW0KaG6Qe5Z4OUnlmzosvvih3CEmXbfecbffL86l9z9HyAe2Zk+KG9sZli2y752y7XyD97pkSZYwm8n+sXP9RyBWznH8Jsu2e5Yw5m+6ZEmWM6C9Rcs6dqGy7Z0qUyTk3I4cHxWvevHnQaDSjHnP06FHRsVpjoXPT49p0bnpcW+pzfT7fiEHqlCgJIWQM9OhNCCFjoERJCCFjoERJCCFjyNiZOenO6/UKiw17PB60tLRkzRRMi8WCNWvWZMX9ut1ucBwnTLvV6/UyRyQdjuOEtWI5joPBYBgx3ThlJXnQO4mR1WqNeK/VamWMJnm6urp4AKPucZQpXC4XbzKZeJ4f3ONJrVbLHJG0hv43zfO8cO/pgB69U5DX641YOMBgMMDr9Y7Y4iITDa1dZTqz2Qyr1QpgcCEXl8slc0TS6uzslDuEcaNEmYK0Wi1aWlqEz+EFh4dvb5FpnE4nDAaD3GEkBcdxwhqpXq8XwWAw4/+BUKlU0Ol0wiN4WVmZ3CHFjBJlihqaMDo7O6HX6zO6zS4YDGb0/Q3n9XqhUqngdDqhVqths9ngdDrlDktS4RW8NBoNHA5HWv2jSJ05KS4YDMLpdEYsQpyJ7HY7TCaT3GEkjd/vB8dxwj+AJpMJubm5o+5UmO7cbjesVis4joPZbAYANDc3yxxVbKhGmeIsFgtcLldG17bcbjeWLFkidxhJpVarwbKs8P9r+E+v1ytfUBLiOA4ejwd6vR4mkwk+nw92uz1t2t2pRpnCmpqaYLFYoFarhXbKTE2YdrtdeM9xHDZt2oSqqqqM3V4409sjh/N6vSgpKRE+q9VqrFmzRnTDv5Qjd7c7ic7hcPAul4vneZ4PBAJ8c3OzzBElDwDe5/PJHYbk9Hq9MAwq04cH+Xw+vqGhIeK74Z9TGS2KkYI4jhuxmhHLsggEAjJFlBzBYBA2mw0WiwUmkwlmszlja5TA4P1aLBbodDp0dXUJTw+Zyu12w+v1Ck9Fer0+be6XEiUhhIyBOnMIIWQMlCgJIWQMlCgJIWQMlCgJIWQMlCgJIWQMlCgJIWQMlCgJIWQMlChJyggPNg+vpON2u2Gz2ZI+H9jtdkOn08Fms03od3Q6XcavCJQtaK43SQllZWWorKwUFrIFBucHl5WVwefzJTUWvV6PqqqqCf+O1WpFcXFxxHfZtpxcpqAaJZFdU1MTAIxYZk2r1ab10mvD1xDlOC5i8Q+SPihREtlt2rRJWJ9wuMrKyiRHI52htWWSXihREllxHDfqNghDF07wer1oamqC0+lEU1OT0HY5tE0x3K5ZWVkJjuOEc4YmYrfbDY1GA7PZDKfTCZvNBrPZPOqSX263W7i2xWIBMLh1hU6nExYwaWpqgkajQVNTE7xeb0Q7p9vtxoEDB+ByuWCz2fDkk09Co9EI/xAEg0FoNBrht0lqoTZKkhY4jhMWMQ7T6XTYt28f9Ho99Ho9XC6XsN2Aw+GA0+lEQ0MDtFotNBqN0D6o1+thMBgwc+ZMYTsCp9OJysrKqBt8ha8dXmXe7/ejqakJDQ0N0Ov1KC0tFX67q6tLeNwe2s4ZjlGj0QjNCZMnTxZ+k2VZYdUkknqoRklkFa4tivVsh2uczc3NI5ZcU6vVQpvfzJkzIxaGZVk2opbKsiz8fn/E+UPbDw0GA9xud9RaZXNzM1QqFdxud8Re6+HfaGlpgU6nQ3FxcVwdNSaTSYif47gRHT8kdVCiJLJraGgQ3Ttl6PqF8Upk77JWqxVqhSaTSai5hq+j1Wrj2o41nJBNJhNsNhu8Xm9Gr72Z7ihREtmFOzmGj1scWrurqqoSanNhXq931L12xtpmYGi50+kc0UsdLo927fDnYDAIt9sNh8MBjuNGjJsUiyF8/tC9vUnqojZKkhJcLpewR5BGoxH2MA+3IWq1WlitVjQ1NUGtVsPj8cDhcAj7Yodrc3q9XujECZ/ndrvBcRysViusVquQDH0+n/C4Hf49ACN+L3xti8UiPN7r9XrYbDZYrVaho6ikpAR1dXXC7oqdnZ1QqVQwGAxQq9VCUrTZbNDr9QAGmw+0Wm1abd2ajWiFc5KVwgk5FTpPnE4nJcoURzVKQmRgNptRWVkJlUpFbZNpgNooSdYJ9147HA7Z9tGurKxEMBiE3+9Pmw22shk9ehNCyBioRkkIIWOgREkIIWOgREkIIWP4f4bjdHxqOKZMAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 350x262.5 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "with plt.style.context('science'):\n",
    "    fig, ax = plt.subplots()\n",
    "    ax.plot(complexities, rel_bic_scores, color='black', linewidth=2.5//2)\n",
    "    ax.scatter(complexities, rel_bic_scores, marker='o', c='black', s=25)\n",
    "    plt.xlim(min(complexities)-0.3,max(complexities)+0.3)\n",
    "    plt.ylabel(r\"Relative BIC ($\\times 10^{3}$)\")\n",
    "    plt.xlabel(\"Complexity\")\n",
    "    plt.tight_layout()\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1ddb549f-abde-4cee-a577-53348ad71b35",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:py3.7]",
   "language": "python",
   "name": "conda-env-py3.7-py"
  },
  "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.7.12"
  },
  "widgets": {
   "application/vnd.jupyter.widget-state+json": {
    "state": {},
    "version_major": 2,
    "version_minor": 0
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}