figure_3_alt_notebook_colab.ipynb

{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/github/pachterlab/GVFP_2021/blob/dev/figure_3_alt_notebook_colab.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "-U4Z5kp35WaY"
      },
      "source": [
        "# \"Interpretable and tractable models of transcriptional noise for the rational design of single-molecule quantification experiments\"\n",
        "## by Gennady Gorin<sup>*</sup>, John Vastola<sup>*</sup>, Meichen Fang, and Lior Pachter\n",
        "### <sup>*</sup> These authors contributed equally to this work.\n",
        "# Companion notebook for Figure 3 (alternate version)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "YQ9yXp055Wad"
      },
      "source": [
        "This notebook is self-contained, and generates the panels used in the alternative version Figure 3a presented in the SI (albeit without axis labels). "
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "JXFHzB7B5Wad"
      },
      "outputs": [],
      "source": [
        "# Scientific computing imports\n",
        "import numpy as np\n",
        "from numpy.fft import irfftn\n",
        "from numpy.random import choice\n",
        "from scipy.stats import rv_discrete, poisson, nbinom, gaussian_kde\n",
        "import matplotlib.pyplot as plt\n",
        "\n",
        "# PyMC3-related imports\n",
        "import pymc3 as pm\n",
        "import theano.tensor as tt\n",
        "\n",
        "# Plotting\n",
        "import matplotlib.pyplot as plt\n",
        "import arviz as az\n",
        "import matplotlib as mpl\n",
        "\n",
        "# Timer\n",
        "import time as ti"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "W-_Wk2PC5Waf",
        "outputId": "46537a2e-bd1e-4c4c-f722-93e19ac19174"
      },
      "outputs": [
        {
          "data": {
            "application/javascript": [
              "IPython.OutputArea.prototype._should_scroll = function(lines) {\n",
              "    return false;\n",
              "}\n"
            ],
            "text/plain": [
              "<IPython.core.display.Javascript object>"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        }
      ],
      "source": [
        "%%javascript\n",
        "IPython.OutputArea.prototype._should_scroll = function(lines) {\n",
        "    return false;\n",
        "}"
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "!mkdir -p results"
      ],
      "metadata": {
        "id": "6XHw0NrF5a3C"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "5XmUK5O-5Wag"
      },
      "source": [
        "------------------------------------------------"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "evZmWYb85Wag"
      },
      "source": [
        "## Functions"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "UKcGopKu5Wag"
      },
      "source": [
        "### General-purpose functions"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "tnZLSVS95Wah"
      },
      "outputs": [],
      "source": [
        "# Utility function that converts between (x, y) qualitative regime coordinates to original parameters.\n",
        "def convert_xy_to_params(x, y, beta, gamma, K_avg):\n",
        "    kappa = (beta + gamma)*(x/(1-x))\n",
        "    a_over_th = 1/y - 1\n",
        "    theta = np.sqrt(kappa*K_avg/a_over_th)\n",
        "    a = K_avg*kappa/theta\n",
        "    return a, kappa, theta\n",
        "\n",
        "\n",
        "# Sample from flattened 2D probability distribution p. \n",
        "# mx = [mx0, mx1] is the shape of the 2D domain we're sampling.\n",
        "def sample_from_p(mx, num_data_points, p_flat):\n",
        "    x_N, x_M = np.arange(mx[0]), np.arange(mx[1])\n",
        "    X_N, X_M = np.meshgrid(x_N, x_M, indexing='ij')    # get grid of possible values\n",
        "    x_choices = np.arange(np.prod(mx))                 # represent each grid point as a point in a 1D vector\n",
        "    \n",
        "    samples = choice(x_choices, size=num_data_points, replace=True, p=p_flat)\n",
        "    d_N, d_M = X_N.flatten()[samples], X_M.flatten()[samples]    \n",
        "    return d_N, d_M\n",
        "\n",
        "\n",
        "# Get maximum of a 2D array\n",
        "def get_2D_max(array):\n",
        "    return np.unravel_index(array.argmax(), array.shape)\n",
        "\n",
        "\n",
        "# Get KDE for smooth-looking heatmaps.\n",
        "def get_2D_KDE(x_stats, y_stats, x_min=0, x_max=1, y_min=0, y_max=1):\n",
        "    num_pts = 100                          # hyperparameter\n",
        "    x_arg = np.linspace(x_min, x_max, num_pts)\n",
        "    y_arg = np.linspace(y_min, y_max, num_pts)\n",
        "\n",
        "    X_arg, Y_arg = np.meshgrid(x_arg, y_arg, indexing='ij')      # grid of points      X0, X1\n",
        "\n",
        " \n",
        "    positions = np.vstack([X_arg.ravel(), Y_arg.ravel()])\n",
        "    kernel = gaussian_kde([x_stats, y_stats])\n",
        "    KDE = np.reshape(kernel(positions).T, X_arg.shape)\n",
        "    return KDE, x_arg, y_arg\n",
        "    "
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "ReAEZDuU5Wai"
      },
      "source": [
        "### Trivial model (constitutive/Poisson and mixture/NB) likelihood functions"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "9-BVbez-5Wai"
      },
      "outputs": [],
      "source": [
        "# Constitutive model likelihood function.\n",
        "def get_Poiss_2sp(mx, params):\n",
        "    beta, gamma, a, kappa, theta = params \n",
        "    K_avg = (a*theta)/kappa\n",
        "    mu_N, mu_M = K_avg/beta, K_avg/gamma\n",
        "    \n",
        "    x_N, x_M = np.arange(mx[0]), np.arange(mx[1])\n",
        "    X_N, X_M = np.meshgrid(x_N, x_M, indexing='ij')\n",
        "    \n",
        "    return poisson.pmf(X_N, mu_N)*poisson.pmf(X_M, mu_M)\n",
        "\n",
        "\n",
        "# Mixture model likelihood function.\n",
        "def get_NB_2sp(mx, params):\n",
        "    beta, gamma, a, kappa, theta = params \n",
        "    \n",
        "    # Get generating function argument\n",
        "    u = []\n",
        "    half = mx[:]\n",
        "    half[-1] = mx[-1]//2 + 1\n",
        "    for i in range(len(mx)):\n",
        "        l = np.arange(half[i])\n",
        "        u_ = np.exp(-2j*np.pi*l/mx[i])-1\n",
        "        u.append(u_)\n",
        "    g = np.meshgrid(*[u_ for u_ in u], indexing='ij')\n",
        "    for i in range(len(mx)):\n",
        "        g[i] = g[i].flatten()\n",
        "    \n",
        "    # Get generating function\n",
        "    gf = np.exp(- (a/kappa)*np.log(1 - theta*(g[0]/beta + g[1]/gamma)))\n",
        "    gf = gf.reshape(tuple(half))\n",
        "                              \n",
        "    Pss = irfftn(gf, s=mx)                        # Get Pss by inverse fast Fourier transform\n",
        "    Pss = np.abs(Pss)/np.sum(np.abs(Pss))           # Normalize\n",
        "    return Pss"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Sfbe7iav5Waj"
      },
      "source": [
        "### $\\Gamma$-OU likelihood functions"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "m1f-IW2d5Waj"
      },
      "outputs": [],
      "source": [
        "# Get 2 species GOU generating function using ODE method\n",
        "def get_gf_GOU_2sp_ODE(g0, g1, params):\n",
        "    beta_0, beta_1, a, kappa, theta = params     # get parameters\n",
        "    \n",
        "    c0 = (g0) + (beta_0/(beta_1 - beta_0))*(g1)       #  relevant linear combinations of g_i\n",
        "    c1 = - (beta_0/(beta_1 - beta_0))*(g1)   \n",
        "    \n",
        "    min_fudge, max_fudge = 0.25, 10                                     # Determine integration time scale / length\n",
        "    dt = np.min([1/kappa, 1/theta, 1/a, 1/beta_0, 1/beta_1])*min_fudge\n",
        "    t_max = np.max([1/kappa, 1/theta, 1/a, 1/beta_0, 1/beta_1])*max_fudge\n",
        "    num_tsteps = int(np.ceil(t_max/dt))\n",
        "    t_array = np.linspace(0, t_max, num_tsteps+1)\n",
        "    t_array = t_array.reshape((1, num_tsteps + 1))\n",
        "    \n",
        "    q = np.zeros((g0.shape[0], num_tsteps + 1), dtype=np.complex64)    # initialize array to store ODE\n",
        "    c0 = c0.reshape((c0.shape[0],1))\n",
        "    c1 = c1.reshape((c1.shape[0],1))\n",
        "    q0 = theta*c0*(np.exp(-beta_0*t_array) - np.exp(-kappa*t_array))/(kappa - beta_0) \n",
        "    q1 = theta*c1*(np.exp(-beta_1*t_array) - np.exp(-kappa*t_array))/(kappa - beta_1)\n",
        "    q = q0 + q1\n",
        "    \n",
        "\n",
        "    integrand = q/(1-q)\n",
        "    integral = np.trapz(integrand, dx=dt, axis=1)     # integrate ODE solution\n",
        "    gf = np.exp( a*integral)               # get generating function\n",
        "    return gf\n",
        "\n",
        "\n",
        "# Get Pss for 2 species GOU model via ODE method\n",
        "def get_GOU_2sp(mx, params):\n",
        "    # Get generating function argument\n",
        "    u = []\n",
        "    half = mx[:]\n",
        "    half[-1] = mx[-1]//2 + 1\n",
        "    for i in range(len(mx)):\n",
        "        l = np.arange(half[i])\n",
        "        u_ = np.exp(-2j*np.pi*l/mx[i])-1\n",
        "        u.append(u_)\n",
        "    g = np.meshgrid(*[u_ for u_ in u], indexing='ij')\n",
        "    for i in range(len(mx)):\n",
        "        g[i] = g[i].flatten()\n",
        "    \n",
        "    # Get generating function\n",
        "    gf = get_gf_GOU_2sp_ODE(g[0], g[1], params)                    \n",
        "    gf = gf.reshape(tuple(half))\n",
        "    \n",
        "    Pss = irfftn(gf, s=mx)                        # Get Pss by inverse fast Fourier transform\n",
        "    Pss = np.abs(Pss)/np.sum(np.abs(Pss))           # Normalize\n",
        "    return Pss\n",
        "\n",
        "\n",
        "# Log likelihood of GOU model given data. Uses (x,y) as input instead of (kappa, theta).\n",
        "def ll_GOU(phi, const, mx, data):\n",
        "    # Get parameters\n",
        "    x, y = phi\n",
        "    beta, gamma, K_avg = const\n",
        "    \n",
        "    # Convert from (x, y) to original parameters\n",
        "    a, kappa, theta = convert_xy_to_params(x, y, beta, gamma, K_avg)\n",
        "    params = [beta, gamma, a, kappa, theta]\n",
        "    \n",
        "    Pss = get_GOU_2sp(mx, params)    # Compute Pss\n",
        "\n",
        "    lp = np.log(Pss)\n",
        "    result = np.sum(lp[tuple(data)])\n",
        "    return result\n",
        "\n",
        "\n",
        "# Combines the above functions into one to reduce overhead associated with Python function calls.\n",
        "# Helpful when doing expensive posterior sampling (since many likelihood function calls are required).\n",
        "def ll_GOU2(phi, const, mx, data):\n",
        "    # Get parameters\n",
        "    x, y = phi\n",
        "    beta, gamma, K_avg = const\n",
        "    \n",
        "    \n",
        "    # Convert from (x, y) to original parameters\n",
        "    a, kappa, theta = convert_xy_to_params(x, y, beta, gamma, K_avg)\n",
        "    params = [beta, gamma, a, kappa, theta]\n",
        "    \n",
        "    \n",
        "    # Get generating function argument\n",
        "    u = []\n",
        "    half = mx[:]\n",
        "    half[-1] = mx[-1]//2 + 1\n",
        "    for i in range(len(mx)):\n",
        "        l = np.arange(half[i])\n",
        "        u_ = np.exp(-2j*np.pi*l/mx[i])-1\n",
        "        u.append(u_)\n",
        "    g = np.meshgrid(*[u_ for u_ in u], indexing='ij')\n",
        "    for i in range(len(mx)):\n",
        "        g[i] = g[i].flatten()\n",
        "    \n",
        "    # Get generating function\n",
        "    beta_0 = beta\n",
        "    beta_1 = gamma\n",
        "    \n",
        "    c0 = (g[0]) + (beta_0/(beta_1 - beta_0))*(g[1])       #  relevant linear combinations of g_i\n",
        "    c1 = - (beta_0/(beta_1 - beta_0))*(g[1])   \n",
        "    \n",
        "    min_fudge, max_fudge = 0.5, 10                                     # Determine integration time scale / length\n",
        "    dt = np.min([1/kappa, 1/theta, 1/beta_0, 1/beta_1])*min_fudge\n",
        "    t_max = np.max([1/kappa, 1/theta, 1/beta_0, 1/beta_1])*max_fudge\n",
        "    num_tsteps = int(np.ceil(t_max/dt))\n",
        "    \n",
        "    t_array = np.linspace(0, t_max, num_tsteps+1)\n",
        "    t_array = t_array.reshape((1, num_tsteps + 1))\n",
        "    \n",
        "    q = np.zeros((c0.shape[0], num_tsteps + 1), dtype=np.complex64)    # initialize array to store ODE\n",
        "    c0 = c0.reshape((c0.shape[0],1))\n",
        "    c1 = c1.reshape((c1.shape[0],1))\n",
        "    q0 = theta*c0*(np.exp(-beta_0*t_array) - np.exp(-kappa*t_array))/(kappa - beta_0) \n",
        "    q1 = theta*c1*(np.exp(-beta_1*t_array) - np.exp(-kappa*t_array))/(kappa - beta_1)\n",
        "    q = q0 + q1\n",
        "\n",
        "    integrand = q/(1-q)\n",
        "    integral = np.trapz(integrand, dx=dt, axis=1)     # integrate ODE solution\n",
        "    gf = np.exp( a*integral)               # get generating function\n",
        "    \n",
        "    gf = gf.reshape(tuple(half))\n",
        "    \n",
        "    Pss = irfftn(gf, s=mx)                        # Get Pss by inverse fast Fourier transform\n",
        "    Pss = np.abs(Pss)/np.sum(np.abs(Pss))           # Normalize\n",
        "\n",
        "    lp = np.log(Pss)\n",
        "    result = np.sum(lp[tuple(data)])\n",
        "    return result"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "TdUuF7965Wak"
      },
      "source": [
        "### CIR likelihood functions"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "nybLIuPt5Wak"
      },
      "outputs": [],
      "source": [
        "# ODE for 2 species CIR model\n",
        "def f_2sp(q, t, c0, c1, params):\n",
        "    beta_0, beta_1, a, kappa, theta = params\n",
        "    result = - kappa*q + theta*q*q + kappa*( c0*np.exp(-beta_0*t) + c1*np.exp(-beta_1*t)  ) \n",
        "    return result\n",
        "\n",
        "\n",
        "# Vectorized RK4 implementation for 2 species CIR model\n",
        "def RK4_2sp(q, f, t, c0, c1, step_size, param):\n",
        "    j1 = f(q, t, c0, c1, param)\n",
        "    j2 = f(q + (step_size/2)*j1, t + (step_size/2), c0, c1, param)   \n",
        "    j3 = f(q + (step_size/2)*j2, t + (step_size/2), c0, c1, param)   \n",
        "    j4 = f(q + (step_size)*j3, t + (step_size), c0, c1, param)  \n",
        "    \n",
        "    q_new = q + (step_size/6)*(j1 + 2*j2 + 2*j3 + j4)\n",
        "    return q_new\n",
        "\n",
        "\n",
        "# Get 2 species CIR generating function using ODE method\n",
        "def get_gf_CIR_2sp(g0, g1, params):\n",
        "    beta_0, beta_1, a, kappa, theta = params     # get parameters\n",
        "    \n",
        "    c0 = (g0) + (beta_0/(beta_1 - beta_0))*(g1)       #  relevant linear combinations of g_i\n",
        "    c1 = - (beta_0/(beta_1 - beta_0))*(g1)   \n",
        "    \n",
        "    min_fudge, max_fudge = 0.25, 10                                     # Determine integration time scale / length\n",
        "    dt = np.min([1/kappa, 1/theta, 1/a, 1/beta_0, 1/beta_1])*min_fudge\n",
        "    t_max = np.max([1/kappa, 1/theta, 1/a, 1/beta_0, 1/beta_1])*max_fudge\n",
        "    num_tsteps = int(np.ceil(t_max/dt))\n",
        "     \n",
        "    q = np.zeros((g0.shape[0], num_tsteps + 1), dtype=np.complex64)    # initialize array to store ODE\n",
        "    \n",
        "    # Solve ODE using RK4 method \n",
        "    for i in range(0, num_tsteps):\n",
        "        t = i*dt\n",
        "        q[:,i+1] = RK4_2sp(q[:,i], f_2sp, t, c0, c1, dt, params)\n",
        "        \n",
        "    integral = np.trapz(q, dx=dt, axis=1)     # integrate ODE solution\n",
        "    gf = np.exp((a*theta/kappa)*integral)               # get generating function\n",
        "    return gf\n",
        "\n",
        "\n",
        "# Get Pss for 2 species CIR model via ODE method\n",
        "def get_CIR_2sp(mx, params):\n",
        "    # Get generating function argument\n",
        "    u = []\n",
        "    half = mx[:]\n",
        "    half[-1] = mx[-1]//2 + 1\n",
        "    for i in range(len(mx)):\n",
        "        l = np.arange(half[i])\n",
        "        u_ = np.exp(-2j*np.pi*l/mx[i])-1\n",
        "        u.append(u_)\n",
        "    g = np.meshgrid(*[u_ for u_ in u], indexing='ij')\n",
        "    for i in range(len(mx)):\n",
        "        g[i] = g[i].flatten()\n",
        "    \n",
        "    # Get generating function\n",
        "    gf = get_gf_CIR_2sp(g[0], g[1], params)                    \n",
        "    gf = gf.reshape(tuple(half))\n",
        "                              \n",
        "    Pss = irfftn(gf, s=mx)                        # Get Pss by inverse fast Fourier transform\n",
        "    Pss = np.abs(Pss)/np.sum(np.abs(Pss))           # Normalize\n",
        "    return Pss\n",
        "\n",
        "\n",
        "# Log likelihood of CIR model given data. Uses (x,y) as input instead of (kappa, theta).\n",
        "def ll_CIR(phi, const, mx, data):\n",
        "    # Get parameters\n",
        "    x, y = phi\n",
        "    beta, gamma, K_avg = const\n",
        "    \n",
        "    # Convert from (x, y) to original parameters\n",
        "    a, kappa, theta = convert_xy_to_params(x, y, beta, gamma, K_avg)\n",
        "    params = [beta, gamma, a, kappa, theta]\n",
        "    \n",
        "    Pss = get_CIR_2sp(mx, params)    # Compute Pss\n",
        "\n",
        "    lp = np.log(Pss)\n",
        "    result = np.sum(lp[data])\n",
        "    return result"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "NnMSyjJN5Wal"
      },
      "source": [
        "### Bayes factor functions"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "9rpg-6z_5Wal"
      },
      "outputs": [],
      "source": [
        "# Generates synthetic data (either CIR or GOU) and computes log Bayes factor (e.g. P(CIR)/P(GOU))).\n",
        "# Does this averaged over many trials (num_trials).\n",
        "# Does this assuming different numbers of data points (e.g. num_data_points = [100, 200, 1000]). \n",
        "def log_bayes_factor_avg(model, num_data_points, num_trials, mx, params):\n",
        "    \n",
        "    # Initialize log BF arrays \n",
        "    num_kinds = len(num_data_points) \n",
        "    log_bf_joint = np.zeros((num_kinds, num_trials))\n",
        "    log_bf_nascent = np.zeros((num_kinds, num_trials))\n",
        "    log_bf_mature = np.zeros((num_kinds, num_trials))\n",
        "    \n",
        "    # Get Pss for each model given this parameter set\n",
        "    if model=='CIR':\n",
        "        pss = get_CIR_2sp(mx, params)\n",
        "        pss_other = get_GOU_2sp(mx, params)\n",
        "    elif model=='GOU':\n",
        "        pss = get_GOU_2sp(mx, params)\n",
        "        pss_other = get_CIR_2sp(mx, params)\n",
        "    p_flat = pss.flatten()\n",
        "    \n",
        "    \n",
        "    \n",
        "\n",
        "    # Log-likelihood functions\n",
        "    lp = np.log10(pss)    # shape: x_N_domain by x_M_domain\n",
        "    lp_nascent = np.log10(np.sum(pss, axis=1))\n",
        "    lp_mature = np.log10(np.sum(pss, axis=0))\n",
        "\n",
        "    lp_other = np.log10(pss_other)\n",
        "    lp_other_nascent = np.log10(np.sum(pss_other, axis=1))\n",
        "    lp_other_mature = np.log10(np.sum(pss_other, axis=0))\n",
        "    \n",
        "    \n",
        "\n",
        "    \n",
        "    # Get fake data\n",
        "    for r in range(0, num_kinds):\n",
        "            \n",
        "        # Sample nascent and mature counts\n",
        "        d_N, d_M = sample_from_p(mx, (num_data_points[r], num_trials), p_flat)\n",
        "\n",
        "        # Log-likelihoods given data\n",
        "        ll_joint = np.sum(lp[d_N, d_M], axis=0)\n",
        "        ll_nascent = np.sum(lp_nascent[d_N], axis=0)\n",
        "        ll_mature = np.sum(lp_mature[d_M], axis=0)\n",
        "\n",
        "        ll_other_joint = np.sum(lp_other[d_N, d_M], axis=0)\n",
        "        ll_other_nascent = np.sum(lp_other_nascent[d_N], axis=0) \n",
        "        ll_other_mature = np.sum(lp_other_mature[d_M], axis=0)\n",
        "\n",
        "        # Bayes factor\n",
        "        log_bf_joint[r,:] = ll_joint - ll_other_joint\n",
        "        log_bf_nascent[r,:] = ll_nascent - ll_other_nascent\n",
        "        log_bf_mature[r,:] = ll_mature - ll_other_mature\n",
        "\n",
        "                          \n",
        "    log_bf_joint_avg = np.mean(log_bf_joint, axis=1)\n",
        "    log_bf_nascent_avg = np.mean(log_bf_nascent, axis=1)\n",
        "    log_bf_mature_avg = np.mean(log_bf_mature, axis=1)\n",
        "    \n",
        "    return log_bf_joint_avg, log_bf_nascent_avg, log_bf_mature_avg"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "QLT682Ij5Wal"
      },
      "source": [
        "------------------------------------------"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "M54VXOLq5Wam"
      },
      "source": [
        "## Log Bayes factors across space of qualitative regimes"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "oDIS2TqY5Wam"
      },
      "source": [
        "Define a function that does everything. Will run this multiple times."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "scrolled": true,
        "id": "R5-WvgWn5Wam"
      },
      "outputs": [],
      "source": [
        "def get_BF_plot(name, model, K_avg, beta, gamma, mx):\n",
        "\n",
        "\n",
        "    # Hyperparameters\n",
        "    num_data_points = np.array([100, 200, 1000, 5000])     # do everything assuming we have _ cells (e.g. 100)\n",
        "    num_kinds = len(num_data_points)\n",
        "    num_trials = 100    # for a given parameter set, how many fake data sets to average over\n",
        "\n",
        "\n",
        "    # Set up qualitative regimes grid\n",
        "    x_min, x_max, num_x = 0.1, 0.9, 10\n",
        "    y_min, y_max, num_y = 0.1, 0.9, 10\n",
        "    x, y = np.linspace(x_min, x_max, num_x), np.linspace(y_min, y_max, num_y)\n",
        "    X, Y = np.meshgrid(x, y, indexing='ij')\n",
        "\n",
        "    a, kappa, theta = convert_xy_to_params(X, Y, beta, gamma, K_avg)   # values of other params at each grid point\n",
        "\n",
        "\n",
        "    # Output arrays\n",
        "    log_bf_joint = np.zeros((num_kinds, num_x, num_y))\n",
        "    log_bf_nascent = np.zeros((num_kinds, num_x, num_y))\n",
        "    log_bf_mature = np.zeros((num_kinds, num_x, num_y))\n",
        "    \n",
        "    \n",
        "    # Compute log Bayes factors\n",
        "    start_time = ti.time()\n",
        "    print(\"Computing log Bayes factors on a \", num_x, \" x \", num_y, \" grid...\")\n",
        "    for i in range(0, num_x):\n",
        "        print(\"Row \", i+1, \"/\", num_x, \"...\")\n",
        "        for j in range(0, num_y):\n",
        "            params = [beta, gamma, a[i,j], kappa[i,j], theta[i,j]]   # parameters for this grid point\n",
        "\n",
        "            log_bf_joint[:,i, j], log_bf_nascent[:,i, j], log_bf_mature[:,i, j] = (\n",
        "                                log_bayes_factor_avg(model, num_data_points, num_trials, mx, params) )\n",
        "\n",
        "\n",
        "    print(\"--- %s seconds ---\" % (ti.time() - start_time))\n",
        "    \n",
        "    \n",
        "    \n",
        "    # Plot results\n",
        "    \n",
        "    # Clip arrays before plotting. Lower bound: 0 (no evidence either way), upper bound: 2 (decisive evidence).\n",
        "    bf_min, bf_max = 0, 2\n",
        "\n",
        "    for r in range(0, num_kinds):\n",
        "        gs_kw = dict(hspace=-0.37, wspace=0, width_ratios=[1,1,1])\n",
        "        fig, axs = plt.subplots(1,3, sharey=True, sharex=True, squeeze=True, gridspec_kw=gs_kw)\n",
        "\n",
        "        (ax1, ax2, ax3) = axs\n",
        "\n",
        "        ax1.imshow(np.transpose(log_bf_joint[r,:,:]), origin='lower', extent=[x_min, x_max, y_min, y_max], vmin = bf_min, vmax = bf_max)\n",
        "        ax2.imshow(np.transpose(log_bf_nascent[r,:,:]), origin='lower', extent=[x_min, x_max, y_min, y_max], vmin = bf_min, vmax = bf_max)\n",
        "        ax3.imshow(np.transpose(log_bf_mature[r,:,:]), origin='lower', extent=[x_min, x_max, y_min, y_max], vmin = bf_min, vmax = bf_max)\n",
        "\n",
        "        plt.savefig('results/'+name+str(num_data_points[r])+'c.png', bbox_inches='tight')\n",
        "        plt.savefig('results/'+name+str(num_data_points[r])+'c.pdf', bbox_inches='tight')\n",
        "        plt.show()\n",
        "\n",
        "\n",
        "    #####\n",
        "    \n",
        "    \n",
        "    return"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "TKsAxlU-5Wam"
      },
      "source": [
        "----------------------------------------------------"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "Ce_1Eodj5Wan",
        "outputId": "2cd3c92e-82bf-40b4-d564-348feb278780"
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Computing log Bayes factors on a  10  x  10  grid...\n",
            "Row  1 / 10 ...\n",
            "Row  2 / 10 ...\n",
            "Row  3 / 10 ...\n",
            "Row  4 / 10 ...\n",
            "Row  5 / 10 ...\n",
            "Row  6 / 10 ...\n",
            "Row  7 / 10 ...\n",
            "Row  8 / 10 ...\n",
            "Row  9 / 10 ...\n",
            "Row  10 / 10 ...\n",
            "--- 10.5771005153656 seconds ---\n"
          ]
        },
        {
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAACPCAYAAADqQSA9AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAP1klEQVR4nO3dbYhc93XH8e9vZne1q0dLlh9AK7uSo6ZRik3bVCWuS52ENLLBiEKhTkLbBAchqPvwoiGGQt7kTdtQKLQOQjUmFEL1JmnqpnIcp08pVezICbYcOZG9lh1ppTh6tB5X0u7M6Ysd1ePtSnv+0oxm9ur3gcUzd87cObr/8dk7d+f8/4oIzMxs/qv1OgEzM+sMF3Qzs4pwQTczqwgXdDOzinBBNzOrCBd0M7OKGOjVCw8OLYrhkeWpWN0+md7vbUOnivIY0VQ6dkj1dOzFaKRjTzaH07FvnV2WjgWon1U+9mL+K6wXTh9nZHBpKnZy6WB6v438oZhWz+c8MJAfk6nJ/FhT8M3f2oX8eAAMTOR3rkYu9vy544wMLEkmkC8RjeGCY1aq4Bir4KvY2WMGoMn8+weARj4+Gs107Ok4fjQibpntsZ4V9OGR5fzSfX+ci/3sofR+//SObxflsX7oWDr2joHF6dj9U2fSsTvOvDcd+1e7PpaOBVi+c0E6dvGh/Btw9zN/zQfv/INU7M8+fFt6v2+/r6wvork0/wt55a35X/ZH38r9sgJgKv9Bd/FY2f9yN+/Jn8wsOHEhFfe9F/6We2/53VRs85ab0q9/el3ZyQYFv9tUUEtrk/n30IITF9OxAz87mU8C4EQ+vnnmbDr22fNf+cnlHvMlFzOzinBBNzOrCBd0M7OKcEE3M6sIF3Qzs4pwQTczqwgXdDOzinBBNzOriFRBl7RR0l5JY5Iem+XxZZL+RdJLkvZI+nTnUzUzsyuZs21NUh14HPgoMA7skvRURLzSFvaHwCsR8ZCkW4C9kr4SEZdtw6pNTLJo98FUksf//o5UHMCf3flIOhZg4hcn0rFDr42kYwfyuyUKPieN/jjfGQkwfCSfyOB4vmuWySnirSOp0Nv+I98SuPjQzfkcgLO356cVWPjWTenYnyvoNiTybdsLDhd2Gxa8N2r7D6fiYmqKxpGjqdh6LZ/AsiNvp2MBGqtWpmPrb53I71gFLai1fGwUdHMCNE/lu8WjYJqAK8mM1gZgLCL2tQr0dmDTzHyAJZIELAaOA2WVx8zMrkmmoK8CDrTdH29ta/d3wPuAQ8DLwJ9EFJy2mJnZNcvMFDTbZ5KZn0c/BrwIfBi4C3hW0n9HxLtmQ5K0GdgMMFxPzvZmfePA6d0cOP0yABfjfI+zsRIHJl5hfGL6Kulk5CbxsvknU9DHgdVt90eZPhNv92ngLyIigDFJbwC/AHyvPSgitgHbAJYN3VY2rZ713Oold7N6yd0A7Nz/Dz3OxkqsHlnP6pH1AOw8sr3H2Vi3ZC657ALWSVojaQh4GHhqRsx+4CMAkm4D3gvs62SiZmZ2ZXOeoUfElKRHgWeAOvBkROyRtKX1+FbgC8CXJb3M9CWaz0VE7s/oZmbWEanZ9iNiB7BjxratbbcPAb/V2dTMzKyEO0XNzCrCBd3MrCJc0M3MKqJni0TTaNA8dToVumws33KrWFSUxvDR/DLzI8fz7bnNggXQ6yWL2h7LLxoMMHAkd4wB4uy5fGyzSfNMrrV54Gy+52DhG/mFnAFGfpp/C5es2j65Ij/NQ30i3xRdunJ87djb6djmuez4CZQ7l2sczU8HUVuSX0QdoPbmT/PBA/lxjgv5hZ9jIj81RvNC2ff3NZCfloLm9Wv9NzOzecAF3cysIlzQzcwqwgXdzKwiXNDNzCrCBd3MrCJc0M3MKsIF3cysIlzQzcwqwgXdzKwiXNDNzCqid3O5ADRz60hrIj9/yZI38/ORAEQ9/zutOZSPHTiVn/ehZH6P2un83BMATOTX/myeyc+ZUyIK5sDQkeNF+64N5ufLiCUL07FDB0+mYzWZn8sl3i6bq6ZkrfW4WDbPT0qzYKXIqfxxANLzyQA0J3LzBk3Hdme926K5WYCYzM8p0yk+QzczqwgXdDOzinBBNzOriFRBl7RR0l5JY5Ieu0zM/ZJelLRH0n91Nk0zM5vLnH8UlVQHHgc+CowDuyQ9FRGvtMXcBHwJ2BgR+yXd2qV8zczsMjJn6BuAsYjYFxEXge3AphkxnwC+FhH7ASLicGfTNDOzuWQK+irgQNv98da2dj8PLJf0n5K+L+n3O5WgmZnlZL6Hrlm2zfxy6gDwK8BHgBHgu5Kei4hX37UjaTOwGWBYZWt/Wu+NT73G+NRrAExStr6i9dZ4Y4zx5hgAk9Gd72lb72UK+jiwuu3+KHBolpijEXEWOCvpO8A9wLsKekRsA7YBLKuvLOhYsH4wOrCO0YF1AHx34l97nI2VGK2/h9H6ewB4bvKbPc7GuiVzyWUXsE7SGklDwMPAUzNi/hn4DUkDkhYCvwb8qLOpmpnZlcx5hh4RU5IeBZ4B6sCTEbFH0pbW41sj4keSvgnsBprAExHxwyvut9lMt+jWC9rBa7UV6VgAnSu4dLBgKB0ag/V0bO3wiXwOC0fysUCcL/h4XStoSxConvs3xrn8dAVaVHYprnky30pfK5jmgYI2+uap0+lYjQzncwCi5P1Zm+3q6Gw7jXRbugYL3vOlUw8k3z9QNn1EreAYl+yXgmkYAFByPACiMxcsUnO5RMQOYMeMbVtn3P8i8MWOZGVmZsXcKWpmVhEu6GZmFeGCbmZWES7oZmYV4YJuZlYRLuhmZhXhgm5mVhEu6GZmFeGCbmZWEalO0a6Q0GDnX752PN+GDZStVD5R0CZ8Or9KOStuSofG2XP5/VLWjh0Fq9cTEI1GKlTNfMt081S+lR/y0w8ARME0AV3L4ULZSvAlY5Idj6LXL9hn6esXHbdmvjU+zp4tSKKgPb9Uh9r5S/gM3cysIlzQzcwqwgXdzKwiXNDNzCrCBd3MrCJc0M3MKsIF3cysIlzQzcwqwgXdzKwiUgVd0kZJeyWNSXrsCnG/Kqkh6Xc6l6KZmWXM2XsvqQ48DnwUGAd2SXoqIl6ZJe4vgWdSrxyRbmuOkjb6hSP5WECDg+nYxtFj6dj68pvSsXEyP11BTEykYwGikW+7L24dT7Y2N4tWbC8bPwqmFWi8fTIdW1+6NB1b1B5fMBXD9BMKVpovic22vDcL3hPFbfT51v+iPGoF+y1RkkOPZM7QNwBjEbEvIi4C24FNs8T9EfBV4HAH8zMzs6RMQV8FHGi7P97a9n8krQJ+G9jaudTMzKxEZrrD2T5Hzfys/TfA5yKioSt87JK0GdgMMMzCZIrWL8abr3MwXgdgkoKZJ63nxmMfB2Mf4LGrskxBHwdWt90fBQ7NiPkAsL1VzFcCD0qaioivtwdFxDZgG8BSrbj+c0vaNRmt3cUodwHwfONbPc7GSoxqLaNaC8DzzW/3OBvrlkxB3wWsk7QGOAg8DHyiPSAi1ly6LenLwDdmFnMzM+uuOQt6RExJepTpb6/UgScjYo+kLa3Hfd3czKwPpJYMiogdwI4Z22Yt5BHxqWtPy8zMSrlT1MysIlzQzcwqwgXdzKwiUtfQuybbqlyyanxBe/d0CvlvT9YXL0rHNo6dSMdqeEE6tqSVHyAm86vMa6Dg7dAg3+qt/HlD83zhd6QL2rG1IH+cG6fz0zHUCvYbU4Wt/yUKjnN6RfqSNvrC1viiY1E8rUDSPGjnL+EzdDOzinBBNzOrCBd0M7OKcEE3M6sIF3Qzs4pwQTczqwgXdDOzinBBNzOrCBd0M7OKcEE3M6sIF3Qzs4ro7VwuSTE11ZVYKJvfo3mhYJ6RWn7uiWbJvCHDw/kcACI/xKXHLj0fCAXzz6T3WS4u5ue1KdE8fz4dq8Ghon1Ho2CukW7MS1Kyz9L5VkrGumTf2TmiSnXz39chPkM3M6sIF3Qzs4pwQTczq4hUQZe0UdJeSWOSHpvl8U9K2t362Snpns6namZmVzJnQZdUBx4HHgDWAx+XtH5G2BvAb0bE3cAXgG2dTtTMzK4sc4a+ARiLiH0RcRHYDmxqD4iInRFxaYme54DRzqZpZmZzyRT0VcCBtvvjrW2X8wjw9LUkZWZm5TJfUp7ty5ezfsFS0oeYLuj3XebxzcBmgGEWJlO0fjEe+zjIGwBMUrj2p/WUx+7GkCno48DqtvujwKGZQZLuBp4AHoiIY7PtKCK20bq+vlQrrv+37u2ajGoto6wF4Pn4tx5nYyU8djeGzCWXXcA6SWskDQEPA0+1B0i6A/ga8HsR8Wrn0zQzs7nMeYYeEVOSHgWeAerAkxGxR9KW1uNbgc8DNwNf0nR77FREfGDOV0+2xha1pBe25xa1g3erlbcg56LpB6An7cd9mQP0xfjFZHemH+gL3RznfngP9UMOc0hN9BERO4AdM7Ztbbv9GeAznU3NzMxKuFPUzKwiXNDNzCrCBd3MrCJc0M3MKsIF3cysIlzQzcwqwgXdzKwiXNDNzCrCBd3MrCIUPWpnlXQE+EmHd7sSONrhfV6tqufyy8APOrzPqh+zq+GxK1f1XO6MiFtme6BnBb0bJL2QmkPmOnAu5fopz37JpV/ymEs/5Xkj5+JLLmZmFeGCbmZWEVUr6P20OLVzKddPefZLLv2Sx1z6Kc8bNpdKXUM3M7uRVe0M3czshjVvCrqkjZL2ShqT9Ngsj39S0u7Wz05J97Q99qaklyW9KOmFLudxv6STrdd6UdLns8/tQi6fbcvjh5Iakla0HuvYMelAntdl7JK5ePzK8/T49cv4RUTf/zC99N3rwFpgCHgJWD8j5l5geev2A8DzbY+9Cay8TnncD3zjap7b6VxmxD8E/Hunj8l8GTuPn8fvRhi/+XKGvgEYi4h9EXER2A5sag+IiJ0RcaJ19zlgtBd5dOm5ndjfx4F/vIbXu1r9MnapXLr03E7sz+Pn8ZvTfCnoq4ADbffHW9su5xHg6bb7AXxL0vclbb4OeXxQ0kuSnpb0/sLndjoXJC0ENgJfbdvcqWMyl34Zu5JcPH7v8PhdWy7XdfxSi0T3gdmWVZ/16zmSPsT0m+q+ts2/HhGHJN0KPCvpxxHxnS7l8QOmW3PPSHoQ+DqwLvncTudyyUPA/0TE8bZtnTomc+mXscvm4vF7N4/f1edyyXUbv/lyhj4OrG67Pwocmhkk6W7gCWBTRBy7tD0iDrX+exj4J6Y/LnUlj4g4FRFnWrd3AIOSVmb/DZ3Mpc3DzPi418FjMpd+GbtULh6//8fjd5W5tLl+49fpi/Ld+GH6k8Q+YA3v/AHi/TNi7gDGgHtnbF8ELGm7vRPY2MU8bued7/dvAPYz/dt8zud2OpdW3DLgOLCoG8dkvoydx8/jdyOMX8ffAN36AR4EXmX6L8t/3tq2BdjSuv0EcAJ4sfXzQmv72tbBfgnYc+m5Xczj0dbrvMT0H4juvdJzu5lL6/6ngO0zntfRYzJfxs7j5/Gr+vi5U9TMrCLmyzV0MzObgwu6mVlFuKCbmVWEC7qZWUW4oJuZVYQLuplZRbigm5lVhAu6mVlF/C+g+2uqaKpTEQAAAABJRU5ErkJggg==\n",
            "text/plain": [
              "<Figure size 432x288 with 3 Axes>"
            ]
          },
          "metadata": {
            "needs_background": "light"
          },
          "output_type": "display_data"
        },
        {
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAACPCAYAAADqQSA9AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAP8UlEQVR4nO3db4xc5XXH8e+Z2V3v+u9iMInqXRobaIppoUrAVShVoYjGICGrUl84iZomAlluS9SqahVLVSNVedMKteqLGlkWQumLqE7V/KmbmrhJ/5HWQAwEA7YxWpYYD3aw1//Wa7ze3ZnTFzsuw3S9ex77jmf28e8jjTxz75k7Z+9z9/jO3fs8j7k7IiIy/5XanYCIiBRDBV1EJBMq6CIimVBBFxHJhAq6iEgmVNBFRDLR1a4PLi9c5N39y0Oxv7DieIuzydfhyYXh2NOji8Kx/t5Jehf0h2KrP1MLb/fm3pFwLECPlZPio857POfRWm84dmwqHgvw/rkF4diu8Vjc5KkT9HUvC8VOLYqXiGraj0bKDdPliYTtJhwS5fGELMzisUBpIn4M2YX4Dzg6NTLi7itmWte2gt7dv5yPPfZHodgf/d6TLc4mX3949K5w7L/84O5w7PE//yvW/tLvhmLP/tlYeLvfuP3vwrEAN3UtToqPOjDxfjh217k14dj/PnlLUh4/3hOPX/56rOAMff0JPnXLo6HYE5+InXQBnLo9HAqAJ9THxYfiwZNL4ttdfnAqHFvrSivoi96JH0Plt94Nx+4a2XboUut0yUVEJBMq6CIimVBBFxHJhAq6iEgmVNBFRDKhgi4ikgkVdBGRTKigi4hkIlTQzWydmR00syEz2zzD+mVm9s9mttfM9pnZF4tPVUREZjNnT1EzKwNbgAeBCrDHzHa4+/6GsN8H9rv7I2a2AjhoZl9390v2Z+05W2PlD8+Hkvz4L34+FAfw15/8h3AsQDfVcOzdC84kbTtq89EHwrHPHkrradi1J95t7ub/OBuOHZmYovtQbEiG8lOD4e3ef+8fh2MBvCveddt74l2xu5fGu2JX3+sLxy4YSftSPLA33pNx8evvheKGJqfgyLFQ7JIV8eEg+kbShmEoX4i3hyf00uw5GRwDAfBSfLtdx0fDsQB+Oh5fOxv/3ZtN5OhaCwy5+3C9QG8H1jfFOLDEzAxYDJwE4keiiIhcsUhBXwkcbnhdqS9r9LfAbcAR4DXgD9wTRjcSEZErFhmca6bvJM3fcz8NvAL8OnAz8H0z+6G7f+g7h5ltBDYCLFgQG+1NOkfl2ItUjr0EwGQtdrlMOsPh0Vc5PPoqABMevyQh80ukoFeAxougA0yfiTf6IvAX7u7AkJm9Dfw88KPGIHffBmwDWLpkZcromdIBBm68i4Ebp0dvfOHHW9qcjaQYXHoHg0vvAGD3obQRLWX+iFxy2QPcamarzKwH2ADsaIp5B3gAwMw+AnwcGC4yURERmd2cZ+juPmVmjwO7gDLwtLvvM7NN9fVbga8CXzOz15i+RPNld0+bqUBERK5IaIILd98J7GxatrXh+RHgN4pNTUREUqinqIhIJlTQRUQyoYIuIpKJtk0SbVWn60TsXuZlz1wX3u6XRn4nKY+U7uD9N8a75469Ec+52hu/g/Ojz4VDAegaj3fYLY8kdG12xycmQ6FL3zgdz+FCWv+E0Y/FD+FaV7xrevdYymz38e7jy/dfCMcCdJ9JuGd8MtbWXq1SPRNr654DlfDH9/QvDccC2Hh8X9SWxYcgKB07lZRHOIez8cnOpxOJny97rZi7uHWGLiKSCRV0EZFMqKCLiGRCBV1EJBMq6CIimVBBFxHJhAq6iEgmVNBFRDKhgi4ikgkVdBGRTKigi4hkom1juTBVpXQ6NjbK9S/HN1ua6k9KY2JxfBdYLT4+S9/i+PgeSyrxcRxqiS3W99P4WCA+9n58w9UqfjbWflarhje7sCvtHKPr/MJ4cMJwGd4Vb7+e4JhEALUFaQ1Yen8iHOvng3mYYd2xPDw45gtAyeL7DICu+L4onUgYZyhhu7WEn49q/DgG8AsJ4/Yk/I7MRmfoIiKZUEEXEcmECrqISCZCBd3M1pnZQTMbMrPNl4i5z8xeMbN9ZvZfxaYpIiJzmfOvB2ZWBrYADwIVYI+Z7XD3/Q0x/cCTwDp3f8fMbmxRviIicgmRM/S1wJC7D7v7BLAdWN8U81ngW+7+DoC7Hys2TRERmUukoK8EDje8rtSXNfo54Doz+08ze8nMPl9UgiIiEhO5YXOmm0ub7+jtAj4JPAD0Ac+Z2fPu/uaHNmS2EdgI0Ftekp6ttNXh8/upnJ++0jZB2tyY0l6Hxw9QGT8AwKQnzFMq80qkoFeAwYbXA8CRGWJG3P0ccM7MngXuBD5U0N19G7ANYFnPR4qZFVWumsG+NQz2rQHguePfaHM2kmKw9zYGe28DYPeJf2xzNtIqkUsue4BbzWyVmfUAG4AdTTH/BPyqmXWZ2ULgl4EDxaYqIiKzmfMM3d2nzOxxYBdQBp52931mtqm+fqu7HzCz7wGvAjXgKXd/fdYN16r4aLDr+KK+UBxA3/GpcCxA76l4d+XJhfHb9hcfjX8Bsal4bNd4Whfh8qmE7vwJXZUdx6di+7rU0xPP4Vza5YDekTPhWO+PX+arLornXDo/GY8dS7tUZecShhWYCObhjk8Gf0/K5fjnn0s41gArx3+fvFpL2nYruKddVPDa1b8IERr0wN13Ajublm1tev0E8ERxqYmISAr1FBURyYQKuohIJlTQRUQyoYIuIpIJFXQRkUyooIuIZEIFXUQkEyroIiKZUEEXEclE4hzyxfFqjerYuVBsuRbv9tv31omkPGpL4rPGL1gQ7wbtCTOglybj3fnLR0+GYwFI2He1lFnNnXDX//Bs9AApsQD9y+Kx0e7uQNfx+BAE3p3wa3TiVDyWhO78kDYrfXCW+ZTO6342NpTHRdaVsN8sfu6ZNKRA8BieL3SGLiKSCRV0EZFMqKCLiGRCBV1EJBMq6CIimVBBFxHJhAq6iEgmVNBFRDKhgi4ikolQQTezdWZ20MyGzGzzLHF3m1nVzH6ruBRFRCRizr63ZlYGtgAPAhVgj5ntcPf9M8T9JbAr9MlmWHBGcTsbGyIAgJTuxEA5YWZu743PBO+lhK7KlaPx7XbHcwDwc/F95xMTSdumFGs/H0+b6T6FHxsJx5au649vOKEbvY/EZ7u3hCEhIG3f+VTCMAHBtkthiccmpYR9kdIeU/HhLpK6/qfuM4/nUZRI1VkLDLn7sLtPANuB9TPEfQn4JnCswPxERCQoUtBXAocbXlfqy/6Pma0EfhPYWlxqIiKSInJ9YqbvRc3XKf4G+LK7V2f7SmlmG4GNAL3ERzmUzlCpDVGpvQXAJK27jCLFq/gw7/owoLbLWaSgV4DBhtcDwJGmmLuA7fVifgPwsJlNuft3GoPcfRuwDWBp6fqUkTmlAwyUbmGgdAsAz0/G/lQinWHAVjNgqwF4ofaDNmcjrRIp6HuAW81sFfAusAH4bGOAu6+6+NzMvgZ8t7mYi4hIa81Z0N19ysweZ/rulTLwtLvvM7NN9fW6bi4i0gFC9/i5+05gZ9OyGQu5u3/hytMSEZFU6ikqIpIJFXQRkUyooIuIZCKtn3yR3PHJWFfz2ukz4c2WVtyQlkdCl2I7n3D/7mjCDOgJXaZrp9JmjU+SMLM6EJ45vjYe38fl/mWJOSQM3TA+Ht/s6Fg4ttTXG46tjo6GYwEsZSiLlPYLtl2S4FAeF6UMNREdJiQ5j5Su/63YZwXTGbqISCZU0EVEMqGCLiKSCRV0EZFMqKCLiGRCBV1EJBMq6CIimVBBFxHJhAq6iEgmVNBFRDKhgi4ikon2jeXSItWjP02Kt76+cGzS2BM98fFZrJwwRkTieBkpY9VEx9ZJVornXD2TOtZJdzjWxxLG7UnYz0nHRcK4PQA+NZkQ3N5ZHZNyhaR8PWXMlYRjPjc6QxcRyYQKuohIJlTQRUQyESroZrbOzA6a2ZCZbZ5h/efM7NX6Y7eZ3Vl8qiIiMps5C7qZlYEtwEPAGuAzZramKext4Nfc/Q7gq8C2ohMVEZHZRc7Q1wJD7j7s7hPAdmB9Y4C773b3i1PpPA8MFJumiIjMJVLQVwKHG15X6ssu5VHgmStJSkRE0kXuQ7cZls14A6mZ3c90Qb/3Eus3AhsBelkYTFE6RcWHeZe3AZgkYX5VaTu13bUhUtArwGDD6wHgSHOQmd0BPAU85O4nZtqQu2+jfn19qS1vby8ISTZgqxlgNQAv+L+1ORtJoba7NkQuuewBbjWzVWbWA2wAdjQGmNlNwLeA33b3N4tPU0RE5jLnGbq7T5nZ48AuoAw87e77zGxTff1W4CvA9cCTZgYw5e53FZVk7ULrviL62FhCcMKXioTux7WUbs0J3egB8FpafCvUWtcVu1XDFXjCfvPJFn7ZtJmueHao1KEHUn62lG23eQiEdgqN5eLuO4GdTcu2Njx/DHis2NRERCSFeoqKiGRCBV1EJBMq6CIimVBBFxHJhAq6iEgmVNBFRDKhgi4ikgkVdBGRTKigi4hkwrxN3WTN7DhwqODN3gCMFLzNy5V7Lp8AXi54m7nvs8uhtkuXey4/6+4rZlrRtoLeCmb2YpFjyFwJ5ZKuk/LslFw6JY+5dFKe13IuuuQiIpIJFXQRkUzkVtA7aXJq5ZKuk/LslFw6JY+5dFKe12wuWV1DFxG5luV2hi4ics2aNwXdzNaZ2UEzGzKzzTOs/5yZvVp/7DazOxvW/cTMXjOzV8zsxRbncZ+Znal/1itm9pXoe1uQy5805PG6mVXNbHl9XWH7pIA8r0rbBXNR+6XnqfbrlPZz945/MD313VvAaqAH2AusaYq5B7iu/vwh4IWGdT8BbrhKedwHfPdy3lt0Lk3xjwD/XvQ+mS9tp/ZT+10L7TdfztDXAkPuPuzuE8B2YH1jgLvvdvdT9ZfPAwPtyKNF7y1ie58B/v4KPu9ydUrbhXJp0XuL2J7aT+03p/lS0FcChxteV+rLLuVR4JmG1w78q5m9ZGYbr0IenzKzvWb2jJndnvjeonPBzBYC64BvNiwuap/MpVPaLiUXtd8H1H5XlstVbb/QJNEdYKbpwWe8PcfM7mf6oLq3YfGvuPsRM7sR+L6ZveHuz7Yoj5eZ7po7ZmYPA98Bbg2+t+hcLnoE+B93P9mwrKh9MpdOabtoLmq/D1P7XX4uF1219psvZ+gVYLDh9QBwpDnIzO4AngLWu/uJi8vd/Uj932PAt5n+utSSPNx91N3H6s93At1mdkP0ZygylwYbaPq6V+A+mUuntF0oF7Xf/6P2u8xcGly99iv6onwrHkx/kxgGVvHBHyBub4q5CRgC7mlavghY0vB8N7CuhXl8lA/u718LvMP0/+ZzvrfoXOpxy4CTwKJW7JP50nZqP7XftdB+hR8ArXoADwNvMv2X5T+tL9sEbKo/fwo4BbxSf7xYX766vrP3AvsuvreFeTxe/5y9TP+B6J7Z3tvKXOqvvwBsb3pfoftkvrSd2k/tl3v7qaeoiEgm5ss1dBERmYMKuohIJlTQRUQyoYIuIpIJFXQRkUyooIuIZEIFXUQkEyroIiKZ+F+SDVIz3ZCjBgAAAABJRU5ErkJggg==\n",
            "text/plain": [
              "<Figure size 432x288 with 3 Axes>"
            ]
          },
          "metadata": {
            "needs_background": "light"
          },
          "output_type": "display_data"
        },
        {
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAACPCAYAAADqQSA9AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAP1UlEQVR4nO3db3Bc5XXH8e/RWsK2LGyEbVIkAXZxWogDbUjNnyQtSUpjaBxPp7wgydAmE0Zxp+50+iINM51Jp5O+aJo3aSdQj8tQpmlatxkINYwpJWTatBiDCTHGxrUr5D9am/qf/FfYlrR7+kJLWDaS9lz7rld6/PvMaNi99+yzR/dZju9e3ed5zN0REZHpr6XZCYiISD5U0EVEEqGCLiKSCBV0EZFEqKCLiCRCBV1EJBEzmvXG8zsLfl1Pa7PeXi7Qtr3QcXV7KPba1qEGZyNZ/KSvQOu8zlDs0vmHG5xNunadnRuOLR+K18Kh4/uPuPuC8fY1raBf19PKy8/2NOvt5QJd9YlBfvs794RiH+7a1OBsJIvZt5zlmt/7o1Dsyw/8TYOzSdendnw6HHvmr68Ox77w5B/vnWifLrmIiCRCBV1EJBEq6CIiiVBBFxFJhAq6iEgiVNBFRBKhgi4ikggVdBGRRIQKupktN7OdZtZnZg+Os3+umT1lZq+Z2XYz+2L+qYqIyGTqjhQ1swLwEHAXUAQ2m9l6d3+jKuz3gTfcfYWZLQB2mtl33X14onYPl1pZc7wrlOSqeftDcXLxnB6czQv/+KFY8Fc0UvRiuGdnbORu69DjvO+l0VDsXR9ZEX7/5254KhzbSP9yOj7k/m8HfjUcO/BCd6Y8Fr5aCsfO2XUsU9sTiZyhLwP63L2/UqDXAStrYhzoMDMD5gCDQOwTIyIiuYgU9C5goOp5sbKt2reBG4ADwOvAH7p7OZcMRUQkJDI5l42zrXZl6U8BW4BPAD8PPGdm/+XuJ9/TkFkv0AvQefVlmZOV5lr7nRM88g9jXVp6e3aTs5EsBp7aRvGpbQCMDE94JVSmuUhBLwLV0yJ2M3YmXu2LwF+4uwN9ZrYb+EXg5eogd18LrAW4dmlH7T8KMsX13j+X3vvHrk+2/9LZJmcjWfSsWErPiqUAvLDy8SZnI40SueSyGVhiZovMrA24D1hfE7MP+CSAmV0F/ALQn2eiIiIyubpn6O4+amargWeBAvCou283s1WV/WuArwOPmdnrjF2i+aq7H2lg3iIiUiO0wIW7bwA21GxbU/X4APAb+aYmIiJZaKSoiEgiVNBFRBKhgi4ikoimLRJ9cHAef/XPtQNOx7fqyw83OBsBWPzEl8OxM858i843Yvczr95/a7jdb3e9FI69FNz75q+HY/c8f10oruXMCO1b3wrFvvnf14Tf//2HfjccC7C0q/bu54lt2Z1hQfnjbeHQ9v3xc9qFO7MNfp/zv8fDsb43n+lNdIYuIpIIFXQRkUSooIuIJEIFXUQkESroIiKJUEEXEUmECrqISCJU0EVEEqGCLiKSCBV0EZFEqKCLiCSiaXO5FM7BvF2xdaTXD8XXr/xM+9vnm1KSssyjcvmuQjj22LkSs3YfC8X+cP0t4Xb/9N7T4ViAP1uwPVN8Izx/Jn7c/u7gxzK1/eqr14djr95RCsW9VSrjp06FYq97OhYHcHJxezgWoP+qJeHYy66ItxutKwCzD54Lx7YNnoknAdhQPL5civVdPTpDFxFJhAq6iEgiVNBFRBIRKuhmttzMdppZn5k9OEHMnWa2xcy2m9l/5pumiIjUU/ePomZWAB4C7gKKwGYzW+/ub1TFzAMeBpa7+z4zW9igfEVEZAKRM/RlQJ+797v7MLAOqF1q6HPAE+6+D8DdD+WbpoiI1BMp6F3AQNXzYmVbtfcDV5jZf5jZj83sd/JKUEREYiL3ods423ycdm4BPgnMAl40s03uvus9DZn1Ar0AbbMz3FgqU8Lg6y8yuO1FALyk+/2nk4P9mzi4e2y91pJnu59apo9IQS8C1Su0dgO1q7sWgSPuPgQMmdmPgJuB9xR0d18LrAVov7Kn9h8FmeI6P3g7nR+8HYC9a/6yydlIFlctvo2rFt8GwPbvf6PJ2UijRC65bAaWmNkiM2sD7gPW18T8K/AxM5thZrOBW4Ed+aYqIiKTqXuG7u6jZrYaeBYoAI+6+3YzW1XZv8bdd5jZvwFbgTLwiLtvm6zdwrkyHXtjX/3+fNdvhuIAPvPL3wvHTkef3nV3pvjtO3rqB1Vct2skHLu3VIJjJ0OxHXviNz199wfZhsa/efuCcOz+obnh2Fkz4sfirZOXh2PPbe4MxwJ0DGWIffN4MNKhHPuCXJoZnx2k7XS24esWzAFgXl98OH/ryXjftR04Ho7l7WyXqjxDvNl4V7azC/WWu28ANtRsW1Pz/JvAN3PJSkREMtNIURGRRKigi4gkQgVdRCQRKugiIolQQRcRSYQKuohIIlTQRUQSoYIuIpIIFXQRkUTEx/XmzIZHad17OBR7ZE98+Pptdm+mPJYt3BuOXTIrPs37iMdXgn+s79Zw7PDWeeFYgK6t8SHTs/YdD8d6qUT5+IlQ7OzDo+F22ZLtI/mTozeGY0fmxIeae4Y0Zh+ID9ueeyjeHwAdu+PDx1tOBOcJKDs+PBwKbT10Kvz+hbcvC8cCtB+Lz2tQurIjHNty+mw8CY9/Json48ciq/LZDDlPQmfoIiKJUEEXEUmECrqISCJU0EVEEqGCLiKSCBV0EZFEqKCLiCRCBV1EJBEq6CIiiQgVdDNbbmY7zazPzB6cJO5XzKxklnG4poiIXLC6A5zNrAA8BNwFFIHNZrbe3d8YJ+4bwLOhd/Yyfu5cKLTr+VAYAEeWxleYB3h6/pXh2BkLM6z63d8eDi0Mx4eOX747PlQZoL2YYeXxoWyrlFtrbHz8rP7BeLuj88KxAC2jreHYUlv8OM84Gz/OrUPx1e5bRrP134zB+PB4Dw5NdxwfjU3H0HIiPty95eixcOxYHnGFUvwY+9lYXQGgJX6RwjPkAECWeIt/Nic7cJHfZhnQ5+797j4MrANWjhP3B8DjQHzCExERyU2koHcBA1XPi5VtP2VmXcBvAWvyS01ERLKIfGce77tA7Un/t4CvunvJJvnqYGa9QC/AzJY5wRRlqhg4tZWBU68DMOz5zA4nF8fA2R0Uz+4AYMQzXJKQaSVS0ItA9fy13cCBmpgPA+sqxXw+cI+Zjbr7k9VB7r4WWAswt3VBtouJ0nQ9HTfR03ETABv3/X2Ts5EsembeQM/MGwDYePR7Tc5GGiVS0DcDS8xsEbAfuA/4XHWAuy9657GZPQY8XVvMRUSkseoWdHcfNbPVjN29UgAedfftZraqsl/XzUVEpoDQfWfuvgHYULNt3ELu7l+48LRERCQrjRQVEUmECrqISCJU0EVEEpFtifUmad8XH/5s5dmZ2j7TWQjHusWH87edjt+V2VKKrwQ/Z0/8WAC0DMXvOfbTGYaZu+MjweHjZ+I5zOo/Go4FuOxIvL9HO+Kr0reci/1umWOPnw7HAjCaYch7OfiZ83js6MH4wO9C5xXhWMg2RL98KsNxyzTkPsM5rcf/PwXC0yvkSWfoIiKJUEEXEUmECrqISCJU0EVEEqGCLiKSCBV0EZFEqKCLiCRCBV1EJBEq6CIiiVBBFxFJhAq6iEgimjaXi5dKlI+fCMW2DMTnW5nZ9nOZ8ph5eOI1UGtZhjlXCoPxeVHKl8+Kt3v0VDgWwN8+E8/jZIb5MtzxkeFYu4ePhJttWTA/ngNgZ2I5ALT9X3yeGM/Q19YRn+OHDO0CeIY5TEqnMnw2ysH5TiZZI/hnmswwF9BY0/G2ic5TQ4Y5bQBrydYfmWT5/bKY5NfTGbqISCJU0EVEEqGCLiKSiFBBN7PlZrbTzPrM7MFx9n/ezLZWfjaa2c35pyoiIpOpW9DNrAA8BNwN3Ah81sxurAnbDfyau98EfB1Ym3eiIiIyucgZ+jKgz9373X0YWAesrA5w943ufqzydBPQnW+aIiJST6SgdwEDVc+LlW0T+RLwzIUkJSIi2UXuQx/vZspx74Q0s48zVtA/OsH+XqAXYCbZ1v6U5iuW+iiW+wAYIb4epDRf0fvZ7/2A+i5lkYJeBHqqnncDB2qDzOwm4BHgbncfdwSHu6+lcn398pbO+N3/MiV0F66nu3A9AJuG9SVsOum2xXTbYgBeKv+gydlIo0QuuWwGlpjZIjNrA+4D1lcHmNk1wBPA/e6+K/80RUSknrpn6O4+amargWeBAvCou283s1WV/WuArwFXAg9XhvOOuvuH80oyy/D1ls07MrXdMm9uONZmxKcg8LPxr7W2Pz58vTQ8Eo4FsEJ8qEF0KP9PtcSORzlDzuWBYqYUChn6rzwyGo71cxn6L8OQex+N55CZN+BLr2X4/Axn+/x4hnxtRoZZSjw+nN9H07pQEDpK7r4B2FCzbU3V4weAB/JNTUREstBIURGRRKigi4gkQgVdRCQRKugiIolQQRcRSYQKuohIIlTQRUQSoYIuIpIIFXQRkURYluG3ub6x2WFgb87NzgfiS8w3Vuq5fAh4Nec2Uz9m50N9l13quVzr7gvG29G0gt4IZvZKnnPIXAjlkt1UynOq5DJV8qhnKuV5KeeiSy4iIolQQRcRSURqBX0qLU6tXLKbSnlOlVymSh71TKU8L9lckrqGLiJyKUvtDF1E5JI1bQq6mS03s51m1mdmD46z//NmtrXys9HMbq7at8fMXjezLWb2SoPzuNPMTlTea4uZfS362gbk8pWqPLaZWcnMOiv7cjsmOeR5UfoumIv6L3ue6r+p0n/uPuV/GFv67k1gMdAGvAbcWBNzB3BF5fHdwEtV+/YA8y9SHncCT5/Pa/POpSZ+BfDDvI/JdOk79Z/671Lov+lyhr4M6HP3fncfBtYBK6sD3H2jux+rPN0EdDcjjwa9No/2Pgv80wW83/maKn0XyqVBr82jPfWf+q+u6VLQu4CBqufFyraJfAl4puq5A/9uZj82s96LkMftZvaamT1jZh/I+Nq8c8HMZgPLgcerNud1TOqZKn2XJRf137vUfxeWy0XtvwxLaTeVjbNt3NtzzOzjjH2oPlq1+SPufsDMFgLPmdn/uPuPGpTHq4wNzT1tZvcATwJLgq/NO5d3rABecPfBqm15HZN6pkrfRXNR/72X+u/8c3nHReu/6XKGXgR6qp53Awdqg8zsJuARYKW7H31nu7sfqPz3EPB9xr4uNSQPdz/p7qcrjzcArWY2P/o75JlLlfuo+bqX4zGpZ6r0XSgX9d/PUP+dZy5VLl7/5X1RvhE/jH2T6AcW8e4fID5QE3MN0AfcUbO9HeioerwRWN7APN7Hu/f3LwP2Mfaved3X5p1LJW4uMAi0N+KYTJe+U/+p/y6F/sv9A9CoH+AeYBdjf1n+k8q2VcCqyuNHgGPAlsrPK5XtiysH+zVg+zuvbWAeqyvv8xpjfyC6Y7LXNjKXyvMvAOtqXpfrMZkufaf+U/+l3n8aKSoikojpcg1dRETqUEEXEUmECrqISCJU0EVEEqGCLiKSCBV0EZFEqKCLiCRCBV1EJBH/D5YEMKUnQtf3AAAAAElFTkSuQmCC\n",
            "text/plain": [
              "<Figure size 432x288 with 3 Axes>"
            ]
          },
          "metadata": {
            "needs_background": "light"
          },
          "output_type": "display_data"
        },
        {
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAACPCAYAAADqQSA9AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAOUUlEQVR4nO3da4xc5X3H8e9vxzbGGN8xgfVCnOIkhRar2HG5NUBSgkGlbtq0ci60oKCVpTrtm1ZBapUXzas2rdSqJbEsiyQKTaw2IdSNljqJKieoxsUsxVews/Ea7+KC75e1wevd+ffFjsMw3d15Zj3jnX38+0grZs75zzl/5ll+nDl7njOKCMzMbOJrGe8GzMysPhzoZmaZcKCbmWXCgW5mlgkHuplZJhzoZmaZmDReO543pxDvb5s8Xru3i7S/5zwev4nJYzexdW4/dyQirhlu3bgF+vvbJvPixrbx2r1dpGUP9Hj8JiiP3cRWuK7r9ZHW+ZSLmVkmHOhmZplwoJuZZcKBbmaWCQe6mVkmHOhmZplwoJuZZcKBbmaWiaRAl7Rc0h5JXZKeGGb9TEn/LmmbpF2SHqt/q2ZmNpqqM0UlFYAngfuBXmCrpA0Rsbus7I+B3RHxsKRrgD2S/jki+hvStZmNWRCcj8Gk2skqNLgbA7jp26tqqP6zEdekHKEvA7oiYl8poNcDKypqArhakoDpwDFgoIYOzczsIqUEeivQU/a8t7Ss3D8BvwwcBHYAfxoRxbp0aGZmSVJuzqVhllV+s/QDwCvAx4BfAn4k6fmIOPWeDUntQDvADa3jdl8wG6O13zrJuqeHhvTw0bSP7NYcPHaXh5Qj9F6g/NZsCxg6Ei/3GPBMDOkCuoEPV24oItZGxNKIWHrNXJ+bm2jaH5nJixvbeHFjGx6/icVjd3lICfStwCJJCyVNAVYCGypqDgAfB5B0LfAhYF89GzUzs9FVPe8REQOSVgMbgQLwVETskrSqtH4N8GXgG5J2MHSK5osRcaSBfZuZWYWkE9kR0QF0VCxbU/b4IPCJ+rZmZma18ExRM7NMONDNzDLhQDczy4QvBje7zOw8Op+bn16dVPuzR77W4G7ydfPmzyXXtm5Kn4fZPco6H6GbmWXCgW5mlgkHuplZJhzoZmaZcKCbmWXCgW5mlgkHuplZJhzoZmaZcKCbmWXCgW5mlgkHuplZJnwvFxuTnX1z+eBP/iipdu8932xwN1aLwjmYuXe8u5h47t7+uzXVq3NGcu203hM1djM8H6GbmWXCgW5mlgkHuplZJpICXdJySXskdUl6YoSaeyW9ImmXpJ/Ut00zM6um6h9FJRWAJ4H7gV5gq6QNEbG7rGYW8FVgeUQckDS/Qf2amdkIUo7QlwFdEbEvIvqB9cCKiprPAM9ExAGAiDhU3zbNzKyalEBvBXrKnveWlpX7IDBb0iZJnZL+sF4NmplZmpTr0DXMshhmO0uAjwNXAi9I2hIR77naVVI70A5wQ6svgZ9o1n7rJOuePgXA4Kmp49yN1aJ87AbemTbO3VijpKRqL9BW9nwBcHCYmiMRcQY4I+mnwGLgPYEeEWuBtQBLF0+t/J+CNbn2R2bS/shMAKbddXacu7FalI/d9F99Z5y7sUZJOeWyFVgkaaGkKcBKYENFzb8BvyFpkqRpwK8Dr9a3VTMzG03VI/SIGJC0GtgIFICnImKXpFWl9Wsi4lVJ/wFsB4rAuojYOdp2+yLY8s5gUpO3Ty0k1dml03K6has3JX50v6exvdiQB179raS6SW//C3N39CXV3rvzd5L3v+lXnk2ubaR/PH5jcu0/bLsvuXb6ltpOVc3rHkiu1dlzNW17JEknsiOiA+ioWLam4vlXgK/UpSszM6uZZ4qamWXCgW5mlgkHuplZJhzoZmaZcKCbmWXCgW5mlgkHuplZJhzoZmaZcKCbmWVi3G552N03j8+98HhSbdd9X29wNwawpPMPkmsL577J1QfSpjZ3n0+bZg6wcPL05Fp7r/2b26oXARosUuhLm2p+vOP65P3/2tsrk2sB7rq+O7m288iC5No3e+ck17b0pd9WZOrRYnItwJUHzyTXFl/vrWnbI/ERuplZJhzoZmaZcKCbmWXCgW5mlgkHuplZJhzoZmaZcKCbmWXCgW5mlgkHuplZJpICXdJySXskdUl6YpS6j0galPSp+rVoZmYpqk79l1QAngTuB3qBrZI2RMTuYer+GtiYsmP1t9Dy+tTaO7aanCy+nVx7duu85FoVg8l9aVP/f/P5LyRv9+cfy/s2DwcG0m+DAPCpHY8l187ZHUl1h4pFdCbt92LG64PJ+z//r7OSawE2LViSXNtyPn27s9J/5Zm5P33DUw+/k75hoHAsfayLk+pzF5aUI/RlQFdE7IuIfmA9sGKYui8A3wMO1aUzMzOrSUqgtwI9Zc97S8t+QVIr8ElgTf1aMzOzWqQc52uYZZWf7f4e+GJEDErDlZc2JLUD7QCTZs1ObNGaxYnOFzjZ+cLQk/Ppd5Kz8XeoawuHurYAUBys4ZyETSgpgd4LlN+XcwFwsKJmKbC+FObzgIckDUTEs+VFEbEWWAswdUFb2gk/axqzltzBrCV3APC/f/e349yN1WL+Tbcz/6bbAXjtu38zzt1Yo6QE+lZgkaSFwBvASuAz5QURsfDCY0nfAH5QGeZmZtZYVQM9IgYkrWbo6pUC8FRE7JK0qrTe583NzJpA0rUyEdEBdFQsGzbII+LRi2/LzMxq5ZmiZmaZcKCbmWXCgW5mlon6zDcdgykni9zYkXY97Ibfn5a83d++6uxYW8rSR55flVx7fWf6NOg3zw0wpTttUvDcH9+QvN37r3s4uRbgtjk91YtK/mr+1uTaP3njo8m1+06n3zLhrdPTk2sBtCl9vsb7utOmmsf5AQbfTBu7Kw/NTd5/XHdFci3A7D3ptxWYeqw/uVb9xeTawulzybUtJ2u7bUPx2PH02rdru63ASHyEbmaWCQe6mVkmHOhmZplwoJuZZcKBbmaWCQe6mVkmHOhmZplwoJuZZcKBbmaWCQe6mVkmHOhmZpkYt3u5qH+AKT1Hk2r/8muPJm/327/3Wk19fGj6W8m1n5zxcnJtx+lbk2vX77stufbsz2Yl1wJc8z/p3/R31atvpm94cJDi8RNJpXN2zUne7JHvLEjvAfjhFW3Vi0qeve6O5NopJ0b+btxKLem3wGHW/oH0YmByDfcaKXTXMH7FtN+LSXvfSN7kjMMz0vcPUEy/5wqjfFfxxYiD6f/9R6FQ27bPpY8dUcN7MQofoZuZZcKBbmaWCQe6mVkmkgJd0nJJeyR1SXpimPWflbS99LNZ0uL6t2pmZqOpGuiSCsCTwIPAzcCnJd1cUdYN3BMRtwJfBtbWu1EzMxtdyhH6MqArIvZFRD+wHlhRXhARmyPiwtdzbAFqu1TBzMwuWkqgtwLl3/PVW1o2ks8Dz11MU2ZmVruU69CHuwB02ItYJd3HUKDfPcL6dqAdYGrh6sQWrVn0nNpGz6ntAPRHfb4D0S6NnrO76D27G4DzHrtspQR6L1A+e2MBcLCySNKtwDrgwYgYdsZQRKyldH595hXXps94sabQNmMxbTOG/t69ufvr49yN1aJt2i20TbsFgM1vrR/nbqxRUk65bAUWSVooaQqwEthQXiDpBuAZ4JGI2Fv/Ns3MrJqqR+gRMSBpNbARKABPRcQuSatK69cAXwLmAl/V0BTdgYhYOuqGi0Wi70xSk60/PpZUB3C088bkWoDvL/lwcu1zxz6aXKvB9B5mH08vbu05kb5hoOXY6eTa4tH09zmimDy1uWXP68nbvfbwrORagMHZ6afuWvrTp90PzJyaXDvpZA2nMI4cr15Trobp49Gfeg+CSJ5qXjxxMnn/OpP23/Mv6qddmV4c6R/oi6f60jd7vj+5VpOnJNcCxEBtt3moh6R7uUREB9BRsWxN2ePHgcfr25qZmdXCM0XNzDLhQDczy4QD3cwsEw50M7NMONDNzDLhQDczy4QD3cwsEw50M7NMONDNzDKhqGFKbV13LB0G0ueEp5kHHKnzNscq915uA16u8zZzf8/GwmNXu9x7uTEirhluxbgFeiNIeqnqPWQuEfdSu2bqs1l6aZY+qmmmPi/nXnzKxcwsEw50M7NM5BbozfTl1O6lds3UZ7P00ix9VNNMfV62vWR1Dt3M7HKW2xG6mdlla8IEuqTlkvZI6pL0xDDrPytpe+lns6TFZev2S9oh6RVJLzW4j3slnSzt6xVJX0p9bQN6+fOyPnZKGpQ0p7Subu9JHfq8JGOX2IvHr/Y+PX7NMn4R0fQ/DH313c+BDwBTgG3AzRU1dwKzS48fBP67bN1+YN4l6uNe4AdjeW29e6mofxj4z3q/JxNl7Dx+Hr/LYfwmyhH6MqArIvZFRD+wHlhRXhARmyPiwhc2bgEWjEcfDXptPbb3aeA7F7G/sWqWsUvqpUGvrcf2PH4ev6omSqC3Aj1lz3tLy0byeeC5sucB/FBSp6T2S9DHHZK2SXpO0i01vrbevSBpGrAc+F7Z4nq9J9U0y9jV0ovH710ev4vr5ZKOX9KXRDcBDbNs2MtzJN3H0C/V3WWL74qIg5LmAz+S9FpE/LRBfbzM0NTcPkkPAc8CixJfW+9eLngY+K+IOFa2rF7vSTXNMnapvXj83svjN/ZeLrhk4zdRjtB7gbay5wuAg5VFkm4F1gErIuLoheURcbD0z0PA9xn6uNSQPiLiVET0lR53AJMlzUv9d6hnL2VWUvFxr47vSTXNMnZJvXj8/h+P3xh7KXPpxq/eJ+Ub8cPQJ4l9wELe/QPELRU1NwBdwJ0Vy68Cri57vBlY3sA+3se71/cvAw4w9H/zqq+tdy+lupnAMeCqRrwnE2XsPH4ev8th/Or+C9CoH+AhYC9Df1n+i9KyVcCq0uN1wHHgldLPS6XlHyi92duAXRde28A+Vpf2s42hPxDdOdprG9lL6fmjwPqK19X1PZkoY+fx8/jlPn6eKWpmlomJcg7dzMyqcKCbmWXCgW5mlgkHuplZJhzoZmaZcKCbmWXCgW5mlgkHuplZJv4PPt2jF72j5bwAAAAASUVORK5CYII=\n",
            "text/plain": [
              "<Figure size 432x288 with 3 Axes>"
            ]
          },
          "metadata": {
            "needs_background": "light"
          },
          "output_type": "display_data"
        }
      ],
      "source": [
        "name = \"same_but_GOU_data\"\n",
        "model = \"GOU\"\n",
        "mx = [30, 30]    # get prob dist for this domain ; bigger = more accurate\n",
        "K_avg = 5\n",
        "beta, gamma = 1, 1.5\n",
        "\n",
        "get_BF_plot(name, model, K_avg, beta, gamma, mx)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "jJUv-aiZ5Wan",
        "outputId": "5ae82787-edba-4877-c4bc-ea35dfcb1843"
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Computing log Bayes factors on a  10  x  10  grid...\n",
            "Row  1 / 10 ...\n"
          ]
        },
        {
          "name": "stderr",
          "output_type": "stream",
          "text": [
            "<ipython-input-7-0aa11ad1c9c5>:29: RuntimeWarning: divide by zero encountered in log10\n",
            "  lp_other = np.log10(pss_other)\n"
          ]
        },
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Row  2 / 10 ...\n",
            "Row  3 / 10 ...\n",
            "Row  4 / 10 ...\n",
            "Row  5 / 10 ...\n",
            "Row  6 / 10 ...\n",
            "Row  7 / 10 ...\n",
            "Row  8 / 10 ...\n",
            "Row  9 / 10 ...\n",
            "Row  10 / 10 ...\n",
            "--- 92.30604004859924 seconds ---\n"
          ]
        },
        {
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAACPCAYAAADqQSA9AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAP40lEQVR4nO3da4yc1X3H8e9v9uK18YWbDY3XoeYWShqQGnBbmqqkEa1BQlalvjBBrRIRWZZK1VdVkCrlTd60SlVFaoksCyFaqarfJE1d5BRSRW0qXKghMhhDjRbj2GsXfMHYu17be5l/X+xQhul693/sGc/48e8jrTzzzH/O/Pc5M/995vFzzlFEYGZmV75atxMwM7P2cEE3M6sIF3Qzs4pwQTczqwgXdDOzinBBNzOriP5uvfBg/5JYPHhtKvb89X3pdiMfOpvH0FQ6VuQv8Ryo1dOxE1MD6dh6XelYAE3m/2b3n8n/fpOnTrC4f3kqdmrFYLrd0v6rLyq47Lao7YJ2p/P7WDMlOYDybyP6J3I5nx/P993MknzflR4e1gv6o2Q/REEetemSHMou8a6dL+jsyXwdOl0/cTwiVs71WNcK+uLBa/m1259Ixe5/7Lp0u1MrCnoeuOVz76dj+wuK9E2Lx9Kxe47+Qjp2YmJROhag9t7idOxNu/JvwD3P/yUP3LQxFfvB+lvS7U6uKPuDNXZn/hOpJfnYWl/+wxvH8n0ycKqs6vVP5GNX7p5Mxe3+yV/xwGceT8We/uJn0q8/M1jWd2dvyO+LwbF8f0xdk89h8fH8Z3pgoqy2LDlwOh0bB0bTsS+O/+3PL/SYT7mYmVWEC7qZWUW4oJuZVYQLuplZRbigm5lVhAu6mVlFuKCbmVWEC7qZWUWkCrqk9ZL2SRqR9NQcj6+Q9M+SXpe0V9LX25+qmZnNZ8GRopL6gKeBh4BRYJek7RHxVlPYHwFvRcSjklYC+yT9fURceOja1BT6n6OpJG/ZMZSKAzjxy/mRkQDHDq1Ox9byo3M5VhA7cDY/Cm7leNnw46mC3TF09Hw+uD5DjI2nQlftPJ5udvzO/KhggMXH8uPHazP5L6Qlw8f7z+X7ZOjEuXzDwLnr89NCLMq2PT1D/diJVOiyPfnXn165LB0LsPRgPrY2lR+lqfMF03mcK4idLpy3YeJsOrQ+XTAHwTwyb9t1wEhE7G8U6G3AhpaYAJZJErAU+BBoT4ZmZpaSKeirgUNN90cb25r9DfBLwBFgD/AnEVE28YGZmV2SzORcc8240/od83eB3cBvA7cBP5b0HxHxqdlpJG0CNgEM1ZYWJ2vddfj9XRx+/1UApuplpw6su0aPvsrosdcAmAr3XVVlCvoosKbp/jCzR+LNvg78eUQEMCLpPeAu4L+agyJiK7AVYMXAyrKTwdZ1q2++n9U33w/Arpe/291krMjwqvsYXnUfAK+89tddzsY6JXPKZRdwh6S1kgaBjcD2lpiDwFcAJN0EfA7Y385EzcxsfgseoUfEtKQngReYXSLg2YjYK2lz4/EtwLeB5yTtYfYUzTcjIn9pg5mZXbLUAhcRsQPY0bJtS9PtI8DvtDc1MzMr4ZGiZmYV4YJuZlYRLuhmZhXRtUWikdBAblhx/4dn0s0uP1D2K/VP5IfznrwrPwVBX8GlvovG8jn0FQwzB1g+kh9+3Dd6LN9wwOxVqgurfZRfMHvJaMEq88DSU/lVlM/cNedC6XMqGWqO8osj950pmBMCWP5+ft+RnEqjyAf5axsGZsrGEkatYFHpgfznWmcLprCYKhjQXi/7/eofnSqKbwcfoZuZVYQLuplZRbigm5lVhAu6mVlFuKCbmVWEC7qZWUW4oJuZVYQLuplZRbigm5lVhAu6mVlFuKCbmVVE1+ZyielpZo6fSMXWrluebnfo8HhRHhO/mG/7urfz86JEf/5vZW0yP5dLfVFfOhagdjqfM4vy86hEvU59LDfPiJYsTrfb98FH6ViAqVvy87MsOXB64aCGGMjv59pEfu6Q+rL8fEAAmiiYFGhpbp3eOFmnfjb3vuhbld+/cbJs7hItL1hXuGB+lpjIv+c1tCjf7pn8vEEA9OXfQ3EmP1/VfHyEbmZWES7oZmYV4YJuZlYRqYIuab2kfZJGJD11gZgHJe2WtFfSv7c3TTMzW8iC/ykqqQ94GngIGAV2SdoeEW81xVwLfA9YHxEHJa3qUL5mZnYBmSP0dcBIROyPiElgG7ChJearwA8i4iBARHRg6RQzM5tPpqCvBg413R9tbGt2J3CdpH+T9JqkP2xXgmZmlpO5Dn2uhf9aF5PsB74IfAVYDPynpJcj4p1PNSRtAjYBDLGkPFvrqkNjb3BobA8AUxSs22hdd2h8D4fG3wTcd1WWKeijwJqm+8PAkTlijkfEGeCMpJ8C9wKfKugRsRXYCrC8dn3ZasfWdWuW3cOaZfcA8NKB57qbjBVZs/QLrFn6BQBeOvR3Xc7GOiVzymUXcIektZIGgY3A9paYfwJ+U1K/pCXArwJvtzdVMzObz4JH6BExLelJ4AWgD3g2IvZK2tx4fEtEvC3pX4A3gDrwTES8OX/DEPXkQfr7x3NxgG6+MR0LsOjD/NfP2vnpfB7n88P5da7gK/BA2WwNOj+Zjo3xwuHHmuts3BztnsoPude1K4pS6D+ab5tawXQMZwqmTCjQd6zwdEdJ/50rmCYgcp+9+kf54fy1JWWnUeN02TQd6XYnC/ZZcgoEgJipl+UxPZUPTn6WZhu+8EOp6hARO4AdLdu2tNz/DvCdfFZmZtZOHilqZlYRLuhmZhXhgm5mVhEu6GZmFeGCbmZWES7oZmYV4YJuZlYRLuhmZhXhgm5mVhFl48jbLZJDafsK/u5M5YfnAwwczE8rECvyq5TrVH5Yc8mQbQ0MpGMB6iXD7vs79HYoGHJfOryaYyfSoVq2LJ9HwXQMMTaWz2Hx4nQslA1jL9nPadnpOSgbRg+gwcHSbHJ5TOaH3Guw4PM0XVZbstMrtJOP0M3MKsIF3cysIlzQzcwqwgXdzKwiXNDNzCrCBd3MrCJc0M3MKsIF3cysIlzQzcwqIlXQJa2XtE/SiKSn5om7X9KMpN9vX4pmZpax4FhvSX3A08BDwCiwS9L2iHhrjri/AF5IvbKE+nPDbuNswdD4dGTDNfmVynX6TL7dgmG/UbKye0EslA3nr09MFLWd/R2Lhq/39RWloIJpIUqG83O+ILZkaoOSdqFoNfj6eH66iWy7MTOTbzM7lcdFKJkSouQ9Ufye73GZ33wdMBIR+yNiEtgGbJgj7o+B7wNH25ifmZklZQr6auBQ0/3Rxrb/I2k18HvAlvalZmZmJTLfx+f6btb6Xfu7wDcjYkbzfJWTtAnYBDBE/lSH9YbR+gij9XcBmKLw1IF11Wjs53DsB9x3VZYp6KPAmqb7w8CRlpj7gG2NYn4j8Iik6Yj4YXNQRGwFtgIsr91w+eeWtEsyXLud4drtALw8lfuvEusNw7qVYd0KwCv1f+1yNtYpmYK+C7hD0lrgMLAR+GpzQESs/fi2pOeA51uLuZmZddaCBT0ipiU9yezVK33AsxGxV9LmxuM+b25m1gNS17RFxA5gR8u2OQt5RHzt0tMyM7NSHilqZlYRLuhmZhXhgm5mVhEdWuY9IYKYTq7OXSsY0D9eMDy/kwpWCC8Zqlxq5nTBcPAODd2OkpXjz5QNxa6VTN1AfgqCesl0EyUrxxcqyaNolflsbOSH/gdl0zbEZPLzD6hgSoiSaQKK9lnBNAzd4iN0M7OKcEE3M6sIF3Qzs4pwQTczqwgXdDOzinBBNzOrCBd0M7OKcEE3M6sIF3Qzs4pwQTczqwgXdDOziujeXC6QnkchzufXQNSyZWUpnM/P70E9P0dEyTwVJWImP7cGlM6BUZpNUj3fcG1oqKztgv1RP9eZtTSL+rpwvhwVzGOUnkulU/1cquB9ESX7rWR+lhKdareNfIRuZlYRLuhmZhXhgm5mVhGpgi5pvaR9kkYkPTXH449LeqPxs1PSve1P1czM5rNgQZfUBzwNPAzcDTwm6e6WsPeA34qIe4BvA1vbnaiZmc0vc4S+DhiJiP0RMQlsAzY0B0TEzog42bj7MjDc3jTNzGwhmYK+GjjUdH+0se1CngB+dClJmZlZucx16HNdCDvnBZmSvsxsQf/SBR7fBGwCGCK/FqT1htH6uxyOdwGYojPXdFtnuO+uDpmCPgqsabo/DBxpDZJ0D/AM8HBEnJiroYjYSuP8+nJd3/tX6dunDNduY5jbAHhl5sUuZ2Ml3HdXh8wpl13AHZLWShoENgLbmwMkfRb4AfAHEfFO+9M0M7OFLHiEHhHTkp4EXgD6gGcjYq+kzY3HtwDfAm4AvicJYDoi7lvw1ZUc1lww5LY+NpaOLcoBUP9AOjamCqYUKMihdPhxFAyv7gX1c+c613in9nMtP71CyXD3K04nfzcVDJnp2BwWvS81l0tE7AB2tGzb0nT7G8A32puamZmV8EhRM7OKcEE3M6sIF3Qzs4pwQTczqwgXdDOzinBBNzOrCBd0M7OKcEE3M6sIF3Qzs4pQdGkla0nHgJ+3udkbgeNtbvNiVT2XXwF+1uY2q77PLob7rlzVc7klIlbO9UDXCnonSHo1NYfMZeBcyvVSnr2SS6/ksZBeyvNqzsWnXMzMKsIF3cysIqpW0HtpcWrnUq6X8uyVXHolj4X0Up5XbS6VOoduZnY1q9oRupnZVeuKKeiS1kvaJ2lE0lNzPP64pDcaPzsl3dv02AFJeyTtlvRqh/N4UNKpxmvtlvSt7HM7kMufNuXxpqQZSdc3HmvbPmlDnpel75K5uP/K83T/9Ur/RUTP/zC79N27wK3AIPA6cHdLzAPAdY3bDwOvND12ALjxMuXxIPD8xTy33bm0xD8K/KTd++RK6Tv3n/vvaui/K+UIfR0wEhH7I2IS2AZsaA6IiJ0RcbJx92VguBt5dOi57WjvMeAfLuH1Llav9F0qlw49tx3tuf/cfwu6Ugr6auBQ0/3RxrYLeQL4UdP9AF6U9JqkTZchj1+X9LqkH0n6fOFz250LkpYA64HvN21u1z5ZSK/0XUku7r9PuP8uLZfL2n+pRaJ7wFzLtc95eY6kLzP7pvpS0+bfiIgjklYBP5b03xHx0w7l8TNmh+aOS3oE+CFwR/K57c7lY48CL0XEh03b2rVPFtIrfZfNxf33ae6/i8/lY5et/66UI/RRYE3T/WHgSGuQpHuAZ4ANEXHi4+0RcaTx71HgH5n9utSRPCLidESMN27vAAYk3Zj9HdqZS5ONtHzda+M+WUiv9F0qF/ff/+P+u8hcmly+/mv3SflO/DD7TWI/sJZP/gPi8y0xnwVGgAdatl8DLGu6vRNY38E8buaT6/vXAQeZ/Wu+4HPbnUsjbgXwIXBNJ/bJldJ37j/339XQf21/A3TqB3gEeIfZ/1n+s8a2zcDmxu1ngJPA7sbPq43ttzZ29uvA3o+f28E8nmy8zuvM/gfRA/M9t5O5NO5/DdjW8ry27pMrpe/cf+6/qvefR4qamVXElXIO3czMFuCCbmZWES7oZmYV4YJuZlYRLuhmZhXhgm5mVhEu6GZmFeGCbmZWEf8LgEdT/P3EmpMAAAAASUVORK5CYII=\n",
            "text/plain": [
              "<Figure size 432x288 with 3 Axes>"
            ]
          },
          "metadata": {
            "needs_background": "light"
          },
          "output_type": "display_data"
        },
        {
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAACPCAYAAADqQSA9AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAQZ0lEQVR4nO3db4xc5XXH8e/Z2b82uzbGdgneJZiYpjiqkUhwCU1U0oTUICGrUisRUCqiJJalkPKqAqlV+iKq1Aip6ouCLIvSqFVVt80f6kambqKqhMZ1aqAYMMRkMbZ32IC92Ov12t6/c/pix2GYrnfPY9/x7D7+faSVd+6cuXP2PrNn74zveR5zd0REZPFraXYCIiJSDBV0EZFMqKCLiGRCBV1EJBMq6CIimVBBFxHJRGuznritY6l3LF0Riu350Onwfq9rHbvYlLJ0dHJpOPb0REc4dvrIKTo7l4dix1eFd8uqrtF4MHBt63hSfLOdSbxM+Mi52O8IwPREKRTnb5+gq21ZKHZ8RUKJ6KzEY4H21ql4bMt0OHY04XXMWMI5beIV3u2j8QfYZPznGxl7Z8jdZ/2talpB71i6gl+/6+FQ7Of/5Lnwfv901WsXm1KW/nDwtnDss+V14dgjD/w1t936tVDsW1vDu+WrG/4rHgw8cs3Pk+Kb7YXxiaT4r77yxXDsycFYkR5++C/45A0PhmIP3b86/Py2Pn7iBfDha06EY9csORWOfe7wjeFY+3n8hKdlysKxAGuejZ9cth8/E47d/eqfHbnQffrIRUQkEyroIiKZUEEXEcmECrqISCZU0EVEMqGCLiKSCRV0EZFMqKCLiGQiVNDNbJOZHTSzfjN7dJb7l5nZv5rZfjM7YGZfKj5VERGZy7ydomZWAh4H7gLKwD4z2+nutS2ZXwNec/d7zWwVcNDM/t7dL9gWVxqbpudgrPvrO393ZygO4OlPbwjHAjy54W/Dsf80vDEcu7bjeDj2R+/dHI49OBTv3AMYfW9JOPZDP4o3Dtv4JO39vwjFrt55Q3i/Tx36XDgW4PsfvyUce2yoJxx73erhcOzb5Xh7vp2Nteef1z4cfxN9/b5Y+/jw1DSciP3uXb873kk58kY8FuDwR68Kxx6LvdQAWJYwG0TbmXh7fs+htGkpbHwyHjyQ8APOIfJq2Qj0u/uhaoHeAWyui3Gg28wMuAo4AcQnahARkUsWKehrgIGa2+Xqtlp/BdwMDAKvAA+7e9pMPSIickki77Fnm5Gm/n3K7wAvAb8NfAT4oZk95+4jH9iR2RZgC0BnW/ztrywMx/r3cuzNvQBMV841ORtJ8c7hvbxz5KcATGnsshUp6GWgr+Z2LzNn4rW+BPy5uzvQb2ZvAb8G/E9tkLtvB7YDLFtyXeJklNJsq9fdzup1twPw+ne+1eRsJMW1N9zOtTfMjN3Lux5rcjbSKJGPXPYBN5nZWjNrB+4DdtbFHAU+C2BmvwJ8FDhUZKIiIjK3ec/Q3X3KzB4CdgMl4Cl3P2BmW6v3bwO+CXzbzF5h5iOaR9x9qIF5i4hIndB1au6+C9hVt21bzfeDwOeLTU1ERFKoU1REJBMq6CIimVBBFxHJRNMWiaZSwc7FFsxd8bN40+npM1cnpfH7R78ejm1JWCG8siShryph8dmOk2l/g3tOxmOXvBtf1BZ3fCI2ft1vxRfA7TyZsGI7MDQcnwphWcLl1yNd14ZjexJ6olsS+6fbR+JX93aVo8fZoRKbJqDtnfjizKmdJV1DbeHYttPxNvrJ7oT9jsT3WzoePxYAJCwIXtQ13DpDFxHJhAq6iEgmVNBFRDKhgi4ikgkVdBGRTKigi4hkQgVdRCQTKugiIplQQRcRyYQKuohIJlTQRUQy0by5XCan8MF3Q6Fdy5eGd3t21VVJaXS/WQrHTnbH99s2GN+vJ/xZ7T4am4PjvEpbfJ6Y9mPxOVd8aprpk7G5LUqjK8P7nVyeNpfL8v745ChnV8bHpHswfpzHu+P77TmcMF8OMLayPRxbOnk6FJc0duFnh9autLFrmYzPd2QT8XEujY7HczgVf80znbbufXSuI4DKubTXxYXoDF1EJBMq6CIimVBBFxHJRKigm9kmMztoZv1m9ugFYu40s5fM7ICZPVtsmiIiMp95/1PUzErA48BdQBnYZ2Y73f21mpjlwBPAJnc/ambxVQdERKQQkTP0jUC/ux9y9wlgB7C5LuZ+4HvufhTA3Y8Vm6aIiMwnUtDXAAM1t8vVbbV+FbjazP7TzF4wsz8oKkEREYmJXIc+24XM9UvgtQIfBz4LdAH/bWZ73f2ND+zIbAuwBaDT4teWy8IwcOJFyif+F4BJ4tf6SvMNjLzMwOlXAI1dziIFvQz01dzuBQZniRly9zPAGTP7MXAL8IGC7u7bge0Ay0ori1oXVS6TvhW30rfiVgD2vPpEk7ORFH09G+jr2QDAT976myZnI40S+chlH3CTma01s3bgPmBnXcy/AJ82s1YzWwL8BvB6samKiMhc5j1Dd/cpM3sI2M1MJ/BT7n7AzLZW79/m7q+b2b8BLwMV4El3f3XOHbcY1hlrFbaRc6E4gJ4j8VZpgFMficd3DMf3m9LO33Eq3lLcMRxvgQZofzfe2mwnYu3gANbSQkt0/E6Nhvfb6Wlv3M7esDwc2/12/Nh5fMYElh2Pvz5bJtKmbujZfzIc69H2cTOsLTbrh09Mhp/fhuK5ApS6OpPio7y9LR57Nt5y72OJ7fnT8bH2yfg0AXMJjaq77wJ21W3bVnf7MeCxQrISEZFk6hQVEcmECrqISCZU0EVEMqGCLiKSCRV0EZFMqKCLiGRCBV1EJBMq6CIimVBBFxHJRKz/txEqHm5VbhlLWMU7tbX6SEJr83S8LX26I/63smMoof24Ne1vsE0mTBWQ0HbvlQqV4LjY2XhrPIkrx7ePxMevZSweS8IC7y1n46/PSndiu3vK+I0nzKJYKX5uPLOE+RIAEn6vfTphQKJTIAB++nQ4tpLY+m8dCa/lllI8do4SpzN0EZFMqKCLiGRCBV1EJBMq6CIimVBBFxHJhAq6iEgmVNBFRDKhgi4ikgkVdBGRTIQKupltMrODZtZvZo/OEXebmU2b2e8Vl6KIiETM2/pvZiXgceAuoAzsM7Od7v7aLHHfAnaHnz3aKpzQUtz+dtrK497ZHg9OyGNqWVc4tnQ6oaU4sWXbRs+GY5PaqxNWjidhFXaGR+KxQFvCz0drfKYLXxofPxuPr9heSsmXxJXmS8H2cXc8uCK9TcanS/CUKR6ASkKLvqW8hloa9MFDSns+4CnTNlTSpiy5kMhPvhHod/dD7j4B7AA2zxL3deC7wLFCMhMRkSSRgr4GGKi5Xa5u+yUzWwP8LrCtuNRERCRF5D3obJ8z1L/v/0vgEXefnmvGNTPbAmwB6LSlwRRloRg49xrlczOftE162sxz0lwDY69THnsdgEkSZmWURSVS0MtAX83tXmCwLuYTwI5qMV8J3GNmU+7+dG2Qu28HtgMsK60sfv5Oaai+rvX0da0HYM/xf2xyNpKir/Nm+jpvBmDP0D83ORtplEhB3wfcZGZrgbeB+4D7awPcfe35783s28AP6ou5iIg01rwF3d2nzOwhZq5eKQFPufsBM9tavV+fm4uILACh67jcfRewq27brIXc3R+89LRERCSVOkVFRDKhgi4ikgkVdBGRTMR7oQvm7vhErK248k68+dS6u5PySFqnPKH1v/XYe/HdtsXbmj1hpXQAT2iZnh4aStix4xOxlnc/ldDOn3AsgHi7O2DBdncAL/8inkN7fPqI6DE7z7o6w7HT752I79hj0zykrHRvKVNHAFaKn0/6VEIbfco4J8Qm88t/ZbbO0EVEMqGCLiKSCRV0EZFMqKCLiGRCBV1EJBMq6CIimVBBFxHJhAq6iEgmVNBFRDKhgi4ikgkVdBGRTDRtLhfc8anYXC7WviS828pIwrwhQEtwTgsA6+iI7zhhXovp4fgcKpYwbwgA4/G5X6w1YR6VCcJzVaTMw+EJ+QK0pMzFsTRhHduE/U4nzFXTkjC3DoCPngnHWnRemwqNmWck4XcJoDIRP8bhnw3wSsLPlnIcEuZyahadoYuIZEIFXUQkEyroIiKZCBV0M9tkZgfNrN/MHp3l/gfM7OXq1x4zu6X4VEVEZC7zFnQzKwGPA3cD64EvmNn6urC3gN9y9w3AN4HtRScqIiJzi5yhbwT63f2Qu08AO4DNtQHuvsfdT1Zv7gV6i01TRETmEynoa4CBmtvl6rYL+TLwzKUkJSIi6SLXoc928eWsF2+a2WeYKeifusD9W4AtAJ3Ery2XhaE83U+50g/AJGnXi0tzlSv9lCtvAhq7nEUKehnoq7ndCwzWB5nZBuBJ4G53n3WFZHffTvXz9R5bcflXUJVL0ltaR29pHQB7J/QmbDHpbVlHb0t17CZ3NzkbaZTIRy77gJvMbK2ZtQP3ATtrA8zseuB7wBfd/Y3i0xQRkfnMe4bu7lNm9hCwGygBT7n7ATPbWr1/G/AN4BrgCZtpj51y908UlWRKO7i1ps1m4OMT4djK8Kl4Hikt+hZvB6iMNfDtcmLrdni3Ke35Ddx3ZSg+xULKmKS0pVcmYtNd/FKDxmRBqMTHzhNiG6YR0yUULFT93H0XsKtu27aa778CfKXY1EREJIU6RUVEMqGCLiKSCRV0EZFMqKCLiGRCBV1EJBMq6CIimVBBFxHJhAq6iEgmVNBFRDJh3qR2VjM7DhwpeLcrgYT+7obKPZdbgRcL3mfux+xiaOzS5Z7Lh9191Wx3NK2gN4KZPV/kHDKXQrmkW0h5LpRcFkoe81lIeV7JuegjFxGRTKigi4hkIreCvpAWp1Yu6RZSngsll4WSx3wWUp5XbC5ZfYYuInIly+0MXUTkirVoCrqZbTKzg2bWb2aPznL/A2b2cvVrj5ndUnPfYTN7xcxeMrPnG5zHnWZ2qvpcL5nZN6KPbUAuf1STx6tmNm1mK6r3FXZMCsjzsoxdMBeNX3qeGr+FMn7uvuC/mFn67k3gRqAd2A+sr4u5A7i6+v3dwE9r7jsMrLxMedwJ/OBiHlt0LnXx9wL/UfQxWSxjp/HT+F0J47dYztA3Av3ufsjdJ4AdwObaAHff4+4nqzf3Ar3NyKNBjy1if18A/uESnu9iLZSxC+XSoMcWsT+Nn8ZvXouloK8BBmpul6vbLuTLwDM1tx34dzN7wcy2XIY8Pmlm+83sGTP7WOJji84FM1sCbAK+W7O5qGMyn4Uydim5aPzep/G7tFwu6/iFFoleAGyWbbNenmNmn2HmRfWpms2/6e6DZrYa+KGZ/czdf9ygPF5kpjV31MzuAZ4Gbgo+tuhczrsX+Im7n6jZVtQxmc9CGbtoLhq/D9L4XXwu51228VssZ+hloK/mdi8wWB9kZhuAJ4HN7v7e+e3uPlj99xjwfWbeLjUkD3cfcffR6ve7gDYzWxn9GYrMpcZ91L3dK/CYzGehjF0oF43f/6Pxu8hcaly+8Sv6Q/lGfDHzTuIQsJb3/wPiY3Ux1wP9wB1125cC3TXf7wE2NTCPa3n/+v6NwFFm/prP+9iic6nGLQNOAEsbcUwWy9hp/DR+V8L4Ff4CaNQXcA/wBjP/s/zH1W1bga3V758ETgIvVb+er26/sXqw9wMHzj+2gXk8VH2e/cz8B9Edcz22kblUbz8I7Kh7XKHHZLGMncZP45f7+KlTVEQkE4vlM3QREZmHCrqISCZU0EVEMqGCLiKSCRV0EZFMqKCLiGRCBV1EJBMq6CIimfg/e/Fuan2YJucAAAAASUVORK5CYII=\n",
            "text/plain": [
              "<Figure size 432x288 with 3 Axes>"
            ]
          },
          "metadata": {
            "needs_background": "light"
          },
          "output_type": "display_data"
        },
        {
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAACPCAYAAADqQSA9AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAPzElEQVR4nO3dfYxU53XH8e/ZV5aXXYwxtsNubOxQu44EaYuJ7aapXcsptmqhKlFkJ06bNClFqqv+UyuWqkaV8ldlqWrlOEHYtaKmUZFaHBenuDRylSYNBmMcwLx48ZoYGOPYLLAsLAvL7Jz+sYM9THZ3zoUZZvbZ30caMffeM8+cvc/s4c7d+zzX3B0REZn6muqdgIiIVIcKuohIIlTQRUQSoYIuIpIIFXQRkUSooIuIJKKlXm88f16z39jTWq+3l8u066Ax87o5odiPtQ/WOBvJYufbTbRf2xmKvbVjoLbJCAC7h+aFY8+99W6/u18z3ra6FfQbe1p5ZVNPvd5eLtO8u0/y6ac/H4p9fvGmGmcjWXR++hRLvvXHodj/W/JcjbMRgNs2PxKO7f3s3x6caJtOuYiIJEIFXUQkESroIiKJUEEXEUmECrqISCJU0EVEEqGCLiKSCBV0EZFEhAq6ma0ws14z6zOzx8fZ3mVmL5jZTjPbY2ZfqX6qIiIymYojRc2sGXgKuA/IAdvMbIO77y0J+3Ngr7s/aGbXAL1m9n13H5mo3X3DV3Hnzs+Gknx56fpQnFw5IyfaeXv9zbHgXzkEkHoqDLQy9MJ1seAltc0lZTf9+5+FY2/4z9FwbO8k2yJH6MuBPnc/UCzQ64CVZTEOzDEzA2YDx4F8OEMREblskYK+EDhcspwrriv1LeDXgSPA68BfunuhKhmKiEhIZHIuG2dd+Z2lfx/YAfwecDPwIzP7qbtfNM2ema0CVgG0LYjN9iaNY+33TvLMv4x1aX54Zp2zkSzUd9ND5Ag9B5ROi9jN2JF4qa8Az/mYPuAXwK3lDbn7Wndf5u7LWrr0oZpqVn2pi1c29fDKph5aOmbVOx3JQH03PUQK+jZgsZktMrM24CFgQ1nMIeBeADO7FrgFOFDNREVEZHIVT7m4e97MHgU2Ac3As+6+x8xWF7evAb4JfNfMXmfsFM3X3b2/hnmLiEiZ0A0u3H0jsLFs3ZqS50eAz1Q3NRERyUIjRUVEEqGCLiKSCBV0EZFE1O0m0X6ihfPrF8SCl9Y2Fxnz5IkbwrEtw28wf9fZUOwTx4NTBACPzXsrHCsX+4P994fiWk/9G9f/+Fgo9uaPrQ6//3ce/KdwLMBnZp7PFB+1ZqB83OPEfvDub4Rj39wbbxeg+6Xy4ToTm7n/aKa2J6IjdBGRRKigi4gkQgVdRCQRKugiIolQQRcRSYQKuohIIlTQRUQSoYIuIpIIFXQRkUSooIuIJEIFXUQkEXWby8UK0DIcm+sgy9wMq+e+c6kpTXv/uPOecKzl99LafyYUu3ZjfKr8xx75Tjh2Osjy2e/demMs0IHz+VDo7IPxY76/2v25cCzArfPfD8eeybeFYw8PzA3Hjrx2VTi2Kzb9zQfaTsbmOhpLpDrz2ugIXUQkESroIiKJUEEXEUlEqKCb2Qoz6zWzPjN7fIKYu81sh5ntMbP/rW6aIiJSScU/ippZM/AUcB+QA7aZ2QZ331sSMxf4NrDC3Q+ZWfDOFSIiUi2RI/TlQJ+7H3D3EWAdsLIs5gvAc+5+CMDd43++FhGRqogU9IXA4ZLlXHFdqV8DrjKzH5vZdjP7o2olKCIiMZHr0G2cdeUXkLcAvwXcC3QAL5vZFnfff1FDZquAVQBts+LXf0pjGHxpG6de2gZAU36oztlIFoNbXmZw6xYAmkZj4wdk6okU9BzQU7LcDRwZJ6bf3YeAITP7CWO3dr6ooLv7WmAtwKyre+J3UJWG0Hnv7XTeezsAx/7kqTpnI1l03nEnnXfcCUD/3zxR52ykViKnXLYBi81skZm1AQ8BG8pi/gP4HTNrMbOZwCeBfdVNVUREJlPxCN3d82b2KLAJaAaedfc9Zra6uH2Nu+8zs/8CdgEF4Bl33z1Zu015Z8bx0VCS3zv0yVAcwOq5z4Vjp4O/P35TOLZ998x4wwXHgsOVm4fHO2s3vqWvPBzPAXh66T+HY5e3t4ZjzxRGwrF9+UI4tpV4LMCT++4Ox3YeiMX1FwrY0HAodu5bsSkCAPo7sp1G3UM8/uyC+H6bdTg+vKazP36ioOut2D67oOlsfDh/YfBUprYnEprLxd03AhvL1q0pW34C0Hc5EZE60UhREZFEqKCLiCRCBV1EJBEq6CIiiVBBFxFJhAq6iEgiVNBFRBKhgi4ikggVdBGRRIRGitZC0/kCHe/Ehru+u/3acLu35z+fKY8NS54Nx17fMjtT21HfP3V1OPaF/qWZ2t66f1E4dl6GYdAURrHB06HQBT+P/3zHh7MNH/9T4jM13/GRt8OxB0/PC8fmMtxlfuhohukVgI5cfLqCGSeCw+MLo/ipWN81nY9/Jq7dHp8uAeBcV3M4tnlfPI+RWfEc5hw6F47NMpQfoPmd/nDs6Ei2fTcRHaGLiCRCBV1EJBEq6CIiiVBBFxFJhAq6iEgiVNBFRBKhgi4ikggVdBGRRKigi4gkIlTQzWyFmfWaWZ+ZPT5J3O1mNmpmn6teiiIiElFx6L+ZNQNPAfcBOWCbmW1w973jxP0dsCn0zvk8TUcHQqFdb8aHgx/riA/ZBni47ZFw7NBIWzi2kGEU/ay2+JDid/bFp0EAuKrXwrHX/Dw2HBzgQMHxs2dDsTPeiw+vnplhODjA2Z/GPxs/a43HnpsX78CmDCPCu96L9wdAx9H43e67dhwNxXnBKQzH+q6jLz583Wd3hGMBRtvjU2m0nsqHY+dsPRLPofuacGzT4HA4FqAwdCYc6+fjP99kIkfoy4E+dz/g7iPAOmDlOHF/AawH3q9KZiIikkmkoC8EDpcs54rrPmBmC4E/BNZULzUREckiMtvieN8Ry7+P/gPwdXcfNZv4K6WZrQJWAcxors3MhVI7ufdfJff+dgDOe7avn1Jfhwd2kBvYCcB5j51ukaknUtBzQE/JcjdQfpJqGbCuWMznAw+YWd7dny8Ncve1wFqArrYFGc4ySyPoXrCM7gXLANi6/ck6ZyNZ9Mz9BD1zPwHA5v1P1zcZqZlIQd8GLDazRcA7wEPAF0oD3P2DSbfN7LvAD8uLuYiI1FbFgu7ueTN7lLGrV5qBZ919j5mtLm7XeXMRkQYQumORu28ENpatG7eQu/uXLz8tERHJSiNFRUQSoYIuIpIIFXQRkUSEzqHXhDt+PjZmuv1kfPhzZ1+2oePvDXwkHHvu2tFwbNNwfIj36MH4/6tdGe7CDjBjIMMQ9oGheMOFAn42NqS/5Vh8SoE5uWwfybbT8fiz8+L7eUZ8xDttp+L7uG0w2xDv1jPxeDsTu77czLDW4H4bGAy/f1OW+S6A2XsyDI1vjvdd4eq54dimvlw4NiubmWEqhNPx35HJ6AhdRCQRKugiIolQQRcRSYQKuohIIlTQRUQSoYIuIpIIFXQRkUSooIuIJEIFXUQkESroIiKJUEEXEUlE/eZyaW7GOueEQufsPhpv9uzVmdJoH4zP/dK2K8O8KOfi88+Yx+eImfHLDPOtAOTjeXD0WDy2pZmmuV2hUD8+EG62PTrHSFHrYFs4ds6bsbmDALw1w5xATfHjouZjp+LtAt4UnxOocDI474oZ1t4ebDOeb5Nnm8vFZs+KB5+Kf+4tH5//xrPEBucu+vAFGX73Mu67iegIXUQkESroIiKJUEEXEUlEqKCb2Qoz6zWzPjN7fJztXzSzXcXHZjNbWv1URURkMhULupk1A08B9wO3AQ+b2W1lYb8AftfdlwDfBNZWO1EREZlc5Ah9OdDn7gfcfQRYB6wsDXD3ze5+ori4BeiubpoiIlJJpKAvBA6XLOeK6ybyVeDFy0lKRESyi1z0O96FsONeNGlm9zBW0D81wfZVwCqAGS2xa9ClcRw+9wa5c70AjBC7f6U0htxIL7mRC32X8XpqmTIiBT0H9JQsdwNHyoPMbAnwDHC/u487QsXd11I8v94147rqXEkvV0xP+630tN8KwMunN9Q5G8miu+0WuttuAWDL0At1zkZqJXLKZRuw2MwWmVkb8BBw0W+zmX0UeA74krvvr36aIiJSScUjdHfPm9mjwCagGXjW3feY2eri9jXAN4CrgW+bGUDe3ZdN2vDIeQpHfhlK0jo6QnEAszIMEQaY2RU/9TPaFc+DDEO2m49nyDnDMHoAa4kPpc8Png7HeiHPaH9sqoAsOXBuJB4LNGVpOwO7/pp47PGT4Vg/M5wpj8KZM/G2g8PY3UcZPXGiciCA1WDqgQtOxz/3TXNmZ2s7yIfj/eGj8Sk6AJq7OsOxowPxz9BkQr8N7r4R2Fi2bk3J868BX6tKRiIickk0UlREJBEq6CIiiVBBFxFJhAq6iEgiVNBFRBKhgi4ikggVdBGRRKigi4gkQgVdRCQR5lW623TmNzY7ChyscrPzgf4qt3mpUs/lN4HXqtxm6vvsUqjvsks9lxvcfdy5KepW0GvBzF6tOIfMFaJcsmukPBsll0bJo5JGynM656JTLiIiiVBBFxFJRGoFvZFuTq1csmukPBsll0bJo5JGynPa5pLUOXQRkekstSN0EZFpa8oUdDNbYWa9ZtZnZo+Ps/2LZrar+NhsZktLtr1tZq+b2Q4ze7XGedxtZieL77XDzL4RfW0NcnmsJI/dZjZqZvOK26q2T6qQ5xXpu2Au6r/sear/GqX/3L3hH4zd+u4t4CagDdgJ3FYWcxdwVfH5/cDWkm1vA/OvUB53Az+8lNdWO5ey+AeB/6n2Ppkqfaf+U/9Nh/6bKkfoy4E+dz/g7iPAOmBlaYC7b3b3CzdK3AJ01yOPGr22Gu09DPzrZbzfpWqUvgvlUqPXVqM99Z/6r6KpUtAXAodLlnPFdRP5KvBiybID/21m281s1RXI404z22lmL5rZxzO+ttq5YGYzgRXA+pLV1donlTRK32XJRf33IfXf5eVyRfuvNrdMr77xbj0+7uU5ZnYPYx+qT5Ws/m13P2JmC4Afmdkb7v6TGuXxGmNDc0+b2QPA88Di4GurncsFDwI/c/fjJeuqtU8qaZS+i+ai/ruY+u/Sc7ngivXfVDlCzwE9JcvdwJHyIDNbAjwDrHT3YxfWu/uR4r/vAz9g7OtSTfJw90F3P118vhFoNbP50Z+hmrmUeIiyr3tV3CeVNErfhXJR//0K9d8l5lLiyvVftU/K1+LB2DeJA8AiPvwDxMfLYj4K9AF3la2fBcwpeb4ZWFHDPK7jw+v7lwOHGPvfvOJrq51LMa4LOA7MqsU+mSp9p/5T/02H/qv6B6BWD+ABYD9jf1n+6+K61cDq4vNngBPAjuLj1eL6m4o7eyew58Jra5jHo8X32cnYH4jumuy1tcyluPxlYF3Z66q6T6ZK36n/1H+p959GioqIJGKqnEMXEZEKVNBFRBKhgi4ikggVdBGRRKigi4gkQgVdRCQRKugiIolQQRcRScT/A60dEUOvo51lAAAAAElFTkSuQmCC\n",
            "text/plain": [
              "<Figure size 432x288 with 3 Axes>"
            ]
          },
          "metadata": {
            "needs_background": "light"
          },
          "output_type": "display_data"
        },
        {
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAACPCAYAAADqQSA9AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAANwklEQVR4nO3dfWxV933H8ffHBsJDEgJxSBMwDelYM7rBlBISSKcmrdqQTBWqtEokXaZErTy2sU6bVAVtVf/JtGmqNk3bUjGEokpNW7oqWUMjMpo2XZsOCBDCQ4CQOiYYlw7C8zPG5rs/fFlubm3f3zX3cq9//rwki3vO+d5zvtyf9fHx8XlQRGBmZsNfU70bMDOz6nCgm5llwoFuZpYJB7qZWSYc6GZmmXCgm5llYlS9NtwyuTluax1dr83bFXpn/0U8fsOTx254e237hcMRcVN/y+oW6Le1jmbj2tZ6bd6u0LwH9nv8himP3fDWfEv7voGW+ZCLmVkmHOhmZplwoJuZZcKBbmaWCQe6mVkmHOhmZplwoJuZZcKBbmaWiaRAl7RQ0h5J7ZKW9bN8oqQfSNomaaekx6vfqpmZDaZsoEtqBp4CHgRmAQ9LmlVS9mfAroiYA9wH/KOkMVXu1czMBpGyhz4PaI+IjojoBlYBi0pqArhOkoBrgaNAT1U7NTOzQaUE+lRgf9F0V2FesX8Dfgs4AOwA/iIiLlWlQzMzS5Jycy71M6/0ydIPAFuBTwAfAl6S9EpEnHzfiqQ2oA1g+tS63RfMhmjFN0+w8pm+IX33SG+du7FKeOxGhpQ99C6g+NZs0+jbEy/2OPBc9GkH9gJ3lK4oIlZExNyImHvTjc1D7dnqpO3RiWxc28rGta14/IYXj93IkBLom4CZkmYU/tC5GFhdUtMJfBJA0s3Ah4GOajZqZmaDK3vcIyJ6JC0F1gLNwNMRsVPSksLy5cCTwDck7aDvEM0TEXG4hn2bmVmJpAPZEbEGWFMyb3nR6wPAp6vbmpmZVcJXipqZZcKBbmaWCQe6mVkmfDK42QjTeXECXzpwV1Ltv9y6qcbdGMCcjQ9XUP3kgEu8h25mlgkHuplZJhzoZmaZcKCbmWXCgW5mlgkHuplZJhzoZmaZcKCbmWXCgW5mlgkHuplZJhzoZmaZ8L1cbEg6uq9j8d5PJNWumvFyjbuxSpw4N44f7JidVOt7uVwd+vGkqqzHe+hmZplwoJuZZcKBbmaWiaRAl7RQ0h5J7ZKWDVBzn6StknZK+ml12zQzs3LK/lFUUjPwFPApoAvYJGl1ROwqqrkB+DqwMCI6JU2pUb9mZjaAlD30eUB7RHRERDewClhUUvMI8FxEdAJExKHqtmlmZuWkBPpUYH/RdFdhXrHfBCZJ+m9Jr0n6o2o1aGZmaVLOQ1c/86Kf9XwU+CQwDlgvaUNEvPW+FUltQBvA9Kk+BX64WfHNE6x85iQAF46NrnM3Vonises9Na7O3VitpKRqF9BaND0NONBPzeGIOAOckfQzYA7wvkCPiBXACoC5c8aW/lCwBtf26ETaHp0IQMv9x+vbjFWkeOzGzT9X526sVlIOuWwCZkqaIWkMsBhYXVLzPPB7kkZJGg/cDeyubqtmZjaYsnvoEdEjaSmwFmgGno6InZKWFJYvj4jdkv4L2A5cAlZGxBu1bNzq68yZsWzcMjOt2Jf+NxR1i2s6r6l3G8POr3pOV1R/70++lFw7/e2LlbbTr6QD2RGxBlhTMm95yfTXgK9VpSszM6uYrxQ1M8uEA93MLBMOdDOzTDjQzcwy4UA3M8uEA93MLBMOdDOzTDjQzcwy4UA3M8tE3W55uO/iBP64a35S7b9PW1/jbqxSTRdh/C+b692GFTl96XxS3ahzMHnXpaTarxz6neTt/+2UHcm1jeI/Tk9Mrv3rTY9XtO7Jr6TfXmH8pvaK1j0Q76GbmWXCgW5mlgkHuplZJhzoZmaZcKCbmWXCgW5mlgkHuplZJhzoZmaZcKCbmWUiKdAlLZS0R1K7pGWD1N0lqVfSH1SvRTMzS1H20n9JzcBTwKeALmCTpNURsaufun8A1qZs+OSp8fzo53PSulzsS/+vht3dZ5Nrm883MWlPb1LtD8+OTl7vp8dX5+nnI9G9m9MuTW8+v5KJe04m1X7r1XuSt7955vTkWoA/bf1Jcu2HRh9Jrv3WsbuTa7/70wXJtdccqeyARnN3JNdemjYlfcWHBl6U0uE8oD0iOiKiG1gFLOqn7s+BZwffnJmZ1UpKoE8F9hdNdxXm/T9JU4HPAsur15qZmVUi5W6L6mde6e8S/ww8ERG9Un/lhRVJbUAbQPOkSYktWqP43rfP8L1v9x2WuXhhQp27sUocXbuFo2u39E30pB9as+ElJdC7gNai6WnAgZKaucCqQpi3AA9J6omI7xcXRcQKYAXANa2t6QeYrCF87pEJfO6RviC/e45PkBpOJj9wJ5MfuBOAzsUr69yN1UpKoG8CZkqaAfwSWAw8UlwQETMuv5b0DeCF0jA3M7PaKhvoEdEjaSl9Z680A09HxE5JSwrLfdzczKwBJD2xKCLWAGtK5vUb5BHx2JW3ZWZmlfKBUDOzTDjQzcwy4UA3M8tE0jH0WmjqgbHv+udJI/n95/8qubap558Y++6FpNo/Wf+Hyev9+cf/NbkW4JZR11ZUX29vXzxdUf2X9302uVY/Tru2I5rEpTHNSbU3bkyPiF+cbS1fVOQrLzyWXHvqjvRbQozrTL/VREtn+tnT4w5XdluK8fvSbq8AEHv3ly9K4EQ1M8uEA93MLBMOdDOzTDjQzcwy4UA3M8uEA93MLBMOdDOzTDjQzcwy4UA3M8uEA93MLBMOdDOzTNTtXi5jjvfwwecPJ9XO+mj6vUB2LXhmqC0NC08c/N2K6p99+Z7k2klvDfw82FLHzl5g1Ou/SKr9wOrfTl7vgnN/mVwL0HLrieTa2S2lT04c2F3X702uff309OTaH+25I7kW4Pr145Jrrz3Ym1Z49jx6fU9S6ZTTtyVvf+zxyp4TfHZK+n1Ubn0pfd9zzKkK7rkS6T2MfyP9+weg91BavgE0XV/BPYnODLKe9LWYmVkjc6CbmWXCgW5mlomkQJe0UNIeSe2SlvWz/POSthe+1kmaU/1WzcxsMGUDXVIz8BTwIDALeFjSrJKyvcDHI2I28CSwotqNmpnZ4FL20OcB7RHRERHdwCpgUXFBRKyLiGOFyQ3AtOq2aWZm5aQE+lSg+PlIXYV5A/kC8OKVNGVmZpVLOQ+9v5OT+z15U9L99AX6xwZY3ga0AYwdfX1ii9YoDu9az5HdGwCIOF/nbqwSBzs2cHDvqwD0euyylRLoXUDx01+nAb92hr2k2cBK4MGIONLfiiJiBYXj6xPH3ZJ+Rr81hJZZ82mZNR+AjhV/V+durBI3334PN9/ed5HZG8/+fZ27sVpJOeSyCZgpaYakMcBiYHVxgaTpwHPAoxHxVvXbNDOzcsruoUdEj6SlwFqgGXg6InZKWlJYvhz4KnAj8HVJAD0RMXfQFUvEqLTT4Cd/Z0JSHcCMQ23JtQC/8eFfJdfu/d+W5FrtH5tc21TBlcqTd1f2i82US+m1N7x2MLm2o0loXNr/8YZX3kle7+gKLqMHOHvzjcm125Ve++bpj1TUR6rbjvRUVN908Vxy7eg3u5LqJKExo9Nqz11I3v71bx5PrgW4tnNMcm3T8UGudy+h02eTay+dPJVe25t4a4WCponXpRd3VxACg0i6l0tErAHWlMxbXvT6i8AXq9KRmZkNia8UNTPLhAPdzCwTDnQzs0w40M3MMuFANzPLhAPdzCwTDnQzs0w40M3MMuFANzPLhKKCp15XdcPSu8C+Kq+2BUh/1HZt5d7LncCWKq8z989sKDx2lcu9lw9GxE39LahboNeCpM1l7yFzlbiXyjVSn43SS6P0UU4j9TmSe/EhFzOzTDjQzcwykVugN9LDqd1L5Rqpz0bppVH6KKeR+hyxvWR1DN3MbCTLbQ/dzGzEGjaBLmmhpD2S2iUt62f55yVtL3ytkzSnaNk7knZI2ippc437uE/SicK2tkr6aup7a9DLl4v6eENSr6TJhWVV+0yq0OdVGbvEXjx+lffp8WuU8YuIhv+i79F3bwO3A2OAbcCskpoFwKTC6weBV4uWvQO0XKU+7gNeGMp7q91LSf1ngJer/ZkMl7Hz+Hn8RsL4DZc99HlAe0R0REQ3sApYVFwQEesi4lhhcgMwrR591Oi91Vjfw8B3rmB7Q9UoY5fUS43eW431efw8fmUNl0CfCuwvmu4qzBvIF4AXi6YD+KGk1yRV9hTpofUxX9I2SS9Kuvy04Ur/D9XqBUnjgYXAs0Wzq/WZlNMoY1dJLx6/93j8rqyXqzp+SQ+JbgDqZ16/p+dIup++b6qPFc2+NyIOSJoCvCTpzYj4WY362ELfpbmnJT0EfB+Ymfjeavdy2WeA/4mIo0XzqvWZlNMoY5fai8fv/Tx+Q+/lsqs2fsNlD70LaC2angYcKC2SNBtYCSyKiCOX50fEgcK/h4D/pO/XpZr0EREnI+J04fUaYLSkltT/QzV7KbKYkl/3qviZlNMoY5fUi8fv13j8hthLkas3ftU+KF+LL/p+k+gAZvDeHyA+UlIzHWgHFpTMnwBcV/R6HbCwhn18gPfO758HdNL307zse6vdS6FuInAUmFCLz2S4jJ3Hz+M3Esav6t8AtfoCHgLeou8vy39TmLcEWFJ4vRI4BmwtfG0uzL+98GFvA3Zefm8N+1ha2M42+v5AtGCw99ayl8L0Y8CqkvdV9TMZLmPn8fP45T5+vlLUzCwTw+UYupmZleFANzPLhAPdzCwTDnQzs0w40M3MMuFANzPLhAPdzCwTDnQzs0z8H4LGRcutPdEkAAAAAElFTkSuQmCC\n",
            "text/plain": [
              "<Figure size 432x288 with 3 Axes>"
            ]
          },
          "metadata": {
            "needs_background": "light"
          },
          "output_type": "display_data"
        }
      ],
      "source": [
        "name = \"same_but_larger_mean\"\n",
        "model = \"CIR\"\n",
        "mx = [120, 120]    # get prob dist for this domain ; bigger = more accurate\n",
        "K_avg = 20\n",
        "beta, gamma = 1, 1.5\n",
        "\n",
        "get_BF_plot(name, model, K_avg, beta, gamma, mx)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "SSCjROlC5Wan",
        "outputId": "95d1dd90-dc51-4d90-e1c6-7c6a891ad225"
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Computing log Bayes factors on a  10  x  10  grid...\n",
            "Row  1 / 10 ...\n",
            "Row  2 / 10 ...\n",
            "Row  3 / 10 ...\n",
            "Row  4 / 10 ...\n",
            "Row  5 / 10 ...\n",
            "Row  6 / 10 ...\n",
            "Row  7 / 10 ...\n",
            "Row  8 / 10 ...\n",
            "Row  9 / 10 ...\n",
            "Row  10 / 10 ...\n",
            "--- 11.113316774368286 seconds ---\n"
          ]
        },
        {
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAACPCAYAAADqQSA9AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAPi0lEQVR4nO3dXYxc9XnH8e9vZ3a9fieElyKvU5lgSgkFqU2dlqYtaUpjkJBVqRckqFWiRK6lUvWqClKl3ESqUuWmFxBZFkKRqqruRVLqpiYkbdWmCpAaIt7sYFgMwYOT4jewvV57X+bpxQ4wbNfe52/PeHf//n2klXfOPHPm2fM/+/jM2fP8jyICMzNb+gYWOgEzM+sNF3Qzs0q4oJuZVcIF3cysEi7oZmaVcEE3M6tEc6HeeGhwZQwvuyIVO7Em///O8JqzRXmsbubjlysf21D+ctCfTaxNx46PLUvHziSSz2PwhNKxE8cPM6yVqVgN53NuD5UdY0wtz+fcbhStOm1guiB2smzdUZBzc6ydijszfozlzdw+116eLxElY1Gs4Opq5TYDAANT+djGmYKBBkL57aHJ/I5xYvLwkYi4eq7nFqygDy+7gk/c8qep2IOfWZ1e742//2pRHr/z4VfSsbct/2k6dvXAmXTs1w7enY594akb0rEAU1fk99jr/j1fPfb941/zG4ObU7G6KZ/z6ZFV6ViAI7+S34Un1pZUhXzo0Dv54FWtsr6Ps2vz677mmbFU3J49D3L7uvtSsWM3X5t+/8O3DaZjSzUmCmLH87ErjuSr/5qX3s6vGIhl+e3RaB1Ox3730IPnLEQ+5WJmVgkXdDOzSrigm5lVwgXdzKwSLuhmZpVwQTczq4QLuplZJVzQzcwqkSrokjZL2i9pVNIDczy/VtK/SHpO0l5JX+h9qmZmdj7zttlJagAPAXcCLWCPpF0Rsa8r7M+AfRFxj6Srgf2S/j4iztnfpYkpmoeOpZJc+2quxRzg9Xc+mo4F2DtyfTp2oGBWgenhfOyyY/luwDVHyzoNlx/Lfwhb++Qb+RVHENO5VuiBN/NdcCvPlvXGD53Md5ZODec7YaeX57ebCtrHm+Nl7eNrXst3Mjb3vZ4LnJqiffhoKnR5a0X6/a8bz8cCTC/Lj0dzLL+RBybz27hxIt/RrZOn07EAcSZfMKZPnCha97lk9tpNwGhEHOgU6J3AllkxAayWJGAVcAwo2M3NzOxiZQr6OuBg1+NWZ1m3B4FfBg4BLwB/EREFU+SYmdnFysxsNNf5gNmf+z8DPAv8HvBR4PuS/jsiPvA5QtJWYCvAcCM/4ZYtDgdPvcjBsRcBmKRsVktbWAfPvkTr7H4AJiJ/msGWlkxBbwHrux6PMHMk3u0LwNciIoBRSa8BNwH/0x0UETuAHQBrh64tOxlsC279qltYv+oWAH7Y+rsFzsZKrF92E+uX3QTAk+88urDJWN9kTrnsATZK2iBpCLgX2DUr5g3g0wCSrgV+CTjQy0TNzOz85j1Cj4gpSfcDjwMN4JGI2CtpW+f57cBXgW9KeoGZUzRfjogjfczbzMxmSd0dICJ2A7tnLdve9f0h4A96m5qZmZVwp6iZWSVc0M3MKuGCbmZWiQW7STQRxNncnV/X7j+ZX63Krm8fznVAAzB4On+l5fSyfDt/YyLfg9Vult1Zfc2+4/ngdmEvWDvXYh0n8+M3MLysKIXB6XzOg4P53b09nI/VdH6/aA/l290Bmj/Lj18MDhWtO2PgaL4lvVlwU2SA5sn8nZ8bPy/YDqvzUxDoeEHLfaNs7Npvv5OOjcmCu2Cfh4/Qzcwq4YJuZlYJF3Qzs0q4oJuZVcIF3cysEi7oZmaVcEE3M6uEC7qZWSVc0M3MKuGCbmZWCRd0M7NKLOhcLiTnLxg4MZ5e7Yr/HS7LY6BgzpXx3NwlM+vNh2oqPxdI88ip/IoBGvlE2ifyc64gULP3u0+cGiuKV7Ngfo2T+XUPrMrPB0IUjN9Yfl8ulZ0zJ9pt2qdPp2K1Ir8dBt4uGzveLpijKTlvEAAF+3E/74OpgnmJPJeLmZl9gAu6mVklXNDNzCqRKuiSNkvaL2lU0gPniLlD0rOS9kr6r96maWZm85n3r1qSGsBDwJ1AC9gjaVdE7OuKuQL4BrA5It6QdE2f8jUzs3PIHKFvAkYj4kBETAA7gS2zYj4HfDsi3gCIiLd6m6aZmc0nU9DXAQe7Hrc6y7rdCHxI0n9KekbSn/QqQTMzy8lcSDzXhdqzL99sAr8GfBpYDjwp6amIePkDK5K2AlsBhgdWlmdrC6o1+TKtyZkhnYyzC5yNlWhNvUJr6hUAJvHY1SpT0FvA+q7HI8ChOWKORMQYMCbpB8BtwAcKekTsAHYArG1e3c9r+q0PRgZvZGTwRgCeHNu1wNlYiZHmRkaaGwF4cvxfFzgb65fMKZc9wEZJGyQNAfcCs3+b/xn4bUlNSSuATwA/6W2qZmZ2PvMeoUfElKT7gceBBvBIROyVtK3z/PaI+Imk7wLPA23g4Yh48bwrbreJM7mPfgOnz6TiAJa/9PN0LMD0dVemY9vN/GX7jbH+fKzVibLW/5gqaJmeLogNiGx8o6A9f7ywNb5g6gYNDeVjC6YJKPn5YizXcv+ewfz0ClEwBUF2uoLsdAIAnMn/ngJoxfJ8bDM/djE5lY8t2ufLTipkp1fopdTeEhG7gd2zlm2f9fjrwNd7l5qZmZVwp6iZWSVc0M3MKuGCbmZWCRd0M7NKuKCbmVXCBd3MrBIu6GZmlXBBNzOrhAu6mVklen/b9j5oHz2WjtXKgru1AwPH8630AwXt4CVt5kzlW5XbBXc0L9UubN1Ot0KXtFc3y3bJGM/n3D5VMH4FecTEZH69jcJjqIJ9IyYK7hyv3JQJJfuEptv59y/UTk4TAkAU5FEyLUXJfrxAfIRuZlYJF3Qzs0q4oJuZVcIF3cysEi7oZmaVcEE3M6uEC7qZWSVc0M3MKuGCbmZWiVRBl7RZ0n5Jo5IeOE/cr0ualvRHvUvRzMwy5u1vltQAHgLuBFrAHkm7ImLfHHF/AzyeeeOIoH02187bWLsmFddJJB8L8PaJdGhJm3lEfkoBDQ3mYwvuAg/QHhvPr7uk7T7f7V50Z/Wiu7ADKmjdjoI2+iIlreaFH4pLWt7VTO5HJTMEDBZMYTFQ9rvXHjudX3XB70j7bMFUDKrrJEXmp9kEjEbEgYiYAHYCW+aI+3PgW8BbPczPzMySMgV9HXCw63Grs+w9ktYBfwhs711qZmZWIvMZe67PUbOn2ftb4MsRMa3znPKQtBXYCjBM2ayItvBa7VFa7VcBmKRg9jtbcK3pUVrtUcBjV7NMQW8B67sejwCHZsV8HNjZKeZXAXdLmoqIR7uDImIHsANgja5Mzr1qi8XIwA2MDNwAwFOTqT+V2CIx0riBkUZn7CYeW+BsrF8yBX0PsFHSBuBN4F7gc90BEbHh3e8lfRP4zuxibmZm/TVvQY+IKUn3M3P1SgN4JCL2StrWed7nzc3MFoHUdWoRsRvYPWvZnIU8Ij5/8WmZmVmpui7CNDO7jLmgm5lVwgXdzKwSZX3kvZa8a/z0iXwbfWPVysIU8ldPlrTdx2TJ3drzrcpFdz+nMOeJPt3VvGQbl0w/AES7ZPzybewl41eiZKyL110ybUJyTGIyP09A0TQBlE3b0O7Tdiv5+ZYCH6GbmVXCBd3MrBIu6GZmlXBBNzOrhAu6mVklXNDNzCrhgm5mVgkXdDOzSrigm5lVwgXdzKwSLuhmZpVY2LlcznP/0QvVHj9TFF8y/0W/5nIh2vkcCua/KM5jESiaj6R43SXB+TFBBcdFJetdYornRSn5/S/axpfv3S19hG5mVgkXdDOzSrigm5lVIlXQJW2WtF/SqKQH5nj+PknPd76ekHRb71M1M7PzmbegS2oADwF3ATcDn5V086yw14DfjYhbga8CO3qdqJmZnV/mCH0TMBoRByJiAtgJbOkOiIgnIuJ45+FTwEhv0zQzs/lkCvo64GDX41Zn2bl8EXjsYpIyM7NymQur57pYdM4LPSV9ipmC/slzPL8V2AowzIpkirZYtOIAb8YBACYpu7epLSyP3eUhU9BbwPquxyPAodlBkm4FHgbuioijc60oInbQOb++Rldevlf/L1Ejup4RXQ/Aj9r/tsDZWAmP3eUhc8plD7BR0gZJQ8C9wK7uAEkfAb4N/HFEvNz7NM3MbD7zHqFHxJSk+4HHgQbwSETslbSt8/x24CvAh4FvaKaddyoiPj7vu2dbdAt6tqOwtVrNwfy6Jwpam0vaj/swBcJ72v1rpe+LfrZtF/X+l+hjO/9SamPv535c8ZQJvZSanCQidgO7Zy3b3vX9l4Av9TY1MzMr4U5RM7NKuKCbmVXCBd3MrBIu6GZmlXBBNzOrhAu6mVklXNDNzCrhgm5mVgkXdDOzSigWqLVY0mHgpz1e7VXAkR6v80LVnsuvAj/u8Tpr32YXwmNXrvZcfjEirp7riQUr6P0g6enUHDKXgHMpt5jyXCy5LJY85rOY8rycc/EpFzOzSrigm5lVoraCvphuTu1cyi2mPBdLLoslj/kspjwv21yqOoduZnY5q+0I3czssrVkCrqkzZL2SxqV9MAcz98n6fnO1xOSbut67nVJL0h6VtLTfc7jDknvdN7rWUlfyb62D7n8ZVceL0qalnRl57mebZMe5HlJxi6Zi8evPE+P32IZv4hY9F/M3PruVeB6YAh4Drh5VsztwIc6398F/KjrudeBqy5RHncA37mQ1/Y6l1nx9wD/0ettslTGzuPn8bscxm+pHKFvAkYj4kBETAA7gS3dARHxREQc7zx8ChhZiDz69NperO+zwD9cxPtdqMUydqlc+vTaXqzP4+fxm9dSKejrgINdj1udZefyReCxrscBfE/SM5K2XoI8flPSc5Iek/Sxwtf2OhckrQA2A9/qWtyrbTKfxTJ2Jbl4/N7n8bu4XC7p+KVuEr0IzHU78Tkvz5H0KWZ2qk92Lf6tiDgk6Rrg+5Jeiogf9CmPHzPTmntK0t3Ao8DG5Gt7ncu77gF+GBHHupb1apvMZ7GMXTYXj98HefwuPJd3XbLxWypH6C1gfdfjEeDQ7CBJtwIPA1si4ui7yyPiUOfft4B/YubjUl/yiIgTEXGq8/1uYFDSVdmfoZe5dLmXWR/3erhN5rNYxi6Vi8fv//H4XWAuXS7d+PX6pHw/vpj5JHEA2MD7f4D42KyYjwCjwO2zlq8EVnd9/wSwuY95/ALvX9+/CXiDmf/N531tr3PpxK0FjgEr+7FNlsrYefw8fpfD+PV8B+jXF3A38DIzf1n+q86ybcC2zvcPA8eBZztfT3eWX9/Z2M8Be999bR/zuL/zPs8x8wei28/32n7m0nn8eWDnrNf1dJsslbHz+Hn8ah8/d4qamVViqZxDNzOzebigm5lVwgXdzKwSLuhmZpVwQTczq4QLuplZJVzQzcwq4YJuZlaJ/wPnATbL3VFEdwAAAABJRU5ErkJggg==\n",
            "text/plain": [
              "<Figure size 432x288 with 3 Axes>"
            ]
          },
          "metadata": {
            "needs_background": "light"
          },
          "output_type": "display_data"
        },
        {
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAACPCAYAAADqQSA9AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAQR0lEQVR4nO3dfYxc5XXH8e+Z2XdiG4whTrwGmZeEQAVqGlxKU5U0ojVUyEKqKpKoVVCQZbVUVVtFQYqUP5I/kjR/tFJL6loIRVXbOJXy5qYmJlHVUMUFTKhtsoRFawfweGvj9frdXmZ37+kfOw7DZO05j33HO378+0gr79w5c+fsfWaP79yd5zzm7oiIyKWvstAJiIhIOVTQRUQyoYIuIpIJFXQRkUyooIuIZEIFXUQkEz0L9sQDV3jfoqWh2Fvfe7DD2eRryotw7O6D7w7Hzh6YYLBnUSh2ekl/fL9D4VAA+gfq4djlfUfDsYM2G46d8mo49ngxEI4FOHT6inBs72ELxdWPTzLYd2Uodmp5/JzvPUPx4wtgFv/I9HQRP8aT9fiLaHY6vl+bjh3fM3pOx2Orp2fCscem9k+4+zXzPmf8KcvVt2gptzz4F6HY5z//Dx3OJl+vTp8Mxz74j38Zjh3/6y9y97I/DMUe+P0bwvud/GD8PyCA932gFo799PXfD8fe2hsvTqPTi8OxPzpxSzgW4J9HVodj3/ONvlDcS1v/hrve/0godvSvBsPP/9k7t4RjAQYq0+HYWj128gfwb6/9ajh28kB87Pr3xY7vGVePxF/LS3ZOhGO3vvKl1892ny65iIhkQgVdRCQTKugiIplQQRcRyYQKuohIJlTQRUQyoYIuIpIJFXQRkUyECrqZrTGzUTMbM7PH5rl/iZn9u5ntNLMRM3u4/FRFRORc2s4UNbMq8DhwL1ADtpvZZnd/uSnsT4GX3f0BM7sGGDWzf3H3s87L7pkquGp0KpTkTV9fH4oDeO+vHAjHAnz+pu+GY+8ZjM/8+qdjy8KxRcIbpS/uXBOOBfDd8anj79kVn348PltQHD8Rir32mfiY9J6+NhwLMPG/14Vj/2TVunBsfWl86n/1ZHz8Bg+kvSlefCw+PX7g4KlY4MwslYnYTNjrNr0r/Pwbf/hgOBbg+Mr4VPrZwfhxWDoSj722N57D0hcnw7EAlRPB8QCKQ4eT9n3W5wzErAbG3H1Po0BvAta2xDiwyMwMeBcwCcSrg4iIXLBIQV8B7G26XWtsa/b3wAeAceAl4M/dE7pCiYjIBYs055rvPUnre5rfA3YAvwPcCPzAzP7b3Y+9Y0dm64B1AP39S5KTlYW1//Vn2f/6cwDMeOxymXSHffu3s2//CwBMFwltAOWSEinoNWBl0+1h5s7Emz0MfMndHRgzs58DtwDPNwe5+0ZgI8DiRSviF7qkKyy//i6WX38XALv+48sLnI2kWLH8TlYsvxOA7c//3QJnI50SueSyHbjZzFaZWR/wELC5JeYN4KMAZvZu4P3AnjITFRGRc2t7hu7uM2b2KLAVqAJPuvuIma1v3L8B+ALwNTN7iblLNJ9x93iDXxERuWChBS7cfQuwpWXbhqbvx4HfLTc1ERFJoZmiIiKZUEEXEcmECrqISCYWbJFomynoPRRbwHjZi/FV4/dV4ivXAzw8Hm87s+iq+FTe4+OLwrE9x+Irjw++mbby+NCb8fldg/vjPx/uMBubHm/1+GLAV+46Es8BmLk6vogxNhAO7d8RP87m8U/g9p6IHwuAviNvhWN79h+JBbrjwTEZOBB/TVgxFI4FqEwnvO4Pxo9bdSretqF3MuHnOx0fi2RFOfMwdYYuIpIJFXQRkUyooIuIZEIFXUQkEyroIiKZUEEXEcmECrqISCZU0EVEMqGCLiKSCRV0EZFMqKCLiGRiwXq5UJ/GX6uFQhcviffrOLniiqQ0fH+8T0zlVDz2qpl4f49KPRxKz1Raz4crR0+EY6u1g+FYd6d4K9jb4uCh8H5ZlDZ+1ZPxHh9LRxKOXcICiZXpeO+Qoi/evwSgOnE8HOuTR2JxszMUhw/Hnj9hPHoHesOxAH2T8XVpq0fj66D6YF88iYQ+PEyl9XKZPRhf48eLclbk1Bm6iEgmVNBFRDKhgi4ikolQQTezNWY2amZjZvbYWWLuMbMdZjZiZj8qN00REWmn7R9FzawKPA7cC9SA7Wa22d1fboq5EvgqsMbd3zCzazuUr4iInEXkDH01MObue9y9DmwC1rbEfBz4lru/AeDub5abpoiItBMp6CuAvU23a41tzd4HXGVm/2VmPzGzPy4rQRERiYl8Dn2+xRVbPzTZA/wa8FFgEPgfM3vW3V99x47M1gHrAAYs7fPGsvD2nhph76kRAKbp4PqKUrq99Veo1UcBmHaNXa4iBb0GrGy6PQyMzxMz4e4ngZNm9gxwB/COgu7uG4GNAEuqy8r5JL1cNCuHbmPl0G0A/Pj//nWBs5EUK/tuYWXfLQBsO/btBc5GOiVyyWU7cLOZrTKzPuAhYHNLzHeB3zKzHjMbAn4d+Fm5qYqIyLm0PUN39xkzexTYClSBJ919xMzWN+7f4O4/M7PvA7uAAnjC3X/a9tmrsWnQ1VPxufHXvJjWzaAyG3+jUF8S33dlOmHq/1vxKel9h+PTpQGqh0+GY31mJmnfnWBvtL75axM/OBCOrfYkvDZ64lP0vS8+5b06m9a6oUhom2DRn8/jU839ePz10zsbb4EAhH//5xKJ/z7ZyXibABJe88WRo/H9ApWhoXDs7LFjSfs+m9ArwN23AFtatm1ouf0V4CulZCUiIsk0U1REJBMq6CIimVBBFxHJhAq6iEgmVNBFRDKhgi4ikgkVdBGRTKigi4hkQgVdRCQTafPkF4idjE9377f5mkOeXTEQPwT9CStzV4/F2xWk5FBJaIMAaVO3kwWnY/t051oKWMp084H+eOxb8eNsCbEkvj4p4q0CvB7Lw8yoRNsVvJXQmTGhXQKA1xNem5WEVgxT8XoRbpcAeGrbhuPH48Epr4tz/NrpDF1EJBMq6CIimVBBFxHJhAq6iEgmVNBFRDKhgi4ikgkVdBGRTKigi4hkQgVdRCQToYJuZmvMbNTMxszssXPE3Wlms2b2B+WlKCIiEW3nvZpZFXgcuBeoAdvNbLO7vzxP3JeBrZEn9qKgOHEilGQ1YVqscXU4FqBnImGK/pWL4nkkTEmv1BOmxh+cjMcCfiq+Anp06ngqn06YRl9Ne9OYshJ7JWHl+BTWG5xGDxRH01Z3TxkTD7amcHeK4JR+S5jubqmvn4TjZpbwO5LSDqK/Lx7raVP/rZrQrmCmnPYYkd+e1cCYu+9x9zqwCVg7T9yfAd8E3iwlMxERSRIp6CuAvU23a41tv2BmK4AHgQ3lpSYiIikircbmu97R+t7ub4HPuPusnePyiJmtA9YBDDAUTFG6Ra0Yo1bsBmCahC58suBqxW72ucYud5GCXgNWNt0eBsZbYj4EbGoU82XA/WY24+7faQ5y943ARoDFtrQzFzSlY4YrNzFcuQmAZ6dDfyqRLjFcuZFhbgTgudmnFzgb6ZRIQd8O3Gxmq4B9wEPAx5sD3H3Vme/N7GvA91qLuYiIdFbbgu7uM2b2KHOfXqkCT7r7iJmtb9yv6+YiIl0gtFyHu28BtrRsm7eQu/snLzwtERFJpZmiIiKZUEEXEcmECrqISCbiS153gsX+P/GUqbwHJtJySJj6W5mIT89NWSHcElZLL07HVzQH8OmUnBOOc4KUldWL+nTavnsTVm1POHZJq8GnHOMOtVcAoEgYv2gbhITp7sVU2s+WcjZZJLRt8GBbAwBLGbtga4VfSBmPkugMXUQkEyroIiKZUEEXEcmECrqISCZU0EVEMqGCLiKSCRV0EZFMqKCLiGRCBV1EJBMq6CIimVBBFxHJxML2cglK6c2QyhJ6uRQnT8V3XMR7YDDQH99tYi+XSsK+O9XLpVP7BSBh30VC3w4S1t2s9PWGY1OPhfXFX5/hY5HQksRnUo5ZGk/oz0JKH5VzrGv8SznMJPQOSsl3gegMXUQkEyroIiKZUEEXEclEqKCb2RozGzWzMTN7bJ77P2Fmuxpf28zsjvJTFRGRc2lb0M2sCjwO3AfcCnzMzG5tCfs58NvufjvwBWBj2YmKiMi5Rc7QVwNj7r7H3evAJmBtc4C7b3P3w42bzwLD5aYpIiLtRAr6CmBv0+1aY9vZfAp46kKSEhGRdJHPoc/3oc55P5BpZh9hrqB/+Cz3rwPWAQwwFExRukWt2M0+3w3AdMLntGXh1XwP+3wPoLHLWaSg14CVTbeHgfHWIDO7HXgCuM/dD823I3ffSOP6+mJb2v2f0pd3GK7cyDA3AvDc7NMLnI2kGLYbGLYbAHiu+OECZyOdErnksh242cxWmVkf8BCwuTnAzK4DvgX8kbu/Wn6aIiLSTtszdHefMbNHga1AFXjS3UfMbH3j/g3A54Crga/a3LTbGXf/UNtn99j0eJ9JOJmvVOOxQHHkaDjW+hOm0dfjU4qLUwktBRKnHxenTyfFd0QHp0wnTaVPySNh+ngxldaOIYUntStYYAnHDMDr9Q4lkpJEXhcKQr1c3H0LsKVl24am7x8BHik3NRERSaGZoiIimVBBFxHJhAq6iEgmVNBFRDKhgi4ikgkVdBGRTKigi4hkQgVdRCQTKugiIpmwpJW3y3xis4PA6yXvdhkwUfI+z1fuuXwQeLHkfeZ+zM6Hxi5d7rlc7+7XzHfHghX0TjCzF0I9ZC4C5ZKum/Lslly6JY92uinPyzkXXXIREcmECrqISCZyK+jdtDi1cknXTXl2Sy7dkkc73ZTnZZtLVtfQRUQuZ7mdoYuIXLYumYJuZmvMbNTMxszssXnu/4SZ7Wp8bTOzO5rue83MXjKzHWb2QofzuMfMjjaea4eZfS762A7k8ummPH5qZrNmtrRxX2nHpIQ8L8rYBXPR+KXnqfHrlvFz967/Ym7pu93ADUAfsBO4tSXmbuCqxvf3Ac813fcasOwi5XEP8L3zeWzZubTEPwD8Z9nH5FIZO42fxu9yGL9L5Qx9NTDm7nvcvQ5sAtY2B7j7Nnc/3Lj5LDC8EHl06LFl7O9jwNcv4PnOV7eMXSiXDj22jP1p/DR+bV0qBX0FsLfpdq2x7Ww+BTzVdNuBp83sJ2a27iLk8RtmttPMnjKz2xIfW3YumNkQsAb4ZtPmso5JO90ydim5aPzepvG7sFwu6viFFonuAvMtJz7vx3PM7CPMvag+3LT5N9193MyuBX5gZq+4+zMdyuNF5qbmnjCz+4HvADcHH1t2Lmc8APzY3SebtpV1TNrplrGL5qLxeyeN3/nncsZFG79L5Qy9Bqxsuj0MjLcGmdntwBPAWnc/dGa7u483/n0T+DZzb5c6koe7H3P3E43vtwC9ZrYs+jOUmUuTh2h5u1fiMWmnW8YulIvG75do/M4zlyYXb/zKvijfiS/m3knsAVbx9h8gbmuJuQ4YA+5u2X4FsKjp+23Amg7msZy3P9+/GniDuf/N2z627FwacUuASeCKThyTS2XsNH4av8th/Ep/AXTqC7gfeJW5vyx/trFtPbC+8f0TwGFgR+Prhcb2GxoHeycwcuaxHczj0cbz7GTuD0R3n+uxncylcfuTwKaWx5V6TC6VsdP4afxyHz/NFBURycSlcg1dRETaUEEXEcmECrqISCZU0EVEMqGCLiKSCRV0EZFMqKCLiGRCBV1EJBP/D4VKdr315U9XAAAAAElFTkSuQmCC\n",
            "text/plain": [
              "<Figure size 432x288 with 3 Axes>"
            ]
          },
          "metadata": {
            "needs_background": "light"
          },
          "output_type": "display_data"
        },
        {
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAACPCAYAAADqQSA9AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAP30lEQVR4nO3df2xd9XnH8ffjazuJndiJY35UiWnDRAmBQbV1oTCqBVC6gESjTRvjx5hKQVakMVWqNDXrpP7RTtWmalWnDZRGiLXdGP6jpTRrQwNSW7oRoIEuIQnEqetAbBzyEzs/Hdv3PvvDN+Vy69jPSe7Ntb/5vKQr7j3nud/z+HwvT84993y/x9wdERGZ+epqnYCIiFSGCrqISCJU0EVEEqGCLiKSCBV0EZFEqKCLiCSivlYbbm/L+Uc6Gmq1eTlP23blmNXcFopd2nGwytlIFr19IyxanAvFzrJYnJyf7UPt4diRt9855O6XTLSuZgX9Ix0N/GJTR602L+epZdkwv7vyc6HYF//lm1XORrK44VNv8cyPYgVkScPcKmcjAEuefTgcu/eza98+2zqdchERSYQKuohIIlTQRUQSoYIuIpIIFXQRkUSooIuIJEIFXUQkESroIiKJCBV0M1tlZt1m1mNmaydY32pm/21m28xsp5k9WPlURURkMlOOFDWzHPAosBLoB7aY2QZ3f6Mk7K+BN9z9LjO7BOg2syfdfeRs7e4bm81XD10dSvKL7d2hOLlwcsN5WnYfDcVe+fxnw+32rnziXFO66P3o5OxQXM+By/n0v34+FLv984+dT0pJ+YdDSzPF//tPVoRjP/xcPhy7d5J1kSP05UCPu/cWC3QXsLosxoF5ZmbAXOAIMBbOUEREzlukoC8C+kpe9xeXlfo34BpgANgOfM7dCxXJUEREQiKTc9kEy8rvLP3HwFbgNuB3gOfN7H/c/QPfyc2sE+gEaP3QnMzJSm2t/48hHv/P8S4dGYt9vZfp4cdPHeG5riMA5E/p/71URQp6P1A6LeJixo/ESz0I/KO7O9BjZnuApcAvSoPcfT2wHmDRtfPL/1GQaa7zgVY6H2gFYMGVJ2qcjWSx6t42Vt07Pt3x3TcP1jYZqZrIKZctwFVmtsTMGoF7gA1lMXuB2wHM7DLgaqC3komKiMjkpjxCd/cxM3sE2ATkgCfcfaeZrSmuXwd8BfiWmW1n/BTNF9z9UBXzFhGRMqEbXLj7RmBj2bJ1Jc8HgE9VNjUREclCI0VFRBKhgi4ikggVdBGRRNTsJtGHB+fx7R/cFor94kMa+n8h3N17ezx47GnqDrwXCm3/6fxws1fP/qt4DkD3J7+TKX6mefLYwnDsl797dyiufvjrtO2KDeS+b8+t4e3/15KfhmOni798a0U49rVNyzK13d4bvzK76a2hTG2fjY7QRUQSoYIuIpIIFXQRkUSooIuIJEIFXUQkESroIiKJUEEXEUmECrqISCJU0EVEEqGCLiKSCBV0EZFE1Gwul9wwtL2pu9BV289Oxf/N3vq/H403nC9QOHosFDp7MB9u9tSrc+M5ALe0/mk4dvXibeHYBovnvGzWO+HYdQMrwrEAW7dfGY5d+OtY3P5RZ87+U6HYLS8uDW//1pMt4ViABzteDMd27Vsejt2979Jw7KztTeHYxtgu+42574zEg/vfzdb4WegIXUQkESroIiKJUEEXEUlEqKCb2Soz6zazHjNbe5aYFWa21cx2mtkLlU1TRESmMuWPomaWAx4FVgL9wBYz2+Dub5TEzAceA1a5+14zi/8qISIiFRE5Ql8O9Lh7r7uPAF3A6rKY+4Cn3X0vgLsfqGyaIiIylUhBXwT0lbzuLy4r9VFggZn9zMxeM7Ns9xETEZHzFrkO3SZYVn4BeT3w+8DtwBzgJTN72d13f6Ahs06gE6CxaUH2bKWmjr7yEkdfeRkA8+EaZyNZHNr1Eoe6x/uuMHaixtlItUQKej/QUfJ6MTAwQcwhdz8BnDCznwM3AB8o6O6+HlgP0LywQ6OKZpiWG2+i5cabADj4d1+tcTaSRfvSm2hfOt53v3rqn2ucjVRL5JTLFuAqM1tiZo3APcCGspgfAJ80s3ozawJuBN6sbKoiIjKZKY/Q3X3MzB4BNgE54Al332lma4rr17n7m2b2Y+B1oAA87u47Jms3d7pAS0/sq991L98figPY8Yknw7Ez0brB8p8vJveNHbeFYxe8MXXMGQcBPPYla+7uwXC7hfpsp+Lee+HycOw3F8Yvvsq3xof+z7vkeDh29P+y/X2th+OxLW+fDsXZ6Ci5/oOh2LYd8akY9jZ+KBwL8OWtfx4Prot/oc+dmugs8cSaB+Ltzt03Go4FaDySca6ACgjN5eLuG4GNZcvWlb3+GvC1yqUmIiJZaKSoiEgiVNBFRBKhgi4ikggVdBGRRKigi4gkQgVdRCQRKugiIolQQRcRSYQKuohIIkIjRavB8gVyx2NDlf3V9nC7dy64M1Men74sfif4v5j3q3Ds86fiw6C73s1wR/ND2e4dMueFeeHYlrfjQ5W9UKAwHOu/+iND4Xbn7mkIxwK4xYemnzyaizds8WOdfGN8OP+ck9nmpGvrju1jgIaDJ2OB7vhI7I70rT3BNgHLzwnHAuQb47Gn2+L90bS/EI5t2ROfNTQ3lG2GUTt4JBxbOFWZ2Ut1hC4ikggVdBGRRKigi4gkQgVdRCQRKugiIolQQRcRSYQKuohIIlTQRUQSoYIuIpKIUEE3s1Vm1m1mPWa2dpK4PzCzvJn9WeVSFBGRiCmH/ptZDngUWAn0A1vMbIO7vzFB3D8Bm0JbLjh2KjasuWVPfCjv7teuCMcCfH1xfFqBb7d+Ihx7eKg5nkRvPLa1J94sQP1wfN817DuarfEgPx0bZg5QdyzbEOiG403h2Ob4rsAyxDYeHQvH1p3OxxsGcsMZ2j48GAt0IB/Lo/7gsfD2m5qyzSRSl49Pg5DbORqOHZsXn1Mgy99HhqH8QHgfA3iG2MlEjtCXAz3u3uvuI0AXsHqCuL8BvgccqEhmIiKSSaSgLwL6Sl73F5f9hpktAv4EWFe51EREJIvIdySbYFn5d6VvAF9w97zZROHFhsw6gU6A2fXxWQBleugb3Er/4PjslKPEZwGU2us7sYO+kzsBGPH4rJoys0QKej/QUfJ6MTBQFvNxoKtYzNuBO81szN2fKQ1y9/XAeoDWWZdnm0dUaq5j/sfomP8xADZ3r69tMpJJR/N1dDRfB8Dm/V01zkaqJVLQtwBXmdkS4B3gHuC+0gB3X3LmuZl9C/hheTEXEZHqmrKgu/uYmT3C+NUrOeAJd99pZmuK63XeXERkGghdZ+TuG4GNZcsmLOTu/pnzT0tERLLSSFERkUSooIuIJEIFXUQkEdnG6lZSIY8fjQ27nbc3flf1Qn3GO4+/GR92f2J+/A7zs85+Of5vaX43Ps68eSA+BBqg/mSG+CODmdqmEBuuXDh+ItxkblaGW8EDzVvjbRcujX+OvD5+rGMZhvNnuRM8gF/WFo4N7+dCAR+OjSOoG45P2zB73/FwLAAjGT6b9bl46MH4FBaF/QfDsdaUrbZk+dzjGeaamISO0EVEEqGCLiKSCBV0EZFEqKCLiCRCBV1EJBEq6CIiiVBBFxFJhAq6iEgiVNBFRBKhgi4ikggVdBGRRNRsLhfPFygcjc390DDwXrjdtmPDmfIYbWsKx9YPxe+j6Q0Z5gIZjc/jUHcy2708bSg2Xw6M90m2xmMT1ng+PtdJ4cTJbCnkMuznvfF9Z40N8djgfgDwQrY7L9o7BzIEB/vDnUJ0HpXBofj2D2ebp6auJX5f4fx7g+FYq4+XtSyxhSz7AvCxsXhwXXyumkmbqUgrIiJScyroIiKJUEEXEUlEqKCb2Soz6zazHjNbO8H6+83s9eJjs5ndUPlURURkMlMWdDPLAY8CdwDLgHvNbFlZ2B7gj9z9euArwPpKJyoiIpOLHKEvB3rcvdfdR4AuYHVpgLtvdvczl6K8DCyubJoiIjKVSEFfBPSVvO4vLjubh4BnzycpERHJLnIR5kQXt054Ma2Z3cp4Qb/lLOs7gU6A2cSv/5bpoW/4TfpP7wJglGzXw0tt9Y3son+kG1DfpSxS0PuBjpLXi4GB8iAzux54HLjD3Q9P1JC7r6d4fr2lbmG2ERZScx2zr6Fj9jUAbD7y3RpnI1l0NC6lo3EpAJuHvl/jbKRaIqdctgBXmdkSM2sE7gE2lAaY2RXA08AD7r678mmKiMhUpjxCd/cxM3sE2ATkgCfcfaeZrSmuXwd8CVgIPFYcBj3m7h+vVJKFA4fCsXY826mcxv3x4co2Z3a84dMj4VAfDQ7DBsg4PD9/LMPQ/yxDlbPweM5Zh/5nkWuZG471k6fisRmmNvDRbPvYxzJ8NjzDl95gn2Tpj7oM0yUA5A/Hp/TIso/JMtVExqkYMskwJQSFDH/fJEITGbj7RmBj2bJ1Jc8fBh6uSEYiInJONFJURCQRKugiIolQQRcRSYQKuohIIlTQRUQSoYIuIpIIFXQRkUSooIuIJEIFXUQkEeZZhgtXcsNmB4G3K9xsOxCfJ6C6Us/l94BfVrjN1PfZuVDfZZd6Lh9290smWlGzgl4NZvZqJeeQOR/KJbvplOd0yWW65DGV6ZTnxZyLTrmIiCRCBV1EJBGpFfTpdHNq5ZLddMpzuuQyXfKYynTK86LNJalz6CIiF7PUjtBFRC5aM6agm9kqM+s2sx4zWzvB+vvN7PXiY7OZ3VCy7i0z225mW83s1SrnscLMhorb2mpmX4q+twq5/G1JHjvMLG9mbcV1FdsnFcjzgvRdMBf1X/Y81X/Tpf/cfdo/GL/13a+BK4FGYBuwrCzmZmBB8fkdwCsl694C2i9QHiuAH57LeyudS1n8XcBPKr1PZkrfqf/UfxdD/82UI/TlQI+797r7CNAFrC4NcPfN7n7mJoUvA4trkUeV3luJ9u4FnjqP7Z2r6dJ3oVyq9N5KtKf+U/9NaaYU9EVAX8nr/uKys3kIeLbktQPPmdlrZtZ5AfK4ycy2mdmzZnZtxvdWOhfMrAlYBXyvZHGl9slUpkvfZclF/fc+9d/55XJB+y90k+hpYKLbZ094eY6Z3cr4h+qWksV/6O4DZnYp8LyZ7XL3n1cpj18yPjT3uJndCTwDXBV8b6VzOeMu4EV3P1KyrFL7ZCrTpe+iuaj/Pkj9d+65nHHB+m+mHKH3Ax0lrxcDA+VBZnY98Diw2t0Pn1nu7gPF/x4Avs/416Wq5OHuR939ePH5RqDBzNqjf0MlcylxD2Vf9yq4T6YyXfoulIv677eo/84xlxIXrv8qfVK+Gg/Gv0n0Akt4/weIa8tirgB6gJvLljcD80qebwZWVTGPy3n/+v7lwF7G/zWf8r2VzqUY1wocAZqrsU9mSt+p/9R/F0P/VfwDUK0HcCewm/Fflv++uGwNsKb4/HHgPWBr8fFqcfmVxZ29Ddh55r1VzOOR4na2Mf4D0c2TvbeauRRffwboKntfRffJTOk79Z/6L/X+00hREZFEzJRz6CIiMgUVdBGRRKigi4gkQgVdRCQRKugiIolQQRcRSYQKuohIIlTQRUQS8f9SGRmy28HTSAAAAABJRU5ErkJggg==\n",
            "text/plain": [
              "<Figure size 432x288 with 3 Axes>"
            ]
          },
          "metadata": {
            "needs_background": "light"
          },
          "output_type": "display_data"
        },
        {
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAACPCAYAAADqQSA9AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAORElEQVR4nO3df2xd9X3G8fdjJyGEYCckQCExIzDaDcYP8SMI2gpaRglMKJu6IQpiakVrZRrV/qgqkDZVqipNmti0ShtVlmUMqVSLokHbtAtLEdNgbSAYWAgEMHIDjd1khUCc4ITEsf3ZH76Uy+2N7/c45+Zef/O8JCv3nPO553xyv9ePj4/vOUcRgZmZzXwdrW7AzMzK4UA3M8uEA93MLBMOdDOzTDjQzcwy4UA3M8vErFZtePFpnXFuz+xWbd6O0ZuDR/D4zUweu5nt+W2H90TE6fWWtSzQz+2ZzbObelq1eTtGy28a9PjNUB67ma3zrIFfHG2ZD7mYmWXCgW5mlgkHuplZJhzoZmaZcKCbmWXCgW5mlgkHuplZJhzoZmaZSAp0SSsk9UsakHRfneXdkn4k6UVJ2yV9qfxWzcxsKg3PFJXUCTwA3AgMAX2SNkTEK1Vlfw68EhG3Sjod6Jf0vYgYbUrX1nKvHlzIVS/cllTbd/n6JndjZpC2h74cGIiIHZWAXgesrKkJ4FRJAuYD7wJjpXZqZmZTSgn0JcBg1fRQZV61fwR+F9gFvAT8RURMlNKhmZklSbk4l+rMq72z9E3AVuCzwPnA45L+JyL2f2RFUi/QC3DOkpZdF8ymac1397H24ckhHdt3Uou7sSKqx+7td8Zb3I01S8oe+hBQfWm2pUzuiVf7EvBoTBoA3gB+p3ZFEbEmIq6MiCtPX9Q53Z6tRXrv6ubZTT08u6mHWd3zWt2OFVA9dv7ey1dKoPcBF0haJmkOcDuwoaZmJ3ADgKQzgU8AO8ps1MzMptbwuEdEjEm6B9gEdAIPRsR2Sasqy1cD3wIekvQSk4do7o2IPU3s28zMaiQdyI6IjcDGmnmrqx7vAj5XbmtmZlaEzxQ1M8uEA93MLBMOdDOzTPjD4DY9e2fRsX5RWu3lzW3FbKZb9h9fKVB971GXeA/dzCwTDnQzs0w40M3MMuFANzPLhAPdzCwTDnQzs0w40M3MMuFANzPLhAPdzCwTDnQzs0w40M3MMuFrudi0aAJmH0y7D/jBidHk9c7rmDPdlizRqwcXctULtyXV9l2+vsnd5Ov5w+nv+7MfT78t4M4plnkP3cwsEw50M7NMONDNzDKRFOiSVkjqlzQg6b6j1Fwvaauk7ZKeLLdNMzNrpOEfRSV1Ag8ANwJDQJ+kDRHxSlXNAuA7wIqI2CnpjCb1a2ZmR5Gyh74cGIiIHRExCqwDVtbU3AE8GhE7ASLirXLbNDOzRlICfQkwWDU9VJlX7ePAQkn/Lel5SX9aVoNmZpYm5XPoqjMv6qznCuAG4GTgaUnPRMTrH1mR1Av0ApyzxB+Bn2nWfHcfax/eD8CRwye3uBsronrsxvad1OJurFlSUnUI6KmaXgrsqlOzJyIOAAckPQVcCnwk0CNiDbAG4MpL59b+ULA213tXN713dQPQdeGhFndjRVSP3amfHmlxN9YsKYdc+oALJC2TNAe4HdhQU/ND4NOSZkmaB1wNvFpuq2ZmNpWGe+gRMSbpHmAT0Ak8GBHbJa2qLF8dEa9K+k9gGzABrI2Il6da796JTh4Z6Upq8vPz9yfV2fHTeWCU7r7aX9Tq++y2O5LX+8xl/z7dlk5420bTfmsaPzCLkb7FaSu9/BgaOsH9yY++mlx77vBYKdtMOpAdERuBjTXzVtdM3w/cX0pXZmZWmM8UNTPLhAPdzCwTDnQzs0w40M3MMuFANzPLhAPdzCwTDnQzs0w40M3MMuFANzPLRMsuebj7UBd/3b8iqfbzV/jO48fDN9++ML04nobRI0mlw309jYsqbuu6Ib0HYP15TxSqz9nKn6Sdaj7r0D+w8LWJpNqv7U4/9//vznohuXYm+r1n7ixUv+CV9P3luTvKuYWE99DNzDLhQDczy4QD3cwsEw50M7NMONDNzDLhQDczy4QD3cwsEw50M7NMONDNzDKRFOiSVkjqlzQg6b4p6q6SNC7pj8tr0czMUjQ89V9SJ/AAcCMwBPRJ2hARr9Sp+xtgU8qGY2QWh3+aeOfxK9LK7Ng89Own04vjaWJsPKn0jP9Nv6P5Nj6R3gPwldlpd7oH+OeenxVadzPsHBspVH/Tlj9Lrl2wdXZS3d7D43T9PK2P72+/LHn71546kFwLcNO89NPdd6S/hXjyQPp76O/7fj+5dsGWOelNAGds2Z9cq5GDhdZ9NCl76MuBgYjYERGjwDpgZZ26rwKPAOVclMDMzApJCfQlwGDV9FBl3q9JWgL8EbC6vNbMzKyIlKstqs68qJn+NnBvRIxL9corK5J6gV6AWV0LE1u0djHy5BZGntoCQMfE+y3uxop4Z/vTvLv9aQBirJxf7639pAT6EFB9/dOlwK6amiuBdZUwXwzcImksIn5QXRQRa4A1ACef1VP7Q8Ha3Pzrrmb+dVcDMLzqb1vcjRWx6KJrWHTRNQC88a/3t7gba5aUQO8DLpC0DPglcDtwR3VBRCz74LGkh4Af14a5mZk1V8NAj4gxSfcw+emVTuDBiNguaVVluY+bm5m1gaQ7FkXERmBjzby6QR4RXzz2tszMrCifKWpmlgkHuplZJhzoZmaZSDqG3pQNHwjO7DucVHv3zk8lr/dfzvnpdFvK0sVb7mhcVHHGU+lvh73j40wM70uqnf/6cPJ6x+aellwL8MSCi5Nrl+85O7n24OH007xnd6ZdAgFg3755ybUAXc/NTa49rT/t++nNsXE6304bu8VPLE3e/tcpdgmnb3aln8swum1Bcm3HkaOfC1Orazi5lK43j6QXA51730uuTf1easR76GZmmXCgm5llwoFuZpYJB7qZWSYc6GZmmXCgm5llwoFuZpYJB7qZWSYc6GZmmXCgm5llwoFuZpaJll3LpePQEeb2706q7Vt/SfJ6b/yDYvcq/fb565Nrf3t2+sv15Pvp1+z4q/4/TK59p39Rci1A10D6z+xFz76VvuIIYjztGiYxmDbOAN1zir0lT9l9UnLt3o8vTu9jeCK5tnM0/W6Kp8xOv84IwClDI8m1s4beSSscG2diz7tJpYu2zk/e/pwDXcm1AEfmpY9d93sFxuNweu3ctw8l13bs/FVyLUCMjaXXRjl35PQeuplZJhzoZmaZcKCbmWUiKdAlrZDUL2lA0n11lt8paVvla7OkS8tv1czMptIw0CV1Ag8ANwMXAl+QdGFN2RvAdRFxCfAtYE3ZjZqZ2dRS9tCXAwMRsSMiRoF1wMrqgojYHBF7K5PPAOm3OTEzs1KkBPoSYLBqeqgy72juBh47lqbMzKy4lA/91vvgbN0PTUr6DJOBXvcmoJJ6gV6AuZ2nJrZo7WJweCtDw1sBOELa/SutPQyOvMTgyMsAjEb6Z69tZkkJ9CGgp2p6KbCrtkjSJcBa4OaIqHuGQ0SsoXJ8vXvOmeV8kt6Om54Fl9Gz4DIANr/2T61txgrpmX8xPfMnb6i9+ZcPt7gba5aUQy59wAWSlkmaA9wObKgukHQO8ChwV0S8Xn6bZmbWSMM99IgYk3QPsAnoBB6MiO2SVlWWrwa+ASwCviMJYCwirpxyxePjTAzvS2ry7Kf2J9UB7N811eH933THx76WXNuRfiYv83elnRYPkH4CNJy/+/0C1TD7/4aTayd2Fzu1OVW8n95zx+49hdY951fpp1KcubPAK32owCGluQXWO5b+vgCI995Lrh0fOZC2zghi9EhSbcdQ+nuie0/a9/OvFbiUBh0FTpk5UuCU+wLvzThY7HtvosC6KenU/6RXNCI2Ahtr5q2uevxl4MuldGRmZtPiM0XNzDLhQDczy4QD3cwsEw50M7NMONDNzDLhQDczy4QD3cwsEw50M7NMONDNzDKhsu42XXjD0tvAL0pe7WKg2LnjzZN7L5cDL5S8ztxfs+nw2BWXey+/FRGn11vQskBvBknPNbyGzHHiXoprpz7bpZd26aORdurzRO7Fh1zMzDLhQDczy0Rugd5ON6d2L8W1U5/t0ku79NFIO/V5wvaS1TF0M7MTWW576GZmJ6wZE+iSVkjqlzQg6b46y++UtK3ytVnSpVXL3pT0kqStkp5rch/XS9pX2dZWSd9IfW4Tevl6VR8vSxqXdFplWWmvSQl9HpexS+zF41e8T49fu4xfRLT9F5O3vvs5cB4wB3gRuLCm5lpgYeXxzcCWqmVvAouPUx/XAz+eznPL7qWm/lbgv8p+TWbK2Hn8PH4nwvjNlD305cBAROyIiFFgHbCyuiAiNkfE3srkM8DSVvTRpOeWsb4vAP92DNubrnYZu6RemvTcMtbn8fP4NTRTAn0JMFg1PVSZdzR3A49VTQfwE0nPS+o9Dn1cI+lFSY9Juqjgc8vuBUnzgBXAI1Wzy3pNGmmXsSvSi8fvQx6/Y+vluI5fgdtut5TqzKv78RxJn2HyTfWpqtmfjIhdks4AHpf0WkQ81aQ+XmDy1NwRSbcAPwAuSHxu2b184FbgZxHxbtW8sl6TRtpl7FJ78fh9lMdv+r184LiN30zZQx8CeqqmlwK7aoskXQKsBVZGxDsfzI+IXZV/3wK+z+SvS03pIyL2R8RI5fFGYLakxan/hzJ7qXI7Nb/ulfiaNNIuY5fUi8fvN3j8ptlLleM3fmUflG/GF5O/SewAlvHhHyAuqqk5BxgArq2ZfwpwatXjzcCKJvbxMT78fP9yYCeTP80bPrfsXip13cC7wCnNeE1myth5/Dx+J8L4lf4GaNYXcAvwOpN/Wf7LyrxVwKrK47XAXmBr5eu5yvzzKi/2i8D2D57bxD7uqWznRSb/QHTtVM9tZi+V6S8C62qeV+prMlPGzuPn8ct9/HymqJlZJmbKMXQzM2vAgW5mlgkHuplZJhzoZmaZcKCbmWXCgW5mlgkHuplZJhzoZmaZ+H8ni5KYlhJVLwAAAABJRU5ErkJggg==\n",
            "text/plain": [
              "<Figure size 432x288 with 3 Axes>"
            ]
          },
          "metadata": {
            "needs_background": "light"
          },
          "output_type": "display_data"
        }
      ],
      "source": [
        "name = \"same_but_beta_bigger\"\n",
        "model = \"CIR\"\n",
        "mx = [30, 30]    # get prob dist for this domain ; bigger = more accurate\n",
        "K_avg = 5\n",
        "beta, gamma = 1.5, 1\n",
        "\n",
        "get_BF_plot(name, model, K_avg, beta, gamma, mx)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "FqQ62WoP5Wan"
      },
      "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.8.5"
    },
    "colab": {
      "provenance": [],
      "collapsed_sections": [],
      "include_colab_link": true
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}