2DIsing.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 2D Ising Model\n",
    "The 2D Ising model using the metropolis-based Monte Carlo method, as described in Fig. 8.8 Schroeder (an undergraduate statistical mechanics/thermodynamics textbook). The magnetization M, the average energy U, the heat capacity C, and the magnetic susceptibility χ are all plotted as functions of temperature. The signature discontinuity in C and χ mark the critical temperature and phase transition of the 2D spin lattice."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "\n",
       "    <div class=\"bk-root\">\n",
       "        <a href=\"https://bokeh.pydata.org\" target=\"_blank\" class=\"bk-logo bk-logo-small bk-logo-notebook\"></a>\n",
       "        <span id=\"1001\">Loading BokehJS ...</span>\n",
       "    </div>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/javascript": [
       "\n",
       "(function(root) {\n",
       "  function now() {\n",
       "    return new Date();\n",
       "  }\n",
       "\n",
       "  var force = true;\n",
       "\n",
       "  if (typeof root._bokeh_onload_callbacks === \"undefined\" || force === true) {\n",
       "    root._bokeh_onload_callbacks = [];\n",
       "    root._bokeh_is_loading = undefined;\n",
       "  }\n",
       "\n",
       "  var JS_MIME_TYPE = 'application/javascript';\n",
       "  var HTML_MIME_TYPE = 'text/html';\n",
       "  var EXEC_MIME_TYPE = 'application/vnd.bokehjs_exec.v0+json';\n",
       "  var CLASS_NAME = 'output_bokeh rendered_html';\n",
       "\n",
       "  /**\n",
       "   * Render data to the DOM node\n",
       "   */\n",
       "  function render(props, node) {\n",
       "    var script = document.createElement(\"script\");\n",
       "    node.appendChild(script);\n",
       "  }\n",
       "\n",
       "  /**\n",
       "   * Handle when an output is cleared or removed\n",
       "   */\n",
       "  function handleClearOutput(event, handle) {\n",
       "    var cell = handle.cell;\n",
       "\n",
       "    var id = cell.output_area._bokeh_element_id;\n",
       "    var server_id = cell.output_area._bokeh_server_id;\n",
       "    // Clean up Bokeh references\n",
       "    if (id != null && id in Bokeh.index) {\n",
       "      Bokeh.index[id].model.document.clear();\n",
       "      delete Bokeh.index[id];\n",
       "    }\n",
       "\n",
       "    if (server_id !== undefined) {\n",
       "      // Clean up Bokeh references\n",
       "      var cmd = \"from bokeh.io.state import curstate; print(curstate().uuid_to_server['\" + server_id + \"'].get_sessions()[0].document.roots[0]._id)\";\n",
       "      cell.notebook.kernel.execute(cmd, {\n",
       "        iopub: {\n",
       "          output: function(msg) {\n",
       "            var id = msg.content.text.trim();\n",
       "            if (id in Bokeh.index) {\n",
       "              Bokeh.index[id].model.document.clear();\n",
       "              delete Bokeh.index[id];\n",
       "            }\n",
       "          }\n",
       "        }\n",
       "      });\n",
       "      // Destroy server and session\n",
       "      var cmd = \"import bokeh.io.notebook as ion; ion.destroy_server('\" + server_id + \"')\";\n",
       "      cell.notebook.kernel.execute(cmd);\n",
       "    }\n",
       "  }\n",
       "\n",
       "  /**\n",
       "   * Handle when a new output is added\n",
       "   */\n",
       "  function handleAddOutput(event, handle) {\n",
       "    var output_area = handle.output_area;\n",
       "    var output = handle.output;\n",
       "\n",
       "    // limit handleAddOutput to display_data with EXEC_MIME_TYPE content only\n",
       "    if ((output.output_type != \"display_data\") || (!output.data.hasOwnProperty(EXEC_MIME_TYPE))) {\n",
       "      return\n",
       "    }\n",
       "\n",
       "    var toinsert = output_area.element.find(\".\" + CLASS_NAME.split(' ')[0]);\n",
       "\n",
       "    if (output.metadata[EXEC_MIME_TYPE][\"id\"] !== undefined) {\n",
       "      toinsert[toinsert.length - 1].firstChild.textContent = output.data[JS_MIME_TYPE];\n",
       "      // store reference to embed id on output_area\n",
       "      output_area._bokeh_element_id = output.metadata[EXEC_MIME_TYPE][\"id\"];\n",
       "    }\n",
       "    if (output.metadata[EXEC_MIME_TYPE][\"server_id\"] !== undefined) {\n",
       "      var bk_div = document.createElement(\"div\");\n",
       "      bk_div.innerHTML = output.data[HTML_MIME_TYPE];\n",
       "      var script_attrs = bk_div.children[0].attributes;\n",
       "      for (var i = 0; i < script_attrs.length; i++) {\n",
       "        toinsert[toinsert.length - 1].firstChild.setAttribute(script_attrs[i].name, script_attrs[i].value);\n",
       "      }\n",
       "      // store reference to server id on output_area\n",
       "      output_area._bokeh_server_id = output.metadata[EXEC_MIME_TYPE][\"server_id\"];\n",
       "    }\n",
       "  }\n",
       "\n",
       "  function register_renderer(events, OutputArea) {\n",
       "\n",
       "    function append_mime(data, metadata, element) {\n",
       "      // create a DOM node to render to\n",
       "      var toinsert = this.create_output_subarea(\n",
       "        metadata,\n",
       "        CLASS_NAME,\n",
       "        EXEC_MIME_TYPE\n",
       "      );\n",
       "      this.keyboard_manager.register_events(toinsert);\n",
       "      // Render to node\n",
       "      var props = {data: data, metadata: metadata[EXEC_MIME_TYPE]};\n",
       "      render(props, toinsert[toinsert.length - 1]);\n",
       "      element.append(toinsert);\n",
       "      return toinsert\n",
       "    }\n",
       "\n",
       "    /* Handle when an output is cleared or removed */\n",
       "    events.on('clear_output.CodeCell', handleClearOutput);\n",
       "    events.on('delete.Cell', handleClearOutput);\n",
       "\n",
       "    /* Handle when a new output is added */\n",
       "    events.on('output_added.OutputArea', handleAddOutput);\n",
       "\n",
       "    /**\n",
       "     * Register the mime type and append_mime function with output_area\n",
       "     */\n",
       "    OutputArea.prototype.register_mime_type(EXEC_MIME_TYPE, append_mime, {\n",
       "      /* Is output safe? */\n",
       "      safe: true,\n",
       "      /* Index of renderer in `output_area.display_order` */\n",
       "      index: 0\n",
       "    });\n",
       "  }\n",
       "\n",
       "  // register the mime type if in Jupyter Notebook environment and previously unregistered\n",
       "  if (root.Jupyter !== undefined) {\n",
       "    var events = require('base/js/events');\n",
       "    var OutputArea = require('notebook/js/outputarea').OutputArea;\n",
       "\n",
       "    if (OutputArea.prototype.mime_types().indexOf(EXEC_MIME_TYPE) == -1) {\n",
       "      register_renderer(events, OutputArea);\n",
       "    }\n",
       "  }\n",
       "\n",
       "  \n",
       "  if (typeof (root._bokeh_timeout) === \"undefined\" || force === true) {\n",
       "    root._bokeh_timeout = Date.now() + 5000;\n",
       "    root._bokeh_failed_load = false;\n",
       "  }\n",
       "\n",
       "  var NB_LOAD_WARNING = {'data': {'text/html':\n",
       "     \"<div style='background-color: #fdd'>\\n\"+\n",
       "     \"<p>\\n\"+\n",
       "     \"BokehJS does not appear to have successfully loaded. If loading BokehJS from CDN, this \\n\"+\n",
       "     \"may be due to a slow or bad network connection. Possible fixes:\\n\"+\n",
       "     \"</p>\\n\"+\n",
       "     \"<ul>\\n\"+\n",
       "     \"<li>re-rerun `output_notebook()` to attempt to load from CDN again, or</li>\\n\"+\n",
       "     \"<li>use INLINE resources instead, as so:</li>\\n\"+\n",
       "     \"</ul>\\n\"+\n",
       "     \"<code>\\n\"+\n",
       "     \"from bokeh.resources import INLINE\\n\"+\n",
       "     \"output_notebook(resources=INLINE)\\n\"+\n",
       "     \"</code>\\n\"+\n",
       "     \"</div>\"}};\n",
       "\n",
       "  function display_loaded() {\n",
       "    var el = document.getElementById(\"1001\");\n",
       "    if (el != null) {\n",
       "      el.textContent = \"BokehJS is loading...\";\n",
       "    }\n",
       "    if (root.Bokeh !== undefined) {\n",
       "      if (el != null) {\n",
       "        el.textContent = \"BokehJS \" + root.Bokeh.version + \" successfully loaded.\";\n",
       "      }\n",
       "    } else if (Date.now() < root._bokeh_timeout) {\n",
       "      setTimeout(display_loaded, 100)\n",
       "    }\n",
       "  }\n",
       "\n",
       "\n",
       "  function run_callbacks() {\n",
       "    try {\n",
       "      root._bokeh_onload_callbacks.forEach(function(callback) {\n",
       "        if (callback != null)\n",
       "          callback();\n",
       "      });\n",
       "    } finally {\n",
       "      delete root._bokeh_onload_callbacks\n",
       "    }\n",
       "    console.debug(\"Bokeh: all callbacks have finished\");\n",
       "  }\n",
       "\n",
       "  function load_libs(css_urls, js_urls, callback) {\n",
       "    if (css_urls == null) css_urls = [];\n",
       "    if (js_urls == null) js_urls = [];\n",
       "\n",
       "    root._bokeh_onload_callbacks.push(callback);\n",
       "    if (root._bokeh_is_loading > 0) {\n",
       "      console.debug(\"Bokeh: BokehJS is being loaded, scheduling callback at\", now());\n",
       "      return null;\n",
       "    }\n",
       "    if (js_urls == null || js_urls.length === 0) {\n",
       "      run_callbacks();\n",
       "      return null;\n",
       "    }\n",
       "    console.debug(\"Bokeh: BokehJS not loaded, scheduling load and callback at\", now());\n",
       "    root._bokeh_is_loading = css_urls.length + js_urls.length;\n",
       "\n",
       "    function on_load() {\n",
       "      root._bokeh_is_loading--;\n",
       "      if (root._bokeh_is_loading === 0) {\n",
       "        console.debug(\"Bokeh: all BokehJS libraries/stylesheets loaded\");\n",
       "        run_callbacks()\n",
       "      }\n",
       "    }\n",
       "\n",
       "    function on_error() {\n",
       "      console.error(\"failed to load \" + url);\n",
       "    }\n",
       "\n",
       "    for (var i = 0; i < css_urls.length; i++) {\n",
       "      var url = css_urls[i];\n",
       "      const element = document.createElement(\"link\");\n",
       "      element.onload = on_load;\n",
       "      element.onerror = on_error;\n",
       "      element.rel = \"stylesheet\";\n",
       "      element.type = \"text/css\";\n",
       "      element.href = url;\n",
       "      console.debug(\"Bokeh: injecting link tag for BokehJS stylesheet: \", url);\n",
       "      document.body.appendChild(element);\n",
       "    }\n",
       "\n",
       "    for (var i = 0; i < js_urls.length; i++) {\n",
       "      var url = js_urls[i];\n",
       "      var element = document.createElement('script');\n",
       "      element.onload = on_load;\n",
       "      element.onerror = on_error;\n",
       "      element.async = false;\n",
       "      element.src = url;\n",
       "      console.debug(\"Bokeh: injecting script tag for BokehJS library: \", url);\n",
       "      document.head.appendChild(element);\n",
       "    }\n",
       "  };var element = document.getElementById(\"1001\");\n",
       "  if (element == null) {\n",
       "    console.error(\"Bokeh: ERROR: autoload.js configured with elementid '1001' but no matching script tag was found. \")\n",
       "    return false;\n",
       "  }\n",
       "\n",
       "  function inject_raw_css(css) {\n",
       "    const element = document.createElement(\"style\");\n",
       "    element.appendChild(document.createTextNode(css));\n",
       "    document.body.appendChild(element);\n",
       "  }\n",
       "\n",
       "  var js_urls = [\"https://cdn.pydata.org/bokeh/release/bokeh-1.3.1.min.js\", \"https://cdn.pydata.org/bokeh/release/bokeh-widgets-1.3.1.min.js\", \"https://cdn.pydata.org/bokeh/release/bokeh-tables-1.3.1.min.js\", \"https://cdn.pydata.org/bokeh/release/bokeh-gl-1.3.1.min.js\"];\n",
       "  var css_urls = [];\n",
       "\n",
       "  var inline_js = [\n",
       "    function(Bokeh) {\n",
       "      Bokeh.set_log_level(\"info\");\n",
       "    },\n",
       "    \n",
       "    function(Bokeh) {\n",
       "      \n",
       "    },\n",
       "    function(Bokeh) {} // ensure no trailing comma for IE\n",
       "  ];\n",
       "\n",
       "  function run_inline_js() {\n",
       "    \n",
       "    if ((root.Bokeh !== undefined) || (force === true)) {\n",
       "      for (var i = 0; i < inline_js.length; i++) {\n",
       "        inline_js[i].call(root, root.Bokeh);\n",
       "      }if (force === true) {\n",
       "        display_loaded();\n",
       "      }} else if (Date.now() < root._bokeh_timeout) {\n",
       "      setTimeout(run_inline_js, 100);\n",
       "    } else if (!root._bokeh_failed_load) {\n",
       "      console.log(\"Bokeh: BokehJS failed to load within specified timeout.\");\n",
       "      root._bokeh_failed_load = true;\n",
       "    } else if (force !== true) {\n",
       "      var cell = $(document.getElementById(\"1001\")).parents('.cell').data().cell;\n",
       "      cell.output_area.append_execute_result(NB_LOAD_WARNING)\n",
       "    }\n",
       "\n",
       "  }\n",
       "\n",
       "  if (root._bokeh_is_loading === 0) {\n",
       "    console.debug(\"Bokeh: BokehJS loaded, going straight to plotting\");\n",
       "    run_inline_js();\n",
       "  } else {\n",
       "    load_libs(css_urls, js_urls, function() {\n",
       "      console.debug(\"Bokeh: BokehJS plotting callback run at\", now());\n",
       "      run_inline_js();\n",
       "    });\n",
       "  }\n",
       "}(window));"
      ],
      "application/vnd.bokehjs_load.v0+json": "\n(function(root) {\n  function now() {\n    return new Date();\n  }\n\n  var force = true;\n\n  if (typeof root._bokeh_onload_callbacks === \"undefined\" || force === true) {\n    root._bokeh_onload_callbacks = [];\n    root._bokeh_is_loading = undefined;\n  }\n\n  \n\n  \n  if (typeof (root._bokeh_timeout) === \"undefined\" || force === true) {\n    root._bokeh_timeout = Date.now() + 5000;\n    root._bokeh_failed_load = false;\n  }\n\n  var NB_LOAD_WARNING = {'data': {'text/html':\n     \"<div style='background-color: #fdd'>\\n\"+\n     \"<p>\\n\"+\n     \"BokehJS does not appear to have successfully loaded. If loading BokehJS from CDN, this \\n\"+\n     \"may be due to a slow or bad network connection. Possible fixes:\\n\"+\n     \"</p>\\n\"+\n     \"<ul>\\n\"+\n     \"<li>re-rerun `output_notebook()` to attempt to load from CDN again, or</li>\\n\"+\n     \"<li>use INLINE resources instead, as so:</li>\\n\"+\n     \"</ul>\\n\"+\n     \"<code>\\n\"+\n     \"from bokeh.resources import INLINE\\n\"+\n     \"output_notebook(resources=INLINE)\\n\"+\n     \"</code>\\n\"+\n     \"</div>\"}};\n\n  function display_loaded() {\n    var el = document.getElementById(\"1001\");\n    if (el != null) {\n      el.textContent = \"BokehJS is loading...\";\n    }\n    if (root.Bokeh !== undefined) {\n      if (el != null) {\n        el.textContent = \"BokehJS \" + root.Bokeh.version + \" successfully loaded.\";\n      }\n    } else if (Date.now() < root._bokeh_timeout) {\n      setTimeout(display_loaded, 100)\n    }\n  }\n\n\n  function run_callbacks() {\n    try {\n      root._bokeh_onload_callbacks.forEach(function(callback) {\n        if (callback != null)\n          callback();\n      });\n    } finally {\n      delete root._bokeh_onload_callbacks\n    }\n    console.debug(\"Bokeh: all callbacks have finished\");\n  }\n\n  function load_libs(css_urls, js_urls, callback) {\n    if (css_urls == null) css_urls = [];\n    if (js_urls == null) js_urls = [];\n\n    root._bokeh_onload_callbacks.push(callback);\n    if (root._bokeh_is_loading > 0) {\n      console.debug(\"Bokeh: BokehJS is being loaded, scheduling callback at\", now());\n      return null;\n    }\n    if (js_urls == null || js_urls.length === 0) {\n      run_callbacks();\n      return null;\n    }\n    console.debug(\"Bokeh: BokehJS not loaded, scheduling load and callback at\", now());\n    root._bokeh_is_loading = css_urls.length + js_urls.length;\n\n    function on_load() {\n      root._bokeh_is_loading--;\n      if (root._bokeh_is_loading === 0) {\n        console.debug(\"Bokeh: all BokehJS libraries/stylesheets loaded\");\n        run_callbacks()\n      }\n    }\n\n    function on_error() {\n      console.error(\"failed to load \" + url);\n    }\n\n    for (var i = 0; i < css_urls.length; i++) {\n      var url = css_urls[i];\n      const element = document.createElement(\"link\");\n      element.onload = on_load;\n      element.onerror = on_error;\n      element.rel = \"stylesheet\";\n      element.type = \"text/css\";\n      element.href = url;\n      console.debug(\"Bokeh: injecting link tag for BokehJS stylesheet: \", url);\n      document.body.appendChild(element);\n    }\n\n    for (var i = 0; i < js_urls.length; i++) {\n      var url = js_urls[i];\n      var element = document.createElement('script');\n      element.onload = on_load;\n      element.onerror = on_error;\n      element.async = false;\n      element.src = url;\n      console.debug(\"Bokeh: injecting script tag for BokehJS library: \", url);\n      document.head.appendChild(element);\n    }\n  };var element = document.getElementById(\"1001\");\n  if (element == null) {\n    console.error(\"Bokeh: ERROR: autoload.js configured with elementid '1001' but no matching script tag was found. \")\n    return false;\n  }\n\n  function inject_raw_css(css) {\n    const element = document.createElement(\"style\");\n    element.appendChild(document.createTextNode(css));\n    document.body.appendChild(element);\n  }\n\n  var js_urls = [\"https://cdn.pydata.org/bokeh/release/bokeh-1.3.1.min.js\", \"https://cdn.pydata.org/bokeh/release/bokeh-widgets-1.3.1.min.js\", \"https://cdn.pydata.org/bokeh/release/bokeh-tables-1.3.1.min.js\", \"https://cdn.pydata.org/bokeh/release/bokeh-gl-1.3.1.min.js\"];\n  var css_urls = [];\n\n  var inline_js = [\n    function(Bokeh) {\n      Bokeh.set_log_level(\"info\");\n    },\n    \n    function(Bokeh) {\n      \n    },\n    function(Bokeh) {} // ensure no trailing comma for IE\n  ];\n\n  function run_inline_js() {\n    \n    if ((root.Bokeh !== undefined) || (force === true)) {\n      for (var i = 0; i < inline_js.length; i++) {\n        inline_js[i].call(root, root.Bokeh);\n      }if (force === true) {\n        display_loaded();\n      }} else if (Date.now() < root._bokeh_timeout) {\n      setTimeout(run_inline_js, 100);\n    } else if (!root._bokeh_failed_load) {\n      console.log(\"Bokeh: BokehJS failed to load within specified timeout.\");\n      root._bokeh_failed_load = true;\n    } else if (force !== true) {\n      var cell = $(document.getElementById(\"1001\")).parents('.cell').data().cell;\n      cell.output_area.append_execute_result(NB_LOAD_WARNING)\n    }\n\n  }\n\n  if (root._bokeh_is_loading === 0) {\n    console.debug(\"Bokeh: BokehJS loaded, going straight to plotting\");\n    run_inline_js();\n  } else {\n    load_libs(css_urls, js_urls, function() {\n      console.debug(\"Bokeh: BokehJS plotting callback run at\", now());\n      run_inline_js();\n    });\n  }\n}(window));"
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "from bokeh.io import output_notebook, show\n",
    "from bokeh.plotting import figure\n",
    "output_notebook()\n",
    "\n",
    "from numba import jit, njit\n",
    "import numpy as np\n",
    "\n",
    "import random\n",
    "import time"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Helper functions\n",
    "\n",
    "# Calculate the energy of the lattice.\n",
    "# The energy is the sum of interactions between spins divided by the total number of spins\n",
    "@njit\n",
    "def calcE(s):\n",
    "    E = 0\n",
    "    for i in range(width):\n",
    "        for j in range(width):\n",
    "            E += -deltaU(s,i,j)/2\n",
    "    return E/width**2\n",
    "\n",
    "# Calculate magnetization. \n",
    "# Magnetization is the sum of all spins divided by the total number of spins\n",
    "@njit\n",
    "def calcM(s):\n",
    "    m = np.abs(s.sum())\n",
    "    return m/width**2\n",
    "\n",
    "# Calculate interaction energy between spins. Assume periodic boundaries\n",
    "# Interaction energy will be the difference in energy due to flipping spin i,j \n",
    "# (i.e. 2*spin_value*neighboring_spins)\n",
    "@njit\n",
    "def deltaU(s,i,j):\n",
    "    #top\n",
    "    if i == 0:\n",
    "        t = s[width-1,j]  # periodic boundary\n",
    "    else:\n",
    "        t = s[i-1,j]\n",
    "    #bottom\n",
    "    if i == width-1:\n",
    "        b = s[0,j]  # periodic boundary\n",
    "    else:\n",
    "        b = s[i+1,j]\n",
    "    #left\n",
    "    if j == 0:\n",
    "        l = s[i,width-1]  # periodic boundary\n",
    "    else:\n",
    "        l = s[i,j-1]\n",
    "    #right\n",
    "    if j == width-1:\n",
    "        r = s[i,0]  # periodic boundary\n",
    "    else:\n",
    "        r = s[i,j+1]\n",
    "    return 2*s[i,j]*(t+b+r+l)  # difference in energy if i,j is flipped\n",
    "\n",
    "# Monte carlo cycle\n",
    "@njit\n",
    "def montecc(s,Temp,cycles):\n",
    "    for m in range(cycles):\n",
    "        i = random.randrange(width)  # choose random row\n",
    "        j = random.randrange(width)  # choose random column\n",
    "        ediff = deltaU(s,i,j)\n",
    "        if ediff <= 0:\n",
    "            s[i,j] = -s[i,j]  # flip spin\n",
    "        elif random.random() < np.exp(-ediff/Temp):\n",
    "            s[i,j] = -s[i,j]\n",
    "    return s\n",
    "\n",
    "# Compute physical quantities\n",
    "@njit\n",
    "def computePQ(s,T,cycles):\n",
    "    En = 0\n",
    "    En_sq = 0\n",
    "    Mg = 0\n",
    "    Mg_sq = 0\n",
    "    for p in range(cycles):\n",
    "        s = montecc(s,T,1)\n",
    "        E = calcE(s)\n",
    "        M = calcM(s)\n",
    "        En += E\n",
    "        Mg += M\n",
    "        En_sq += E*E\n",
    "        Mg_sq += M*M\n",
    "    aveEn = En/cycles\n",
    "    mag = Mg/cycles\n",
    "    sus = (Mg_sq/cycles-(Mg/cycles)**2)/T\n",
    "    heatCap = (En_sq/cycles-(En/cycles)**2)/(T**2)\n",
    "    return aveEn, mag, sus, heatCap"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Elapsed = 153.1240358352661\n"
     ]
    }
   ],
   "source": [
    "# Initial parameters\n",
    "\n",
    "width = 32  # width of sqaure spin lattice; 64 points may take hours\n",
    "s = np.random.choice([1,-1],size=(width,width))  # initial spin sites randomly (+1 or -1)\n",
    "cycles = 1000 * width**2\n",
    "\n",
    "# Initlaize temperature range (crosses critical temperature)\n",
    "temps = np.arange(1.60,3.15,0.02)\n",
    "\n",
    "# Inititalize magnetization, average energy, heat capacity, and susceptibility\n",
    "mag = np.zeros(len(temps))\n",
    "aveEn = np.zeros(len(temps))\n",
    "heatCap = np.zeros(len(temps))\n",
    "sus = np.zeros(len(temps))\n",
    "\n",
    "start = time.time()\n",
    "# Compute physical quanities (magnetization, etc.)\n",
    "for ind, T in enumerate(temps):\n",
    "    # Cycle spins\n",
    "    s = montecc(s,T,cycles)\n",
    "    # Compute physical quanitites with 1000 cycles per spin\n",
    "    aveEn[ind], mag[ind], sus[ind], heatCap[ind] = computePQ(s,T,cycles)\n",
    "end = time.time()\n",
    "print(\"Elapsed = %s\" % (end - start))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY4AAAEGCAYAAABy53LJAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAZgklEQVR4nO3de7RkZXnn8e9P8ILICNp4AxscBS8LEeWIjjiMSAeRiRAVL9EI3oZx1NjJrOUQxfGeLCe6MpJlVBBM1IUjxoiQSIBuoxJHAbsZBBQCBG3ptC5A8a5gwzN/1D56ONQ5Z1efU1W7Tn0/a53Ve1fts+ux7c1T7/u8+9mpKiRJause4w5AkjRZTBySpIGYOCRJAzFxSJIGYuKQJA1k53EHMAxr1qypfffdd9xhSNLE2Lx58y1VtWebY1dl4th3333ZtGnTuMOQpImRZEvbY52qkiQNxMQhSRrIWBJHkhck+WaSO5PMLHLcd5JcmeTyJM49SVIHjKvGcRXwPODUFsceXlW3DDkeSVJLY0kcVXU1QJJxfLwkaRm6XuMo4MIkm5OcuNiBSU5MsinJpptvvnlE4UnS9BnaiCPJRuAhfd46uarOaXmaQ6tqW5IHARuSXFNVF/U7sKpOA04DmJmZseWvJA3J0BJHVa1bgXNsa/68KcnZwCFA38QhSRqNzk5VJdk1yW6z28CR9IrqkjR1Nm+5lePPuITNW25ttT9M41qO+9wkW4H/AHw+yQXN6w9Lcl5z2IOBryT5BnAp8PmqOn8c8UrSuJ2y8Vouuu4WTtl4bav9YRrXqqqzgbP7vL4NOLrZvgF4wohDk6RO2LzlVk7ZeC3r1+3Pwfvswfp1+wO0/nOYshofHTszM1P2qpI0SeYniuPPuISLrruFw/Zbw8df9ZShf36SzVW14A3Zc63KJoeSNGlmp5oAPv6qp4x0BDEoE4ckjcFSU1EH77PHSEYaO8LEIUljMH+E0eVEMV9nl+NK0iRbarns+nX7c9h+azo5FbUURxySNATzRxSTPMKYz8QhSUMwzuWyw+ZyXEnaAfOL2/P3J80gy3GtcUjSDhjnndvj5lSVJLUw6J3cq5lTVZLUwqjv5B41p6okaZlW0/LZlWbikDSVlrrPYn7NYnb57CQWvleaNQ5JU2mp+yymqWYxKBOHpKm0VHF7km/QGzaL45KmxqTfazFMFsclqY9putdimEwcklYtV0YNh4lD0qrhyqjRMHFImlhLJQpHGMPhqipJE2upJbSujBoOE4ekiWWiGA+nqiRNjPlTU9YsxsPEIWliuJy2G0wckjrL5bTdZOKQ1Bkup50MFscldYaNBieDiUNSZ7hKajKYOCR1holiMljjkDQySz08SZPBxCFpaJYqdru8djI5VSVpaJYqdlv8nkw+yEnS0PjgpMnR+Qc5JXlvkmuSXJHk7CS7L3DcUUn+Jcn1Sf5k1HFKGowtQabDuGocG4ADqupA4FrgTfMPSLIT8FfAs4HHAb+f5HEjjVLSQKxZTIexJI6qurCqtje7FwN79znsEOD6qrqhqm4HPgUcO6oYJS3NliDTqQvF8VcCZ/V5fS/gxjn7WwEXeEsdMr/47X0Y02FoiSPJRuAhfd46uarOaY45GdgOnNnvFH1eW7CSn+RE4ESAtWvXDhyvpMG5Kmo6DS1xVNW6xd5PcgLwu8AR1X9p11bg4XP29wa2LfJ5pwGnQW9V1cABSxqYI4zpNK5VVUcBJwHHVNUvFjjs68B+SR6R5F7Ai4FzRxWjpLvzTm/B+FZVfQDYDdiQ5PIkHwZI8rAk5wE0xfPXAxcAVwOfrqpvjileSbhqSj1jKY5X1aMWeH0bcPSc/fOA80YVl6TFWdMQ2KtK0hxLNSH0hj6BiUPSHDYhVBtduI9DUkfYhFBt2ORQmmI2IdSszjc5lNQNTkVpR5g4pClibymtBGsc0hSxt5RWgolDmiIWu7USTBzSFHGEoZVgjUNa5ewvpZVm4pBWOVdOaaU5VSWtctY1tNJMHNIqZ11DK82pKmmVsaahYTNxSBNmqQ621jQ0bE5VSRNm/k188/etaWjYlkwcSV4PnFlVjnulDliqg601DQ3bkt1xk7yb3vO+LwM+ClxQHW+pa3dcSRrMinbHraq3APsBZwAvB65L8mdJHrmsKCVJE6lVcbwZYXy/+dkO7AF8JsmfDzE2SVIHtalxvAE4AbgFOB14Y1X9Osk9gOuA/zHcECVJXdJmVdUa4HlVtWXui1V1Z5LfHU5YkqSuapM43g+Q5AFzXvtpVf26qq4eTliSpK5qU+O4DLgZuJbe1NTNwLeTXJbk4GEGJ0nqnjaJ43zg6KpaU1UPBJ4NfBp4LfDBYQYnyRYi6p42iWOmqi6Y3amqC4HDqupi4N5Di0yaUrYQUde1qXH8MMlJwKea/RcBtybZCbhzaJFJU8oWIuq6NonjJcDbgM81+19pXtsJeOGQ4pKmli1E1HWLJo5mVHFSVf3hAodcv/IhSdPNRKGuW7TGUVV3AK6ckiT9Rpvi+P9Lcm6SlyV53uzP0COTpoSrpjRp2iSOBwA/AJ4JPKf58Y5xaQe5akqTbsnieFW9YhSBSKvV5i23csrGa1m/bn8O3mcPV01p4rVpcrg/8CHgwVV1QJIDgWOq6t1Dj05aBZZKFBbDNWnaLMf9CPBG4FSAqroiySeBHU4cSd5Lb8rrduBfgVdU1Y/6HPcd4KfAHcD2tg8ZkbrERKHVpk2N475Vdem817Yv83M3AAdU1YH0emC9aZFjD6+qg0wamlSzieLgffYYdyjSimiTOG5pnvZXAEmOA763nA+tqgurajb5XAzsvZzzSZJGp81U1euA04DHJPk34NvAH6xgDK8EzlrgvQIuTFLAqVV12kInSXIicCLA2rVrVzA8SdJcbVZV3QCsS7IrcI+q+mmbEyfZCDykz1snV9U5zTEn05v2OnOB0xxaVduSPAjYkOSaqrpogThPo5fgmJmZqTYxSpIG12ZV1b2B5wP7AjsnAaCq3rnY71XVuiXOewK9+0GOaJ5p3u8c25o/b0pyNnAI0DdxSJJGo02N4xzgWHojg5/P+dlhSY4CTqK3rPcXCxyza5LdZreBI4GrlvO5kqTla1Pj2Luqjlrhz/0AvWd5bGhGMBdX1WuSPAw4vaqOBh4MnN28vzPwyao6f4XjkCQNqE3i+GqSx1fVlSv1oVX1qAVe3wYc3WzfADxhpT5TkrQy2iSOpwMvT/Jt4DYgQDX3YEiSpkybxPHsoUchSZoYCxbHkzwToKq20FuGu2X2B5/RIS3INula7RZbVfW+Odt/N++9twwhFmki2SZd02axqaossN1vX5patknXtFkscdQC2/32pall91tNm8USx79Pci690cXsNs3+I4YemTQhTBSaNosljmPnbL9v3nvz9yVJU2LBxFFVXx5lIJKkydCmV5UkSb9h4pAG5H0amnatE0fToVaaet6noWm3ZOJI8rQk3wKubvafkOSDQ49M6qj16/bnsP3WeJ+GplabXlX/G3gWcC5AVX0jyWFDjUrqMJffatq1mqqqqhvnvXTHEGKRJE2ANiOOG5M8Dagk9wLeQDNtJUmaPm1GHK8BXgfsBWwFDmr2JUlTaMkRR1XdArx0BLFIkibAkokjyV/2efnHwKaqOmflQ5IkdVmbqar70Jueuq75ORB4APCqJO8fYmySpA5qUxx/FPDMqtoOkORDwIXA7wBXDjE2SVIHtRlx7AXMvWt8V+BhVXUHcNtQopIkdVabEcefA5cn+RK9Z3EcBvxZ04Jk4xBjkyR1UJtVVWckOQ84hF7ieHNVbWvefuMwg5MkdU/bJoe/Ar4H/BB4lC1HNE3shivdVZvluK8G1gN7A5cDTwW+BjxzuKFJ3TDbDRewR5VEuxHHeuDJwJaqOhx4InDzUKOSOsRuuNJdtSmO/6qqfpWEJPeuqmuSPHrokUkdYTdc6a7ajDi2Jtkd+BywIck5wLYlfkeaGNYwpMG0WVX13Gbz7Um+CNwfOH+oUUkjZA1DGsyiiSPJPYArquoAgKr68kiikkZotnZhDUNqZ9HEUVV3JvlGkrVV9d1RBSWNkjUMaTBtiuMPBb6Z5FLg57MvVtUxQ4tKktRZbRLHO4YehSRpYiy5qqqpa3wHuGez/XXgsuV+cJJ3JbkiyeVJLkzysAWOOyHJdc3PCcv9XEnS8iyZOJL8F+AzwKnNS3vRW5q7XO+tqgOr6iDgH4C39vnsBwBvA55Cr1fW25LssQKfLS3I5bnS4trcx/E64FDgJwBVdR3woOV+cFX9ZM7urkD1OexZwIaq+mFV3QpsAI5a7mdLi5ldnnvKxmvHHYrUSW1qHLdV1e1JAEiyM/3/Iz+wJH8KHE/vUbSH9zlkL+DGOftbm9f6netE4ESAtWvXrkR4mlIuz5UW12bE8eUkbwZ2SfI7wN8Cf9/m5Ek2Jrmqz8+xAFV1clU9HDgTeH2/U/R5rW/SqqrTqmqmqmb23HPPNuFJfc0uzz14H2dFpX7ajDj+BHgVvcfE/lfgPOD0NievqnUt4/gk8Hl69Yy5tgLPmLO/N/CllueUJA1Bm8RxLPDxqvrISn5wkv2aegnAMcA1fQ67gN7TBme/+h0JvGkl45AkDabNVNUxwLVJPpHkPzc1jpXwnmba6gp6CWE9QJKZJKcDVNUPgXfRWwL8deCdzWuSpDFJ1dJ17iT3BJ4NvAh4Or2VTq8ecmw7bGZmpjZt2jTuMCRpYiTZXFUzbY5tNXqoql8n+Ud6held6E1fdTZxSJKGp80NgEcl+RvgeuA4eoXxhw45LklSR7Wpcbyc3p3i+1fVCVV1XlVtH25Y0vB4Z7i0PG16Vb24qj5XVbcBJDk0yV8NPzRpZcxPFN4ZLi1PqxpHkoOAlwAvBL4NfHaYQUkraf4T/rwzXFqeBRNHkv2BFwO/D/wAOIveKqx+rUGkzpqfKHxwk7Q8i404rgH+GXhOVV0PkOSPRxKVtIJMFNLKWqzG8Xzg+8AXk3wkyRH07x0ldY4FcGl4FkwcVXV2Vb0IeAy9/lB/DDw4yYeSHDmi+KQdYgFcGp4li+NV9XN63WvPbB6s9AJ6jQ8vHHJs0g6zAC4NT6uWI5PGliOSNJhBWo60uQFQ6jxrGtLomDg0kbypTxqflWqRLo2UN/VJ42Pi0ETYvOVWTtl4LevX7c/B++zhTX3SGJk4NBHmjzBMFNL4mDg0EZyKkrrDxKGJ4AhD6g5XVUmSBmLikCQNxMQhSRqIiUOSNBAThzrJFiJSd5k41Em2EJG6y8ShTpg/wli/bn8O22+N921IHWTi0Fgs1aRw9r6Ng/fZY5xhSurDGwA1FjYplCaXiUNjYZNCaXKZODQWJgppclnjkCQNxMShkfC+DGn1MHFoKHy0q7R6WePQULhqSlq9xpI4krwLOBa4E7gJeHlVbetz3B3Alc3ud6vqmNFFqeVw1ZS0eqWqRv+hyb+rqp80228AHldVr+lz3M+q6n6Dnn9mZqY2bdq0ApFK0nRIsrmqZtocO5Yax2zSaOwKjD57aUVZ/Jamx9iK40n+NMmNwEuBty5w2H2SbEpycZLfW+J8JzbHbrr55ptXPF4tzuK3ND2GljiSbExyVZ+fYwGq6uSqejhwJvD6BU6zthk6vQR4f5JHLvR5VXVaVc1U1cyee+654v97dFc2JZSm19CK41W1ruWhnwQ+D7ytzzm2NX/ekORLwBOBf12pGLXj5q+asvgtTY+xTFUl2W/O7jHANX2O2SPJvZvtNcChwLdGE6GW4ghDml7jqnG8p5m2ugI4ElgPkGQmyenNMY8FNiX5BvBF4D1VZeIYk/lTU7Y9l6bXWO7jqKrnL/D6JuDVzfZXgcePMi791uYtt3LKxmtZv25/Dt5nj7tNTUmaXt45rr6881vSQkwcAu4+wvDOb0kLMXEIcJWUpPZMHALu3ltKkhZiW/UpMX9VlKukJO0oE8eUmN8SxBYhknaUU1VTYv5UlFNTknbUWNqqD5tt1e++SkqSFtP5tuoaPqeiJA2LU1WrlFNRkobFEccEW2yllKukJA2LiWOCuVJK0jg4VTVBlmoL4vSUpFFwVdUEOf6MS7jouls4bL81tgORtKJcVbVK+HhWSV1k4uiQ+Ylifs3CgrekLjBxjNFSicIRhqQusjg+Qks9Vc9nYEiaBCaOITJRSFqNTBwryEQhaRqYOFaQiULSNLA4PoClHoY0v5jtKihJq5EjjgHMH1H4nG5J08jEsQhbfEjS3dlyZI75icIWH5KmxSAtRxxxzLFUcVuSZOK4C1dBSdLSTBxzmCgkaWkux5UkDcTEIUkaiIlDkjQQE4ckaSAmDknSQEwckqSBmDgkSQNZlS1HktwMbNnBX18D3LKC4ay0LsdnbDuuy/F1OTbodnxdjg3uGt8+VbVnm19alYljOZJsatuvZRy6HJ+x7bgux9fl2KDb8XU5Ntjx+JyqkiQNxMQhSRqIiePuTht3AEvocnzGtuO6HF+XY4Nux9fl2GAH47PGIUkaiCMOSdJATBySpIFMZeJI8tEkNyW5apFjnpHk8iTfTPLlLsWX5P5J/j7JN5r4XjHC2B6e5ItJrm4+e32fY5LkL5Ncn+SKJE/qUGwvbWK6IslXkzxhFLG1jW/OsU9OckeS47oU27iui5b/347lukhynySXzvncd/Q55t5JzmquiUuS7Nuh2P57km8118QXkuyz5Imraup+gMOAJwFXLfD+7sC3gLXN/oM6Ft+bgf/VbO8J/BC414hieyjwpGZ7N+Ba4HHzjjka+EcgwFOBSzoU29OAPZrtZ48qtrbxNe/tBPwTcB5wXFdiG+d10TK+sVwXzb/z+zXb9wQuAZ4675jXAh9utl8MnDWiv7c2sR0O3LfZ/m9tYpvKEUdVXUTvH9VCXgJ8tqq+2xx/00gCa7SIr4DdkgS4X3Ps9hHF9r2quqzZ/ilwNbDXvMOOBT5ePRcDuyd5aBdiq6qvVtWtze7FwN7DjmuQ+Bp/CPwdMLJ/dy1jG9t10TK+sVwXzb/znzW792x+5q86Ohb4WLP9GeCIJs6xx1ZVX6yqXzS7ra6JqUwcLewP7JHkS0k2Jzl+3AHN8wHgscA24EpgfVXdOeogmuH2E+l9i5lrL+DGOftb6f8fyKFZJLa5XkVvZDRyC8WXZC/gucCHRx/Vb2LYl/5/d524LhaJb2zXRZKdklxOL9lvqKoFr4mq2g78GHhgR2Kbq9U14TPH+9sZOBg4AtgF+FqSi6vq2vGG9RvPAi4Hngk8EtiQ5J+r6iejCiDJ/eh9K/6jPp/b75vUyNZ9LxHb7DGH07tInj6quOZ89mLxvR84qaruGMEX0kFjG/t1sUR8Y7suquoO4KAkuwNnJzmgqubWKMd2TbSIDYAkfwDMAP9pqXM64uhvK3B+Vf28qm4BLgJGVkRt4RX0pgyqqq4Hvg08ZlQfnuSe9C7eM6vqs30O2Qo8fM7+3vS+BXYhNpIcCJwOHFtVPxhFXAPENwN8Ksl3gOOADyb5vY7ENtbrokV8Y70uAKrqR8CXgKPmvfWbayLJzsD9WXw6epSxkWQdcDJwTFXdttS5TBz9nQP8xyQ7J7kv8BR6c6pd8V163/pI8mDg0cANo/jgZl72DODqqvqLBQ47Fzi+WV31VODHVfW9LsSWZC3wWeBlox5Btomvqh5RVftW1b705sJfW1Wf60JsjPG6aBnfWK6LJHs23+ZJsguwDrhm3mHnAic028cB/1RNNXrcsSV5InAqvaTRqm41lVNVSf4P8AxgTZKtwNvoFY2oqg9X1dVJzgeuAO4ETu83tBtXfMC7gL9JciW9IfBJzTfAUTgUeBlwZTNvCr3VLGvnxHcevZVV1wO/oPdNsCuxvZXe3PIHm6mg7TW67qVt4huXJWMb83XR5u9uXNfFQ4GPJdmJ3pfxT1fVPyR5J7Cpqs6ll/Q+keR6eiONF48grraxvZfeYoK/ba6J71bVMYud1JYjkqSBOFUlSRqIiUOSNBAThyRpICYOSdJATBySpIFM5XJcTackDwS+0Ow+BLgDuLnZP6Sqbh9LYItI8krgvKr6/rhjkWa5HFdTKcnbgZ9V1fs6EMtOTVuIfu99BXh9VV3e7/0Ffmfnph+SNBROVUlAkhOa5xZcnuSDSe7R3CH9oyTvTXJZkguSPCXJl5PckOTo5ndfneTs5v1/SfKWlud9d5JLgUOSvCPJ15NcleTDzV33LwIOAs5qfv9eSbbOuRP4qUk2NtvvTnJqkg3AXzef8RfNZ1+R5NWj/1vVamXi0NRLcgC9jrRPq6qD6E3hzt7Ze3/gwqp6EnA78HZ6bS1eALxzzmkOaX7nScBLkhzU4ryXVdUhVfU14JSqejLw+Oa9o6rqLHpN+15UVQe1mEp7IvCcqnoZcCJwU1UdAjwZeF3TbkVaNmscUq9/z5OBTU3LhV34bVv4X1bVhmb7Snp9t7Y3bS32nXOOC2af85Hkc/S67u68yHlvB86e8/tHJHkjcB9gDbCZwVu+n1NVv2q2jwQem2RuotqPXj8naVlMHFKvr9FHq+p/3uXFXhfTud/y7wRum7M99/qZXyysJc77y9kmd03DwA/Qe8LdvyV5N70E0s92fjtTMP+Yn8/73/TaqvoC0gpzqkqCjcALk6yB3uqrHZjWOTLJ7k0SOBb4vwOcdxd6ieiWJLsBz5/z3k/pPSp11nfoPRODecfNdwHw2iZJkeTRTXdUadkccWjqVdWVSd4BbExyD+DXwGsY7BkiXwE+Se8BQp+YXQXV5rxV9YMkHwOuArZw1yfb/TVwepJf0qujvB34SJLvA5cuEs+p9DrHXt5Mk91EL6FJy+ZyXGmZmhVLB1TVH407FmkUnKqSJA3EEYckaSCOOCRJAzFxSJIGYuKQJA3ExCFJGoiJQ5I0kP8P/u5Y/gbJaoYAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZsAAAEGCAYAAACzYDhlAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAeu0lEQVR4nO3df7RdZX3n8ffHRBAqvyRXiwmQrCFoo1WBY6DaMmpQgrVEZ3CMjpJSnKwqlNhOW2FsFdGuVZYdFQZQGYIGKk2YiHrbZYkJqNQRQm4oikkM3ILIFToEiYgoYOJ3/tjPhcPh3HP2uec8Z5977+e11lnZZ+9n7/09N3ff73l+7GcrIjAzM8vpOVUHYGZm05+TjZmZZedkY2Zm2TnZmJlZdk42ZmaW3eyqAxgEc+bMifnz51cdhpnZlLJ169aHImKoTFknG2D+/PmMjIxUHYaZ2ZQi6d6yZd2MZmZm2TnZmJlZdk42ZmaWnZONmZll52RjZmbZOdmYmVl2TjZmZpadk42ZmWXnZGNmZtk52ZiZWXZONmZmlp2TjZmZZedkY2Zm2TnZmJlZdk42ZmaWnZONmZll52RjZmbZOdmYmVl2lSYbSUsl7ZQ0KuncJtv3lbQubd8saX5af6ikb0j6uaRLGvY5TtIdaZ+LJak/n8bMzCZSWbKRNAu4FDgFWAS8U9KihmJnArsj4ijgU8CFaf3jwF8Df97k0J8BVgIL02tp76M3M7NOVFmzWQyMRsTdEfEksBZY1lBmGbAmLa8HlkhSRDwWEd+mSDpPkXQYcGBE3BwRAVwFvDXrpzAzs7aqTDZzgfvq3o+ldU3LRMQe4BHg0DbHHGtzTAAkrZQ0Imlk165dHYZuZmadqDLZNOtLiUmUmVT5iLg8ImoRURsaGmpxSDMz61aVyWYMOLzu/Tzg/onKSJoNHAQ83OaY89oc08zM+qzKZLMFWChpgaR9gOXAcEOZYWBFWj4NuDH1xTQVEQ8Aj0o6IY1COx34au9DNzOzTsyu6sQRsUfS2cAGYBZwZURsk3QBMBIRw8Bq4GpJoxQ1muXj+0v6IXAgsI+ktwJviojtwPuALwD7Af+cXmZmViG1qCjMGLVaLUZGRqoOw8xsSpG0NSJqZcp6BgEzM8vOycbMzLJzsjEzs+ycbMzMLDsnGzMzy87JxszMsnOyMTOz7JxszMwsOycbMzPLzsnGzMyyc7IxM7PsnGzMzCw7JxszM8vOycbMzLJzsjEzs+ycbMzMLDsnGzMzy87JxszMsnOyMTOz7JxszMwsOycbMzPLrtJkI2mppJ2SRiWd22T7vpLWpe2bJc2v23ZeWr9T0sl16/9U0jZJ35f0D5Ke159PY2ZmE6ks2UiaBVwKnAIsAt4paVFDsTOB3RFxFPAp4MK07yJgOfAyYClwmaRZkuYC5wC1iHg5MCuVMzOzClVZs1kMjEbE3RHxJLAWWNZQZhmwJi2vB5ZIUlq/NiKeiIh7gNF0PIDZwH6SZgP7A/dn/hxmZtZGlclmLnBf3fuxtK5pmYjYAzwCHDrRvhHxY+DvgB8BDwCPRMTXs0RvZmalVZls1GRdlCzTdL2kQyhqPQuAFwO/IendTU8urZQ0Imlk165dHYRtZmadqjLZjAGH172fx7ObvJ4qk5rFDgIebrHvScA9EbErIn4FXAe8ptnJI+LyiKhFRG1oaKgHH8fMzCZSZbLZAiyUtEDSPhQd+cMNZYaBFWn5NODGiIi0fnkarbYAWAjcStF8doKk/VPfzhJgRx8+i5mZtTC7qhNHxB5JZwMbKEaNXRkR2yRdAIxExDCwGrha0ihFjWZ52nebpGuB7cAe4KyI2AtslrQeuC2t/1fg8n5/NjMzeyYVFYWZrVarxcjISNVhmJlNKZK2RkStTFnPIGA2CVvv3c3pqzez9d7dVYdiNiU42ZhNwkWb7uSmux7iok13Vh2K2ZRQWZ+N2VS26qSjn/GvmbXmZGM2CccdeQhXnXl81WGYTRluRjMzs+ycbMzMLDsnGzMzy87JxszMsnOyMTOz7JxszMwsOycbMzPLzsnGzMyyc7IxM7PsnGzMzCw7JxszM8vOycbMzLJzsjEzs+zaJhtJI5LOknRIPwIyM7Ppp0zNZjnwYmCLpLWSTpakzHGZmdk00jbZRMRoRHwIOBq4BrgS+JGkj0p6Qe4Azcxs6ivVZyPpFcD/BD4BfAk4DfgZcGO+0MzMbLpo+6ROSVuBnwKrgXMj4om0abOk1+YMzszMpocyNZu3R8SSiLhmPNFIWgAQEf+pm5NLWippp6RRSec22b6vpHVp+2ZJ8+u2nZfW75R0ct36gyWtl/QDSTsk/U43MZqZWffKJJv1Jdd1RNIs4FLgFGAR8E5JixqKnQnsjoijgE8BF6Z9F1EMXHgZsBS4LB0P4CLg+oh4KfBKYEe3sZqZWXcmbEaT9FKKP+YHSaqvwRwIPK8H514MjEbE3el8a4FlwPa6MsuA89PyeuCSNBJuGbA21bTukTQKLJa0DTgR+EOAiHgSeLIHsZqZWRda9dm8BHgLcDDwB3XrHwX+Ww/OPRe4r+79GHD8RGUiYo+kR4BD0/pbGvadC/wS2AV8XtIrga3Aqoh4rPHkklYCKwGOOOKIHnwcMzObyITNaBHx1Yg4A3hLRJxR9zonIr7Tg3M3u1cnSpaZaP1s4FjgMxFxDPAY8Ky+IICIuDwiahFRGxoaKh+12SRsvXc3p6/ezNZ7d1cdilklJkw2kv4yLb5L0sWNrx6ceww4vO79POD+icpImg0cBDzcYt8xYCwiNqf16ymSj1mlLtp0Jzfd9RAXbbrzWduciGwmaNWMNt6xPpLp3FuAhWlk248pOvzf1VBmGFgB3Exxb8+NERGShoFrJH2SYnaDhcCtEbFX0n2SXhIRO4ElPLMPyKwSq046+hn/1htPRABXndnYkmw2PUyYbCLiH9O/a3KcOPXBnA1sAGYBV0bENkkXACMRMUxxb8/VaQDAwxQJiVTuWopEsgc4KyL2pkP/CfBFSfsAdwNn5IjfrBPHHXnIhImkVSIymy4U0dhN0lBA2khxr81P0/tDKEaCndxyxymkVqvFyEiuCpyZ2fQkaWtE1MqULXOfzdB4ogGIiN3ACycbnJl1x308NhWVSTZ7JT01NljSkTx71JiZ9UmrwQZmg6rt3GjAh4BvS/pWen8i6f4UM2tu6727uWjTnaw66WiOO7K3j4JyH49NRWUeMXA9xfDhdcC1wHERsSF3YGbdqLqpKWftY3ywQa+TmFlOZWo2AHuBBymmqVkkiYi4KV9YZt2pejixax9mz1TmEQPvBVZR3Dh5O3ACxX0vb8gbmtnkVf3HvtVQZ7OZqMwAgVXAq4F7I+L1wDEU84+ZDSw3NZkNljLJ5vGIeByK58tExA8oJuk0m7Kq7tMxm2nK9NmMSToY+AqwUdJunj2HmdmUUnWfjtlM0zbZRMTb0uL5kr5BMRnm9VmjMsus6j4ds5mmTDMako6VdA7wCopZlf1AMpvSquzTcROezURtk42kDwNrKB5aNofiwWR/lTsws+nKMwDYTFSmz+adwDF1gwT+FrgN+HjOwMymKzfh2UxUJtn8kOJmzsfT+32Bf8sVkNl053twbCYqk2yeALalRw0E8EaKudIuBoiIczLGZ2Zm00CZZPPl9Br3zTyhmJnZdFVm6HOWJ3Wa9VPOWZhz63XsU/lnYVNXmdFoCyWtl7Rd0t3jr34EZ9YrU3kEWK9jn8o/C5u6ytxn83ngM8Ae4PXAVcDVOYMy67VVJx3NiQvnZBsBlvPemU5jbxdL7p+FWTOKaP3QzfSM6eMk3RERv53W/UtE/F5fIuyDWq0WIyMjVYdhU9jpqzdz010PceLCOZWPNBukWGx6S/mhVqZsqYk4JT0HuEvS2ZLeBrywqwjNBlynNZVe1ha6rSW55mKDqEyy+QCwP3AOcBzwbmBFL04uaamknZJGJZ3bZPu+ktal7Zslza/bdl5av1PSyQ37zZL0r5L+qRdx2szTab9GL6e/6bZPxY9XsEE04Wg0Sc8DDoiILWnVz4EzJL0IeKTbE0uaBVxKcd/OGLBF0nBEbK8rdiawOyKOkrQcuBB4h6RFwHLgZcCLgU2Sjo6IvWm/VcAO4MBu47SZqcq7/D3DgE1HrWo2FwPN+mVOAj7Vg3MvBkYj4u40sedaYFlDmWUU87IBrAeWSFJavzYinoiIe4DRdDwkzQN+H7iiBzHaDFVl7cA1E5uOWiWb342I6xpXRsQXgRN7cO65wH1178fSuqZlImIPRY3q0Db7fhr4S+DXrU4uaaWkEUkju3b5waNTjWdONptaWiUbTXK/spodv3Fo3ERlmq6X9BbgwYjY2u7kEXF5RNQiojY0NNQ+Whso/b5XZJCT2yDHZjauVdJ4UNLixpWSXg30oiowBhxe934ez34C6FNlJM2meHDbwy32fS1wqqQfUjTLvUHS3/cgVhsw/R5xNcg3Qg5ybGbjWk1X8xfAtZK+AIzXFGrA6RSd893aAiyUtAD4cTrmuxrKDFOMfLsZOA24MSJC0jBwjaRPUgwQWAjcGhE3A+cBSHod8OcR8e4exGoDpt8zJw9yp/0gx2Y2bsJkExG3pprNWcAfptXbgOMj4sFuTxwReySdDWwAZgFXRsQ2SRcAIxExDKwGrpY0SlGjWZ723SbpWmA7xcwGZ9WNRDPruUF+LMAgx2Y2ru0MAjOBZxAwM+tcr2cQMJv23MlulpeTjRnTu5PdidQGQZlHDLy9zDqznHL/wZzO84l1m0idrKwXytRsziu5zmzS2v1By13zmM537XebSKdzrc/6p9XcaKcAbwbmSrq4btOBFCPAzHpm/A8a0HRklYf3Tl63o9U6/dn7SaDWTKv7bO4HRoBTefo+G4BHgT/NGZTNPO3+oHl4b3U6/dm3++JgM1Or+2y+C3xX0jUR8as+xmQzkJPJ9OFaqDVTps9mvqT1krZLunv8lT0ysw64E3twTOf+L5u8Msnm88BnKPppXg9cBVydMyizTrkTe/K6TdRO9FZGmWSzX0TcQDHbwL0RcT7whrxhmXVmOg9dzq3TRN2YXJzorYxWAwTGPS7pOcBdaS6zHwMvzBuWWWfc5zN5nfaxNA4AcB+NldF2brT0SIEdwMHAxyiGPn8iIm7JH15/eG40s/I8tNnGdTI3WumJOCX9RkQ81lVkA8rJpv/8B8ts6uvpRJySfkfSdoraDZJeKemyLmO0Gcbt/GYzW5kBAp8GTgZ+Ak/df3NizqBs+mlMLu069D3CqTr+2VsOpWZ9joj7Glb5QWXWkcbk0u5eDNd8qtPrn31j8nIym5nKjEa7T9JrgJC0D3AOqUnNrKxOR4t5hFN1ev2zbxy95ulsZqYyyeaPgYuAucAY8HWKR0WbZeOhzNXp9c++MXk1/tvNYBEPNJk6/FhoPBptshovdF/4Nhmnr97MTXc9xIkL53Sc5LrZ17rXyWi0Vo8Y+F/AhJkoIs6ZRGw2jbRrHnHysTLqazqd/s64uXXqaDVAYITi0QJbefoxA/Uvm2EaO3YbO/0b37uT38qoHyzS6e9M40CTdoMPPDihOq0eMbBmfFnSB+rf28zQ+C2zsebS2Lbf+N7fOq1T3f7OtKtde3BCdUoNfaZFc1o3JC2VtFPSqKRzm2zfV9K6tH2zpPl1285L63dKOjmtO1zSNyTtkLRN0qoccc8Und4b08hTzVun2tVU2tVM2tWuPWFrdcqMRstC0izgUuCNFKPctkgajojtdcXOBHZHxFGSlgMXAu+QtAhYDrwMeDGwSdLRFI9B+O8RcZukA4CtkjY2HNNKavyW6RFi1m+dDptuV7v273B1Wg0QeJSnazT7S/rZ+CYgIuLALs+9GBiNiLvT+dYCy4D6xLAMOD8trwcukaS0fm1EPAHcI2kUWBwRNwMPUAT4qKQdFEO2nWwmwRemVa3dsOl2uv0d9iCX3pmwGS0iDoiIA9Nrdt3yAT1INFAkgfqZCcbSuqZlImIP8AhwaJl9U5PbMcDmZieXtFLSiKSRXbt2TfpDTGW+s9sGXWOzWr+bZj3IpXfK9tnkoCbrGvuGJirTcl9Jzwe+BHwgIn7WpCwRcXlE1CKiNjQ0VDLk6aXxQvKFZTNNp31A7fb3F7aJVdZnQ1EbObzu/Tzg/gnKjEmaDRwEPNxqX0nPpUg0X4yI6/KEPjU1Ngl020RhNtV12gfUbn+PdptYlTWbLcBCSQvSnGvLgeGGMsPAirR8GnBjFFMeDAPL02i1BcBC4NbUn7Ma2BERn+zLpxhg7ab1r7qJwqxqjTWXTu/TaXevWTszqSZUWbJJfTBnAxsoJva8NiK2SbpA0qmp2Grg0DQA4M+Ac9O+24BrKTr+rwfOioi9wGuB9wBvkHR7er25rx9sgHjYp1lrjV+w2jUld/qFrV0z20xquq6yGY2I+BrwtYZ1H65bfhx4+wT7/g3wNw3rvk3z/pwZycM+zTrTrim506bmds1svZyUdNB5Ik6mz0Sc0/kX1Wwq6nSy2lYTiw7i9d3Tx0JbdTptz51JVXKzqaDTftH6pu7p1uTmZDPAGn+52rX/uk/GbGprNSlpt4MZOt3ea042A6zdPE/tOivNbOpq9yj1dl9GO92eW6UDBKa7dm2s7dpz283z5PtizKavdgN6Gq//doMP2m3PzQMEyDdAoN1TBBu3+6mDZjZZnX657YWePKnTymn1H9huWKNrKmbWK+1qQlXf+uCaDd3VbBprI62Sj2suZjaduGbTR+3aRVuVNTObKVyzobd9NoN445WZWQ6u2VSo6nZRM7NB5PtszMwsOycbMzPLzsnGzMyyc7IxM7PsnGzMzCw7JxszM8vOycbMzLJzsjEzs+ycbMzMLDsnGzMzy67SZCNpqaSdkkYlndtk+76S1qXtmyXNr9t2Xlq/U9LJZY9pZmb9V1mykTQLuBQ4BVgEvFPSooZiZwK7I+Io4FPAhWnfRcBy4GXAUuAySbNKHtPMzPqsyprNYmA0Iu6OiCeBtcCyhjLLgDVpeT2wRJLS+rUR8URE3AOMpuOVOaaZmfVZlclmLnBf3fuxtK5pmYjYAzwCHNpi3zLHBEDSSkkjkkZ27drVxccwM7N2qkw2arKu8eE6E5XpdP2zV0ZcHhG1iKgNDQ21DNTMzLpTZbIZAw6vez8PuH+iMpJmAwcBD7fYt8wxzcysz6pMNluAhZIWSNqHosN/uKHMMLAiLZ8G3BjFo0WHgeVptNoCYCFwa8ljmplZn1X2pM6I2CPpbGADMAu4MiK2SboAGImIYWA1cLWkUYoazfK07zZJ1wLbgT3AWRGxF6DZMfv92czM7JlUVBRmtlqtFiMjI1WHYWY2pUjaGhG1MmU9g4CZmWXnZGNmZtk52ZiZWXZONmZmlp2TjZmZZedkY2Zm2TnZmJlZdk42ZmaWnZONmZll52RjZmbZOdmYmVl2TjZmZpadk42ZmWXnZGNmZtk52ZiZWXZONmZmlp2TjZmZZedkY2Zm2TnZmJlZdk42ZmaWnZONmZllV0mykfQCSRsl3ZX+PWSCcitSmbskrahbf5ykOySNSrpYktL6T0j6gaTvSfqypIP79ZnMzGxiVdVszgVuiIiFwA3p/TNIegHwEeB4YDHwkbqk9BlgJbAwvZam9RuBl0fEK4A7gfNyfggzMyunqmSzDFiTltcAb21S5mRgY0Q8HBG7KRLJUkmHAQdGxM0REcBV4/tHxNcjYk/a/xZgXs4PYWZm5VSVbF4UEQ8ApH9f2KTMXOC+uvdjad3ctNy4vtEfAf/ck2jNzKwrs3MdWNIm4DebbPpQ2UM0WRct1tef+0PAHuCLLeJbSdEUxxFHHFEyJDMzm4xsySYiTppom6T/J+mwiHggNYs92KTYGPC6uvfzgG+m9fMa1t9fd+wVwFuAJamZbaL4LgcuB6jVahOWMzOz7lXVjDYMjI8uWwF8tUmZDcCbJB2SBga8CdiQmt0elXRCGoV2+vj+kpYCHwROjYhf5P4QZmZWTlXJ5m+BN0q6C3hjeo+kmqQrACLiYeBjwJb0uiCtA3gfcAUwCvwbT/fNXAIcAGyUdLukz/bp85iZWQtq0dI0Y9RqtRgZGak6DDOzKUXS1oiolSnrGQTMzCw7JxszM8vOycbMzLJzsjEzs+ycbMzMLDsnGzMzy87JxszMsvN9NoCkXcC9XRxiDvBQj8LpNcc2eYMc3yDHBoMd3yDHBoMdX2NsR0bEUJkdnWx6QNJI2Rub+s2xTd4gxzfIscFgxzfIscFgx9dNbG5GMzOz7JxszMwsOyeb3ri86gBacGyTN8jxDXJsMNjxDXJsMNjxTTo299mYmVl2rtmYmVl2TjZmZpadk01Jkq6U9KCk77co87r00LZtkr41KLFJOkjSP0r6bortjD7Gdrikb0jakc69qkkZSbpY0qik70k6dsDi+68pru9J+o6kVw5KbHVlXy1pr6TTBim2Cq+JMv+vlVwXkp4n6da68360SZl9Ja1L18RmSfP7EVsH8f2ZpO3pmrhB0pFtDxwRfpV4AScCxwLfn2D7wcB24Ij0/oUDFNv/AC5My0PAw8A+fYrtMODYtHwAcCewqKHMmymetirgBGBzH392ZeJ7DXBIWj6lX/GViS1tmwXcCHwNOG1QYqv4migTXyXXRfo9f35afi6wGTihocz7gc+m5eXAuj7+7MrE93pg/7T8vjLxuWZTUkTcRPHLOJF3AddFxI9S+Qf7EhilYgvgAEkCnp/K7ulTbA9ExG1p+VFgBzC3odgy4Koo3AIcLOmwQYkvIr4TEbvT21uAeYMSW/InwJeAfv7OlYmtymuiTHyVXBfp9/zn6e1z06txpNYyYE1aXg8sSXFmVya+iPhGRPwivS11TTjZ9M7RwCGSvilpq6TTqw6oziXAbwH3A3cAqyLi1/0OIjUFHEPxTaneXOC+uvdjNP+jmlWL+OqdSVEL66uJYpM0F3gb8Nl+x1QXw3ya/9wG4ppoEV9l14WkWZJup/iCsDEiJrwmImIP8AhwaD9iKxlfvVLXxOxeBWfMBo4DlgD7ATdLuiUi7qw2LABOBm4H3gD8B2CjpH+JiJ/1KwBJz6f49v2BJudt9o2tr2Py28Q3Xub1FBfW7w5QbJ8GPhgRe/v0xbeT2Cq/JtrEV9l1ERF7gVdJOhj4sqSXR0R9n2ul10SJ+ACQ9G6gBvzHdsd0zaZ3xoDrI+KxiHgIuAnoS0dyCWdQNGdERIwC9wAv7dfJJT2X4oL/YkRc16TIGHB43ft5FN82+6JEfEh6BXAFsCwifjJAsdWAtZJ+CJwGXCbprQMSW6XXRIn4Kr0uACLip8A3gaUNm566JiTNBg6idVN5Fi3iQ9JJwIeAUyPiiXbHcrLpna8CvydptqT9geMp2okHwY8ovl0i6UXAS4C7+3Hi1M68GtgREZ+coNgwcHoalXYC8EhEPDAo8Uk6ArgOeE+fv5W3jS0iFkTE/IiYT9G2//6I+MogxEaF10TJ+Cq5LiQNpRoDkvYDTgJ+0FBsGFiRlk8DbozUGz8I8Uk6BvgcRaIp1RfnZrSSJP0D8DpgjqQx4CMUHWdExGcjYoek64HvAb8GrmhW7awiNuBjwBck3UFRPf9g+qbZD68F3gPckdqAoRgFdERdfF+jGJE2CvyC4htnv5SJ78MU7eWXpaaqPdGfWXnLxFaVtrFVeU2UiY/qrovDgDWSZlF84b82Iv5J0gXASEQMUyTKqyWNUtRolvchrk7i+wTFoIr/k66JH0XEqa0O6ulqzMwsOzejmZlZdk42ZmaWnZONmZll52RjZmbZOdmYmVl2Hvps1oKkQ4Eb0tvfBPYCu9L7xRHxZCWBtSDpj4CvRcS/Vx2L2TgPfTYrSdL5wM8j4u8GIJZZaUqRZtu+DZwdEbc32z7BPrPTHFxmWbgZzWySJK1Iz/24XdJlkp6T7pb/qaRPSLpN0gZJx0v6lqS7Jb057fteSV9O23dK+quSx/24pFuBxZI+KmmLpO9L+myageEdwKuAdWn/fSSN1d0RfoKkTWn545I+J2kj8Pl0jk+mc39P0nv7/1O16crJxmwSJL2cYrbl10TEqyiapMfv8j4I+HpEHAs8CZxPMS3K24EL6g6zOO1zLPAuSa8qcdzbImJxRNwMXBQRrwZ+O21bGhHrKCaXfEdEvKpEM98xwB9ExHuAlcCDEbEYeDVwVpqqx6xr7rMxm5yTKP4gj6TpOvbj6cck/DIiNqblOyjmetuTpkWZX3eMDePPyZH0FYrZpGe3OO6TwJfr9l8i6S+A5wFzgK10/viDr0bE42n5TcBvSapPbgsp5hAz64qTjdnkCLgyIv76GSuLGXrraxO/Bp6oW66/5ho7TKPNcX85PhljmtjyEoqnUf5Y0scpkk4ze3i6FaOxzGMNn+n9EXEDZj3mZjSzydkE/BdJc6AYtTaJJqc3STo4JY5lwP/t4Lj7USSvhyQdAPznum2PUjwKedwPKZ4rQ0O5RhuA96fEhqSXpFl/zbrmmo3ZJETEHZI+CmyS9BzgV8Af09lzeL4NXEPx4K6rx0ePlTluRPxE0hrg+8C9PPMplJ8HrpD0S4p+ofOB/y3p34FbW8TzOYpZkW9PTXgPUiRBs6556LNZBdJIr5dHxAeqjsWsH9yMZmZm2blmY2Zm2blmY2Zm2TnZmJlZdk42ZmaWnZONmZll52RjZmbZ/X+wxn9+kRPApgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYYAAAEGCAYAAABhMDI9AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAYlElEQVR4nO3de5RlZXnn8e9PkHhDYNnt6ADdzbgalTEKdAlEM3jrRCRrYLJCFIygBmUZxeBlsnQmGQeJ/0ycOMGJeIl4wUlQ43ghDoq2QTFGkG7DnYC9WltazAKUgBcUwWf+OLv0nEN1nV1dtc85VfX9rFWrzt77PXs/XbDPc9733e/7pqqQJGnWgyYdgCRpupgYJEkDTAySpAEmBknSABODJGnA3pMOYKHWrFlTGzZsmHQYkrSsbNu27Y6qWtum7LJLDBs2bGDr1q2TDkOSlpUkO9uWtSlJkjSgs8SQ5H1Jbkty3W6OJ8nbk2xPck2SI7uKRZLUXpc1hg8Ax81z/HnAxubnDOCdHcYiSWqps8RQVZcB35+nyInABdVzObB/ksd2FY8kqZ1J9jEcCNzSt72r2SdJmqBJJobMsW/OGf2SnJFka5Ktt99+e8dhSdLqNsnEsAs4uG/7IODWuQpW1XuqaqaqZtaubfUYriRpD00yMVwEnNY8nXQMcFdVfXeC8UiS6PZx1QuBrwKPT7IryelJXpHkFU2Ri4EdwHbgr4BXdhXLrG077+S0869g284759yWJHU48rmqThlxvIBXdXX9uZy75WYu+8YdAFxw+tEP2N62807O3XIzZ20+lE3rD3jAtiStBqtq5PNZmw/l2I1rOGvzoXNuzyaKc7fc/IBtaxuSVotlN1fSYmxafwAXnH70brf7E8bwb2sbklaLVZUYRpkvccyXNMBmKkkrh4mhpcXUNmDhiUOSJsXEsESWOnFYw5A0KSaGMVlo4rCGIWlSTAxTYrE1DElaKukNJ1g+ZmZmyhXcGNnUZI1CUr8k26pqpk3ZVTWOYSWZrWHMfujPNwYDHHchqT2bklaIxXZuS9IsE8MKsdjObROFpFkmhlXCRCGpLRPDKmWikLQ7JgYBC08UklYuE4PmNCpRSFq5fFxVrfQ/HusU5NLKZmLQgo0aMyFpebMpSQs2asyEpOXNKTEkaRVwSgxN1Kg+CPskpOlmYtCSG9UHYZ+ENN3sY9CSG9UHYZ+ENN3sY5CkVcA+BknSHjMxSJIGmBgkSQNMDJKkASYGSdIAE4MkaYCJQZI0wMSgqTNqygyn1JC6ZWLQxA1/0I+aMsMpNaRuOSWGJm542dDhKTOG15t2Sg2pWyYGTdzwB/3wsqLDiWP4uKSlZWLQxI36oLeGII1Xp30MSY5LclOS7UneOMfxdUkuTfJPSa5JcnyX8Wh56l9vei6u9yAtrc4SQ5K9gHcAzwMOA05JcthQsT8BPlpVRwAnA+d1FY9WLtd7kJZWl01JRwHbq2oHQJIPAycCN/SVKeCRzev9gFs7jEcrlOs9SEurs/UYkpwEHFdVL2u2TwWOrqoz+8o8FvgccADwcGBzVW2b41xnAGcArFu3btPOnTs7iVmSVqppWY8hc+wbzkKnAB+oqoOA44EPJXlATFX1nqqaqaqZtWvXdhCqJGlWl4lhF3Bw3/ZBPLCp6HTgowBV9VXgIcCaDmOSJI3QZWK4EtiY5JAk+9DrXL5oqMy3gecAJHkivcRwe4cxSZJG6CwxVNV9wJnAJcCN9J4+uj7JOUlOaIq9Hnh5kquBC4GX1HJbhFqSVphOB7hV1cXAxUP73tT3+gbg6V3GIElaGCfRkyQNMDFo1XFktDQ/E4NWHUdGS/NzEj2tOo6MluZnYtCq47Td0vxsSpIkDTAxSJIGmBgkSQNMDJKkASYGSdIAE4MkaYCJQZI0wMQgSRpgYpAkDTAxSJIGmBi04jmbqrQwJgateM6mKi2Mk+hpxXM2VWlhrDFoxZudTXXT+gPmPG5TkzTIxKBVz6YmaZBNSVr1bGqSBrVKDEmeDpwNrG/eE6Cq6t91F5o0Hi7cIw1qW2M4H3gtsA24v7twJEmT1jYx3FVVn+k0EknSVGibGC5N8lbg48BPZ3dW1dc7iUqSNDFtE8NsA+xM374Cnr204UiSJq1VYqiqZ3UdiCRpOrQax5BkvyRvS7K1+fnzJPt1HZwkafzaDnB7H/AD4PnNz93A+7sKSpI0OW37GB5XVb/Tt/3mJFd1EZAkabLa1hjuSfLrsxvNgLd7uglJkjRJbWsMfwB8sOlXCPB94CVdBSVJmpy2TyVdBTwlySOb7bs7jUqSNDHzJoYkL6qq/5PkdUP7Aaiqt3UYmyRpAkb1MTy8+b3vHD+PGHXyJMcluSnJ9iRv3E2Z5ye5Icn1Sf5mAbFLkjowb42hqt7dvNxSVV/pP9Z0QO9Wkr2AdwC/AewCrkxyUVXd0FdmI/BfgKdX1Z1JHr0H/wZJ0hJq+1TS/265r99RwPaq2lFV9wIfBk4cKvNy4B1VdSdAVd3WMh5JUkdG9TH8GvA0YO1QP8Mjgb1GnPtA4Ja+7V38cs6lWYc21/lKc76zq+qzc8RxBnAGwLp160ZcVpK0GKNqDPvQ60vYm8H+hbuBk0a8N3Psq6HtvYGNwDOBU4D3Jtn/AW+qek9VzVTVzNq1a0dcVpK0GKP6GL4EfCnJB6pq5wLPvQs4uG/7IODWOcpcXlU/A76Z5CZ6ieLKBV5LkrRE2g5w+3GzHsO/Bx4yu7Oq5pt2+0pgY5JDgO8AJwMvHCrzSXo1hQ8kWUOvaWlHy5gkSR1o2/n818A/A4cAbwa+xYhv9VV1H3AmcAlwI/DRqro+yTlJTmiKXQJ8L8kNwKXAH1XV9xb8r5AkLZlUDTf7z1Eo2VZVm5JcU1VPbvZ9qaqe0XmEQ2ZmZmrr1q3jvqwkLWvN5/jM6JLtm5J+1vz+bpLfotdXcNCeBCdJmm5tE8Nbmgn0Xk9v/MIjgdd2FpUkaWLaJoarq+ou4C7gWQBJHtNZVJKkiWnb+fzNJBcmeVjfvou7CEiSNFltE8O1wJeBLyd5XLNvrgFskqRlrm1TUlXVeUmuBv4uyRt44ChmSdIK0DYxBKCqvpLkOcBHgCd0FpUkaWLaJobjZ19U1XeTPJve5HqSpBWm1QpuwCmzq7YNuayTqCRJEzOqxtC/gtsw+xgkaQXqbAU3SdLy1OUKbtKytG3nnZx2/hVs23nnHh2XlrsuV3CTlqVzt9zMZd+4A4ALTh9edHD0cWm5G9XHMLyC26w2K7hJy9JZmw8d+L3Q49Jy13ba7fVVtTPJw6vqR2OIa7ecdluSFm4h02637WP4t81iOjc2F3hKkvP2NEBJ0vRqmxj+Angu8D2AqroaOLaroCRJk9M2MVBVtwztun+JY5EkTYG2U2LckuRpQCXZB/hDmmYlSdLK0rbG8ArgVcCBwC7g8GZbkrTCtKoxVNUdwO91HIskaQq0SgxJ1gIvBzb0v6eqfr+bsCRJk9K2j+FT9FZw24KdzpK0orVNDA+rqjd0GokkaSq07Xz+dJLjRxeTJC13bRPDWfSSwz1J7k7ygyR3dxmYNC2cTVWrTdunkuZaqEdaFZxNVatN26eSjpxj913Azqq6b2lDkqaLs6lqtWnb+XwecCRwbbP9q8DVwKOSvKKqPtdFcNI02LT+AGsKWlXa9jF8CziiqjZV1SZ6I5+vAzYDf9ZRbJKkCWibGJ5QVdfPblTVDfQSxY5uwpIkTUrbpqSbkrwT+HCz/QLg5iS/Avysk8gkSRPRtsbwEmA78BrgtcCOZt/PgGd1EZgkaTLaPq56D/Dnzc+wHy5pRJKkiWpVY0iyMcnHktyQZMfsT4v3HZfkpiTbk7xxnnInJakkrdYjlSR1p21T0vuBdwL30Ws6ugD40HxvSLIX8A7gecBhwClJDpuj3L70Fv65on3YkqSutE0MD62qLwCpqp1VdTbw7BHvOQrYXlU7qupeeh3XJ85R7k/pPfL6k5axSJI61DYx/CTJg4BvJDkzyW8Djx7xngOB/nWidzX7fiHJEcDBVfXp+U6U5IwkW5Nsvf3221uGLEnaE20Tw2uAh9Fr8tkEnAq8eMR7Mse++sXBXqL5X8DrR128qt5TVTNVNbN27dqWIUuS9kTbp5KubF7+EHhpy3PvAg7u2z4IuLVve1/gScAXkwA8BrgoyQlVtbXlNSRJS2zexJDkovmOV9UJ8xy+EtiY5BDgO8DJwAv73nsXsKbvWl8E/rNJQZIma1SN4dfo9RNcSO+pobmah+ZUVfclORO4BNgLeF9VXZ/kHGBrVc2bdCRJkzEqMTwG+A3gFHrf9v8fcGH/vEnzqaqLgYuH9r1pN2Wf2eackqRuzdv5XFX3V9Vnq+rFwDH0psX4YpJXjyU6SdLYjex8bibK+y16tYYNwNuBj3cbliRpUkZ1Pn+Q3pNDnwHeXFXXjSUqSdLEjKoxnAr8CDgU+MPmsVLodUJXVT2yw9gkSRMwb2KoqrYD4CRJK4Qf/JKkASYGSdIAE4MkaYCJQZI0wMQgSRpgYpAkDTAxSJIGmBgkSQNMDJKkASYGaYy27byT086/gm0775x0KNJumRikMTp3y81c9o07OHfLzZMORdotE4O0hEbVCM7afCjHblzDWZsPHXNkUnsj12OQ1N5sjQDggtOPfsDxTesPmHO/NE2sMUhLaLhGYJ+CliMTg7SEZmsEm9YfAIzuUzBxaBrZlCR1aLbmsLs+hVFNT9IkmBikDo3qUxiVOKRJMDFIE2RntKaRfQySpAEmBknSABODJGmAiUGSNMDEIE0xxzloEkwM0hRz0j1Ngo+rSlPMcQ6aBBODNMUc56BJsClJkjTAxCBJGmBikCQN6DQxJDkuyU1Jtid54xzHX5fkhiTXJPlCkvVdxiNJGq2zxJBkL+AdwPOAw4BTkhw2VOyfgJmqejLwMeDPuopH6so4xxo4rkHj0GWN4Shge1XtqKp7gQ8DJ/YXqKpLq+rHzeblwEEdxiN1YpxjDRzXoHHo8nHVA4Fb+rZ3AfM9d3c68Jm5DiQ5AzgDYN26dUsVn7QkxjnWwHENGodUVTcnTn4XeG5VvazZPhU4qqpePUfZFwFnAs+oqp/Od96ZmZnaunVrFyFL0oqVZFtVzbQp22WNYRdwcN/2QcCtw4WSbAb+mBZJQZLUvS77GK4ENiY5JMk+wMnARf0FkhwBvBs4oapu6zAWSVJLnSWGqrqPXvPQJcCNwEer6vok5yQ5oSn2VuARwN8muSrJRbs5nSRpTDqdK6mqLgYuHtr3pr7Xm7u8viRp4Rz5LEkaYGKQVpDFDoBzAJ3AxCCtKAsdADecCBxAJ3A9BmlFWegAuNlEAHDB6Uc7gE5AhwPcuuIAN6m9bTvv5NwtN3PW5kPZtP6ABR9fbHlNj4UMcLMpSVrBhpuGhpuOZleIa/shb1PT6mBTkrSMjfoGP9w0NNx0tFA2Na0OJgZpGRv1QT+8ZvRiP9gXsga1zU7Ll4lBWsYW+kG/kA/2Nub78F9s7USTYx+DtIwttI9goUaNa5ivD+OszYdy7MY1NjstQ9YYJO3WqG/9o/ow+t9j09LyYY1BmiLTNvJ41Lf+4RrLfOVHPSGl6WGNQZoi09Yuv9A+ifnKL/YJKWsc42NikKbISn4cdLFPSE1b0lzJbEqSpkjXncnTZPjfOqppabiZygkDu2NikDQVRo2qHk4kix2FvdgJB1cyE4OkqbDQGsFiH4dd6PVW03QgJgZJU2GhNYJRTVGjthd6vdU0LsPEIGkqLfSDePiDfdT2qOuNSiRL2bQ06lzjbsbyqSRJU2mhj8oOP+U06veo6416CmrU8f7Ha2fL7+5R28Vea6mZGCStCMMf7KO2RxmVSEYd7/8wBwY+2IfHZCz2WkvNhXokaQkMf9jPV2M47fwruOwbd3DsxjUjaxtL9ejyQhbqscYgSUtguLlnuIaykMF9kx7MZ2KQpCWwkOaeUc1akx4Bb1OSJK0CrvksSdpjJgZJ0gATgyRpgIlBkjTAxCBJGmBikCQNMDFIkgaYGCRJA5bdALcktwM7F3GKNcAdI0tNxjTHBtMdn7HtuWmOb5pjg+mObzi29VW1ts0bl11iWKwkW9uO/hu3aY4Npjs+Y9tz0xzfNMcG0x3fYmKzKUmSNMDEIEkasBoTw3smHcA8pjk2mO74jG3PTXN80xwbTHd8exzbqutjkCTNbzXWGCRJ8zAxSJIGrMjEkOR9SW5Lct08ZZ6Z5Kok1yf50rTElmS/JH+X5OomtpeOMbaDk1ya5Mbm2mfNUSZJ3p5ke5Jrkhw5ZfH9XhPXNUn+MclTpiW2vrJPTXJ/kpPGEdtC4pvEfdHyv+sk74uHJPla37XfPEeZX0nykea+uCLJhimK7XVJbmjuiS8kWT/yxFW14n6AY4Ejget2c3x/4AZgXbP96CmK7b8C/6N5vRb4PrDPmGJ7LHBk83pf4GbgsKEyxwOfAQIcA1wxxr9dm/ieBhzQvH7euOJrE1tzbC/g74GLgZOm7G83kfuiZWyTvC8CPKJ5/WDgCuCYoTKvBN7VvD4Z+MgUxfYs4GHN6z9oE9uKrDFU1WX0/sfZnRcCH6+qbzflbxtLYLSKrYB9kwR4RFP2vjHF9t2q+nrz+gfAjcCBQ8VOBC6onsuB/ZM8dlriq6p/rKo7m83LgYOmJbbGq4H/C4zt/7kmpjbxTeS+aBnbJO+LqqofNpsPbn6Gn9o5Efhg8/pjwHOaWCceW1VdWlU/bjZb3RMrMjG0cChwQJIvJtmW5LRJB9TnL4EnArcC1wJnVdXPxx1EUxU+gt43kH4HArf0be9i7g/ATs0TX7/T6dVuxmp3sSU5EPht4F3jjmkojg3M/beb+H0xT2wTvS+S7JXkKnoJ/fNVtdv7oqruA+4CHjUlsfVrdU/svVTBLTN7A5uA5wAPBb6a5PKqunmyYQHwXOAq4NnA44DPJ/lyVd09rgCSPILet9rXzHHdub4FjfWZ5xHxzZZ5Fr2b4NenKLa/AN5QVfeP4cvknEbEN9H7YkRsE70vqup+4PAk+wOfSPKkqurvJ5zYfdEiNgCSvAiYAZ4x6pyrtcawC/hsVf2oqu4ALgPG0knZwkvpVeerqrYD3wSeMK6LJ3kwvZvzr6vq43MU2QUc3Ld9EL1vcWPRIj6SPBl4L3BiVX1vimKbAT6c5FvAScB5Sf7TFMU3sfuiRWwTvS9mVdW/Al8Ejhs69Iv7IsnewH7M32Q8zthIshn4Y+CEqvrpqHOt1sTwKeA/JNk7ycOAo+m1a06Db9P7xkaSfwM8Htgxjgs3baLnAzdW1dt2U+wi4LTm6aRjgLuq6rvTEl+SdcDHgVPHWQNsE1tVHVJVG6pqA7126FdW1SenJT4mdF+0jG2S98Xa5ts4SR4KbAb+eajYRcCLm9cnAX9fTW/vpGNLcgTwbnpJoVW/0YpsSkpyIfBMYE2SXcB/p9cpQ1W9q6puTPJZ4Brg58B756p6TSI24E+BDyS5ll719A3Nt7dxeDpwKnBt02YJvadB1vXFdzG9J5O2Az+m901uXNrE9yZ6bbvnNc0199V4Zr9sE9skjYxvgvdFm7/dJO+LxwIfTLIXvS/TH62qTyc5B9haVRfRS2wfSrKdXk3h5CmK7a30Ouz/trknvl1VJ8x3UqfEkCQNWK1NSZKk3TAxSJIGmBgkSQNMDJKkASYGSdKAFfm4qlanJI8CvtBsPga4H7i92T6qqu6dSGDzSPL7wMVV9S+TjkWa5eOqWpGSnA38sKr+5xTEslczbcFcx/4BOLOqrprr+G7es3czH4/UCZuStCokeXEzb/1VSc5L8qBmhO+/Jnlrkq8nuSTJ0Um+lGRHkuOb974sySea4zcl+ZOW531Lkq8BRyV5c5Irk1yX5F3NyPEXAIcDH2nev0+SXX0jWY9JsqV5/ZYk707yeeD9zTXe1lz7miQvG/9fVSuViUErXpIn0ZvV9GlVdTi9JtTZkan7AZ+rqiOBe4Gz6U298LvAOX2nOap5z5HAC5Mc3uK8X6+qo6rqq8C5VfVU4FebY8dV1UfoTQz3gqo6vEVT1xHAf6yqU4EzgNuq6ijgqcCrmulApEWzj0GrwWZ6H55bmykBHsovpw6/p6o+37y+lt7cT/c1Uy9s6DvHJbPrPCT5JL1ZW/ee57z3Ap/oe/9zkvwR8BBgDbCNhU8J/qmq+knz+jeBJybpT0Qb6c0pJC2KiUGrQYD3VdV/G9jZmwWz/1v6z4Gf9r3uvz+GO+NqxHnvmZ1ErZmQ7i/prVL2nSRvoZcg5nIfv6zJD5f50dC/6ZVV9QWkJWZTklaDLcDzk6yB3tNLe9Ds8ptJ9m8+5E8EvrKA8z6UXqK5I8m+wO/0HfsBveUsZ32L3poIDJUbdgnwyiYJkeTxzeya0qJZY9CKV1XXprdI+pYkDwJ+BryCha0j8Q/A39BbJOZDs08RtTlvVX0vyQeB64CdDK5O9n7gvUnuodePcTbwV0n+BfjaPPG8m97so1c1zVi30UtY0qL5uKo0QvPEz5Oq6jWTjkUaB5uSJEkDrDFIkgZYY5AkDTAxSJIGmBgkSQNMDJKkASYGSdKA/w+R4YJwf11fcgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZIAAAEGCAYAAABPdROvAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAc3ElEQVR4nO3df5xddX3n8debBAL+4EeTWNkESFYSbcQWYQz0hygSIFghVOISRaA0PuIqWHy47ha6xVWq7dJ2sVhRQAIFukooSB0rGgkg+AMCM5ACgSWMkUgAl6REFlHA4Gf/ON8hl+Ode87MmXPvnbnv5+Mxj3t+fM+Zz03umc/9/jjfo4jAzMxsrHbqdABmZjaxOZGYmVklTiRmZlaJE4mZmVXiRGJmZpVM7XQA7TBjxoyYM2dOp8MwM5tQBgcHt0bEzKJyPZFI5syZw8DAQKfDMDObUCRtKlPOTVtmZlaJE4mZmVXiRGJmZpXUmkgkLZb0kKQhSWc12T9N0qq0f62kOWn7kZIGJd2XXt/RcMx30jnXpZ/X1PkezMystdo62yVNAS4EjgQ2A3dJ6o+IBxqKLQe2RcT+kpYB5wEnAluBYyPicUkHAKuBWQ3HnRQR7j03M+sCddZIFgJDEbExIl4ArgaW5MosAa5Iy9cCR0hSRNwTEY+n7euBXSVNqzFWMzMbozoTySzg0Yb1zby8VvGyMhGxHXgamJ4rcwJwT0Q837Dt8tSsdY4kNfvlklZIGpA0sGXLlirvw8zMWqgzkTT7A5+fs75lGUlvJGvu+mDD/pMi4k3AW9PPyc1+eURcEhF9EdE3c2bh/TRmZjZGdSaSzcA+DeuzgcdHKiNpKrAH8FRanw1cD5wSET8cPiAiHkuvzwBfJmtCMzOzDqkzkdwFzJM0V9IuwDKgP1emHzg1LS8Fbo6IkLQn8A3g7Ij4/nBhSVMlzUjLOwPvAu6v8T2YmVmB2hJJ6vM4g2zE1YPANRGxXtK5ko5LxVYC0yUNAR8DhocInwHsD5yTG+Y7DVgt6V5gHfAY8KW63oOZmRVTLzxqt6+vLzzXlpnZ6EgajIi+onK+s93MzCpxIjEzs0qcSMzMrBInEjMzq8SJxMzMKnEiMTOzSpxIzMysEicSMzOrxInEzMwqcSIxM7NKnEjMzKwSJxIzM6vEicTMzCpxIjEzs0qcSMzMrBInEjMzq8SJxMzMKnEiMTOzSpxIzLrM4KZtnLJyLYObtnU6FLNSnEjMuswFazZw28NbuWDNhk6HYlbK1E4HYGYvd+ai+S97Net2TiRmXebg/fbiyuWHdDoMs9LctGVmZpU4kZiZWSVOJGZmVokTiZmZVVJrIpG0WNJDkoYkndVk/zRJq9L+tZLmpO1HShqUdF96fUfDMQen7UOSPidJdb4HMzNrrbZEImkKcCFwDLAAeK+kBbliy4FtEbE/8FngvLR9K3BsRLwJOBW4quGYLwIrgHnpZ3Fd78HMzIrVWSNZCAxFxMaIeAG4GliSK7MEuCItXwscIUkRcU9EPJ62rwd2TbWXvYHdI+L2iAjgSuD4Gt+DmZkVqDORzAIebVjfnLY1LRMR24Gngem5MicA90TE86n85oJzAiBphaQBSQNbtmwZ85swM7PW6kwkzfouYjRlJL2RrLnrg6M4Z7Yx4pKI6IuIvpkzZ5YI18zMxqLORLIZ2KdhfTbw+EhlJE0F9gCeSuuzgeuBUyLihw3lZxec08zM2qjORHIXME/SXEm7AMuA/lyZfrLOdIClwM0REZL2BL4BnB0R3x8uHBFPAM9IOjSN1joF+FqN78HMzArUlkhSn8cZwGrgQeCaiFgv6VxJx6ViK4HpkoaAjwHDQ4TPAPYHzpG0Lv28Ju37EHApMAT8EPhmXe/BzMyKKRv8NLn19fXFwMBAp8MwM5tQJA1GRF9ROd/ZbmZmlTiRmJlZJU4kZmZWiROJmZlV4kRiZmaVOJGYmVklTiRmZlaJE4mZmVXiRGJmZpU4kZiZWSVOJGZmVokTiZmZVeJEYmZmlTiRmJlZJU4kZmZWiROJmZlV4kRiZmaVOJGYTTCDm7Zxysq1DG7a1ulQzAAnErMJ54I1G7jt4a1csGZDp0MxA2BqpwMws9E5c9H8l72adZoTidkEc/B+e3Hl8kM6HYbZS9y0ZWZmlTiRmJlZJU4kZmZWiROJmZlV4kRiZmaVOJGYmVkltSYSSYslPSRpSNJZTfZPk7Qq7V8raU7aPl3SLZJ+JunzuWO+k865Lv28ps73YGZmrdWWSCRNAS4EjgEWAO+VtCBXbDmwLSL2Bz4LnJe2PwecA3x8hNOfFBEHpp8nxz96s/bxlCc20ZVKJJKuk/SHkkaTeBYCQxGxMSJeAK4GluTKLAGuSMvXAkdIUkQ8GxHfI0soZpOapzyxia5sYvgi8D7gYUn/U9IbShwzC3i0YX1z2ta0TERsB54Gppc49+WpWescSWpWQNIKSQOSBrZs2VLilGadceai+Rw2b4anPLEJq1QiiYg1EXEScBDwCHCjpB9IOk3SziMc1uwPfIyhTN5JEfEm4K3p5+QRYr4kIvoiom/mzJkFpzTrnOEpTw7eb69Oh2I2JqWbqiRNB/4Y+ABwD3ABWWK5cYRDNgP7NKzPBh4fqYykqcAewFOt4oiIx9LrM8CXyZrQzMysQ8r2kXwV+C7wCuDYiDguIlZFxEeAV41w2F3APElzJe0CLAP6c2X6gVPT8lLg5ogYsUYiaaqkGWl5Z+BdwP1l3oOZmdWj7Oy/l0bEDY0bJE2LiOcjoq/ZARGxXdIZwGpgCnBZRKyXdC4wEBH9wErgKklDZDWRZQ3nfwTYHdhF0vHAUcAmYHVKIlOANcCXyr9dMzMbb2pRAdhRSLo7Ig4q2tat+vr6YmBgoNNhmJlNKJIGR6osNGpZI5H0WrKRVbtJejM7Osd3J2vmMjOzHlfUtHU0WQf7bOD8hu3PAH9eU0xmZjaBtEwkEXEFcIWkEyLiujbFZGZmE0hR09b7I+KfgDmSPpbfHxHnNznMzMx6SFHT1ivT60hDfM3MrMcVNW1dnF4/1Z5wzMxsoilq2vpcq/0R8afjG46ZmU00RU1bg22JwsxGNLhpGxes2cCZi+Z7Pi7rSmVGbZlZBw1PMw9w5fJDOhyN2a8ratr6+4j4qKSv02RW3og4rrbIzAzgpenlPc28dauipq2r0uvf1R2I2WQx3k1Rw9PMm3WrlrP/RsRger0VuB3YRja54u1pm1nPyz8q1088tF5TavZfSX8IXAT8kGy+rbmSPhgR36wzOLOJIN+H4aYo6zVlp5H/X8DhETEEIOl1wDcAJxLrefnE4aYo6zVlE8mTw0kk2Qg8WUM8ZhOOE4f1uqJRW+9Oi+sl3QBcQzZ66z1kT0A0M7MeV1QjObZh+f8Cb0vLWwDfGWVmZoU3JJ7WrkDMeoXvVLfJpqhp679FxN9I+gea35DoubbMRsl3qttkU9S09WB69QPPzcaJhwfbZFPUtPX1tPjziPjnxn2S3lNbVGaTmEd52WTT8s72BmeX3GZmZj2mqI/kGOCdwKzcs0l2B7bXGZiZmU0MRTWSx8n6R54jezbJ8E8/cHS9oZlNDvm5uMwmm6JJG/8tPZNkf+ArwD3A3cC/RoSvCrMS2j2JoxOXtVvZKVKOBC7GkzaajVq7R2l5eLG1W9nO9vPJJm18e0S8DTgc+GzRQZIWS3pI0pCks5rsnyZpVdq/VtKctH26pFsk/UzS53PHHCzpvnTM5ySp5Hsw64jhUVrtuvnwzEXzOWzeDA8vtrYpm0hGPWmjpCnAhcAxwALgvZIW5IotB7ZFxP5kiem8tP054Bzg401O/UVgBTAv/Swu+R7MekK7E5dZ2USyXtINkv5Y0qnA14G7JL27YWLHvIXAUERsjIgXgKuBJbkyS4Dh58JfCxwhSRHxbER8jyyhvETS3sDuEXF7RARwJXB8yfdgZmY1KJtIdmXHpI1vJ5u08TfIJnV81wjHzAIebVjfnLY1LRMR24Gngekt4piVztPqnABIWiFpQNLAli1bWpzSzMyqKNXZPsbJG5v1XeTn6ypTZkzlI+IS4BKAvr6+Vuc0M7MKStVIJM2XdJOk+9P6b0v6i4LDNgP7NKzPJrsvpWkZSVOBPcieCd/qnLMLzmlmZm1UtmnrS2RTovwSICLuBZYVHHMXME/SXEm7pPL9uTL9wKlpeSlwc+r7aCoingCekXRoGq11CvC1ku/BzMxqUPY+kldExJ25kbYtp0iJiO2SzgBWA1OAyyJivaRzgYGI6AdWAldJGiKribyUnCQ9QjYVyy6SjgeOiogHgA8B/wjsRvbMeN/LYmbWQWUTyVZJryP1R0haCjxRdFBE3ADckNv2iYbl58ge29vs2DkjbB8ADigZt5mZ1axsIjmdrOP6DZIeA34EvL+2qMzMbMIoO2prI7BI0iuBnSLimXrDMjOziaLsqK2/krRnulHwGUl7Sfp03cGZmVn3Kztq65iI+OnwSpr59531hGRmZhNJ2UQyRdK04RVJuwHTWpQ3M7MeUbaz/Z+AmyRdTjZy60/YMUeWmZn1sLKd7X8j6V5gEdk0JX8ZEatrjczMzCaEUokkjdb6dkR8S9LrgddL2jkifllveGZm1u3K9pHcBuwqaRawBjiN7O5yMzPrcWUTiSLi58C7gX+IiD8ie1iVmZn1uNKJRNLvAicB30jbynbUm5nZJFY2kXyUbPbf69PEi/8RuKW+sMzMbKIoO2rrVuDWhvWNwJ/WFZSZmU0cZUdt3UKTJxFGxDvGPSIzM5tQyvZzfLxheVfgBAqeR2JmZr2hbNPWYG7T9yXd2rSwmZn1lLKz//5Gw88MSYuB19Ycm5mNg8FN2zhl5VoGN23rdCg2SZVt2hpkRx/JduARYHkdAZnZ+LpgzQZue3grAFcuP6TD0dhk1DKRSHoL8GhEzE3rp5L1jzwCPFB7dGZW2ZmL5r/s1Wy8FTVtXQy8ACDpMOCvyWb9fZrs0btm1uUO3m8vrlx+CAfvt1enQ7FJqqhpa0pEPJWWTwQuiYjrgOskras3NDMzmwiKaiRTJA0nmyOAmxv2eYoUMzMrTAZfAW6VtBX4BfBdAEn7kzVvmZlZj2tZI4mIzwD/hWzK+D+IiOGRWzsBH6k3NLPO8HBZs9EpvI8kIu6IiOsj4tmGbRsi4u56QzPrjOHhshes2dDpUDrOSdXKcD+HWY6Hy+7ge1CsDCcSs5zh4bLmpGrllH0eyZhIWizpIUlDks5qsn+apFVp/1pJcxr2nZ22PyTp6Ibtj0i6T9I6SQN1xm82GY2mucr3oFgZtdVIJE0BLgSOBDYDd0nqj4jGO+KXA9siYn9Jy4DzgBMlLQCWAW8E/gOwRtL8iHgxHXd4RGytK3azyczNVTbe6qyRLASGImJjRLwAXA0syZVZQnanPMC1wBGSlLZfHRHPR8SPgKF0PjOr6MxF8zls3gw3V9m4qTORzAIebVjfnLY1LRMR28nuTZlecGwA35Y0KGnFSL9c0gpJA5IGtmzZUumNWPfxaKKxc3OVjbc6E4mabMs/ZXGkMq2O/f2IOAg4Bjg9zQH264UjLomIvojomzlzZtmYbYLwEN3OcAK3ZupMJJuBfRrWZwOPj1QmTcWyB/BUq2MjYvj1SeB63OTVk9w80xlO4NZMnYnkLmCepLmSdiHrPO/PlekHTk3LS4Gb093z/cCyNKprLjAPuFPSKyW9GkDSK4GjgPtrfA/Wpdw8U59WtY58AncNxaDGUVsRsV3SGcBqYApwWUSsl3QuMBAR/cBK4CpJQ2Q1kWXp2PWSriF75sl24PSIeFHSbwLXZ/3xTAW+HBHfqus9mPWiVqO68vfYeASYQc03JEbEDcANuW2faFh+DnjPCMd+BvhMbttG4HfGP1Kz8gY3beOCNRs4c9H8SVkjGs1NiL5h0aDmGxLNJqPJ1k+Qb54aTbNhUVk3ffUGT5FiNkqT7Vt4nc1TbvrqDU4kZqM02ebiqjMxFp17sjcT9grteMTI5NXX1xcDA56Wy6zbnLJyLbc9vJXD5s2YVMl5spA0GBF9ReXcR2JmbZPvMxntcGL3uXQnN22ZWdvk+0xGO5zYfS7dyYnEzNqmqM+k6n7rDPeRmFnXcmd8Z7mPxHqa29I7Y7z/3SfbPTuTlZu2bFJyW3pnjPe/u5uyJgYnEpuU/AeoM8b7373qPTtuGmsP95GYTXD+Yzky36dSTdk+EtdIzCY4N+ONzDXT9nBnu00KdXaud3vHvR/ytUOVCSht7JxIbFKoc3RPt48c6uU/lvnEMd7/V93+JaJbOJHYpFDnt3J/4+9e+cRR9H812sTQKjE5yezgPhKbFOqckXeyzfY7meT7QIr+r/L9SUUDFfLnbyxf1DfVS4MgnEjMbMIabZLPJ4aiZNBqLrCijvxeGgThRGI9qZe+LdoO+cTQqsZRVEMpSmK9NGLMfSTWk7q9A93aIz9Qoehz0WpgQ9URY/njq06p384+HCcS60nuQLdmqnwuipJQ0R/2/PFF56u6fzy5act6Ur5Zwk1dBtUGVlTtM8kfX9TsNtr9dfIUKWZ4Kg2rX/4P/Wi/vBR9Ruv4DHsaeZtQRts+XPX8eW7qsrqNtj8mr+gz2snPsJu2rCvkq/3jPXRytMM8zeo22qanos9oJz/DrpFYLYpqGPn1/Lep/HrVGoprHNZtJtPUNrUmEkmLJT0kaUjSWU32T5O0Ku1fK2lOw76z0/aHJB1d9pw2PqoORSwagZJfz19UVZsB8ibTRWvWbWpLJJKmABcCxwALgPdKWpArthzYFhH7A58FzkvHLgCWAW8EFgNfkDSl5Dl7xmi/9bdaH20iKCpfVMMYbQ3BNQqz7lVnjWQhMBQRGyPiBeBqYEmuzBLgirR8LXCEJKXtV0fE8xHxI2Aona/MOcfNeP6hrmN9tN/6W62PNhEUlS+qYYy2huAahVn3qrOzfRbwaMP6ZiDfE/RSmYjYLulpYHrafkfu2FlpueicAEhaAawA2Hfffcf0Boo6gDu9XjTufLSvjcv5jruiqSXcWW3Wu+pMJGqyLX/TykhlRtrerAbV9EaYiLgEuASy+0hGDnNk4/2Herxfi/7Yj3Z9NInAicPMhtV2Q6Kk3wU+GRFHp/WzASLirxvKrE5lbpc0FfgJMBM4q7HscLl0WMtzNuMbEs3MRq8bbki8C5gnaa6kXcg6z/tzZfqBU9PyUuDmyDJbP7AsjeqaC8wD7ix5TjMza6PamrZSn8cZwGpgCnBZRKyXdC4wEBH9wErgKklDwFNkiYFU7hrgAWA7cHpEvAjQ7Jx1vQczMyvmubbMzKypbmjaMjOzHuBEYmZmlTiRmJlZJU4kZmZWiROJmZlV4kRiZmaVOJGYmVklTiRmZlaJE4mZmVXiRGJmZpU4kZiZWSVOJGZmVokTiZmZVeJEYmZmlTiRmJlZJU4kZmZWiROJmZlV4kRiZmaVOJGYmVklTiRmZlaJE4mZmVXiRGJmZpU4kZiZWSVOJGZmVokTiZmZVeJEYmZmlTiRmJlZJYqITsdQO0lbgE1jPHwGsHUcwxlv3RxfN8cG3R2fYxu7bo6vm2ODX49vv4iYWXRQTySSKiQNRERfp+MYSTfH182xQXfH59jGrpvj6+bYYOzxuWnLzMwqcSIxM7NKnEiKXdLpAAp0c3zdHBt0d3yObey6Ob5ujg3GGJ/7SMzMrBLXSMzMrBInEjMzq8SJBJB0maQnJd3foszbJa2TtF7Srd0Un6Q9JH1d0r+l+E5rY2z7SLpF0oPpd5/ZpIwkfU7SkKR7JR3URbGdlGK6V9IPJP1OO2IrG19D2bdIelHS0m6KrVPXRcn/245cF5J2lXRnw+/9VJMy0yStStfEWklzuii2j0l6IF0TN0nar/DEEdHzP8BhwEHA/SPs3xN4ANg3rb+my+L7c+C8tDwTeArYpU2x7Q0clJZfDWwAFuTKvBP4JiDgUGBtF8X2e8BeafmYdsVWNr60bwpwM3ADsLRbYuvkdVEyvo5cF+lz/qq0vDOwFjg0V+bDwEVpeRmwqk3/bmViOxx4RVr+UJnYXCMBIuI2sg/ZSN4HfDUifpzKP9mWwJIS8QXwakkCXpXKbm9TbE9ExN1p+RngQWBWrtgS4MrI3AHsKWnvbogtIn4QEdvS6h3A7LrjGk18yUeA64C2fe5Kxtax66JkfB25LtLn/Gdpdef0kx/VtAS4Ii1fCxyR4ux4bBFxS0T8PK2WuiacSMqZD+wl6TuSBiWd0umAcj4P/BbwOHAfcGZE/KrdQaTq+ZvJvuU0mgU82rC+meZ/MGvTIrZGy8lqTm03UnySZgF/BFzU/qheimEOzf/tuuK6aBFfx64LSVMkrSNL/jdGxIjXRERsB54GpndJbI1KXRNTxyu4SW4qcDBwBLAbcLukOyJiQ2fDesnRwDrgHcDrgBslfTci/l+7ApD0KrJvzR9t8nubfdNq27jzgtiGyxxOdtH8QbviavjdreL7e+DPIuLFNnxhHW1sHb8uCuLr2HURES8CB0raE7he0gER0djH2bFrokRsAEh6P9AHvK3onK6RlLMZ+FZEPBsRW4HbgLZ1ypZwGlkTQ0TEEPAj4A3t+uWSdia7mP93RHy1SZHNwD4N67PJviV2Q2xI+m3gUmBJRPx7O+IaRXx9wNWSHgGWAl+QdHyXxNbR66JEfB29LgAi4qfAd4DFuV0vXROSpgJ70Lr5up2xIWkR8N+B4yLi+aJzOZGU8zXgrZKmSnoFcAhZm2y3+DHZt0Ik/SbwemBjO35xatddCTwYEeePUKwfOCWN3joUeDoinuiG2CTtC3wVOLndNcwy8UXE3IiYExFzyNrSPxwR/9INsdHB66JkfB25LiTNTN/2kbQbsAj4P7li/cCpaXkpcHOk3u1OxybpzcDFZEmkVL+Xm7YASV8B3g7MkLQZ+B9knVBExEUR8aCkbwH3Ar8CLm1WFexUfMBfAv8o6T6yKvOfpW+I7fD7wMnAfandFbLRMvs2xHcD2citIeDnZN8UuyW2T5C1TX8hNR1tj/bNzlomvk4pjK3D10WZf7tOXRd7A1dImkL2Zf2aiPhXSecCAxHRT5YEr5I0RFYTWdaGuMrG9rdkgxP+OV0TP46I41qd1FOkmJlZJW7aMjOzSpxIzMysEicSMzOrxInEzMwqcSIxM7NKPPzXepqk6cBNafW1wIvAlrS+MCJe6EhgLUj6E+CGiPhJp2MxAw//NXuJpE8CP4uIv+uCWKakqSya7fsecEZErGu2f4RjpqY5nczGnZu2zEYg6dT07IZ1kr4gaad0F/dPJf2tpLslrZZ0iKRbJW2U9M507AckXZ/2PyTpL0qe99OS7gQWSvqUpLsk3S/pojQzwInAgcCqdPwukjY33K18qKQ1afnTki6WdCNwefod56fffa+kD7T/X9UmIycSsyYkHUA26+7vRcSBZM3Aw3cf7wF8OyIOAl4APkk2Fcd7gHMbTrMwHXMQ8D5JB5Y4790RsTAibgcuiIi3AG9K+xZHxCqyiQhPjIgDSzS9vRk4NiJOBlYAT0bEQuAtwOlpihizStxHYtbcIrI/tgNpmojd2DEV/i8i4sa0fB/Z3GHb01QccxrOsXr4WSeS/oVsZuGpLc77AnB9w/FHSPqvwK7ADGCQ0U9z/7WIeC4tHwX8lqTGxDWPbE4qszFzIjFrTsBlEXHOyzZmM7U21gJ+BTzfsNx4TeU7IKPgvL8YnrgvTYL4ebKnAD4m6dNkCaWZ7exoXciXeTb3nj4cETdhNo7ctGXW3BrgP0maAdnorjE0Ax0lac+UFJYA3x/FeXcjS0xbJb0aOKFh3zNkj5cd9gjZc0HIlctbDXw4JS0kvT7NAGtWiWskZk1ExH2SPgWskbQT8EvgPzO656h8D/gy2UOVrhoeZVXmvBHx75KuAO4HNvHyp/9dDlwq6Rdk/TCfBL4k6SfAnS3iuZhsdtx1qVntSbIEZ1aJh/+a1SCNiDogIj7a6VjM6uamLTMzq8Q1EjMzq8Q1EjMzq8SJxMzMKnEiMTOzSpxIzMysEicSMzOr5P8Dw3GkRrskI8oAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Plot calculated stuff\n",
    "\n",
    "plt.scatter(temps,aveEn,s=2)\n",
    "plt.xlabel('Temperature')\n",
    "plt.ylabel('Average Energy')\n",
    "plt.show()\n",
    "\n",
    "plt.scatter(temps,heatCap,s=2)\n",
    "plt.xlabel('Temperature')\n",
    "plt.ylabel('Heat Capacity')\n",
    "plt.show()\n",
    "\n",
    "plt.scatter(temps,mag,s=2)\n",
    "plt.xlabel('Temperature')\n",
    "plt.ylabel('Magnetization')\n",
    "plt.show()\n",
    "\n",
    "plt.scatter(temps,sus,s=2)\n",
    "plt.xlabel('Temperature')\n",
    "plt.ylabel('Susceptibility')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Faster Computation\n",
    "The first completed version of the code (with numba) the hyperparameters were\n",
    "\n",
    "```width = 32\n",
    "cycles = 1000 * width**2\n",
    "temps = np.arange(1.60,3.15,0.02)```\n",
    "\n",
    "and took ~ 155 seconds\n",
    "\n",
    "Let's see what optimizations can be made! Before we look at each function, note that `x**2` takes noticably longer than `x*x` and numba offers a fastmasth option that will be incorperated in each function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Initial parameters\n",
    "\n",
    "width = 2**5  # width of sqaure spin lattice; MUST BE POWER OF 2; 64 points may take hours\n",
    "width_log2 = int(np.log2(width))\n",
    "s = np.random.choice([1,-1],size=(width,width))  # initial spin sites randomly (+1 or -1)\n",
    "cycles = 1000 * width**2\n",
    "\n",
    "# Initlaize temperature range (crosses critical temperature)\n",
    "temps = np.arange(1.60,3.15,0.02)\n",
    "\n",
    "# Inititalize magnetization, average energy, heat capacity, and susceptibility\n",
    "mag = np.zeros(len(temps))\n",
    "aveEn = np.zeros(len(temps))\n",
    "heatCap = np.zeros(len(temps))\n",
    "sus = np.zeros(len(temps))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "288 ns ± 17.9 ns per loop (mean ± std. dev. of 7 runs, 1000000 loops each)\n"
     ]
    }
   ],
   "source": [
    "# old way\n",
    "%timeit width**2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "65 ns ± 0.531 ns per loop (mean ± std. dev. of 7 runs, 10000000 loops each)\n"
     ]
    }
   ],
   "source": [
    "# new way\n",
    "%timeit width*width"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Delta u\n",
    "Is there a faster way to implement periodic boundary conditions and add up neighboring spins? An attempt using np.roll and matrix addition/mulitplication was made, but it was x20 slower."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculate interaction energy between spins. Assume periodic boundaries. i and j and be any value\n",
    "# Interaction energy will be the difference in energy due to flipping spin i,j \n",
    "# (i.e. 2*spin_value*neighboring_spins)\n",
    "@njit(fastmath=True)\n",
    "def deltaU_fast(s,i,j):\n",
    "    #top\n",
    "    if i == 0:\n",
    "        t = s[width-1,j]  # periodic boundary\n",
    "    else:\n",
    "        t = s[i-1,j]\n",
    "    #bottom\n",
    "    if i == width-1:\n",
    "        b = s[0,j]  # periodic boundary\n",
    "    else:\n",
    "        b = s[i+1,j]\n",
    "    #left\n",
    "    if j == 0:\n",
    "        l = s[i,width-1]  # periodic boundary\n",
    "    else:\n",
    "        l = s[i,j-1]\n",
    "    #right\n",
    "    if j == width-1:\n",
    "        r = s[i,0]  # periodic boundary\n",
    "    else:\n",
    "        r = s[i,j+1]\n",
    "    return s[i,j] * (t+b+r+l)  # difference in energy if i,j is flipped"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "381 ns ± 2.45 ns per loop (mean ± std. dev. of 7 runs, 1000000 loops each)\n"
     ]
    }
   ],
   "source": [
    "%timeit deltaU(s,15,15)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "383 ns ± 2.83 ns per loop (mean ± std. dev. of 7 runs, 1000000 loops each)\n"
     ]
    }
   ],
   "source": [
    "%timeit deltaU_fast(s,15,15)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "%timeit reveals that `calcE()` takes up the most time. Let's start there for possible improvements.\n",
    "## Calculate Energy\n",
    "Using just matrix manipulations to compute the total energy took x10 longer than the the double for loop with jit. Also, having a function for the bulk spins and a function for the edge spins (in hopes of reducing if logic) actually ran slower as well. I noticed the output of `deltaU()` is multipled by 2, but then immediately divided by two in `calcE()`. Resolving this throughout the functions helped lower runtime. The two other speed ups are from `fastmath=True` and replacing `x**2` with `x*x`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculate the energy of the lattice.\n",
    "# The energy is the sum of interactions between spins divided by the total number of spins\n",
    "@njit(fastmath=True)\n",
    "def calcE_fast(s):\n",
    "    E = 0\n",
    "    for i in range(width):\n",
    "        for j in range(width):\n",
    "            E += -deltaU_fast(s,i,j)\n",
    "    return E / (width*width)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2.87 µs ± 94.3 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)\n"
     ]
    }
   ],
   "source": [
    "# old way\n",
    "%timeit calcE(s)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.79 µs ± 21.7 ns per loop (mean ± std. dev. of 7 runs, 1000000 loops each)\n"
     ]
    }
   ],
   "source": [
    "# new way\n",
    "%timeit calcE_fast(s)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Calculate Magnetization\n",
    "Already fast and not a main time contributor"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculate magnetization. \n",
    "# Magnetization is the sum of all spins divided by the total number of spins\n",
    "@njit(fastmath=True)\n",
    "def calcM_fast(s):\n",
    "    m = np.abs(s.sum()) / (width*width)\n",
    "    return m"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "391 ns ± 4.92 ns per loop (mean ± std. dev. of 7 runs, 1000000 loops each)\n"
     ]
    }
   ],
   "source": [
    "%timeit calcM(s)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "392 ns ± 2.15 ns per loop (mean ± std. dev. of 7 runs, 1000000 loops each)\n"
     ]
    }
   ],
   "source": [
    "%timeit calcM_fast(s)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Monte Carlo Cycle\n",
    "Overall, `montecc()` and `montecc_fast()` take the same time for one cycle, but `montecc_fast()` is noticably faster over many cycles.\n",
    "\n",
    "### Lattice Size\n",
    "If we restrict the lattice size to be a power of 2, the random number generating time deceases by a factor of 10 for that part of the computation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Monte carlo cycle\n",
    "@njit(fastmath=True)\n",
    "def montecc_fast(s,Temp,cycles):\n",
    "    for m in range(cycles):\n",
    "        i = random.getrandbits(width_log2)  # choose random row\n",
    "        j = random.getrandbits(width_log2)  # choose random column\n",
    "        ediff = 2 * deltaU_fast(s,i,j)\n",
    "        if ediff <= 0:\n",
    "            s[i,j] *= -1  # flip spin\n",
    "        elif random.random() < np.exp(-ediff/Temp):\n",
    "            s[i,j] *= -1\n",
    "    return s"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "986 ns ± 3.72 ns per loop (mean ± std. dev. of 7 runs, 1000000 loops each)\n"
     ]
    }
   ],
   "source": [
    "# old way: random number generation\n",
    "%timeit random.randrange(width)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "151 ns ± 0.816 ns per loop (mean ± std. dev. of 7 runs, 10000000 loops each)\n"
     ]
    }
   ],
   "source": [
    "# new way: random number generation\n",
    "%timeit random.getrandbits(width_log2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.47 µs ± 29.5 ns per loop (mean ± std. dev. of 7 runs, 1000000 loops each)\n"
     ]
    }
   ],
   "source": [
    "#generate an s for testing\n",
    "%timeit montecc(s,1.60,20)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The slowest run took 10.30 times longer than the fastest. This could mean that an intermediate result is being cached.\n",
      "3.19 µs ± 3.76 µs per loop (mean ± std. dev. of 7 runs, 1 loop each)\n"
     ]
    }
   ],
   "source": [
    "%timeit montecc_fast(s,1.60,20)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Main computation loop\n",
    "Most of the time comes from calcE function. All previously mentioned speed ups are included here."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compute physical quantities\n",
    "@njit(fastmath=True)\n",
    "def computePQ_fast(s,T,cycles):\n",
    "    En = 0\n",
    "    En_sq = 0\n",
    "    Mg = 0\n",
    "    Mg_sq = 0\n",
    "    for p in range(cycles):\n",
    "        s = montecc_fast(s,T,1)\n",
    "        E = calcE_fast(s)\n",
    "        M = calcM_fast(s)\n",
    "        En += E\n",
    "        Mg += M\n",
    "        En_sq += E*E\n",
    "        Mg_sq += M*M\n",
    "    aveEn = En / cycles\n",
    "    mag = Mg / cycles\n",
    "    sus = (Mg_sq / cycles - (mag * mag))/T\n",
    "    heatCap = (En_sq / cycles - (aveEn * aveEn))/(T*T)\n",
    "    return aveEn, mag, sus, heatCap"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "94 µs ± 1.46 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)\n"
     ]
    }
   ],
   "source": [
    "%timeit computePQ(s,1.60,50)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "79 µs ± 4.62 µs per loop (mean ± std. dev. of 7 runs, 1 loop each)\n"
     ]
    }
   ],
   "source": [
    "%timeit computePQ_fast(s,1.60,50)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Faster Computation\n",
    "Let's use these faster functions and see how well they do! The previous computation was 155s. With the \"fast\" functions and same hyperparameters, it's 126s! An 18.7% faster computation!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Fast helper functions all in one place\n",
    "\n",
    "# Calculate the energy of the lattice.\n",
    "# The energy is the sum of interactions between spins divided by the total number of spins\n",
    "@njit(fastmath=True)\n",
    "def calcE_fast(s):\n",
    "    E = 0\n",
    "    for i in range(width):\n",
    "        for j in range(width):\n",
    "            E += -deltaU_fast(s,i,j)\n",
    "    return E / (width*width)\n",
    "\n",
    "\n",
    "# Calculate magnetization. \n",
    "# Magnetization is the sum of all spins divided by the total number of spins\n",
    "@njit(fastmath=True)\n",
    "def calcM_fast(s):\n",
    "    m = np.abs(s.sum()) / (width*width)\n",
    "    return m\n",
    "\n",
    "# Calculate interaction energy between spins. Assume periodic boundaries. i and j and be any value\n",
    "# Interaction energy will be the difference in energy due to flipping spin i,j \n",
    "# (i.e. 2*spin_value*neighboring_spins)\n",
    "@njit(fastmath=True)\n",
    "def deltaU_fast(s,i,j):\n",
    "    #top\n",
    "    if i == 0:\n",
    "        t = s[width-1,j]  # periodic boundary\n",
    "    else:\n",
    "        t = s[i-1,j]\n",
    "    #bottom\n",
    "    if i == width-1:\n",
    "        b = s[0,j]  # periodic boundary\n",
    "    else:\n",
    "        b = s[i+1,j]\n",
    "    #left\n",
    "    if j == 0:\n",
    "        l = s[i,width-1]  # periodic boundary\n",
    "    else:\n",
    "        l = s[i,j-1]\n",
    "    #right\n",
    "    if j == width-1:\n",
    "        r = s[i,0]  # periodic boundary\n",
    "    else:\n",
    "        r = s[i,j+1]\n",
    "    return s[i,j] * (t+b+r+l)  # difference in energy if i,j is flipped\n",
    "\n",
    "# Monte carlo cycle\n",
    "@njit(fastmath=True)\n",
    "def montecc_fast(s,Temp,cycles):\n",
    "    for m in range(cycles):\n",
    "        i = random.getrandbits(width_log2)  # choose random row\n",
    "        j = random.getrandbits(width_log2)  # choose random column\n",
    "        ediff = 2 * deltaU_fast(s,i,j)\n",
    "        if ediff <= 0:\n",
    "            s[i,j] *= -1  # flip spin\n",
    "        elif random.random() < np.exp(-ediff/Temp):\n",
    "            s[i,j] *= -1\n",
    "    return s\n",
    "\n",
    "\n",
    "# Compute physical quantities\n",
    "@njit(fastmath=True)\n",
    "def computePQ_fast(s,T,cycles):\n",
    "    En = 0\n",
    "    En_sq = 0\n",
    "    Mg = 0\n",
    "    Mg_sq = 0\n",
    "    for p in range(cycles):\n",
    "        s = montecc_fast(s,T,1)\n",
    "        E = calcE_fast(s)\n",
    "        M = calcM_fast(s)\n",
    "        En += E\n",
    "        Mg += M\n",
    "        En_sq += E*E\n",
    "        Mg_sq += M*M\n",
    "    aveEn = En / cycles\n",
    "    mag = Mg / cycles\n",
    "    sus = (Mg_sq / cycles - (mag * mag))/T\n",
    "    heatCap = (En_sq / cycles - (aveEn * aveEn))/(T*T)\n",
    "    return aveEn, mag, sus, heatCap"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Elapsed = 126.50802516937256\n"
     ]
    }
   ],
   "source": [
    "start = time.time()\n",
    "# Compute physical quanities (magnetization, etc.)\n",
    "for ind, T in enumerate(temps):\n",
    "    # Cycle spins\n",
    "    s = montecc_fast(s,T,cycles)\n",
    "    # Compute physical quanitites with 1000 cycles per spin\n",
    "    aveEn[ind], mag[ind], sus[ind], heatCap[ind] = computePQ_fast(s,T,cycles)\n",
    "end = time.time()\n",
    "print(\"Elapsed = %s\" % (end - start))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY4AAAEGCAYAAABy53LJAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAZ60lEQVR4nO3dfbRcdX3v8feHBwGRJaFBeSa2Jj4sxEgOD1e8USBF5Cr4QIViBRUul4qa9q7VRRWftV1ecfWaLi8FBOvDoleqJUIlBRIfoFZ5OIcbAgpNKBCJ0UUiUSkqGPjcP2YfHMY5Z/acc2b2npnPa62zsvfMPnu+xmy+8/t+f/u3ZZuIiIiydqg6gIiIGCxJHBER0ZUkjoiI6EoSR0REdCWJIyIiurJT1QH0wvz5871gwYKqw4iIGBgTExNbbe9d5tihTBwLFixgfHy86jAiIgaGpI1lj02pKiIiupLEERERXUniiIiIrlSSOCT9kaTvS3pS0tg0xz0g6U5JayWlaRERUQNVNcfvAt4IXFLi2GNsb+1xPBERUVIlicP23QCSqvj4iIiYhbr3OAzcIGlC0jnTHSjpHEnjksa3bNnSp/AiIkZPz0YcktYA+7R56wLbV5c8zdG2N0t6DrBa0j22b2p3oO1LgUsBxsbGslZ8RESP9Cxx2F42B+fYXPz5kKSVwBFA28QRETHKJjZuY8Wa9SxftoglB8/r6WfVtlQlaXdJe0xuA8fTaKpHRIyciY3bOOPyW5jYuK3t/oo167lpw1ZWrFnf81iqmo77BkmbgP8CXCvp+uL1/SStKg57LvAdSXcAtwLX2r6uingjIqrWmhha95cvW8TShfNZvmxRz2OpalbVSmBlm9c3AycW2/cBL+1zaBERtdBaeppMCFP9ueTgeXzxrCP7EpuG8ZnjY2NjziKHETHIzrj8Fm7asJWlC+f3JSFImrA95Q3ZzWrb44iIGGadehb9LD11ayiXVY+IqLvJHgXAF8868nf2+1l66lYSR0REBTr1LOosPY6IiB5obW738z6LmUiPIyKiYp2mzw6ylKoiIuZAt9NnB1lKVRERM9CaKPo9fXaupVQVETHHOi3xUefps3MtpaqIiBJap8tWeed21ZI4IiLa6NSzGKVE0SqlqogIOpeiJhNFHafS9ltGHBERdC5FxW8lcUREUO1qs4MmpaqIGEmtpamUospL4oiIkdGcLIbpTu5+S6kqIoZW68yo5j5Gehgzl8QREUNruoZ3ehgzl8QREUMrDe/eSI8jIoZGGt79kcQREQOr00170RspVUXEwMpNe9VI4oiIgZH1o+ohpaqIqK2sH1VPGXFERG2lFFVPSRwRUVspRdVTSlURURuZTjsYkjgiojYynXYwpFQVEbWRHsZgyIgjIvqmtRSV0tRgSuKIiL5pLUWlNDWYUqqKiL5pLUWlNDWYZLvqGObc2NiYx8fHqw4jImJgSJqwPVbm2EpKVZIulHSPpHWSVkrac4rjTpD075LulfSX/Y4zImantYcRw6GqHsdq4BDbhwLrgfe2HiBpR+D/AK8BXgz8saQX9zXKiOhKVqsdDZUkDts32N5e7N4MHNDmsCOAe23fZ/tx4MvAyf2KMSI665Qoli9bxNKF89PDGDJ1aI6/A7iyzev7Aw827W8CplxrQNI5wDkABx100FzGFxFT6LSWVJYIGU49SxyS1gD7tHnrAttXF8dcAGwHrmh3ijavTdnJt30pcCk0muNdBxwRHWVZ84AeJg7by6Z7X9KZwGuB49x+atcm4MCm/QOAzXMXYUR0q3WEkUQxmiopVUk6ATgfeKXtX05x2G3AQknPA34EnAac3qcQI6KN3HcRUF2P4zPALsBqSQA32z5X0n7AZbZPtL1d0ruA64Edgc/Z/n5F8UYEKUVFQyWJw/bzp3h9M3Bi0/4qYFW/4oqIp2vtaURA1qqKiGnkPoxopw7TcSOiptLTiHaSOCJiSulpRDspVUXEU7K2VJSRxBExwrK2VMxESlURI6zTkiER7SRxRIywLBkSM5FSVcSQm65vkWd8x0wkcUQMuea+RZrfMRdSqooYcs3lqNaeRsRMJHFEDLnmvkWa3zEXkjgiRkia3zEX0uOIiIiuJHFEDJk0wKPXkjgihkzu/o5e65g4JL1LUiZ5R9RU6whj+bJFLF04Pw3w6JkyzfF9gNsk3Q58Drh+imeER0QF8hzw6LeOIw7b7wcWApcDbwM2SPprSX/Q49giooSMMKLfSvU4ihHGT4qf7cA84KuSPtnD2CKijdbSVJYNiX4r0+N4j6QJ4JPAvwEvsf2nwBLgTT2OLyJapPkdVSvT45gPvNH2xuYXbT8p6bW9CSsiJk1s3MaKNetZvmwRSw6el7u/o3JlEsenASTt1fTaI7Z/Y/vu3oQVEZPS/I66KdPjuB3YAqwHNhTb90u6XdKSXgYXMYoyvTbqrsyI4zpgpe3rASQdD5wA/CNwEZCvPhFzKCOMqLsyI46xyaQBYPsGYKntm4FdehZZxIjKCCPqrsyI42FJ5wNfLvZPBbZJ2hF4smeRRYyojDCi7sqMOE4HDgC+VvwcWLy2I/Dm3oUWMRqyKGEMmmlHHMWo4nzb757ikHvnPqSI0ZKn8sWgmXbEYfsJGjf6RcQcyaypGHRlehz/T9I1wFeARydftH1Vz6KKGGKZNRWDrkzi2Av4KXBs02sGkjgiZiB3fseg0zCukD42Nubx8fGqw4iIGBiSJmyPlTm2zCKHiyR9Q9Jdxf6hkt4/2yAjImIwlZmO+1ngvcBvAGyvA07rZVAREVFfZXocz7R9q6Tm17bP5kMlXQi8Dngc+A/g7bZ/1ua4B4BHgCeA7WWHURER0TtlRhxbi6f9GUDSKcCPZ/m5q4FDbB9KY/HE905z7DG2FydpRETUQ5kRx3nApcALJf0IuB/4k9l8aLHe1aSbgVNmc76IiOifjonD9n3AMkm7AzvYfmSOY3gHcOVUHw/cIMnAJbYvneokks4BzgE46KCD5jjEiIiY1DFxSNqFxiNiFwA7TfY6bH+0w++tAfZp89YFtq8ujrmARr/kiilOc7TtzZKeA6yWdI/tm9odWCSVS6ExHbfT/66IiJiZMqWqq4GfAxPAY2VPbHvZdO9LOhN4LXCcp7iZxPbm4s+HJK0EjgDaJo6IiOiPMonjANsnzOWHSjoBOB94pe1fTnHMU6WxYvt4YNpRTkRE9F6ZWVXflfSSOf7czwB70Cg/rZV0MYCk/SStKo55LvAdSXcAtwLX2r5ujuOIiIgulRlxvAJ4m6T7aZSqBLiYSjsjtp8/xeubgROL7fuAl870MyLqYmLjNlasWc/yZYtYcvC8qsOJmLUyieM1PY8iYojleRsxbKZMHJKOtf1N2xslPc/2/U3vvRHY2JcIIwZcVsONYTNdj+NTTdv/1PJeFjmMmELrg5omn7eRMlUMi+kSh6bYbrcfEYXJ0tSKNeurDiWiJ6brcXiK7Xb7EVFIaSqG3XSJ4/eLR8aqaZti/3k9jyyipjrNksqjYGPYTZc4Tm7a/lTLe637EUOrNVFkllSMuikTh+0b+xlIRF21JorWUlTu04hRU+Y+joiR1pooWktRGYHEqEniiOigU88izfAYNaUTh6TdbT/ay2AiBlGa4TFqOi5yKOnlkn4A3F3sv1TSRT2PLCIiaqnM6rj/G3g18FMA23cAS3sZVERE1FeZxIHtB1teeqIHsURExAAo0+N4UNLLAUt6BvAeirJVRESMnjIjjnOB84D9gU3A4mI/IiJGUMcRh+2twFv6EEtERAyAjolD0t+2efnnwLjtq+c+pIiIqLMypapdaZSnNhQ/hwJ7AWdJ+nQPY4uIiBoq0xx/PnCs7e0Akv4OuAH4Q+DOHsYWERE1VGbEsT+we9P+7sB+tp8AHutJVBERUVtlRhyfBNZK+jaNZ3EsBf5a0u7Amh7GFhERNVRmVtXlklYBR9BIHO+zvbl4+y96GVxERNRPqTvHgV8DPwYeBp4vKUuORESMqDLTcc8GlgMHAGuBo4DvAcf2NrSIiKijMiOO5cDhwEbbxwAvA7b0NKqIiKitMonj17Z/DSBpF9v3AC/obVgREVFXZWZVbZK0J/A1YLWkbcDmDr8TERFDquOIw/YbbP/M9oeBDwCXA6/vdWARVZnYuI0zLr+FiY3bqg4lopamHXFI2gFYZ/sQANs39iWqiAqtWLOemzZsBcgjYSPamDZx2H5S0h2SDrL9w34FFVGl5csWPe3PiHi6Mj2OfYHvS7oVeHTyRdsn9SyqiAotOXheRhoR0yiTOD7S8ygiImJglGmO3wg8AOxcbN8G3D7bD5b0MUnrJK2VdIOk/aY47kxJG4qfM2f7uRERMTsdE4ek/w58FbikeGl/GlNzZ+tC24faXgx8Hfhgm8/eC/gQcCSNtbI+JGneHHx2RETMUJkbAM8DjgZ+AWB7A/Cc2X6w7V807e4OuM1hrwZW237Y9jZgNXDCbD87IiJmrkzieMz245M7knai/X/kuybpryQ9SOOZ5r8z4qAxunmwaX9T8Vq7c50jaVzS+JYtWRElyst9GxHdKZM4bpT0PmA3SX8IfAX45zInl7RG0l1tfk4GsH2B7QOBK4B3tTtFm9faJi3bl9oesz229957lwkvAvjtfRsr1qyvOpSIgVBmVtVfAmfReEzs/wBWAZeVObntZSXj+AfgWhr9jGabgFc17R8AfLvkOSNKyX0bEd2RPX3VSdIbgFW25/QxsZIWFv0SJL0beKXtU1qO2QuYAA4rXrodWGL74enOPTY25vHx8bkMNyJiqEmasD1W5tgypaqTgPWSviTpvxU9jrnwiaJstQ44nsby7Ugak3QZQJEgPkZjCvBtwEc7JY2ITtLTiJidjiMOAEk7A68BTgVeQWOm09k9jm3GMuKI6Zxx+S3ctGErSxfOzx3iEYVuRhylRg+2fyPpX2g0pncDTgZqmzgippOeRsTslLkB8ARJnwfuBU6h0Rjft8dxRcyZ1tLU5FpUSw7OvaQRM1Gmx/E2GneKL7J9pu1Vtrf3NqyI2WlOFpluGzG3OpaqbJ/WvC/paOB02+f1LKqIWWp+pkZKUxFzq1SPQ9Ji4HTgzcD9wFW9DCpitpqTRZZJj5hbUyYOSYuA04A/Bn4KXEljFtYxfYotYsaSLCJ6Z7oRxz3AvwKvs30vgKQ/70tUERFRW9M1x98E/AT4lqTPSjqO9mtHRUTECJkycdheaftU4IU01of6c+C5kv5O0vF9ii8iImqmzBMAH7V9he3X0lhkcC2NhQ8jImIElbmP4ynFA5UusX1srwKKmImsPxXRP10ljoi6yk1+Ef0zVyvdRlQqN/lF9E8SRwyF3LcR0T8pVcVASk8jojpJHDGQ0tOIqE5KVTGQ0tOIqE5GHDEQ8kyNiPpI4oiBkNJURH0kcUQttY4wli9bxNKF81OaiqiB9DiilpofxDRZksp024h6SOKIWkrzO6K+kjiiljLCiKiv9DgiIqIrSRxRC7kTPGJwJHFELWS6bcTgSI8jaiHN8IjBkRFHVCJ3gkcMriSOqERKUxGDK4kj+iJ3gkcMj/Q4oi9yJ3jE8EjiiL5I8ztieKRUFT2R5nfE8KokcUj6mKR1ktZKukHSflMc90RxzFpJ1/Q7zpi5NL8jhldVpaoLbX8AQNJ7gA8C57Y57le2F/c1spgTKU1FDK9KRhy2f9G0uzvgKuKIuZPSVMToqKzHIemvJD0IvIXGiKOdXSWNS7pZ0us7nO+c4tjxLVu2zHm8Mb2UpiJGh+zefNmXtAbYp81bF9i+uum49wK72v5Qm3PsZ3uzpN8HvgkcZ/s/On322NiYx8fHZxF9dGti4zZWrFnP8mWLMsqIGECSJmyPlTq2V4mjLEkHA9faPqTDcZ8Hvm77q53OmcTRe0kUEcOlm8RR1ayqhU27JwH3tDlmnqRdiu35wNHAD/oTYbRq7WGkNBUxuqqaVfUJSS8AngQ2UsyokjQGnGv7bOBFwCWSnqSR4D5hO4mjIq13fmfWVMToqrxU1QspVc1eaykqpamI4dZNqSpLjkRbWVsqIqaSJUcCyOq1EVFeRhwBZIQREeUlcQSQJUIioryUqkZUlgiJiJlK4hgRuQ8jIuZKSlUjIvdhRMRcSeIYUq33XbQmijS/I2KmkjiGVGZJRUSvJHEMkeZRRkpREdEraY4PsOka3pklFRG9ksQxQDrNjMrd3hHRDylVDZBOM6PSx4iIfkjiqLHMjIqIOkriqLHMjIqIOkqPo0ayQm1EDIIkjgp1anZnZlRE1FFKVRXKMiARMYiSOPooze6IGAYpVfVQSlERMYySOGahNTHkBr2IGAUpVc1Ca48iN+hFxChI4uhCpx5FEkVEjALZrjqGOTc2Nubx8fFZn6c1UZxx+S3ctGErSxfOT0KIiKEiacL2WJlj0+Nokh5FRERnKVU1SY8iIqKzJI4mSRQREZ0lcTRJooiI6Cw9joiI6EoSR0REdCWJIyIiupLEERERXUniiIiIriRxREREV5I4IiKiK0O5VpWkLcDGGf76fGDrHIYz1+ocX2KbuTrHV+fYoN7x1Tk2eHp8B9veu8wvDWXimA1J42UX+qpCneNLbDNX5/jqHBvUO746xwYzjy+lqoiI6EoSR0REdCWJ43ddWnUAHdQ5vsQ2c3WOr86xQb3jq3NsMMP40uOIiIiuZMQRERFdSeKIiIiujGTikPQ5SQ9JumuaY14laa2k70u6sU7xSXq2pH+WdEcR39v7GNuBkr4l6e7is5e3OUaS/lbSvZLWSTqsRrG9pYhpnaTvSnppP2IrG1/TsYdLekLSKXWKrarrouT/t5VcF5J2lXRr0+d+pM0xu0i6srgmbpG0oEax/U9JPyiuiW9IOrjjiW2P3A+wFDgMuGuK9/cEfgAcVOw/p2bxvQ/4X8X23sDDwDP6FNu+wGHF9h7AeuDFLcecCPwLIOAo4JYaxfZyYF6x/Zp+xVY2vuK9HYFvAquAU+oSW5XXRcn4Krkuin/nzyq2dwZuAY5qOeadwMXF9mnAlX36eysT2zHAM4vtPy0T20iOOGzfROMf1VROB66y/cPi+If6ElihRHwG9pAk4FnFsdv7FNuPbd9ebD8C3A3s33LYycAX3XAzsKekfesQm+3v2t5W7N4MHNDruLqJr/Bu4J+Avv27KxlbZddFyfgquS6Kf+f/WezuXPy0zjo6GfhCsf1V4Lgizspjs/0t278sdktdEyOZOEpYBMyT9G1JE5LOqDqgFp8BXgRsBu4Eltt+st9BFMPtl9H4FtNsf+DBpv1NtP8PZM9ME1uzs2iMjPpuqvgk7Q+8Abi4/1E9FcMC2v/d1eK6mCa+yq4LSTtKWksj2a+2PeU1YXs78HPg92oSW7NS10SeOd7eTsAS4DhgN+B7km62vb7asJ7yamAtcCzwB8BqSf9q+xf9CkDSs2h8K/6zNp/b7ptU3+Z9d4ht8phjaFwkr+hXXE2fPV18nwbOt/1EH76Qdhtb5ddFh/gquy5sPwEslrQnsFLSIbabe5SVXRMlYgNA0p8AY8ArO50zI472NgHX2X7U9lbgJqBvTdQS3k6jZGDb9wL3Ay/s14dL2pnGxXuF7avaHLIJOLBp/wAa3wLrEBuSDgUuA062/dN+xNVFfGPAlyU9AJwCXCTp9TWJrdLrokR8lV4XALZ/BnwbOKHlraeuCUk7Ac9m+nJ0P2ND0jLgAuAk2491OlcSR3tXA/9V0k6SngkcSaOmWhc/pPGtD0nPBV4A3NePDy7qspcDd9v+mykOuwY4o5hddRTwc9s/rkNskg4CrgLe2u8RZJn4bD/P9gLbC2jUwt9p+2t1iI0Kr4uS8VVyXUjau/g2j6TdgGXAPS2HXQOcWWyfAnzTRTe66tgkvQy4hEbSKNW3GslSlaT/C7wKmC9pE/AhGk0jbF9s+25J1wHrgCeBy9oN7aqKD/gY8HlJd9IYAp9ffAPsh6OBtwJ3FnVTaMxmOagpvlU0ZlbdC/ySxjfBusT2QRq15YuKUtB292/10jLxVaVjbBVfF2X+7qq6LvYFviBpRxpfxv/R9tclfRQYt30NjaT3JUn30hhpnNaHuMrGdiGNyQRfKa6JH9o+abqTZsmRiIjoSkpVERHRlSSOiIjoShJHRER0JYkjIiK6ksQRERFdGcnpuDGaJP0e8I1idx/gCWBLsX+E7ccrCWwakt4BrLL9k6pjiZiU6bgxkiR9GPhP25+qQSw7FstCtHvvO8C7bK9t9/4Uv7NTsR5SRE+kVBUBSDqzeG7BWkkXSdqhuEP6Z5IulHS7pOslHSnpRkn3STqx+N2zJa0s3v93Se8ved6PS7oVOELSRyTdJukuSRcXd92fCiwGrix+/xmSNjXdCXyUpDXF9sclXSJpNfD3xWf8TfHZ6ySd3f+/1RhWSRwx8iQdQmNF2pfbXkyjhDt5Z++zgRtsHwY8DnyYxrIWfwR8tOk0RxS/cxhwuqTFJc57u+0jbH8PWGH7cOAlxXsn2L6SxqJ9p9peXKKU9jLgdbbfCpwDPGT7COBw4LxiuZWIWUuPI6Kxfs/hwHix5MJu/HZZ+F/ZXl1s30lj3a3txbIWC5rOcf3kcz4kfY3Gqrs7TXPex4GVTb9/nKS/AHYF5gMTdL/k+9W2f11sHw+8SFJzolpIYz2niFlJ4ohorGv0OdsfeNqLjVVMm7/lPwk81rTdfP20Ngvd4by/mlzkrlgw8DM0nnD3I0kfp5FA2tnObysFrcc82vK/6Z22v0HEHEupKgLWAG+WNB8as69mUNY5XtKeRRI4Gfi3Ls67G41EtFXSHsCbmt57hMajUic9QOOZGLQc1+p64J1FkkLSC4rVUSNmLSOOGHm275T0EWCNpB2A3wDn0t0zRL4D/AONBwh9aXIWVJnz2v6ppC8AdwEbefqT7f4euEzSr2j0UT4MfFbST4Bbp4nnEhorx64tymQP0UhoEbOW6bgRs1TMWDrE9p9VHUtEP6RUFRERXcmIIyIiupIRR0REdCWJIyIiupLEERERXUniiIiIriRxREREV/4/CZY7+tecUuMAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZsAAAEGCAYAAACzYDhlAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAei0lEQVR4nO3de5RdZZ3m8e9jIrcWAaFs6QRJHANtdGyBMuClaRWU4KjRbhyjo0TEYbVAg91rpoW2VUBdI2OPCqOIDBeB0Q50RCl7qZGLitoYqCAKAYGSi5TQQ2EQ8QIYfOaPvQtODqfO2XXZ55yqej5r1cq+vHvv34Ha9Tv73e9FtomIiKjTU3odQEREzH1JNhERUbskm4iIqF2STURE1C7JJiIiarew1wH0g912281LlizpdRgREbPKxo0b77c9UKVskg2wZMkShoeHex1GRMSsIumuqmVTjRYREbVLsomIiNol2URERO2SbCIionZJNhERUbskm4iIqF2STURE1C7JJiIiapdkExERtUuyiYiI2iXZRERE7ZJsIiKidkk2ERFRuySbiIioXZJNRETULskmIiJq19NkI2mlpFskjUg6ocX+bSVdVO7fIGlJuX1XSd+S9GtJn246Zj9JN5THnC5J3fk0ERExkZ4lG0kLgM8AhwLLgbdKWt5U7EjgAdvPBT4JnFpufxj4APDfWpz6s8BRwLLyZ+XMRx8REZPRyyebFcCI7dttPwqsBVY1lVkFnF8urwMOkiTbv7H9PYqk8zhJuwNPt321bQMXAG+s9VNERERHvUw2i4C7G9ZHy20ty9jeAjwI7NrhnKMdzgmApKMkDUsaHhsbm2ToERExGb1MNq3epXgKZaZU3vZZtgdtDw4MDLQ5ZURETFcvk80osEfD+mLgnonKSFoI7ARs7nDOxR3OGRERXdbLZHMtsEzSUknbAKuBoaYyQ8Cacvkw4MryXUxLtu8FHpJ0QNkK7XDg0pkPPSIiJmNhry5se4ukY4H1wALgXNubJJ0CDNseAs4BLpQ0QvFEs3r8eEl3Ak8HtpH0RuA1tm8C3gN8Htge+Hr5ExERPaQ2DwrzxuDgoIeHh3sdRkTErCJpo+3BKmUzgkBERNQuySYiImqXZBMREbVLsomIiNol2URERO2SbCIionZJNhERUbskm4iIqF2STURE1C7JJiIiapdkExERtUuyiYiI2iXZRERE7ZJsIiKidkk2ERFRuySbiIioXZJNRETULskmIiJql2QTERG1S7KJiIjaJdlERETteppsJK2UdIukEUkntNi/raSLyv0bJC1p2Hdiuf0WSYc0bP9bSZsk3SjpnyVt151PExERE+lZspG0APgMcCiwHHirpOVNxY4EHrD9XOCTwKnlscuB1cDzgZXAGZIWSFoEHAcM2n4BsKAsFxERPdTLJ5sVwIjt220/CqwFVjWVWQWcXy6vAw6SpHL7WtuP2L4DGCnPB7AQ2F7SQmAH4J6aP0dERHTQy2SzCLi7YX203NayjO0twIPArhMda/vnwD8BPwPuBR60/c1WF5d0lKRhScNjY2Mz8HEiImIivUw2arHNFcu03C5pF4qnnqXAnwB/JOntrS5u+yzbg7YHBwYGJhF2RERMVi+TzSiwR8P6Yp5c5fV4mbJabCdgc5tjDwbusD1m+/fAJcBLa4k+IiIq62WyuRZYJmmppG0oXuQPNZUZAtaUy4cBV9p2uX112VptKbAMuIai+uwASTuU73YOAm7uwmeJiIg2Fvbqwra3SDoWWE/Rauxc25sknQIM2x4CzgEulDRC8USzujx2k6SLgZuALcAxth8DNkhaB1xXbv8hcFa3P1tERGxNxYPC/DY4OOjh4eFehxERMatI2mh7sErZjCAQERG1S7KJiIjaJdlEdNnGux7g8HM2sPGuB3odSkTXJNlEdNlpl9/KVbfdz2mX39rrUCK6pmet0SLmq+MP3murfyPmgySbiC7bb89duODI/XsdRkRXpRotIiJql2QTERG1S7KJiIjaJdlERETtkmwiIqJ2STYRPZZOnjEfJNlE9FhzJ88kn5iL0s8moseaO3mOJx8g/XFizkiyieix5k6eGWEg5qIkm4g+kxEGYi7KO5uIiKhdkk1ERNQuySaiC9LCLOa7jslG0rCkYyTt0o2AIuaizGET812VJ5vVwJ8A10paK+kQSao5rog55fiD9+LAZbulhVnMWx2Tje0R2+8H9gK+CJwL/EzSyZKeUXeAEbNRc7XZeAuz/fZMBUHMT5Xe2Uh6IfC/gI8DXwIOA34FXDmdi0taKekWSSOSTmixf1tJF5X7N0ha0rDvxHL7LZIOadi+s6R1kn4i6WZJL5lOjBFTkWqziK117GcjaSPwS+Ac4ATbj5S7Nkh62VQvLGkB8Bng1cAoRTXdkO2bGoodCTxg+7mSVgOnAm+RtJyieu/5FFV8l0vay/ZjwGnAN2wfJmkbYIepxhgxVemYGbG1Kp0632z79sYNkpbavsP2X07j2iuAkfFzS1oLrAIak80q4KRyeR3w6fJ90SpgbZn47pA0AqyQtAk4EHgngO1HgUenEWPElKRjZsTWqlSjrau4bbIWAXc3rI+W21qWsb0FeBDYtc2xzwHGgPMk/VDS2ZL+qNXFJR1VtrQbHhsbm4GPExERE5kw2Uj6U0l/Bewk6S8bft4JbDcD127Vos0Vy0y0fSGwL/BZ2/sAvwGe9C4IwPZZtgdtDw4MDFSPOqKF9KOJaK9dNdrewOuAnYHXN2x/CPivM3DtUWCPhvXFwD0TlBmVtBDYCdjc5thRYNT2hnL7OiZINhEzKSM1R7Q3YbKxfSlwqaSX2L66hmtfCyyTtBT4OcUL/7c1lRkC1gBXU7SAu9K2JQ0BX5T0CYoGAsuAa2w/JuluSXvbvgU4iK3fAUXUIg0CItqbMNlI+nvb/xN4m6S3Nu+3fdx0Lmx7i6RjgfXAAuBc25sknQIM2x6iaAF3YdkAYDNFQqIsdzFFItkCHFO2RAP4G+ALZUu024EjphNnRBVpEBDRnuzm1yTlDun1tr8qaU2r/bbPrzWyLhocHPTw8HCvw4iImFUkbbQ9WKVsu2q0r5b/zpmkEhERvVFlIM7LJO3csL6LpPX1hhUREXNJlX42A7Z/Ob5i+wHgmfWFFBERc02VZPOYpGePr0jakyf3h4mIiJhQleFq3g98T9J3yvUDgaPqCykiIuaajsnG9jck7QscQNFz/29t3197ZBERMWdUebIBeAy4j2KYmuWSsH1VfWFFRMRcUqU12ruBqyg6X55c/ntSvWFFdFfGNouoV5UGAscDLwbusv1KYB+KkZUj5oxMdhZRryrVaA/bflgSkra1/RNJe9ceWUQXZWyziHpVSTajZafOrwCXSXqAJ4/OHDGrZWyziHp1rEaz/Sbbv7R9EvABisEx31h3YBH9JO90IqanUmu0sunzyyk6c36/nG45Yt7IfDUR01OlNdoHgfMppmPejWLK5X+sO7CIfnL8wXtx4LLd8k4nYoomnGLg8QLSzcA+th8u17cHrrP9vC7E1xWZYmDu2XjXA5x2+a0cf/Be7LfnLr0OJ2JOmswUA1WaPt9J0Zlz3LbAT6cQV0TX1N2UOe9wIianSrJ5BNgk6fOSzgNuBH4t6XRJp9cbXsTUdKr26pQsOu1Pv5yIyanSQODL5c+4b9cTSsTM6dSUudML/0770y8nYnKqDMSZmTpjzumULDrtT7+ciMmp0kBgGfA/gOU0vLux/Zx6Q+ueNBCIiJi8mW4gcB7wWWAL8ErgAuDCqYcXERHzTZVks73tKyiegu4qRxJ41UxcXNJKSbdIGpF0Qov920q6qNy/QdKShn0nlttvkXRI03ELJP1Q0r/ORJwR/Swt42I2qJJsHpb0FOA2ScdKehPwzOleWNIC4DPAoRRVdG+VtLyp2JHAA7afC3wSOLU8djmwGng+sBI4ozzfuOOBm6cbY8RskJZxMRtUSTbvBXYAjgP2A94OrJmBa68ARmzfXg5/sxZY1VRmFcXoBQDrgIMkqdy+1vYjtu8ARsrzIWkx8J+As2cgxoi+l9ENYjaYMNlI2k7SgO1rbf/a9qjtI4Cjgetn4NqLgLsb1kfLbS3L2N4CPEgxbE67Yz8F/D3wh3YXl3SUpGFJw2NjmZ4nZq/xlnFTHSkh1XDRDe2ebE4H/rzF9oMpqrSmSy22NTeNm6hMy+2SXgfcZ3tjp4vbPsv2oO3BgYGBztFGzFGphotuaJdsXm77kuaNtr8AHDgD1x4F9mhYX8yT58l5vIykhcBOwOY2x74MeIOkOymq5V4l6f/OQKwRc0bzk0yq4aIb2iWbVk8PVY6r6lpgmaSlkraheOE/1FRmiCfeDx0GXOmiY9AQsLpsrbYUWAZcY/tE24ttLynPd6Xtt89ArBFzRvOTzHSr4SKqaJc07pO0onmjpBcD037JUb6DORZYT9Fy7GLbmySdIukNZbFzgF0ljQB/B5xQHrsJuBi4CfgGcIztx6YbU8Rc0PzkkieZ6AcTjiBQJpqLgc8D4+9ABoHDgdW2N3QjwG7ICAIxlxx+zgauuu1+Dly2Gxccuf+T1iNmymRGEJhwbDTb15QJ5xjgneXmTcD+tu+bdpQRXTSf5rdpHtctg4ZGP+g4Ntp8kCebuS/f7usznxJ5bG2mx0aL6DuT7Rsym99T9Hs/mDSdjiqqzGcT0Xc6zTfTbDZPCTDZz9ptqaaLKjomG0lvtv0vnbZFdNN8+gPX7591Nify6J4q89lcZ3vfTttms7yziYiYvBlpjSbpUOC1wCJJpzfsejrF3DYRtclL54i5pV0DgXuAYeBhin424z9DwCFtjouYtrx07h+dOolGVNGun82PgB9J+qLt33cxpognvafo9KSTJ6H6NDdQ6PcGC9GfqjR9XiJpnaSbJN0+/lN7ZDGvNH9bbh6vq9OTznx6Eqr7yaLT8DazuRl59E6Vps/nAR+imFbglcARtB+kM2LSOn1b7tQiq99bbM2kup8sms/f3Nosrc9iKqokm+1tXyFJtu8CTpL0XYoEFDEjOiWLTn/g5tMfwLoT60yfP1WcAdWaPn+fYhK1dcCVwM+Bj9neu/7wuiNNn+uRPzIBGSpoLpuRps8N3gvsABwHfJiiKm1N2yMi6P+e79Ed86mKMybWMdnYvhagqEXzEfWHFHNF/sgEzK8qzphYx9Zokl4i6SaKCc6Q9GeSzqg9spj1GluUpW9GxPxWpenzpyg6cf4CHu9/c2CdQcXs1C6hNDdNTvKJmF8qTTFg++6mTZmCOZ6UMNr1dWnumzGf+sXMN/kiEa1UaSBwt6SXApa0DUVDgZvrDStmg+YGAO3e0TTX2+d9ztyVhiHRSpVk89fAacAiYBT4JsVU0THPNSeMybwIzkvjuavOLxJpTj97ZVpo0s8mopuaE0bjOtA2maTPTn+ZqSkG/jcwYSayfdwUYouIea7dwJ7AVvuaE9NkB2ht1i7R5UmpXu2q0Rq/6p9MDcPTSFpJUUW3ADjb9sea9m8LXADsR9Ea7i227yz3nQgcSdFY4Tjb6yXtUZZ/FvAH4Czbp8103BExdc0Jo1W1W3NDEmg9Tttk3w9lBOsest3xB/hhlXKT+aFIMD8FngNsA/wIWN5U5mjgzHJ5NXBRuby8LL8tsLQ8zwJgd2DfssyOwK3N52z1s99++zns4Ts3+x1n/8DDd27udSgRtjv/Tk72d7a5fH7npwcYdsW/+ZWaPtOmOm0aVgAjtm+3/SiwFljVVGYVcH65vA44SJLK7WttP2L7DmAEWGH7XtvXAdh+iKLV3KIaYp8TOjVdThPW6LXmqSY67e800Vtz+U7nj5lTNdnUYRHQ2H9nlCcnhsfL2N4CPAjsWuVYSUuAfYANrS4u6ShJw5KGx8bGpvwhZrPm5NKpL0yST/S75t/Z9OfqHxMmG0kPSfqVpF8BLxxfHt8+A9duNSdO8xPURGXaHivpacCXgPfabhmr7bNsD9oeHBgYqBjy7NLpW15zcmn+lpeOmDHbTHeit5n+QpUvaE9oNy30jjVfexTYo2F9MXDPBGVGJS0EdgI2tztW0lMpEs0XbF9ST+izQ6eXoZOdIyYdMaPfTXeit5luMNB8vvnc+q2X1WjXAsskLS1HJlgNDDWVGeKJ6QwOA64sX0oNAaslbStpKbAMuKZ8n3MOcLPtT3TlU/SxmZ7ON/XbMdc13yOdageadao9mM9V0z3t1CnptRQDfS4AzrX9UUmnULRwGJK0HXAhxbuXzcBq27eXx74feBewhaK67OuSXg58F7iBoukzwD/Y/lq7ONKpMyJaae5E2qlTaaf9zU82s72T6kxPnlabMgl8rWnbBxuWHwbePMGxHwU+2rTte7R+nzMvzOdH9Ig6dOoT1KnTabPJVE3Ptfs5w9Uwd55sZvu3pIjZZqbvucYEM17lNn7uTsmnF8lp1jzZxMzKC/yI7prpe66xQUHzuTs1Xuj30RB62UAgOuj0crJTh7WIqNdM33ONDQo6dUXo1Bih3xofJNn0UKdfhk4d1NLvJWJuaZe8mvc13/+d9jfrdjJKNVoPdXrs7fRyMtVmEfNXp/u/U2OGble7pYEA9TUQ6DSceYY7j4huaW7MMBN/b9JAoE9Mtgd/v7/gi4jZazoz686EJJtpavftYLLVXqkWi4i69Hoq9lSjMb1qtDoeTSMiZoNUo3XRZNvCR0TMR0k205SRkSMiOkuymWG9rheNiOhH6dQZERG1S7KJiIjaJdlERETtkmwiIqJ2STYREVG7JJuIiKhdkk1ERNQuySYiImqXZBMREbXrabKRtFLSLZJGJJ3QYv+2ki4q92+QtKRh34nl9lskHVL1nBER0X09SzaSFgCfAQ4FlgNvlbS8qdiRwAO2nwt8Eji1PHY5sBp4PrASOEPSgornjIiILuvlk80KYMT27bYfBdYCq5rKrALOL5fXAQdJUrl9re1HbN8BjJTnq3LOiIjosl4mm0XA3Q3ro+W2lmVsbwEeBHZtc2yVcwIg6ShJw5KGx8bGpvExIiKik14mG7XY1jyT20RlJrv9yRvts2wP2h4cGBhoG2hERExPL5PNKLBHw/pi4J6JykhaCOwEbG5zbJVzRkREl/Uy2VwLLJO0VNI2FC/8h5rKDAFryuXDgCtdzGM9BKwuW6stBZYB11Q8Z0REdFnPJk+zvUXSscB6YAFwru1Nkk4Bhm0PAecAF0oaoXiiWV0eu0nSxcBNwBbgGNuPAbQ6Z7c/W0REbE3Fg8L8Njg46OHh4V6HERExq0jaaHuwStmMIBAREbVLsomIiNol2URERO2SbCIionZJNhERUbskm4iIqF2STURE1C7JJiIiapdkExERtUuyiYiI2iXZRERE7ZJsIiKidkk2ERFRuySbiIioXZJNRETULskmIiJql2QTERG1S7KJiIjaJdlERETtkmwiIqJ2STYREVG7JJuIiKhdT5KNpGdIukzSbeW/u0xQbk1Z5jZJaxq27yfpBkkjkk6XpHL7xyX9RNKPJX1Z0s7d+kwRETGxXj3ZnABcYXsZcEW5vhVJzwA+BOwPrAA+1JCUPgscBSwrf1aW2y8DXmD7hcCtwIl1foiIiKimV8lmFXB+uXw+8MYWZQ4BLrO92fYDFIlkpaTdgafbvtq2gQvGj7f9TdtbyuN/ACyu80NEREQ1vUo2f2z7XoDy32e2KLMIuLthfbTctqhcbt7e7F3A1ycKQNJRkoYlDY+NjU0y/IiImIyFdZ1Y0uXAs1rsen/VU7TY5jbbG6/9fmAL8IWJTm77LOAsgMHBQU9ULiIipq+2ZGP74In2Sfp/kna3fW9ZLXZfi2KjwCsa1hcD3y63L27afk/DudcArwMOKqvZIiKix3pVjTYEjLcuWwNc2qLMeuA1knYpGwa8BlhfVrs9JOmAshXa4ePHS1oJvA94g+3f1v0hIiKiml4lm48Br5Z0G/Dqch1Jg5LOBrC9GfgwcG35c0q5DeA9wNnACPBTnng382lgR+AySddLOrNLnyciItpQapqKdzbDw8O9DiMiYlaRtNH2YJWyGUEgIiJql2QTERG1S7KJiIjaJdlERETtkmwiIqJ2STYREVG7JJuIiKhdkk1ERNQunToBSWPAXdM4xW7A/TMUzkxLbFPXz/H1c2zQ3/H1c2zQ3/E1x7an7YEqBybZzABJw1V70XZbYpu6fo6vn2OD/o6vn2OD/o5vOrGlGi0iImqXZBMREbVLspkZZ/U6gDYS29T1c3z9HBv0d3z9HBv0d3xTji3vbCIionZ5somIiNol2URERO2SbCqSdK6k+yTd2KbMK8oZQjdJ+k6/xCZpJ0lflfSjMrYjuhjbHpK+Jenm8trHtygjSadLGpH0Y0n79ll8/6WM68eS/k3Sn/VLbA1lXyzpMUmH9VNsPbwnqvx/7cl9IWk7Sdc0XPfkFmW2lXRReU9skLSkG7FNIr6/k3RTeU9cIWnPjie2nZ8KP8CBwL7AjRPs3xm4CXh2uf7MPortH4BTy+UBYDOwTZdi2x3Yt1zeEbgVWN5U5rUUU3sLOADY0MX/dlXieymwS7l8aLfiqxJbuW8BcCXwNeCwfomtx/dElfh6cl+Uv+dPK5efCmwADmgqczRwZrm8Grioi//tqsT3SmCHcvk9VeLLk01Ftq+i+GWcyNuAS2z/rCx/X1cCo1JsBnaUJOBpZdktXYrtXtvXlcsPATcDi5qKrQIucOEHwM6Sdu+X+Gz/m+0HytUfAIv7JbbS3wBfArr5O1cltl7eE1Xi68l9Uf6e/7pcfWr509xSaxVwfrm8DjiojLN2VeKz/S3bvy1XK90TSTYzZy9gF0nflrRR0uG9DqjBp4HnAfcANwDH2/5Dt4MoqwL2ofim1GgRcHfD+iit/6jWqk18jY6keArrqolik7QIeBNwZrdjaohhCa3/u/XFPdEmvp7dF5IWSLqe4gvCZbYnvCdsbwEeBHbtRmwV42tU6Z5YOFPBBQuB/YCDgO2BqyX9wPatvQ0LgEOA64FXAf8BuEzSd23/qlsBSHoaxbfv97a4bqtvbF1tk98hvvEyr6S4sV7eR7F9Cnif7ce69MV3MrH1/J7oEF/P7gvbjwEvkrQz8GVJL7Dd+M61p/dEhfgAkPR2YBD4i07nzJPNzBkFvmH7N7bvB64CuvIiuYIjKKozbHsEuAP4025dXNJTKW74L9i+pEWRUWCPhvXFFN82u6JCfEh6IXA2sMr2L/ootkFgraQ7gcOAMyS9sU9i6+k9USG+nt4XALZ/CXwbWNm06/F7QtJCYCfaV5XXok18SDoYeD/wBtuPdDpXks3MuRT4c0kLJe0A7E9RT9wPfkbx7RJJfwzsDdzejQuX9cznADfb/sQExYaAw8tWaQcAD9q+t1/ik/Rs4BLgHV3+Vt4xNttLbS+xvYSibv9o21/ph9jo4T1RMb6e3BeSBsonBiRtDxwM/KSp2BCwplw+DLjS5dv4fohP0j7A5ygSTaV3calGq0jSPwOvAHaTNAp8iOLFGbbPtH2zpG8APwb+AJzd6rGzF7EBHwY+L+kGisfz95XfNLvhZcA7gBvKOmAoWgE9uyG+r1G0SBsBfkvxjbNbqsT3QYr68jPKqqot7s6ovFVi65WOsfXynqgSH727L3YHzpe0gOIL/8W2/1XSKcCw7SGKRHmhpBGKJ5rVXYhrMvF9nKJRxb+U98TPbL+h3UkzXE1ERNQu1WgREVG7JJuIiKhdkk1ERNQuySYiImqXZBMREbVL0+eINiTtClxRrj4LeAwYK9dX2H60J4G1IeldwNds/3uvY4kYl6bPERVJOgn4te1/6oNYFpRDirTa9z3gWNvXt9o/wTELyzG4ImqRarSIKZK0ppz343pJZ0h6Stlb/peSPi7pOknrJe0v6TuSbpf02vLYd0v6crn/Fkn/WPG8H5F0DbBC0smSrpV0o6QzyxEY3gK8CLioPH4bSaMNPcIPkHR5ufwRSZ+TdBlwXnmNT5TX/rGkd3f/v2rMVUk2EVMg6QUUoy2/1PaLKKqkx3t57wR80/a+wKPASRTDorwZOKXhNCvKY/YF3ibpRRXOe53tFbavBk6z/WLgP5b7Vtq+iGJwybfYflGFar59gNfbfgdwFHCf7RXAi4FjyqF6IqYt72wipuZgij/Iw+VwHdvzxDQJv7N9Wbl8A8VYb1vKYVGWNJxj/fg8OZK+QjGa9MI2530U+HLD8QdJ+u/AdsBuwEYmP/3BpbYfLpdfAzxPUmNyW0YxhljEtCTZREyNgHNtf2CrjcUIvY1PE38AHmlYbrznml+YusN5fzc+GGM5sOWnKWaj/Lmkj1AknVa28EQtRnOZ3zR9pqNtX0HEDEs1WsTUXA78Z0m7QdFqbQpVTq+RtHOZOFYB35/EebenSF73S9oR+KuGfQ9RTIU87k6KeWVoKtdsPXB0mdiQtHc56m/EtOXJJmIKbN8g6WTgcklPAX4P/DWTm4fne8AXKSbuunC89ViV89r+haTzgRuBu9h6FsrzgLMl/Y7ivdBJwP+R9O/ANW3i+RzFqMjXl1V491EkwYhpS9PniB4oW3q9wPZ7ex1LRDekGi0iImqXJ5uIiKhdnmwiIqJ2STYREVG7JJuIiKhdkk1ERNQuySYiImr3/wFmKRFweiQSTAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYYAAAEGCAYAAABhMDI9AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAYq0lEQVR4nO3dfZRkdX3n8fdH0KgRgeOMq8vDDOsBH9YoMC0QzRofJgmSs7A5YRWIoobIMYrBh83R3WQNEv9J3JjgRqJEEHET1Lg+EBcljkExRpAewzNhnIOOjJjDoBPwAUXwu3/UbVNV9nTdnu5bVd39fp3TZ+pW3br1nYFbn/493N9NVSFJ0pyHTLoASdJ0MRgkSQMMBknSAINBkjTAYJAkDdh30gUs1rp162rjxo2TLkOSVpStW7feXVXr2+y74oJh48aNzM7OTroMSVpRkuxou69dSZKkAZ0FQ5KLktyV5KY9vJ4k70iyPckNSY7uqhZJUntdthguBo5f4PUXAIc3P2cCf9FhLZKkljoLhqq6Cvj2ArucBFxSPVcDByR5fFf1SJLameQYw0HAHX3bO5vnJEkTNMlgyDzPzbuiX5Izk8wmmd21a1fHZUnS2jbJYNgJHNK3fTBw53w7VtUFVTVTVTPr17eahitJ2kuTDIbLgNOb2UnHAfdU1TcnWI8kiW6nq14KfBF4YpKdSc5I8sokr2x2uRy4HdgO/CXwqq5qmbN1x25Ov/Aatu7YPe+2JKnDK5+r6tQRrxfw6q4+fz7nbdnGVV+5G4BLzjj2p7a37tjNeVu2cfbmI9i04cCB7bn3z70mSavVilsSYynmvuD39OdCwQG0DhGDQ9JKtqaCYdOGA7nkjGP3uD0qOPofL6X1YXBImmZrKhhGGRUciwmRUcEhSdPKYNhLS2192KKQNK1cXbUjc0Ex9yU/vD3Xgjhvy7Z5t51BJWlSDIYJOXvzETz78HUDLYn+7cUGhyQtF7uSJmSpXVEOdkvqisEwpRYbHA52S1ouBsMKtTdTbSWpjfQuQF45ZmZmyns+L55dTdLalmRrVc202dfB5zXCwWxJbdmVtEYsdTBb0tphMKwRSx3MNiiktcNgELD4oJC0ehkMmtdCQWHrQVrdDAa10h8Up194ja0HaRVzVpIWbXj5Dtd1klYXWwxatOFuJq+6llYXg0FL5lXX0urilc/qnIPV0uR55bOmyvBV1pKmm11J6pxdS9LKYotBnRu+e52zlqTpZjBo7OxakqabwaCxG3UdhKTJMhg0dsNdSy4BLk0XB581cS7YJ00Xg0ETN2plV0njZTBo6gwHhaTxcoxBkjTAYJAkDTAYNPWcpSSNl8GgqecFcdJ4OfisqecsJWm8DAZNPWcpSePVaVdSkuOT3JZke5I3zfP6oUmuTPJPSW5IckKX9UiSRussGJLsA7wTeAHwFODUJE8Z2u33gQ9V1VHAKcD5XdUjSWqnyxbDMcD2qrq9qu4HPgCcNLRPAY9uHu8P3NlhPVqlnLUkLa8ug+Eg4I6+7Z3Nc/3OAV6cZCdwOfCa+Q6U5Mwks0lmd+3a1UWtWkGGg8BZS9Ly6jIYMs9zwzeYPhW4uKoOBk4A3p/kp2qqqguqaqaqZtavX99BqVpJhoPAZbyl5dXlrKSdwCF92wfz011FZwDHA1TVF5M8HFgH3NVhXVrhhqevDs9acnVWaWm6DIZrgcOTHAZ8g97g8mlD+3wdeD5wcZInAw8H7CvSgkZNX/W6B2lpOutKqqoHgLOAK4Bb6c0+ujnJuUlObHZ7A/CKJNcDlwIvq6rh7iZpUbzHtLQ0nV7gVlWX0xtU7n/uzX2PbwGe1WUNkl1L0uJ45bNWPbuWpMUxGLTquaSGtDiuripJGmAwSJIGGAySpAEGgyRpgMEgSRpgMEiSBhgMkqQBBoMkaYDBIEkaYDBIkgYYDFpzXG1VWpjBoDXHW4FKC3MRPa05rrYqLcxg0JrjaqvSwuxKkiQNMBgkSQMMBknSAINBkjTAYJAkDWg1KynJs4BzgA3NewJUVf2H7kqTJE1C2+mqFwKvA7YCD3ZXjiRp0toGwz1V9clOK5EkTYW2wXBlkrcBHwF+OPdkVX25k6okSRPTNhjmLhOd6XuugOctbzmSpElrFQxV9dyuC5EkTYdW01WT7J/k7Ulmm58/SbJ/18VJksav7XUMFwHfAV7Y/NwLvLeroiRJk9N2jOEJVfXrfdtvSXJdFwVJkiarbYvhviS/MLfRXPB2XzclSZImqW2L4beB9zXjCgG+Dbysq6IkSZPTdlbSdcDTkzy62b6306okSROzYDAkeXFV/Z8krx96HoCqenuHtUmSJmDUGMPPNn/uN8/Po0YdPMnxSW5Lsj3Jm/awzwuT3JLk5iR/vYjaJUkdWLDFUFXvbh5uqaov9L/WDEDvUZJ9gHcCvwTsBK5NcllV3dK3z+HAfweeVVW7kzx2L/4OkqRl1HZW0v9u+Vy/Y4DtVXV7Vd0PfAA4aWifVwDvrKrdAFV1V8t6JEkdGTXG8PPAM4H1Q+MMjwb2GXHsg4A7+rZ38m9rLs05ovmcLzTHO6eqPjVPHWcCZwIceuihIz5WkrQUo1oMD6M3lrAvg+ML9wInj3hv5nmuhrb3BQ4HngOcCrwnyQE/9aaqC6pqpqpm1q9fP+JjJUlLMWqM4XPA55JcXFU7FnnsncAhfdsHA3fOs8/VVfUj4KtJbqMXFNcu8rMkScuk7QVu32/ux/AfgYfPPVlVCy27fS1weJLDgG8ApwCnDe3zMXothYuTrKPXtXR7y5okSR1oO/j8V8A/A4cBbwG+xojf6qvqAeAs4ArgVuBDVXVzknOTnNjsdgXwrSS3AFcCv1tV31r030KStGxSNdztP89Oydaq2pTkhqp6WvPc56rqFzuvcMjMzEzNzs6O+2MlaUVrvsdnRu/ZvivpR82f30zyq/TGCg7em+IkSdOtbTC8tVlA7w30rl94NPC6zqqSJE1M22C4vqruAe4BnguQ5HGdVSVJmpi2g89fTXJpkkf2PXd5FwVJkiarbTDcCHwe+HySJzTPzXcBmyRphWvblVRVdX6S64G/TfJGfvoqZknSKtA2GAJQVV9I8nzgg8CTOqtKkjQxbYPhhLkHVfXNJM+jt7ieJGmVaXUHN+DUubu2Dbmqk6qkKbZ1x27O27KNszcfwaYNB066HGnZjWox9N/BbZhjDFqTztuyjau+cjcAl5wxvJK8tPJ1dgc3abU6e/MRA39Kq02Xd3CTVoWtO3Zz+oXXsHXHbgA2bTiQS8441m4krVpd3sFNWhXsOtJaM2qMYfgObnPa3MFNWhXsOtJa03bZ7Q1VtSPJz1bV98ZQ1x657LYkLd5ilt1uO8bw75ub6dzafMDTk5y/twVKkqZX22D4M+BXgG8BVNX1wLO7KkqSNDltg4GqumPoqQeXuRZJ0hRouyTGHUmeCVSShwG/Q9OtJElaXdq2GF4JvBo4CNgJHNlsS5JWmVYthqq6G/iNjmuRJE2BVsGQZD3wCmBj/3uq6je7KUuSNCltxxg+Tu8Obltw0FmSVrW2wfDIqnpjp5VIkqZC28HnTyQ5YfRukqSVrm0wnE0vHO5Lcm+S7yS5t8vCpJViePVVaaVrOytpvhv1SMLVV7X6tJ2VdPQ8T98D7KiqB5a3JGllcfVVrTZtB5/PB44Gbmy2fw64HnhMkldW1d91UZy0EszduEdaLdqOMXwNOKqqNlXVJnpXPt8EbAb+uKPaJEkT0DYYnlRVN89tVNUt9ILi9m7KkiRNStuupNuS/AXwgWb7RcC2JD8D/KiTyiRJE9G2xfAyYDvwWuB1wO3Ncz8CnttFYZKkyWg7XfU+4E+an2HfXdaKpAnbumM3523Zxtmbj2DThgMnXY40dq1aDEkOT/LhJLckuX3up8X7jk9yW5LtSd60wH4nJ6kkre5HKnVp7rqE87Zsm3Qp0kS0HWN4L/AHwJ/S6zp6OZCF3pBkH+CdwC/Ru4fDtUkuawau+/fbj96Nf65ZXOlSN7wuQWtd2zGGR1TVZ4BU1Y6qOgd43oj3HANsr6rbq+p+egPXJ82z3x/Sm/L6g5a1SJ2auy7BbiStVW2D4QdJHgJ8JclZSX4NeOyI9xwE9N8nemfz3E8kOQo4pKo+sdCBkpyZZDbJ7K5du1qWLEnaG22D4bXAI+l1+WwCXgK8dMR75utqqp+82AuaPwXeMOrDq+qCqpqpqpn169e3LFmStDfazkq6tnn4XXrjC23sBA7p2z4YuLNvez/gqcBnkwA8DrgsyYlVNdvyMyRJy2zBYEhy2UKvV9WJC7x8LXB4ksOAbwCnAKf1vfceYF3fZ30W+G+GgiRN1qgWw8/TGye4lN6soQVnIvWrqgeSnAVcAewDXFRVNyc5F5itqgVDR5I0GaOC4XH0ppueSu+3/f8HXNq/btJCqupy4PKh5968h32f0+aYkqRuLTj4XFUPVtWnquqlwHH0lsX4bJLXjKU6SdLYjRx8bhbK+1V6rYaNwDuAj3RbliRpUkYNPr+P3syhTwJvqaqbxlKVJGliRrUYXgJ8DzgC+J1mWin0BqGrqh7dYW2SpAlYMBiqqu0FcJKkVcIvfknSAINBkjTAYJAkDTAYJEkDDAZJ0gCDQZI0wGDQmrd1x25Ov/Aatu7YPelSpKlgMGjNO2/LNq76yt2ct2Vb559lCGklaHWjHmk1O3vzEQN/dmkuhAAuOePYzj9P2hsGg9a8TRsOHNuX9DhDSNpbBoM0RuMMIWlvOcYgSRpgMEiSBhgMkqQBBoMkaYDBIEkaYDBIkgYYDJKkAQaDNEVcMkPTwGCQpsg4122S9sRgkKbI2ZuP4NmHr9vjkhm2KDQOLokhTZFRS2a4CJ/GwWCQVhAX4dM42JUkLaOuu3rmWhSbNhzYyfElMBikZeXgsVYDu5KkZWRXj1YDWwzSMpp0V4+zlrQcDAZpFbErS8vBriRpFbErS8uh0xZDkuOT3JZke5I3zfP665PckuSGJJ9JsqHLeqTVbtJdWVodOguGJPsA7wReADwFODXJU4Z2+ydgpqqeBnwY+OOu6pEktdNli+EYYHtV3V5V9wMfAE7q36Gqrqyq7zebVwMHd1iPJKmFLoPhIOCOvu2dzXN7cgbwyfleSHJmktkks7t27VrGEiVJw7oMhszzXM27Y/JiYAZ423yvV9UFVTVTVTPr169fxhKl6eb0U01Cl8GwEzikb/tg4M7hnZJsBn4POLGqfthhPdKKM2r6qcGhLnQ5XfVa4PAkhwHfAE4BTuvfIclRwLuB46vqrg5rkVakUdNPXW1VXegsGKrqgSRnAVcA+wAXVdXNSc4FZqvqMnpdR48C/iYJwNer6sSuapJWmlHLcHvdgrqQqnm7/afWzMxMzc7OTroMSVpRkmytqpk2+7okhjRB4x4jcExCbRgM0gR1vbbRcBC4lpLaMBikCRp1j+elGg6C4c9b7haELZLVwUX0pAkaNbi8VMOD08Oft9yzmpwltToYDNIqNu5ZTc6SWh2clSRJa4CzkiRJe81gkDq00gZjV1q96obBIHVopU0PXWn1qhsGg9ShrqejLrfFTmddzhaGrZXp4awkqUNdT0ddboudzrqc01Od6jo9DAZJezRq+ulyTk91quv0cLqqJK0BTleV1ImljAMMv3ec4xdaHINBUmtLmbU0/N5Rx5r0DKm1HEyOMUhqbTHjAFt37Oa8Lds4e/MRbNpw4E+9d6njF8PHX25reTDcMQZJnTj9wmu46it38+zD1y3LF+twEIw6/lKDo+vgGTfHGCRN3HJfwzFqCfFR+y/W3NTduVBY7JjISu6KMhgkdWL4i3WphoNg1Bf3qIv1FvtFvtgxkUmPkSyFYwySVoRRFwsOjwmMulhv1PawxY6JLLT/cDfVtHVbGQySVoWlfnGPev9w0Ax/mQ+/3r89vO+oUJp0UDj4LEl7YTGD68P7jmoxjNp/byxm8NkWgyTthcVM3R11i9Xh7eH9xz111haDJE05WwySpAHjXqXX6aqSpAEGgyRpgMEgSRpgMEiSBhgMkqQBBoMkaYDBIEkaYDBIkgasuCufk+wCdizhEOuAu5epnOU2zbXBdNdnbXtvmuub5tpguusbrm1DVa1v88YVFwxLlWS27WXh4zbNtcF012dte2+a65vm2mC661tKbXYlSZIGGAySpAFrMRgumHQBC5jm2mC667O2vTfN9U1zbTDd9e11bWtujEGStLC12GKQJC3AYJAkDViVwZDkoiR3JblpgX2ek+S6JDcn+dy01JZk/yR/m+T6praXj7G2Q5JcmeTW5rPPnmefJHlHku1Jbkhy9JTV9xtNXTck+cckT5+W2vr2fUaSB5OcPI7aFlPfJM6Llv9dJ3lePDzJl/o++y3z7PMzST7YnBfXJNk4RbW9PsktzTnxmSQbRh64qlbdD/Bs4Gjgpj28fgBwC3Bos/3YKartfwB/1DxeD3wbeNiYans8cHTzeD9gG/CUoX1OAD4JBDgOuGaM/3Zt6nsmcGDz+AXjqq9Nbc1r+wB/D1wOnDxl/3YTOS9a1jbJ8yLAo5rHDwWuAY4b2udVwLuax6cAH5yi2p4LPLJ5/NttaluVLYaquore/zh7chrwkar6erP/XWMpjFa1FbBfkgCPavZ9YEy1fbOqvtw8/g5wK3DQ0G4nAZdUz9XAAUkePy31VdU/VtXuZvNq4OBpqa3xGuD/AmP7f66pqU19EzkvWtY2yfOiquq7zeZDm5/hWTsnAe9rHn8YeH5T68Rrq6orq+r7zWarc2JVBkMLRwAHJvlskq1JTp90QX3+HHgycCdwI3B2Vf143EU0TeGj6P0G0u8g4I6+7Z3M/wXYqQXq63cGvdbNWO2ptiQHAb8GvGvcNQ3VsZH5/+0mfl4sUNtEz4sk+yS5jl6gf7qq9nheVNUDwD3AY6aktn6tzol9l6u4FWZfYBPwfOARwBeTXF1V2yZbFgC/AlwHPA94AvDpJJ+vqnvHVUCSR9H7rfa183zufL8FjXXO84j65vZ5Lr2T4BemqLY/A95YVQ+O4ZfJeY2ob6LnxYjaJnpeVNWDwJFJDgA+muSpVdU/Tjix86JFbQAkeTEwA/ziqGOu1RbDTuBTVfW9qrobuAoYyyBlCy+n15yvqtoOfBV40rg+PMlD6Z2cf1VVH5lnl53AIX3bB9P7LW4sWtRHkqcB7wFOqqpvTVFtM8AHknwNOBk4P8l/maL6JnZetKhtoufFnKr6V+CzwPFDL/3kvEiyL7A/C3cZj7M2kmwGfg84sap+OOpYazUYPg78pyT7JnkkcCy9fs1p8HV6v7GR5N8BTwRuH8cHN32iFwK3VtXb97DbZcDpzeyk44B7quqb01JfkkOBjwAvGWcLsE1tVXVYVW2sqo30+qFfVVUfm5b6mNB50bK2SZ4X65vfxknyCGAz8M9Du10GvLR5fDLw99WM9k66tiRHAe+mFwqtxo1WZVdSkkuB5wDrkuwE/oDeoAxV9a6qujXJp4AbgB8D75mv6TWJ2oA/BC5OciO95ukbm9/exuFZwEuAG5s+S+jNBjm0r77L6c1M2g58n95vcuPSpr430+vbPb/prnmgxrP6ZZvaJmlkfRM8L9r8203yvHg88L4k+9D7ZfpDVfWJJOcCs1V1Gb1ge3+S7fRaCqdMUW1vozdg/zfNOfH1qjpxoYO6JIYkacBa7UqSJO2BwSBJGmAwSJIGGAySpAEGgyRpwKqcrqq1KcljgM80m48DHgR2NdvHVNX9EylsAUl+E7i8qv5l0rVIc5yuqlUpyTnAd6vqf01BLfs0yxbM99o/AGdV1XXzvb6H9+zbrMcjdcKuJK0JSV7arFt/XZLzkzykucL3X5O8LcmXk1yR5Ngkn0tye5ITmvf+VpKPNq/fluT3Wx73rUm+BByT5C1Jrk1yU5J3NVeOvwg4Evhg8/6HJdnZdyXrcUm2NI/fmuTdST4NvLf5jLc3n31Dkt8a/7+qViuDQatekqfSW9X0mVV1JL0u1LkrU/cH/q6qjgbuB86ht/TCfwXO7TvMMc17jgZOS3Jki+N+uaqOqaovAudV1TOAn2teO76qPkhvYbgXVdWRLbq6jgL+c1W9BDgTuKuqjgGeAby6WQ5EWjLHGLQWbKb35TnbLAnwCP5t6fD7qurTzeMb6a399ECz9MLGvmNcMXefhyQfo7dq674LHPd+4KN9739+kt8FHg6sA7ay+CXBP15VP2ge/zLw5CT9QXQ4vTWFpCUxGLQWBLioqv7nwJO9VTD7f0v/MfDDvsf958fwYFyNOO59c4uoNQvS/Tm9u5R9I8lb6QXEfB7g31ryw/t8b+jv9Kqq+gzSMrMrSWvBFuCFSdZBb/bSXnS7/HKSA5ov+ZOALyziuI+gFzR3J9kP+PW+175D73aWc75G754IDO037ArgVU0IkeSJzeqa0pLZYtCqV1U3pneT9C1JHgL8CHgli7uPxD8Af03vJjHvn5tF1Oa4VfWtJO8DbgJ2MHh3svcC70lyH71xjHOAv0zyL8CXFqjn3fRWH72u6ca6i15gSUvmdFVphGbGz1Or6rWTrkUaB7uSJEkDbDFIkgbYYpAkDTAYJEkDDAZJ0gCDQZI0wGCQJA34/1O/g5BL/ijsAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZIAAAEGCAYAAABPdROvAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAdzUlEQVR4nO3df5xddX3n8debCQREBJpAZRNCsk2iAu0iDIH+WPxBwGALwRJLFIG18ZGtSsX10d1CW7pK6ba0XW1UFAKBAq4llB91rGgkgFAthsxIBAIFx0hkCC6JRBZRiYOf/eN8B24ud+acO2fO/THzfj4e87jnx/ec+ZyZe+7nnu/3e75HEYGZmdl47dHuAMzMrLs5kZiZWSlOJGZmVooTiZmZleJEYmZmpUxrdwCtMHPmzJg7d267wzAz6yoDAwM7IuKgvHJTIpHMnTuX/v7+dodhZtZVJG0tUs5VW2ZmVooTiZmZleJEYmZmpTiRmJlZKU4kZmZWihOJmZmV4kRiZmalOJGYmVkpTiRmZlZKpYlE0hJJj0oalHRBg/XTJa1N6zdImpuWL5K0Kf18W9I7iu7TzMxaq7JEIqkHuAw4BTgceJekw+uKrQB2RsR84BPApWn5Q0BvRBwFLAGukDSt4D7NzKyFqrwiWQQMRsSWiNgF3AAsrSuzFLg2Td8EnChJEfGTiBhOy/cGRp4HXGSfZmbWQlUmklnAEzXzQ2lZwzIpcTwLzACQdJykzcCDwB+k9UX2Sdp+paR+Sf3bt2+fgMMxM7NGqkwkarAsipaJiA0RcQRwLHChpL0L7pO0/eqI6I2I3oMOyh0F2czMxqnKRDIEHFozPxvYNloZSdOA/YFnagtExCPA88CRBfdpZmYtVGUi2QgskDRP0l7AcqCvrkwfcG6aXgbcGRGRtpkGIOkw4HXA4wX3aWZmLVTZg60iYljSecA6oAe4OiI2S7oY6I+IPmANcL2kQbIrkeVp898CLpD0c+AXwAciYgdAo31WdQxmZpZPEQ2bGCaV3t7e8BMSzcyaI2kgInrzyvnOdjMzK8WJxMzMSnEiMTOzUpxIzMysFCcSMzMrxYnEzMxKcSIxM7NSnEjMzKwUJxKzJg1s3ck5azYwsHVnu0Mx6whOJGZNWrX+Me75zg5WrX+s3aGYdYTKxtoym6zOX7xwt1ezqc6JxKxJxxx2INetOK7dYZh1DFdtmZlZKU4kZmZWihOJmZmV4kRiZmalOJGYmVkpTiRmZlaKE4mZmZXiRGJmZqVUmkgkLZH0qKRBSRc0WD9d0tq0foOkuWn5SZIGJD2YXt9as83X0j43pZ+DqzwGMzMbW2V3tkvqAS4DTgKGgI2S+iLi4ZpiK4CdETFf0nLgUuBMYAdwakRsk3QksA6YVbPdWRHRX1XsZmZWXJVXJIuAwYjYEhG7gBuApXVllgLXpumbgBMlKSLuj4htaflmYG9J0yuM1WzcPBqwTXVVJpJZwBM180PsflWxW5mIGAaeBWbUlTkDuD8iXqhZdk2q1rpIkhr9ckkrJfVL6t++fXuZ4zAbk0cDtqmuykEbG33ARzNlJB1BVt11cs36syLiSUn7ATcDZwPXvWInEauB1QC9vb31v9dswng0YJvqqrwiGQIOrZmfDWwbrYykacD+wDNpfjZwK3BORHx3ZIOIeDK9Pgd8nqwKzaxtRkYDPuawA9sdillbVJlINgILJM2TtBewHOirK9MHnJumlwF3RkRIOgD4EnBhRHxjpLCkaZJmpuk9gd8BHqrwGMzMLEdliSS1eZxH1uPqEeDGiNgs6WJJp6Via4AZkgaBjwAjXYTPA+YDF9V1850OrJP0ALAJeBK4sqpjMDOzfIqY/M0Hvb290d/v3sJmZs2QNBARvXnlfGe7mZmV4kRiZmalOJGYmVkpTiRmZlaKE4mZmZXiRGJmZqU4kZiZWSlOJGZmVooTiZmZleJEYmZmpTiRmJlZKU4kZmZWihOJmZmV4kRiZmalOJGYmVkpTiRmZlaKE4mZmZXiRGJmZqU4kZiZWSlOJGZmVkqliUTSEkmPShqUdEGD9dMlrU3rN0iam5afJGlA0oPp9a012xyTlg9K+qQkVXkMZmY2tsoSiaQe4DLgFOBw4F2SDq8rtgLYGRHzgU8Al6blO4BTI+JXgXOB62u2+SywEliQfpZUdQxm4zGwdSfnrNnAwNad7Q7FrCWqvCJZBAxGxJaI2AXcACytK7MUuDZN3wScKEkRcX9EbEvLNwN7p6uXQ4DXRMS9ERHAdcDpFR6DWdNWrX+Me76zg1XrH2t3KGYtMa3Cfc8CnqiZHwKOG61MRAxLehaYQXZFMuIM4P6IeEHSrLSf2n3OavTLJa0ku3Jhzpw5JQ7DrDnnL16426vZZFdlImnUdhHNlJF0BFl118lN7DNbGLEaWA3Q29vbsIxZFY457ECuW1H/ncls8qqyamsIOLRmfjawbbQykqYB+wPPpPnZwK3AORHx3Zrys3P2aWZmLVRlItkILJA0T9JewHKgr65MH1ljOsAy4M6ICEkHAF8CLoyIb4wUjoingOckHZ96a50DfKHCYzAzsxyVJZKIGAbOA9YBjwA3RsRmSRdLOi0VWwPMkDQIfAQY6SJ8HjAfuEjSpvRzcFr3fuAqYBD4LvDlqo7BzMzyKev8NLn19vZGf39/u8MwA7LuwavWP8b5ixdyzGEHtjscs1FJGoiI3rxyha5IJN0s6bcl+U54s5LcPdgmm6KJ4bPAu4HvSPprSa+vMCazSe38xQs5YcFMdw+2SaOpqi1J+wPvAv6U7P6PK4HPRcTPqwlvYrhqy8yseRNatZV2OAP4L8D7gPuBVcDRwO3jjNHMzCaBQjckSroFeD3ZmFenpm64AGsl+au+mdkUVvTO9qsi4rbaBZKmR8QLRS57zMxs8ipatXVJg2X3TmQgZmbWnca8IpH0WrJBEfeR9EZeHuvqNcCrKo7NzMy6QF7V1tvIGthnAx+vWf4c8CcVxWRmZl1kzEQSEdcC10o6IyJublFMZmbWRfKqtt4TEZ8D5kr6SP36iPh4g83MzGwKyava2je9vrrqQMzMrDvlVW1dkV4/1ppwzMys2+RVbX1yrPUR8aGJDcfMzLpNXtXWQEuiMDOzrlWk15aZmdmo8qq2/j4iPizpi8ArhgmOiNMabGZmZlNIXtXW9en176oOxMzMulNe1dZAer1b0l5kIwAH8GhE7GpBfGZm1uGKDiP/28DlwHfJxtuaJ+m/RsSXqwzOzMw6X9HRf/838JaIeHNEvAl4C/CJvI0kLZH0qKRBSRc0WD9d0tq0foOkuWn5DEl3SfqxpE/XbfO1tM9N6efggsdgZmYVKPo8kqcjYrBmfgvw9FgbSOoBLgNOAoaAjZL6IuLhmmIrgJ0RMV/ScuBS4EzgZ8BFwJHpp95ZEeEHapmZdYC8Xlu/myY3S7oNuJGsjeSdwMacfS8CBiNiS9rXDcBSoDaRLAU+mqZvAj4tSRHxPPB1SfObOBYzM2uDvCuSU2um/y/wpjS9HTgwZ9tZwBM180PAcaOViYhhSc8CM4AdOfu+RtKLwM3AJRHxiq7JZmbWGnm9tt5bYt9qsKz+A79ImXpnRcSTkvYjSyRnA9e94pdLK4GVAHPmzMmP1szMxiWvaut/RMTfSPoUjW9IHGusrSHg0Jr52cC2UcoMSZoG7A88M1ZMEfFken1O0ufJqtBekUgiYjWwGqC3t9dXLGZmFcmr2nokvY6nYXsjsEDSPOBJYDnw7royfcC5ZM9/XwbcOVY1VUo2B0TEDkl7Ar8DrB9HbGZmNkHyqra+mCZ/EhH/VLtO0jtzth2WdB6wDugBro6IzZIuBvojog9YA1wvaZDsSmR5zf4fJ3s2/F6STgdOBrYC61IS6SFLIlcWPVgzM5t4KtJOLelbEXF03rJO1dvbG/397i1sZtYMSQMR0ZtXLq+N5BTg7cCsumeTvAYYLheimZlNBnltJNvI2kdOY/dnkzwH/LeqgjIzs+6R10bybeDbqXeU8KCNZmZWp+gQKScBV+BBG83MrE7RRPJxskEbBwEk/QrwJcCJxMxsiis6+m/TgzaamdnUUPSKpOGgjSODOkbELRXFZ2ZmHa5oItmbVw7a+EtkgzoG4ERiZjZFFUokJQdvNDOzSaxQG4mkhZLukPRQmv81SX9WbWhmZtYNija2XwlcCPwcICIeoGZcLDPrXANbd3LOmg0MbN3Z7lBskiqaSF4VEffVLfMQKWZdYNX6x7jnOztYtf6xdodik1TRxvYd6d6RAJC0DHiqsqjMbMKcv3jhbq9mE61oIvkg2UOiXi/pSeB7wHsqi8rMJswxhx3IdSvqn3JtNnGK9traAiyWtC+wR0Q8V21YZmbWLYr22vpfkg6IiOfTI24PlHRJ1cGZmVnnK9rYfkpE/GhkJiJ2kj2nxMzMpriiiaRH0vSRGUn7ANPHKG9mZlNE0cb2zwF3SLqGrOfW7wPXVhaVmZl1jUJXJBHxN8AlwBuAI4C/SMvMLIdvCLTJrmhj+77AVyPij8i6AU+XtGelkZlNEmVvCHQisk5XtI3kHmBvSbOA9cB7gX/I20jSEkmPShqUdEGD9dMlrU3rN0iam5bPkHSXpB9L+nTdNsdIejBt80lJKngMZm1x/uKFnLBg5rhvCPSd6dbpiiYSRcRPgN8FPhUR7wAOH3MDqQe4DDgllX2XpPptVgA7I2I+8Ang0rT8Z8BFwB812PVngZXAgvSzpOAxmLXFyA2Bxxx24Li2L5uIzKpWOJFI+nXgLLJH7EJ+Q/0iYDAitkTELuAGYGldmaW83Gh/E3CiJKX7Vb5OllBqgzgEeE1E3BsRAVwHnF7wGMy6Ul4ictWXtVvRRPJhstF/b42IzZL+I3BXzjazgCdq5ofSsoZlImIYeBaYkbPPoZx9AiBppaR+Sf3bt2/PCdWse7nqy9qt6BApdwN318xvAT6Us1mjtosYR5lxlY+I1WQdA+jt7R1rn2ZdzYMyWrsVSiSS7qLBB3ZEvHWMzYaAQ2vmZwPbRikzJGkasD/wTM4+Z+fs02xK8aCM1m5Fb0isbfTeGziD/OeRbAQWSJoHPEn2IKx315XpA84F7gWWAXemto+GIuIpSc9JOh7YAJwDfKrgMZiZWQWKVm0N1C36hqS7GxZ+eZthSecB64Ae4OrUvnIx0B8RfcAa4HpJg2RXIi89dVHS48BrgL0knQ6cHBEPA+8n63q8D/Dl9GNmycDWnaxa/xjnL1447p5iZs0oWrX1SzWzewC9wGvztouI24Db6pb9ec30z4B3jrLt3FGW9wNH5gZtNkWNNL4DrvKylihatTXAy20kw8DjZPeAmFmHceO7tdqYiUTSscATETEvzZ9L1j7yOPBw5dGZWdPc+G6tlncfyRXALgBJJwB/RXYD4bOkrrVmZja15VVt9UTESHfcM4HVEXEzcLOkTdWGZmZm3SDviqQn3d8BcCJwZ826ou0rZjYGD3Fi3S4vGfwjcLekHcBPgX8FkDSfrHrLzEpyLyvrdmMmkoj4S0l3AIeQPY9kpOfWHsAfVh2c2VTgXlbW7TTGjeSTRm9vb/T397c7DDOzriJpICJ688oVHf3XzMysIScSM9uNG/+tWU4kZlNcfeLw802sWe7CazbF1fcac+O/NcuJxKzLTPTovvWJw0OsWLOcSMy6zETfd+LEYWW5jcSsy5y/eCEnLJg57qonN6bbRPMViVmXKXsF4TvpbaI5kZhNMW5Mt4nmqi2zKWbkiqZIQ319NZirxawRJxIzG1X9PSW+x8QacdWWmY2qvhrM1WLWSKWDNkpaAqwCeoCrIuKv69ZPB64DjgF+CJwZEY+ndReSPRf+ReBDEbEuLX8ceC4tHy4yoJgHbTQza17RQRsruyKR1ANcBpwEDAEbJfVFRO2z3lcAOyNivqTlwKXAmZIOB5YDRwD/AVgvaWFEvJi2e0tE7KgqdjMzK67KNpJFwGBEbImIXcANwNK6MkvJngEPcBNwoiSl5TdExAsR8T1gMO3PrOXcwGw2tioTySzgiZr5obSsYZmIGCZ76uKMnG0D+KqkAUkrR/vlklZK6pfUv3379lIHYlObG5jNxlZlY7saLKtvkBmtzFjb/mZEbJN0MHC7pH+PiHteUThiNbAasjaS4mGb7c4NzGZjqzKRDAGH1szPBraNUmZI0jRgf+CZsbaNiJHXpyXdSlbl9YpEYjZRPBaV2diqrNraCCyQNE/SXmSN5311ZfqAc9P0MuDO9Fz4PmC5pOmS5gELgPsk7StpPwBJ+wInAw9VeAxmZpajsiuSiBiWdB6wjqz779URsVnSxUB/RPQBa4DrJQ2SXYksT9tulnQj8DAwDHwwIl6U9MvArVl7PNOAz0fEV6o6BjMzy1fpfSSdwveRmFVjop+NYp2l6H0kHiLFrMN0U3dj92gzcCIx6zjd9OHsZ6MYeKwts47TTd2N/WwUAycSs47Tzd2Nm20z6aakaaNzIjGzCdPsFUY3J017mdtIzGzC5LWZuE1kcnIiMbMJk/f0xao7EjhRtYertsysZapuE3HjfXv4isQMf5NtlWaeFw/5/5f69WW7I9v4OJGY0V33bnSTZhN0ffm8/0v9+mYTlU0MJxKbEvI+0PxNthrNJuj68vX/F1+BdCYnEpsS6j+g6j+Q/E22Gs1+0NeXr/+/NHsFUvaKqBlTuXrUje02KdXfGFffyOtG2dZo9j6RvPLNNtY3+38u876Yyu8pJxKblOpP6voPqLE+kDyibedqNjHV/5/z/rdlepWV7ZHWze87DyNvk1KZk/KcNRu45zs7OGHBzCn3zXKya/Z/W+Z9lLdt/fpOfN95GHmb1PLqo8u0ebgBd/LKa7yvN1bbWpltG63v5lEBXLVlXanK+miP/zR51f9v895HY7WtAePe9roVx71ifX1s9Vcs9dvXr29n1ZgTiXUljxprEyHvfVSkbW282+Z9YclLPPXr8xJNldxGYm0x0W/ybm6oNGuk2TaWKtpciraROJFYW+S9yZu9bO/EhkqzdpqIL1cd0dguaYmkRyUNSrqgwfrpktam9Rskza1Zd2Fa/qiktxXdp02MZsc4KjsmUl5DZN4d0m4gN9tdK2+yrSyRSOoBLgNOAQ4H3iXp8LpiK4CdETEf+ARwadr2cGA5cASwBPiMpJ6C+5wy8j7Mm5nP+yBv9oM+r3zeHcv1iSEvUfjOdLP2qfKKZBEwGBFbImIXcAOwtK7MUuDaNH0TcKIkpeU3RMQLEfE9YDDtr8g+J8xEflBXMZ/3Yd7MfN4HebMf9M12bcwbGsOJwqxzVdlraxbwRM38EFBfef1SmYgYlvQsMCMt/2bdtrPSdN4+AZC0ElgJMGfOnHEdQF6viHbP1/fiKPtaO91sj5OyPVTc5dase1WZSNRgWX3L/mhlRlve6AqqYW+BiFgNrIassX30MEc30R/UE/2a92He7PxEftA7MZhNHZX12pL068BHI+Jtaf5CgIj4q5oy61KZeyVNA34AHARcUFt2pFzabMx9NuJeW2ZmzeuEXlsbgQWS5knai6zxvK+uTB9wbppeBtwZWWbrA5anXl3zgAXAfQX3aWZmLVRZ1VZq8zgPWAf0AFdHxGZJFwP9EdEHrAGulzQIPEOWGEjlbgQeBoaBD0bEiwCN9lnVMZiZWT7fkGhmZg11QtWWmZlNAU4kZmZWihOJmZmV4kRiZmalOJGYmVkpTiRmZlaKE4mZmZXiRGJmZqU4kZiZWSlOJGZmVooTiZmZleJEYmZmpTiRmJlZKU4kZmZWihOJmZmV4kRiZmalOJGYmVkpTiRmZlaKE4mZmZXiRGJmZqU4kZiZWSlOJGZmVooTiZmZlaKIaHcMlZO0Hdg6zs1nAjsmMJyJ1snxdXJs0NnxObbx6+T4Ojk2eGV8h0XEQXkbTYlEUoak/ojobXcco+nk+Do5Nujs+Bzb+HVyfJ0cG4w/PldtmZlZKU4kZmZWihNJvtXtDiBHJ8fXybFBZ8fn2Mavk+Pr5NhgnPG5jcTMzErxFYmZmZXiRGJmZqU4kQCSrpb0tKSHxijzZkmbJG2WdHcnxSdpf0lflPTtFN97WxjboZLukvRI+t3nNygjSZ+UNCjpAUlHd1BsZ6WYHpD0b5L+UytiKxpfTdljJb0oaVknxdau86Lg/7Yt54WkvSXdV/N7P9agzHRJa9M5sUHS3A6K7SOSHk7nxB2SDsvdcURM+R/gBOBo4KFR1h8APAzMSfMHd1h8fwJcmqYPAp4B9mpRbIcAR6fp/YDHgMPryrwd+DIg4HhgQwfF9hvAgWn6lFbFVjS+tK4HuBO4DVjWKbG187woGF9bzov0Pn91mt4T2AAcX1fmA8DlaXo5sLZFf7cisb0FeFWafn+R2HxFAkTEPWRvstG8G7glIr6fyj/dksCSAvEFsJ8kAa9OZYdbFNtTEfGtNP0c8Agwq67YUuC6yHwTOEDSIZ0QW0T8W0TsTLPfBGZXHVcz8SV/CNwMtOx9VzC2tp0XBeNry3mR3uc/TrN7pp/6Xk1LgWvT9E3AiSnOtscWEXdFxE/SbKFzwomkmIXAgZK+JmlA0jntDqjOp4E3ANuAB4HzI+IXrQ4iXZ6/kexbTq1ZwBM180M0/sCszBix1VpBduXUcqPFJ2kW8A7g8tZH9VIMc2n8t+uI82KM+Np2XkjqkbSJLPnfHhGjnhMRMQw8C8zokNhqFTonpk1UcJPcNOAY4ERgH+BeSd+MiMfaG9ZL3gZsAt4K/Apwu6R/jYj/16oAJL2a7Fvzhxv83kbftFrW7zwntpEybyE7aX6rVXHV/O6x4vt74I8j4sUWfGFtNra2nxc58bXtvIiIF4GjJB0A3CrpyIiobeNs2zlRIDYAJL0H6AXelLdPX5EUMwR8JSKej4gdwD1AyxplC3gvWRVDRMQg8D3g9a365ZL2JDuZ/09E3NKgyBBwaM38bLJviZ0QG5J+DbgKWBoRP2xFXE3E1wvcIOlxYBnwGUmnd0hsbT0vCsTX1vMCICJ+BHwNWFK36qVzQtI0YH/Grr5uZWxIWgz8KXBaRLyQty8nkmK+APxnSdMkvQo4jqxOtlN8n+xbIZJ+GXgdsKUVvzjV664BHomIj49SrA84J/XeOh54NiKe6oTYJM0BbgHObvUVZpH4ImJeRMyNiLlkdekfiIh/7oTYaON5UTC+tpwXkg5K3/aRtA+wGPj3umJ9wLlpehlwZ6TW7XbHJumNwBVkSaRQu5ertgBJ/wi8GZgpaQj4n2SNUETE5RHxiKSvAA8AvwCuanQp2K74gL8A/kHSg2SXzH+cviG2wm8CZwMPpnpXyHrLzKmJ7zaynluDwE/Ivil2Smx/TlY3/ZlUdTQcrRudtUh87ZIbW5vPiyJ/u3adF4cA10rqIfuyfmNE/Iuki4H+iOgjS4LXSxokuxJZ3oK4isb2t2SdE/4pnRPfj4jTxtqph0gxM7NSXLVlZmalOJGYmVkpTiRmZlaKE4mZmZXiRGJmZqW4+69NaZJmAHek2dcCLwLb0/yiiNjVlsDGIOn3gdsi4gftjsUM3P3X7CWSPgr8OCL+rgNi6UlDWTRa93XgvIjY1Gj9KNtMS2M6mU04V22ZjULSuenZDZskfUbSHuku7h9J+ltJ35K0TtJxku6WtEXS29O275N0a1r/qKQ/K7jfSyTdByyS9DFJGyU9JOnyNDLAmcBRwNq0/V6ShmruVj5e0vo0fYmkKyTdDlyTfsfH0+9+QNL7Wv9XtcnIicSsAUlHko26+xsRcRRZNfDI3cf7A1+NiKOBXcBHyYbieCdwcc1uFqVtjgbeLemoAvv9VkQsioh7gVURcSzwq2ndkohYSzYQ4ZkRcVSBqrc3AqdGxNnASuDpiFgEHAt8MA0RY1aK20jMGltM9mHbn4aJ2IeXh8L/aUTcnqYfJBs7bDgNxTG3Zh/rRp51IumfyUYWnjbGfncBt9Zsf6Kk/w7sDcwEBmh+mPsvRMTP0vTJwBsk1SauBWRjUpmNmxOJWWMCro6Ii3ZbmI3UWnsV8AvghZrp2nOqvgEycvb705GB+9IgiJ8mewrgk5IuIUsojQzzcu1CfZnn647pAxFxB2YTyFVbZo2tB35P0kzIeneNoxroZEkHpKSwFPhGE/vdhywx7ZC0H3BGzbrnyB4vO+JxsueCUFeu3jrgAylpIel1aQRYs1J8RWLWQEQ8KOljwHpJewA/B/6A5p6j8nXg82QPVbp+pJdVkf1GxA8lXQs8BGxl96f/XQNcJemnZO0wHwWulPQD4L4x4rmCbHTcTala7WmyBGdWirv/mlUg9Yg6MiI+3O5YzKrmqi0zMyvFVyRmZlaKr0jMzKwUJxIzMyvFicTMzEpxIjEzs1KcSMzMrJT/DzhBO6wxVEOkAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Plot calculated stuff\n",
    "\n",
    "plt.scatter(temps,aveEn,s=2)\n",
    "plt.xlabel('Temperature')\n",
    "plt.ylabel('Average Energy')\n",
    "plt.show()\n",
    "\n",
    "plt.scatter(temps,heatCap,s=2)\n",
    "plt.xlabel('Temperature')\n",
    "plt.ylabel('Heat Capacity')\n",
    "plt.show()\n",
    "\n",
    "plt.scatter(temps,mag,s=2)\n",
    "plt.xlabel('Temperature')\n",
    "plt.ylabel('Magnetization')\n",
    "plt.show()\n",
    "\n",
    "plt.scatter(temps,sus,s=2)\n",
    "plt.xlabel('Temperature')\n",
    "plt.ylabel('Susceptibility')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Reduced helper functions\n",
    "\n",
    "# Calculate magnetization. \n",
    "# Magnetization is the sum of all spins divided by the total number of spins\n",
    "@njit(fastmath=True)\n",
    "def calcM(s, width):\n",
    "    m = np.abs(s.sum())\n",
    "    return m/(width**2)\n",
    "\n",
    "# Calculate interaction energy between spins. Assume periodic boundaries\n",
    "# Interaction energy will be the difference in energy due to flipping spin i,j \n",
    "# (i.e. 2*spin_value*neighboring_spins)\n",
    "@njit(fastmath=True)\n",
    "def deltaU(s, i, j, width):\n",
    "    #top\n",
    "    if i == 0:\n",
    "        t = s[width-1,j]  # periodic boundary\n",
    "    else:\n",
    "        t = s[i-1,j]\n",
    "    #bottom\n",
    "    if i == width-1:\n",
    "        b = s[0,j]  # periodic boundary\n",
    "    else:\n",
    "        b = s[i+1,j]\n",
    "    #left\n",
    "    if j == 0:\n",
    "        l = s[i,width-1]  # periodic boundary\n",
    "    else:\n",
    "        l = s[i,j-1]\n",
    "    #right\n",
    "    if j == width-1:\n",
    "        r = s[i,0]  # periodic boundary\n",
    "    else:\n",
    "        r = s[i,j+1]\n",
    "    return 2*s[i,j]*(t+b+r+l)  # difference in energy if i,j is flipped\n",
    "\n",
    "# Monte carlo cycle\n",
    "@njit(fastmath=True)\n",
    "def montecc(s, T, cycles, width):\n",
    "    for m in range(cycles):\n",
    "        i = random.randrange(width)  # choose random row\n",
    "        j = random.randrange(width)  # choose random column\n",
    "        ediff = deltaU(s, i, j, width)\n",
    "        if ediff <= 0:\n",
    "            s[i,j] = -s[i,j]  # flip spin\n",
    "        elif random.random() < np.exp(-ediff/T):\n",
    "            s[i,j] = -s[i,j]\n",
    "    return s\n",
    "\n",
    "# Compute physical quantities\n",
    "@njit(fastmath=True)\n",
    "def computePQ(s, T, cycles, width):\n",
    "    Mg = 0\n",
    "    Mg_sq = 0\n",
    "    for p in range(cycles):\n",
    "        s = montecc(s, T, 1, width)\n",
    "        M = calcM(s, width)\n",
    "        Mg += M\n",
    "        Mg_sq += M*M\n",
    "    sus = (Mg_sq/cycles-(Mg/cycles)**2)/T\n",
    "    return sus"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Elapsed = 20.1909019947052\n"
     ]
    }
   ],
   "source": [
    "# Initial parameters\n",
    "\n",
    "width = 16  # width of sqaure spin lattice; 64 points may take hours\n",
    "s = np.random.choice([1,-1], size=(width, width))  # initial spin sites randomly (+1 or -1)\n",
    "cycles = 3000000\n",
    "\n",
    "# Initlaize temperature range (crosses critical temperature)\n",
    "temps = np.arange(1.9, 3.2, 0.04)\n",
    "\n",
    "# Inititalize magnetization, average energy, heat capacity, and susceptibility\n",
    "sus_16 = np.zeros(len(temps))\n",
    "\n",
    "start = time.time()\n",
    "# Compute physical quanities (magnetization, etc.)\n",
    "for ind, T in enumerate(temps):\n",
    "    # Cycle spins\n",
    "    s = montecc(s, T, cycles, width)\n",
    "    sus_16[ind] = computePQ(s, T, cycles, width)\n",
    "end = time.time()\n",
    "print(\"Elapsed = %s\" % (end - start))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Elapsed = 23.322134017944336\n"
     ]
    }
   ],
   "source": [
    "# Initial parameters\n",
    "\n",
    "width = 32  # width of sqaure spin lattice; 64 points may take hours\n",
    "s = np.random.choice([1,-1], size=(width, width))  # initial spin sites randomly (+1 or -1)\n",
    "\n",
    "# Inititalize magnetization, average energy, heat capacity, and susceptibility\n",
    "sus_32 = np.zeros(len(temps))\n",
    "\n",
    "start = time.time()\n",
    "# Compute physical quanities (magnetization, etc.)\n",
    "for ind, T in enumerate(temps):\n",
    "    # Cycle spins\n",
    "    s = montecc(s, T, cycles, width)\n",
    "    sus_32[ind] = computePQ(s, T, cycles, width)\n",
    "end = time.time()\n",
    "print(\"Elapsed = %s\" % (end - start))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Elapsed = 47.70257210731506\n"
     ]
    }
   ],
   "source": [
    "# Initial parameters\n",
    "\n",
    "width = 64  # width of sqaure spin lattice; 64 points may take hours\n",
    "s = np.random.choice([1,-1], size=(width, width))  # initial spin sites randomly (+1 or -1)\n",
    "\n",
    "# Inititalize magnetization, average energy, heat capacity, and susceptibility\n",
    "sus_64 = np.zeros(len(temps))\n",
    "\n",
    "start = time.time()\n",
    "# Compute physical quanities (magnetization, etc.)\n",
    "for ind, T in enumerate(temps):\n",
    "    # Cycle spins\n",
    "    s = montecc(s, T, cycles, width)\n",
    "    sus_64[ind] = computePQ(s, T, cycles, width)\n",
    "end = time.time()\n",
    "print(\"Elapsed = %s\" % (end - start))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "\n",
       "\n",
       "\n",
       "\n",
       "\n",
       "\n",
       "  <div class=\"bk-root\" id=\"f50845e0-0aa6-42fd-bd81-8f8a33b230c4\" data-root-id=\"5847\"></div>\n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/javascript": [
       "(function(root) {\n",
       "  function embed_document(root) {\n",
       "    \n",
       "  var docs_json = {\"54c5b00f-19e2-4e38-9c87-c50a4661194d\":{\"roots\":{\"references\":[{\"attributes\":{\"below\":[{\"id\":\"5856\",\"type\":\"LinearAxis\"}],\"center\":[{\"id\":\"5860\",\"type\":\"Grid\"},{\"id\":\"5865\",\"type\":\"Grid\"}],\"left\":[{\"id\":\"5861\",\"type\":\"LinearAxis\"}],\"plot_height\":350,\"plot_width\":700,\"renderers\":[{\"id\":\"5882\",\"type\":\"GlyphRenderer\"},{\"id\":\"5887\",\"type\":\"GlyphRenderer\"},{\"id\":\"5892\",\"type\":\"GlyphRenderer\"}],\"title\":{\"id\":\"6142\",\"type\":\"Title\"},\"toolbar\":{\"id\":\"5872\",\"type\":\"Toolbar\"},\"x_range\":{\"id\":\"5848\",\"type\":\"DataRange1d\"},\"x_scale\":{\"id\":\"5852\",\"type\":\"LinearScale\"},\"y_range\":{\"id\":\"5850\",\"type\":\"DataRange1d\"},\"y_scale\":{\"id\":\"5854\",\"type\":\"LinearScale\"}},\"id\":\"5847\",\"subtype\":\"Figure\",\"type\":\"Plot\"},{\"attributes\":{\"data_source\":{\"id\":\"5884\",\"type\":\"ColumnDataSource\"},\"glyph\":{\"id\":\"5885\",\"type\":\"Line\"},\"hover_glyph\":null,\"muted_glyph\":null,\"nonselection_glyph\":{\"id\":\"5886\",\"type\":\"Line\"},\"selection_glyph\":null,\"view\":{\"id\":\"5888\",\"type\":\"CDSView\"}},\"id\":\"5887\",\"type\":\"GlyphRenderer\"},{\"attributes\":{\"formatter\":{\"id\":\"6143\",\"type\":\"BasicTickFormatter\"},\"ticker\":{\"id\":\"5862\",\"type\":\"BasicTicker\"}},\"id\":\"5861\",\"type\":\"LinearAxis\"},{\"attributes\":{\"source\":{\"id\":\"5884\",\"type\":\"ColumnDataSource\"}},\"id\":\"5888\",\"type\":\"CDSView\"},{\"attributes\":{},\"id\":\"5862\",\"type\":\"BasicTicker\"},{\"attributes\":{\"line_color\":\"brown\",\"x\":{\"field\":\"x\"},\"y\":{\"field\":\"y\"}},\"id\":\"5890\",\"type\":\"Line\"},{\"attributes\":{\"text\":\"\"},\"id\":\"6142\",\"type\":\"Title\"},{\"attributes\":{},\"id\":\"6143\",\"type\":\"BasicTickFormatter\"},{\"attributes\":{\"dimension\":1,\"ticker\":{\"id\":\"5862\",\"type\":\"BasicTicker\"}},\"id\":\"5865\",\"type\":\"Grid\"},{\"attributes\":{\"line_alpha\":0.1,\"line_color\":\"#1f77b4\",\"x\":{\"field\":\"x\"},\"y\":{\"field\":\"y\"}},\"id\":\"5891\",\"type\":\"Line\"},{\"attributes\":{\"data_source\":{\"id\":\"5879\",\"type\":\"ColumnDataSource\"},\"glyph\":{\"id\":\"5880\",\"type\":\"Line\"},\"hover_glyph\":null,\"muted_glyph\":null,\"nonselection_glyph\":{\"id\":\"5881\",\"type\":\"Line\"},\"selection_glyph\":null,\"view\":{\"id\":\"5883\",\"type\":\"CDSView\"}},\"id\":\"5882\",\"type\":\"GlyphRenderer\"},{\"attributes\":{\"data_source\":{\"id\":\"5889\",\"type\":\"ColumnDataSource\"},\"glyph\":{\"id\":\"5890\",\"type\":\"Line\"},\"hover_glyph\":null,\"muted_glyph\":null,\"nonselection_glyph\":{\"id\":\"5891\",\"type\":\"Line\"},\"selection_glyph\":null,\"view\":{\"id\":\"5893\",\"type\":\"CDSView\"}},\"id\":\"5892\",\"type\":\"GlyphRenderer\"},{\"attributes\":{\"line_alpha\":0.1,\"line_color\":\"#1f77b4\",\"x\":{\"field\":\"x\"},\"y\":{\"field\":\"y\"}},\"id\":\"5881\",\"type\":\"Line\"},{\"attributes\":{\"ticker\":{\"id\":\"5857\",\"type\":\"BasicTicker\"}},\"id\":\"5860\",\"type\":\"Grid\"},{\"attributes\":{\"source\":{\"id\":\"5889\",\"type\":\"ColumnDataSource\"}},\"id\":\"5893\",\"type\":\"CDSView\"},{\"attributes\":{\"active_drag\":\"auto\",\"active_inspect\":\"auto\",\"active_multi\":null,\"active_scroll\":\"auto\",\"active_tap\":\"auto\",\"tools\":[{\"id\":\"5866\",\"type\":\"PanTool\"},{\"id\":\"5867\",\"type\":\"WheelZoomTool\"},{\"id\":\"5868\",\"type\":\"BoxZoomTool\"},{\"id\":\"5869\",\"type\":\"SaveTool\"},{\"id\":\"5870\",\"type\":\"ResetTool\"},{\"id\":\"5871\",\"type\":\"HelpTool\"}]},\"id\":\"5872\",\"type\":\"Toolbar\"},{\"attributes\":{},\"id\":\"5866\",\"type\":\"PanTool\"},{\"attributes\":{},\"id\":\"6145\",\"type\":\"BasicTickFormatter\"},{\"attributes\":{},\"id\":\"5867\",\"type\":\"WheelZoomTool\"},{\"attributes\":{\"bottom_units\":\"screen\",\"fill_alpha\":{\"value\":0.5},\"fill_color\":{\"value\":\"lightgrey\"},\"left_units\":\"screen\",\"level\":\"overlay\",\"line_alpha\":{\"value\":1.0},\"line_color\":{\"value\":\"black\"},\"line_dash\":[4,4],\"line_width\":{\"value\":2},\"render_mode\":\"css\",\"right_units\":\"screen\",\"top_units\":\"screen\"},\"id\":\"6147\",\"type\":\"BoxAnnotation\"},{\"attributes\":{\"overlay\":{\"id\":\"6147\",\"type\":\"BoxAnnotation\"}},\"id\":\"5868\",\"type\":\"BoxZoomTool\"},{\"attributes\":{\"callback\":null,\"data\":{\"x\":{\"__ndarray__\":\"ZXIFMbmCSMD1I91f+FRGwIXVtI43J0TAFoeMvXb5QcBMccjYa5c/wG3UdzbqOzvAjTcnlGjgNsCumtbx5oQywJ77C5/KUizA38FqWsebI8BBEJMriMkVwBBzQokGbvG/cq3jzQklCkA2SjRwi4AeQNpeu/xI9ydATEyuICZXMEAr6f7Cp7I0QAuGT2UpDjlA6iKgB6tpPUDlX/hUluJAQFSuICZXEENAxPxI9xc+RUAzS3HI2GtHQKOZmZmZmUlAE+jBalrHS0CCNuo7G/VNQHlCiQZuEVBAsWkdb04oUUDpkLHXLj9SQCC4RUAPVlNAWN/ZqO9sVECQBm4R0INVQMgtAnqwmlZA\",\"dtype\":\"float64\",\"shape\":[33]},\"y\":{\"__ndarray__\":\"aty7TAq0iD/S0ruqynCRP3UKsRfJL5I/ouIp54immz8brXqlIkOiP2p9DJeA1qs/r4XVx5+IvD+gjHaMlbHBP7wQdwuX8M4/P9LRycKFyz88YC7z70zRP2UcnnItvdM/aCSQHV1g2j/9g8xVrZzYP5Cpc0bgbdc/w0DZ1VDi0z8NQmm+ownTP+tJFm65LM4/emtxYxp6yz8dM6eROqjIP93Z9L4HCMU/ra2HJHRDxD+cxz00qGLAPy6N5J5Ol8A/HzqX8rYgvT+gCgCq4XW4P8/mzEoogLg/g5cW++/gtT9YN0lvnomzP60ClB/6c7Q/Lo/fhnMcsT8wuOQxzkiwP/wUoNmmLa8/\",\"dtype\":\"float64\",\"shape\":[33]}},\"selected\":{\"id\":\"6148\",\"type\":\"Selection\"},\"selection_policy\":{\"id\":\"6149\",\"type\":\"UnionRenderers\"}},\"id\":\"5879\",\"type\":\"ColumnDataSource\"},{\"attributes\":{},\"id\":\"6148\",\"type\":\"Selection\"},{\"attributes\":{\"callback\":null,\"data\":{\"x\":{\"__ndarray__\":\"ZXIFMbmCaMD1I91f+FRmwIXVtI43J2TAFoeMvXb5YcBMccjYa5dfwG3UdzbqO1vAjTcnlGjgVsCumtbx5oRSwJ77C5/KUkzA38FqWsebQ8BBEJMriMk1wBBzQokGbhHAcq3jzQklKkA2SjRwi4A+QNpeu/xI90dATEyuICZXUEAr6f7Cp7JUQAuGT2UpDllA6iKgB6tpXUDlX/hUluJgQFSuICZXEGNAxPxI9xc+ZUAzS3HI2GtnQKOZmZmZmWlAE+jBalrHa0CCNuo7G/VtQHlCiQZuEXBAsWkdb04ocUDpkLHXLj9yQCC4RUAPVnNAWN/ZqO9sdECQBm4R0IN1QMgtAnqwmnZA\",\"dtype\":\"float64\",\"shape\":[33]},\"y\":{\"__ndarray__\":\"+r8tGoMQfj/x71ZrykeBP+MzFM+kJII/V4T0Y9A+lD+Vg9fBvKSYP28KBt/v7JY/g8B5TDx/lj/xgviCB2q+P74tv3ygicw/j0Nlp7yn1z8fbGGyn0/pP0hI1svY+Ng/Q6WgpzaD3j95oEFPLGHdP4/8mse40dU/d/TaFCJHzT+993UK7DXJP3EoJ9/+9cM/2rWviM4Ywj+DqXHGEDq6P/gMB/D/Aro/8s6OTburtj+EmRinbICxP/VWydx6vKs/MCt5eyXPsT/tFIQ93D6qP/qzsUCrM6U/QlDr0Y0Kpj/NAc6/R1OkP4DeGpW71KA/OTZgvj+8nT/3MDcqu3igPzO1YI8X06A/\",\"dtype\":\"float64\",\"shape\":[33]}},\"selected\":{\"id\":\"6150\",\"type\":\"Selection\"},\"selection_policy\":{\"id\":\"6151\",\"type\":\"UnionRenderers\"}},\"id\":\"5884\",\"type\":\"ColumnDataSource\"},{\"attributes\":{},\"id\":\"5869\",\"type\":\"SaveTool\"},{\"attributes\":{},\"id\":\"6149\",\"type\":\"UnionRenderers\"},{\"attributes\":{},\"id\":\"5852\",\"type\":\"LinearScale\"},{\"attributes\":{},\"id\":\"5870\",\"type\":\"ResetTool\"},{\"attributes\":{},\"id\":\"6150\",\"type\":\"Selection\"},{\"attributes\":{\"callback\":null},\"id\":\"5850\",\"type\":\"DataRange1d\"},{\"attributes\":{},\"id\":\"5871\",\"type\":\"HelpTool\"},{\"attributes\":{},\"id\":\"5857\",\"type\":\"BasicTicker\"},{\"attributes\":{},\"id\":\"6151\",\"type\":\"UnionRenderers\"},{\"attributes\":{\"formatter\":{\"id\":\"6145\",\"type\":\"BasicTickFormatter\"},\"ticker\":{\"id\":\"5857\",\"type\":\"BasicTicker\"}},\"id\":\"5856\",\"type\":\"LinearAxis\"},{\"attributes\":{},\"id\":\"6152\",\"type\":\"Selection\"},{\"attributes\":{\"source\":{\"id\":\"5879\",\"type\":\"ColumnDataSource\"}},\"id\":\"5883\",\"type\":\"CDSView\"},{\"attributes\":{},\"id\":\"6153\",\"type\":\"UnionRenderers\"},{\"attributes\":{\"line_color\":\"red\",\"x\":{\"field\":\"x\"},\"y\":{\"field\":\"y\"}},\"id\":\"5880\",\"type\":\"Line\"},{\"attributes\":{\"callback\":null,\"data\":{\"x\":{\"__ndarray__\":\"ZXIFMbmCiMD1I91f+FSGwIXVtI43J4TAFoeMvXb5gcBMccjYa5d/wG3UdzbqO3vAjTcnlGjgdsCumtbx5oRywJ77C5/KUmzA38FqWsebY8BBEJMriMlVwBBzQokGbjHAcq3jzQklSkA2SjRwi4BeQNpeu/xI92dATEyuICZXcEAr6f7Cp7J0QAuGT2UpDnlA6iKgB6tpfUDlX/hUluKAQFSuICZXEINAxPxI9xc+hUAzS3HI2GuHQKOZmZmZmYlAE+jBalrHi0CCNuo7G/WNQHlCiQZuEZBAsWkdb04okUDpkLHXLj+SQCC4RUAPVpNAWN/ZqO9slECQBm4R0IOVQMgtAnqwmpZA\",\"dtype\":\"float64\",\"shape\":[33]},\"y\":{\"__ndarray__\":\"QASYCCpzaT9Bu4OUM3RyP41YTXV3BnI/WVsRn975fT+oaaZljzSEP4hTwzzAJKQ/ky2uuHpQjz+qi5oEo76OPxqwc+ce96A/+f0EDztfsj86KRzvuvq1P+eO+LEwNek/07AzX0nSwT+qFL75sMm5P7zpN0KSZc4/CRJ7063cxD8tSpz95BSqPz9IMRZCFLA/geKm1emvrz9gMLheIF2zPw13zQm32LA/B+BesUd9oT+tOITWAgekPzslVFKFOJo/hFhdxOJcmT/vjecP2iygP1+5QROZc6A/4JZf9APKkz8im5iyHWybP2SIRvOhnJI/wBux3jTpiz+7JE5vRwyQP7Z0kgw0yZQ/\",\"dtype\":\"float64\",\"shape\":[33]}},\"selected\":{\"id\":\"6152\",\"type\":\"Selection\"},\"selection_policy\":{\"id\":\"6153\",\"type\":\"UnionRenderers\"}},\"id\":\"5889\",\"type\":\"ColumnDataSource\"},{\"attributes\":{\"line_color\":\"blue\",\"x\":{\"field\":\"x\"},\"y\":{\"field\":\"y\"}},\"id\":\"5885\",\"type\":\"Line\"},{\"attributes\":{\"callback\":null},\"id\":\"5848\",\"type\":\"DataRange1d\"},{\"attributes\":{},\"id\":\"5854\",\"type\":\"LinearScale\"},{\"attributes\":{\"line_alpha\":0.1,\"line_color\":\"#1f77b4\",\"x\":{\"field\":\"x\"},\"y\":{\"field\":\"y\"}},\"id\":\"5886\",\"type\":\"Line\"}],\"root_ids\":[\"5847\"]},\"title\":\"Bokeh Application\",\"version\":\"1.3.1\"}};\n",
       "  var render_items = [{\"docid\":\"54c5b00f-19e2-4e38-9c87-c50a4661194d\",\"roots\":{\"5847\":\"f50845e0-0aa6-42fd-bd81-8f8a33b230c4\"}}];\n",
       "  root.Bokeh.embed.embed_items_notebook(docs_json, render_items);\n",
       "\n",
       "  }\n",
       "  if (root.Bokeh !== undefined) {\n",
       "    embed_document(root);\n",
       "  } else {\n",
       "    var attempts = 0;\n",
       "    var timer = setInterval(function(root) {\n",
       "      if (root.Bokeh !== undefined) {\n",
       "        embed_document(root);\n",
       "        clearInterval(timer);\n",
       "      }\n",
       "      attempts++;\n",
       "      if (attempts > 100) {\n",
       "        console.log(\"Bokeh: ERROR: Unable to run BokehJS code because BokehJS library is missing\");\n",
       "        clearInterval(timer);\n",
       "      }\n",
       "    }, 10, root)\n",
       "  }\n",
       "})(window);"
      ],
      "application/vnd.bokehjs_exec.v0+json": ""
     },
     "metadata": {
      "application/vnd.bokehjs_exec.v0+json": {
       "id": "5847"
      }
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "nu = 1\n",
    "eta = 2.5\n",
    "tc = 2.35\n",
    "t = (temps-tc)/tc\n",
    "x16 = t * (16**2)**(1/nu)\n",
    "y16 = sus_16 * (16**2)**(eta - 2)\n",
    "x32 = t * (32**2)**(1/nu)\n",
    "y32 = sus_32 * (32**2)**(eta - 2)\n",
    "x64 = t * (64**2)**(1/nu)\n",
    "y64 = sus_64 * (64**2)**(eta - 2)\n",
    "\n",
    "#Interactive plot\n",
    "p = figure(plot_width=700, plot_height=350)\n",
    "p.line(x16, y16, color='red')\n",
    "p.line(x32, y32, color='blue')\n",
    "p.line(x64, y64, color='brown')\n",
    "show(p)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "\n",
       "\n",
       "\n",
       "\n",
       "\n",
       "\n",
       "  <div class=\"bk-root\" id=\"967a6915-46b9-4a15-9adb-5912af1459ad\" data-root-id=\"1140\"></div>\n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/javascript": [
       "(function(root) {\n",
       "  function embed_document(root) {\n",
       "    \n",
       "  var docs_json = {\"9d2f71d0-7344-43ae-aaa8-14c117ef28d4\":{\"roots\":{\"references\":[{\"attributes\":{\"below\":[{\"id\":\"1149\",\"type\":\"LinearAxis\"}],\"center\":[{\"id\":\"1153\",\"type\":\"Grid\"},{\"id\":\"1158\",\"type\":\"Grid\"}],\"left\":[{\"id\":\"1154\",\"type\":\"LinearAxis\"}],\"plot_height\":400,\"renderers\":[{\"id\":\"1175\",\"type\":\"GlyphRenderer\"},{\"id\":\"1180\",\"type\":\"GlyphRenderer\"},{\"id\":\"1185\",\"type\":\"GlyphRenderer\"}],\"title\":{\"id\":\"1201\",\"type\":\"Title\"},\"toolbar\":{\"id\":\"1165\",\"type\":\"Toolbar\"},\"x_range\":{\"id\":\"1141\",\"type\":\"DataRange1d\"},\"x_scale\":{\"id\":\"1145\",\"type\":\"LinearScale\"},\"y_range\":{\"id\":\"1143\",\"type\":\"DataRange1d\"},\"y_scale\":{\"id\":\"1147\",\"type\":\"LinearScale\"}},\"id\":\"1140\",\"subtype\":\"Figure\",\"type\":\"Plot\"},{\"attributes\":{},\"id\":\"1163\",\"type\":\"ResetTool\"},{\"attributes\":{},\"id\":\"1212\",\"type\":\"UnionRenderers\"},{\"attributes\":{\"source\":{\"id\":\"1182\",\"type\":\"ColumnDataSource\"}},\"id\":\"1186\",\"type\":\"CDSView\"},{\"attributes\":{},\"id\":\"1210\",\"type\":\"UnionRenderers\"},{\"attributes\":{},\"id\":\"1208\",\"type\":\"UnionRenderers\"},{\"attributes\":{\"fill_alpha\":{\"value\":0.1},\"fill_color\":{\"value\":\"#1f77b4\"},\"line_alpha\":{\"value\":0.1},\"line_color\":{\"value\":\"#1f77b4\"},\"x\":{\"field\":\"x\"},\"y\":{\"field\":\"y\"}},\"id\":\"1179\",\"type\":\"Scatter\"},{\"attributes\":{\"callback\":null,\"data\":{\"x\":{\"__ndarray__\":\"ZmZmZmZm/j8K16NwPQr/P65H4XoUrv8/KVyPwvUoAEB7FK5H4XoAQM3MzMzMzABAH4XrUbgeAUBxPQrXo3ABQMP1KFyPwgFAFa5H4XoUAkBnZmZmZmYCQLkehetRuAJAC9ejcD0KA0Bdj8L1KFwDQK9H4XoUrgNAAQAAAAAABEBTuB6F61EEQKVwPQrXowRA9yhcj8L1BEBJ4XoUrkcFQJuZmZmZmQVA7VG4HoXrBUA/CtejcD0GQJHC9ShcjwZA43oUrkfhBkA1MzMzMzMHQIfrUbgehQdA2aNwPQrXB0ArXI/C9SgIQH0UrkfheghAz8zMzMzMCEAhhetRuB4JQHM9CtejcAlA\",\"dtype\":\"float64\",\"shape\":[33]},\"y\":{\"__ndarray__\":\"QASYCCpzCT9Bu4OUM3QSP41YTXV3BhI/WVsRn975HT+oaaZljzQkP4hTwzzAJEQ/ky2uuHpQLz+qi5oEo74uPxqwc+ce90A/+f0EDztfUj86KRzvuvpVP+eO+LEwNYk/07AzX0nSYT+qFL75sMlZP7zpN0KSZW4/CRJ7063cZD8tSpz95BRKPz9IMRZCFFA/geKm1emvTz9gMLheIF1TPw13zQm32FA/B+BesUd9QT+tOITWAgdEPzslVFKFODo/hFhdxOJcOT/vjecP2ixAP1+5QROZc0A/4JZf9APKMz8im5iyHWw7P2SIRvOhnDI/wBux3jTpKz+7JE5vRwwwP7Z0kgw0yTQ/\",\"dtype\":\"float64\",\"shape\":[33]}},\"selected\":{\"id\":\"1211\",\"type\":\"Selection\"},\"selection_policy\":{\"id\":\"1212\",\"type\":\"UnionRenderers\"}},\"id\":\"1182\",\"type\":\"ColumnDataSource\"},{\"attributes\":{},\"id\":\"1209\",\"type\":\"Selection\"},{\"attributes\":{},\"id\":\"1164\",\"type\":\"HelpTool\"},{\"attributes\":{},\"id\":\"1204\",\"type\":\"BasicTickFormatter\"},{\"attributes\":{\"callback\":null,\"data\":{\"x\":{\"__ndarray__\":\"ZmZmZmZm/j8K16NwPQr/P65H4XoUrv8/KVyPwvUoAEB7FK5H4XoAQM3MzMzMzABAH4XrUbgeAUBxPQrXo3ABQMP1KFyPwgFAFa5H4XoUAkBnZmZmZmYCQLkehetRuAJAC9ejcD0KA0Bdj8L1KFwDQK9H4XoUrgNAAQAAAAAABEBTuB6F61EEQKVwPQrXowRA9yhcj8L1BEBJ4XoUrkcFQJuZmZmZmQVA7VG4HoXrBUA/CtejcD0GQJHC9ShcjwZA43oUrkfhBkA1MzMzMzMHQIfrUbgehQdA2aNwPQrXB0ArXI/C9SgIQH0UrkfheghAz8zMzMzMCEAhhetRuB4JQHM9CtejcAlA\",\"dtype\":\"float64\",\"shape\":[33]},\"y\":{\"__ndarray__\":\"aty7TAq0SD/S0ruqynBRP3UKsRfJL1I/ouIp54imWz8brXqlIkNiP2p9DJeA1ms/r4XVx5+IfD+gjHaMlbGBP7wQdwuX8I4/P9LRycKFiz88YC7z70yRP2UcnnItvZM/aCSQHV1gmj/9g8xVrZyYP5Cpc0bgbZc/w0DZ1VDikz8NQmm+owmTP+tJFm65LI4/emtxYxp6iz8dM6eROqiIP93Z9L4HCIU/ra2HJHRDhD+cxz00qGKAPy6N5J5Ol4A/HzqX8rYgfT+gCgCq4XV4P8/mzEoogHg/g5cW++/gdT9YN0lvnolzP60ClB/6c3Q/Lo/fhnMccT8wuOQxzkhwP/wUoNmmLW8/\",\"dtype\":\"float64\",\"shape\":[33]}},\"selected\":{\"id\":\"1207\",\"type\":\"Selection\"},\"selection_policy\":{\"id\":\"1208\",\"type\":\"UnionRenderers\"}},\"id\":\"1172\",\"type\":\"ColumnDataSource\"},{\"attributes\":{},\"id\":\"1147\",\"type\":\"LinearScale\"},{\"attributes\":{},\"id\":\"1202\",\"type\":\"BasicTickFormatter\"},{\"attributes\":{},\"id\":\"1207\",\"type\":\"Selection\"},{\"attributes\":{\"callback\":null,\"data\":{\"x\":{\"__ndarray__\":\"ZmZmZmZm/j8K16NwPQr/P65H4XoUrv8/KVyPwvUoAEB7FK5H4XoAQM3MzMzMzABAH4XrUbgeAUBxPQrXo3ABQMP1KFyPwgFAFa5H4XoUAkBnZmZmZmYCQLkehetRuAJAC9ejcD0KA0Bdj8L1KFwDQK9H4XoUrgNAAQAAAAAABEBTuB6F61EEQKVwPQrXowRA9yhcj8L1BEBJ4XoUrkcFQJuZmZmZmQVA7VG4HoXrBUA/CtejcD0GQJHC9ShcjwZA43oUrkfhBkA1MzMzMzMHQIfrUbgehQdA2aNwPQrXB0ArXI/C9SgIQH0UrkfheghAz8zMzMzMCEAhhetRuB4JQHM9CtejcAlA\",\"dtype\":\"float64\",\"shape\":[33]},\"y\":{\"__ndarray__\":\"+r8tGoMQLj/x71ZrykcxP+MzFM+kJDI/V4T0Y9A+RD+Vg9fBvKRIP28KBt/v7EY/g8B5TDx/Rj/xgviCB2puP74tv3ygiXw/j0Nlp7ynhz8fbGGyn0+ZP0hI1svY+Ig/Q6WgpzaDjj95oEFPLGGNP4/8mse40YU/d/TaFCJHfT+993UK7DV5P3EoJ9/+9XM/2rWviM4Ycj+DqXHGEDpqP/gMB/D/Amo/8s6OTburZj+EmRinbIBhP/VWydx6vFs/MCt5eyXPYT/tFIQ93D5aP/qzsUCrM1U/QlDr0Y0KVj/NAc6/R1NUP4DeGpW71FA/OTZgvj+8TT/3MDcqu3hQPzO1YI8X01A/\",\"dtype\":\"float64\",\"shape\":[33]}},\"selected\":{\"id\":\"1209\",\"type\":\"Selection\"},\"selection_policy\":{\"id\":\"1210\",\"type\":\"UnionRenderers\"}},\"id\":\"1177\",\"type\":\"ColumnDataSource\"},{\"attributes\":{},\"id\":\"1211\",\"type\":\"Selection\"},{\"attributes\":{\"data_source\":{\"id\":\"1182\",\"type\":\"ColumnDataSource\"},\"glyph\":{\"id\":\"1183\",\"type\":\"Scatter\"},\"hover_glyph\":null,\"muted_glyph\":null,\"nonselection_glyph\":{\"id\":\"1184\",\"type\":\"Scatter\"},\"selection_glyph\":null,\"view\":{\"id\":\"1186\",\"type\":\"CDSView\"}},\"id\":\"1185\",\"type\":\"GlyphRenderer\"},{\"attributes\":{},\"id\":\"1159\",\"type\":\"PanTool\"},{\"attributes\":{\"bottom_units\":\"screen\",\"fill_alpha\":{\"value\":0.5},\"fill_color\":{\"value\":\"lightgrey\"},\"left_units\":\"screen\",\"level\":\"overlay\",\"line_alpha\":{\"value\":1.0},\"line_color\":{\"value\":\"black\"},\"line_dash\":[4,4],\"line_width\":{\"value\":2},\"render_mode\":\"css\",\"right_units\":\"screen\",\"top_units\":\"screen\"},\"id\":\"1206\",\"type\":\"BoxAnnotation\"},{\"attributes\":{\"callback\":null},\"id\":\"1141\",\"type\":\"DataRange1d\"},{\"attributes\":{},\"id\":\"1150\",\"type\":\"BasicTicker\"},{\"attributes\":{\"source\":{\"id\":\"1172\",\"type\":\"ColumnDataSource\"}},\"id\":\"1176\",\"type\":\"CDSView\"},{\"attributes\":{\"active_drag\":\"auto\",\"active_inspect\":\"auto\",\"active_multi\":null,\"active_scroll\":\"auto\",\"active_tap\":\"auto\",\"tools\":[{\"id\":\"1159\",\"type\":\"PanTool\"},{\"id\":\"1160\",\"type\":\"WheelZoomTool\"},{\"id\":\"1161\",\"type\":\"BoxZoomTool\"},{\"id\":\"1162\",\"type\":\"SaveTool\"},{\"id\":\"1163\",\"type\":\"ResetTool\"},{\"id\":\"1164\",\"type\":\"HelpTool\"}]},\"id\":\"1165\",\"type\":\"Toolbar\"},{\"attributes\":{},\"id\":\"1145\",\"type\":\"LinearScale\"},{\"attributes\":{\"fill_color\":{\"value\":\"brown\"},\"line_color\":{\"value\":\"brown\"},\"x\":{\"field\":\"x\"},\"y\":{\"field\":\"y\"}},\"id\":\"1183\",\"type\":\"Scatter\"},{\"attributes\":{\"fill_alpha\":{\"value\":0.1},\"fill_color\":{\"value\":\"#1f77b4\"},\"line_alpha\":{\"value\":0.1},\"line_color\":{\"value\":\"#1f77b4\"},\"x\":{\"field\":\"x\"},\"y\":{\"field\":\"y\"}},\"id\":\"1174\",\"type\":\"Scatter\"},{\"attributes\":{\"callback\":null},\"id\":\"1143\",\"type\":\"DataRange1d\"},{\"attributes\":{\"data_source\":{\"id\":\"1177\",\"type\":\"ColumnDataSource\"},\"glyph\":{\"id\":\"1178\",\"type\":\"Scatter\"},\"hover_glyph\":null,\"muted_glyph\":null,\"nonselection_glyph\":{\"id\":\"1179\",\"type\":\"Scatter\"},\"selection_glyph\":null,\"view\":{\"id\":\"1181\",\"type\":\"CDSView\"}},\"id\":\"1180\",\"type\":\"GlyphRenderer\"},{\"attributes\":{\"text\":\"\"},\"id\":\"1201\",\"type\":\"Title\"},{\"attributes\":{},\"id\":\"1160\",\"type\":\"WheelZoomTool\"},{\"attributes\":{\"fill_color\":{\"value\":\"red\"},\"line_color\":{\"value\":\"red\"},\"x\":{\"field\":\"x\"},\"y\":{\"field\":\"y\"}},\"id\":\"1173\",\"type\":\"Scatter\"},{\"attributes\":{},\"id\":\"1155\",\"type\":\"BasicTicker\"},{\"attributes\":{},\"id\":\"1162\",\"type\":\"SaveTool\"},{\"attributes\":{\"fill_color\":{\"value\":\"blue\"},\"line_color\":{\"value\":\"blue\"},\"x\":{\"field\":\"x\"},\"y\":{\"field\":\"y\"}},\"id\":\"1178\",\"type\":\"Scatter\"},{\"attributes\":{\"overlay\":{\"id\":\"1206\",\"type\":\"BoxAnnotation\"}},\"id\":\"1161\",\"type\":\"BoxZoomTool\"},{\"attributes\":{\"source\":{\"id\":\"1177\",\"type\":\"ColumnDataSource\"}},\"id\":\"1181\",\"type\":\"CDSView\"},{\"attributes\":{\"fill_alpha\":{\"value\":0.1},\"fill_color\":{\"value\":\"#1f77b4\"},\"line_alpha\":{\"value\":0.1},\"line_color\":{\"value\":\"#1f77b4\"},\"x\":{\"field\":\"x\"},\"y\":{\"field\":\"y\"}},\"id\":\"1184\",\"type\":\"Scatter\"},{\"attributes\":{\"formatter\":{\"id\":\"1204\",\"type\":\"BasicTickFormatter\"},\"ticker\":{\"id\":\"1150\",\"type\":\"BasicTicker\"}},\"id\":\"1149\",\"type\":\"LinearAxis\"},{\"attributes\":{\"formatter\":{\"id\":\"1202\",\"type\":\"BasicTickFormatter\"},\"ticker\":{\"id\":\"1155\",\"type\":\"BasicTicker\"}},\"id\":\"1154\",\"type\":\"LinearAxis\"},{\"attributes\":{\"ticker\":{\"id\":\"1150\",\"type\":\"BasicTicker\"}},\"id\":\"1153\",\"type\":\"Grid\"},{\"attributes\":{\"dimension\":1,\"ticker\":{\"id\":\"1155\",\"type\":\"BasicTicker\"}},\"id\":\"1158\",\"type\":\"Grid\"},{\"attributes\":{\"data_source\":{\"id\":\"1172\",\"type\":\"ColumnDataSource\"},\"glyph\":{\"id\":\"1173\",\"type\":\"Scatter\"},\"hover_glyph\":null,\"muted_glyph\":null,\"nonselection_glyph\":{\"id\":\"1174\",\"type\":\"Scatter\"},\"selection_glyph\":null,\"view\":{\"id\":\"1176\",\"type\":\"CDSView\"}},\"id\":\"1175\",\"type\":\"GlyphRenderer\"}],\"root_ids\":[\"1140\"]},\"title\":\"Bokeh Application\",\"version\":\"1.3.1\"}};\n",
       "  var render_items = [{\"docid\":\"9d2f71d0-7344-43ae-aaa8-14c117ef28d4\",\"roots\":{\"1140\":\"967a6915-46b9-4a15-9adb-5912af1459ad\"}}];\n",
       "  root.Bokeh.embed.embed_items_notebook(docs_json, render_items);\n",
       "\n",
       "  }\n",
       "  if (root.Bokeh !== undefined) {\n",
       "    embed_document(root);\n",
       "  } else {\n",
       "    var attempts = 0;\n",
       "    var timer = setInterval(function(root) {\n",
       "      if (root.Bokeh !== undefined) {\n",
       "        embed_document(root);\n",
       "        clearInterval(timer);\n",
       "      }\n",
       "      attempts++;\n",
       "      if (attempts > 100) {\n",
       "        console.log(\"Bokeh: ERROR: Unable to run BokehJS code because BokehJS library is missing\");\n",
       "        clearInterval(timer);\n",
       "      }\n",
       "    }, 10, root)\n",
       "  }\n",
       "})(window);"
      ],
      "application/vnd.bokehjs_exec.v0+json": ""
     },
     "metadata": {
      "application/vnd.bokehjs_exec.v0+json": {
       "id": "1140"
      }
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Plot calculated stuff\n",
    "\n",
    "p = figure(plot_width=600, plot_height=400)\n",
    "p.scatter(temps, sus_16, color='red')\n",
    "p.scatter(temps, sus_32, color='blue')\n",
    "p.scatter(temps, sus_64, color='brown')\n",
    "\n",
    "show(p)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}