{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.optimize import fsolve\n",
    "from scipy.integrate import quad\n",
    "from scipy.special import eval_laguerre\n",
    "import itertools\n",
    "import copy\n",
    "from image_matrix_helper import compute_master_list, imshow_list, rgb_map, color_to_rgb, list_to_matrix\n",
    "import random\n",
    "import time\n",
    "\n",
    "nb_start = time.time()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Simulation for General Case of Ligand-Receptor Binding \n",
    "\n",
    "In this notebook, we simulate a thermal system of particles of various colors binding onto a grid. We have $R$ different types of particles and particle of type $j$ has $n_j$ copies in the system. Particles can exist both on and off the grid and particles of type $j$ have a binding affinity of $\\gamma_j$ to the grid. Each particle type also has a collection of \"correct\" locations on the grid. A particle of type $j$ binds to this correct location with an additional optimal binding affinity factor of $\\delta_j$ (i.e., its net affinity to such a site is $\\gamma_j \\delta_j$). Here we want to use simulations of this system to affirm analytical calculations of the average number of bound particles and the average number of correctly bound particles as functions of temperature. \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Numerical representations of analytical work"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Equations of Large $N$ approximation \n",
    "\n",
    "In the large $N$ limit, the order parameters for the system can be approximated as \n",
    "\n",
    "\\begin{equation}\n",
    "\\langle k \\rangle =   \\sum_{j=1}^r \\frac{n_j }{\\bar{z} \\gamma_j +1} \\left(\\bar{z} \\gamma_j +  \\frac{\\displaystyle L_{n_j-1} \\left( \\bar{\\phi}_{j}\\right)}{\\displaystyle  L_{n_j} \\left( \\bar{\\phi}_{j} \\right)}\\right) \\qquad \\langle m \\rangle  =   \\sum_{j=1}^r \\frac{n_j \\delta_j}{\\delta_j-1} \\frac{\\displaystyle L_{n_j-1} \\left( \\bar{\\phi}_{j}\\right)}{\\displaystyle  L_{n_j} \\left( \\bar{\\phi}_{j} \\right)},\n",
    "\\end{equation}\n",
    "where $\\bar{z}$ and $\\bar{x}$ are defined as \n",
    "\\begin{equation}\n",
    "\\bar{z} = \\sum_{j=1}^{R} \\frac{n_j}{\\bar{z} \\gamma_j +1} \\left(1- \\frac{\\displaystyle L_{n_j-1} \\left( \\bar{\\phi}_{j}\\right)}{\\displaystyle  L_{n_j} \\left( \\bar{\\phi}_{j} \\right)}\\right), \\qquad \\bar{x} = \\sum_{j=1}^{R} n_j\\left(1- \\frac{\\displaystyle L_{n_j-1} \\left( \\bar{\\phi}_{j}\\right)}{\\displaystyle  L_{n_j} \\left( \\bar{\\phi}_{j} \\right)}\\right)\n",
    "\\end{equation}\n",
    "with\n",
    "\\begin{equation}\n",
    "\\bar{\\phi}_{j}\\equiv \\frac{\\bar{x}}{1-\\delta_j}\\left(1+ \\frac{1}{\\bar{z} \\gamma_j}\\right).\n",
    "\\label{eq:phi_def}\n",
    "\\end{equation}\n",
    "and $L_n(x)$ the $n$th Laguerre polynomial. \n",
    "\n",
    "For these simulations we will take $$\\gamma_j = (\\beta E_V)^{3/2} e^{-\\beta E_j}, \\qquad \\delta_j = e^{\\beta \\Delta_j}$$ where $E_V$ is a volumetric Boltzmann factor associated with free particles (e.g., for a point-particle $E_V \\equiv h^2/2\\pi mV^{2/3}$), and $E_j$ is the binding energy for particles to the grid. We also take where $\\Delta_j$ is the binding energy advantage for particles binding to their correct locations in the grid. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Parameter function definitions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# helper function definitions\n",
    "gamma_func = lambda E0, Ev, T: 4*np.sqrt(2)*np.exp(E0/T)*(Ev/T)**(3/2)\n",
    "delta_func = lambda Del, T: np.exp(Del/T)\n",
    "phi_func = lambda x, z, gamma, delta: x*(1+ 1/(z*gamma))/(1-delta) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Equilibrium equations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def constr_func(X, T, E0s, Dels, Evs, Ns):\n",
    "    \n",
    "    \"\"\"\n",
    "    Equations constraint equations that determine zbar and xbar\n",
    "    \"\"\"\n",
    "    \n",
    "    x = X[0]\n",
    "    z = X[1]\n",
    "    \n",
    "    F = np.ones(2)\n",
    "    \n",
    "    R = len(Ns)\n",
    "    gammas_ = gamma_func(E0s,Evs, T)\n",
    "    deltas_ = delta_func(Dels, T)\n",
    "    phis_ = phi_func(x, z, gammas_, deltas_)\n",
    "    \n",
    "    F[0] = z- np.sum([Ns[j]/(z*gammas_[j]+1)*(1-eval_laguerre(Ns[j]-1, phis_[j])/eval_laguerre(Ns[j], phis_[j])) for j in range(R)])\n",
    "    F[1] = x- np.sum([Ns[j]*(1-eval_laguerre(Ns[j]-1, phis_[j])/eval_laguerre(Ns[j], phis_[j]) ) for j in range(R)])\n",
    "                     \n",
    "    return F\n",
    "    \n",
    "def m_avg(T, E0s, Dels, Evs, Ns):\n",
    "    \n",
    "    \"\"\"\n",
    "    Function that computes m_avg\n",
    "    \"\"\"    \n",
    "    \n",
    "    x, z =  fsolve(constr_func, x0 = (50,500), args = (T, E0s, Dels, Evs, Ns))\n",
    "    \n",
    "    R = len(Ns)\n",
    "    gammas_ = gamma_func(E0s,Evs, T)\n",
    "    deltas_ = delta_func(Dels, T)\n",
    "    phis_ = phi_func(x, z, gammas_, deltas_)\n",
    "    \n",
    "    return np.sum([Ns[j]*deltas_[j]/(deltas_[j]-1)*eval_laguerre(Ns[j]-1, phis_[j])/eval_laguerre(Ns[j], phis_[j]) for j in range(R)] )\n",
    "    \n",
    "def k_avg( T, E0s, Dels, Evs, Ns):\n",
    "    \n",
    "    \"\"\"\n",
    "    Function that computes k_avg\n",
    "    \"\"\"       \n",
    "    \n",
    "    x, z = fsolve(constr_func, x0 = (50,500), args = (T, E0s, Dels, Evs, Ns))\n",
    "\n",
    "    R = len(Ns)\n",
    "    gammas_ = gamma_func(E0s,Evs, T)\n",
    "    deltas_ = delta_func(Dels, T)\n",
    "    phis_ = phi_func(x, z, gammas_, deltas_)    \n",
    "    \n",
    "    return np.sum([Ns[j]/(z*gammas_[j]+1)*(z*gammas_[j] + eval_laguerre(Ns[j]-1, phis_[j])/eval_laguerre(Ns[j], phis_[j])) for j in range(R)])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Important temperatures\n",
    "\n",
    "The temperature $k_BT_{\\text{crit}}$ at which the system settles into its completely correct configuration is defined by \n",
    "\n",
    "\\begin{equation}\n",
    "1 = \\sum_{j=1}^{R}n_j e^{-\\beta_{\\text{crit}} \\Delta_j}\\left( 1+ (\\beta_{\\text{crit}} E_V)^{-3/2} e^{-\\beta_{\\text{crit}} E_j}\\right) + O(e^{-2\\beta_{\\text{crit}} \\Delta_j}). \n",
    "\\label{eq:master_therm_subs}\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# general thermal constraint\n",
    "def master_them_constr(T, E0s, Dels, Evs, Ns):\n",
    "    \n",
    "    F =  1-np.sum(Ns*delta_func(Dels, T)**(-1)*(1+gamma_func(E0s, Evs, T)**(-1)))\n",
    "    \n",
    "    return F\n",
    "\n",
    "# critical temperature\n",
    "kBTcrit_master = lambda E0s, Dels, Evs, Ns: fsolve(master_them_constr, x0 = 100.5, args = (E0s, Dels, Evs, Ns))[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Example analytical plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "<ipython-input-2-a582446d2d66>:4: RuntimeWarning: divide by zero encountered in true_divide\n",
      "  phi_func = lambda x, z, gamma, delta: x*(1+ 1/(z*gamma))/(1-delta)\n",
      "<ipython-input-3-08aba1dbaf99>:17: RuntimeWarning: invalid value encountered in double_scalars\n",
      "  F[0] = z- np.sum([Ns[j]/(z*gammas_[j]+1)*(1-eval_laguerre(Ns[j]-1, phis_[j])/eval_laguerre(Ns[j], phis_[j])) for j in range(R)])\n",
      "<ipython-input-3-08aba1dbaf99>:18: RuntimeWarning: invalid value encountered in double_scalars\n",
      "  F[1] = x- np.sum([Ns[j]*(1-eval_laguerre(Ns[j]-1, phis_[j])/eval_laguerre(Ns[j], phis_[j]) ) for j in range(R)])\n",
      "<ipython-input-2-a582446d2d66>:4: RuntimeWarning: invalid value encountered in multiply\n",
      "  phi_func = lambda x, z, gamma, delta: x*(1+ 1/(z*gamma))/(1-delta)\n",
      "/Users/mobolajiwilliams/opt/miniconda3/lib/python3.8/site-packages/scipy/optimize/minpack.py:175: RuntimeWarning: The iteration is not making good progress, as measured by the \n",
      "  improvement from the last ten iterations.\n",
      "  warnings.warn(msg, RuntimeWarning)\n"
     ]
    }
   ],
   "source": [
    "# temperature vals\n",
    "Tvals = np.linspace(0.1, 3.0, 50)\n",
    "\n",
    "# parameters for the integral\n",
    "np.random.seed(42)\n",
    "\n",
    "R = 50\n",
    "E0_bar, sigma_E = 10.0, 2.0\n",
    "Del_bar, sigma_D = 3.0, 1.0\n",
    "E0s = np.random.randn(R)*sigma_E+E0_bar \n",
    "Dels = np.random.randn(R)*sigma_D+Del_bar\n",
    "Nelems = np.random.randint(1,10,R)\n",
    "Evs = np.ones(R)*0.001\n",
    "\n",
    "# computing analytical values of k and m\n",
    "avg_k_approx_vals = [k_avg(T, E0s, Dels, Evs, Nelems) for T in Tvals]\n",
    "avg_m_approx_vals = [m_avg(T, E0s, Dels, Evs, Nelems) for T in Tvals]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAEYCAYAAABMVQ1yAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAA3iklEQVR4nO3dd3hUZfrG8e+TCqkQEggkgdAhIAQIHZEiLiIIilIsgA2sK6tY13Vd9be69oYgKqKgFLGBgqgoIL2GXqQECEgLEEogCcn7++MEDZEykJmcmZPnc11zzWTmMHO/jN6cnPIeMcaglFLK9/nZHUAppZR7aKErpZRDaKErpZRDaKErpZRDaKErpZRDaKErpZRD2FboXbt2NcAZt507d/7luUu9dRjTwXQY08Ft73cpN3eOx1tuThuT08bjxDE5bTxuGNM52VboBw4c+MtzeXl5NiTxHKeNB5w3JqeNB5w3JqeNBzw3pgCPvKsXeKr9U3ZHUEqpEuXYQr+yxpV2R1BKqRLl2J2iqXtSSd2TancMpZQqMY5dQx/6/VAAZg2aZWsOpZQqKY5dQ1dKqdLmgoUuIqNFZJ+IrDnH6yIib4nIZhFZJSJN3R9TKaXUhbiyhj4G6Hqe168GahfcBgMjih/rEhlj3ZRSqhS64DZ0Y8wcEUk8zyI9gU+MNbH6QhEpJyKVjTG/uytkYTtSZ3Ji7TSCsw8RlHOI4JyDfzwOyj1iZUZAjgOQ92wMRvw4FRBKbmAEuYER5ASG//H4RNlYjoTX5Gh4TY6HxoP4F4z77J9vTMHR/cb8eV/wXF6+IS/fcCrfkJefz/6Mg0RsySU0OIDyoYGUDwkiKjToj/uQIH9ExBN/TUqpUsgdO0XjgJ2Ffk4veO4vhS4ig7HW4omLiyMtLe2M1zMyMi74YVsW/0j7XR9xkHAOmnAOmggOUpkMU5dMwjCAYGgie/HDMMpE408+odkniZTjRHCcSNlNBL9RXo5TQ4798d4nTSBbTBV+M3Gk5tfmy7zLOUrIxf+NnGHfOV8J9BPCy/gTEexPeLA/EQWPy5X1p2JYIJXCgqz78EDCgvy8ovxd+Y58idPGA84bk9PGA8UbU2Ji4jlfK9GjXIwxo4BRACkpKeZswc4XFqD8Tf9m14mnQISyWP9yxF1inkzgSPYRAg9tJujgJgIPbaL6wU3UPbiRXkfn8++QzzlWrw9HG9/BqfI1//hzgiBCwU0QrMf+fkKAn1/BvZC+cwc1qieSlZ3HwawcDh7P4dDxnD8eH87KJfNEDoeO53L4RA4ZWblsPXSSjGM55OTln5E1NMifuPJlqV0pnHqVwqkTG0692HASyofg51eyRX+h78jXOG084LwxOW084JkxuaPQdwEJhX6OL3jOIyLDyhIZduHl5u+cD0CbhDYXWDIU4ioDl5/59O5U/BaNJGLNZ0SsHgO1ukCru6FmZ6u9XRAc4Eegvx+RIX5EhgRSPTrUpT+Xn284cDyb3YdPsvvwCXYfPsGuwyfYeTCL1emZfLfqz19+ygb6Uyc2nCYJ5WhRPYrmiVHEhAe79DlKKWdxR6FPAe4XkQlASyDTU9vPL8aTM58EinEcepVkuG4kdHkWlo6GJR/CuN5QMQn6fQpRNdyWtSg/P6FieBkqhpchOaHcX14/nn2K3/YdY+OeI2zYc5T1vx9hwpIdjJmfBkD16FCaJ5aneWIUl9eOITayjMeyKqW8xwULXUTGAx2AaBFJB/4NBAIYY0YC04BuwGYgC7jNU2FtEVYROjwO7f4Ba7+C7x+Hj7rBgCkQU8eWSKHBASQnlDuj7HPz8lmzK5MlaQdZvO0gM9buZdLSdAAaxUfSpX4lujSoRN1K4V6xLV4p5X6uHOXS/wKvG+A+tyXyVgHB0LgfxDaCT3rCR1fDgK8h9jK7kwEQ6O9Hk6rlaVK1PIPb1yQ/37Bp31Fmrt/Hj+v28uqPm3j1x00kRJXlqqRYul0WS9Oq5bXclXIQx5767zGVkuC26fDJtTCmO9zyJcQ3szvVX/j5CfViI6gXG8F9HWux78hJflq/jx/X7WHsgu18OHcbNWJC6ZOSwPVN4qgYoZtllPJ1eur/pYiuZZV6mUhrbX37fLsTXVDFiDLc1LIqH93WguVPd+GlGxpRITSIF6dvoPWLP3Pnx0uYsXYPuUWOrlFK+Q7HrqG/0fUNz35A+Wpw+/dWoY+9Hvp/BjU7efYz3SQsOIA+KQn0SUlg6/5jfL4snS+WpfPT+n3EhAczqE0it7SsRmRIoN1RlVIXwbFr6MmxySTHJnv2QyKqwKBpUKEmfNYX0pd69vM8oEZMGI91rcf8xzvx4cAU6leO4OUZG2n94kz+M3Ut6Yey7I6olHKRYwv9p60/8dPWnzz/QWExMHAqhETDdw9Dvm9usgjw96Nz/Up8cnsLpj94OV0bxDJ2wXaueHkWfx+/gjW7Mu2OqJS6AMcW+vNznuf5Oc+XzIeFREGX/8DvqZA6rmQ+04PqV47gtb7JzHm0I7e3TeTnDfvo/vZcBn+ylLSDJ+2Op5Q6B8cWeom77EZIaAkzn4WTzlibrVKuLP+8Jol5j3fioS51mL8lg9s/38Kwz1fqphilvJAWuruIwNUvwfEDMPslu9O4VWTZQP7euTZzHu3IDY0qMGXlbjq9Mptnp64j41i23fGUUgW00N2pSjI0vRUWjYT9G+1O43ZRoUHc2zqWWcM6cH3TOMbM30b7l37h3VmbyTnlm/sOlHISLXR36/Q0BIbC90849mIbVcqV5cXejfjhH1fQplY0L32/ka5vzmHubwfsjqZUqebYQn+v+3u81/29kv/gsBhr7pctMymbPrvkP78E1aoYxvsDUvjotubk5Rtu+XAR9322nN8zT9gdTalSybGFXje6LnWj69rz4S3ugui6RC1+BU45fxtzx7oVmTG0PQ93qcNP6/bS+dXZjJy9RTfDKFXCHFvoUzdOZerGqfZ8uH8gdH2BwKM7YcFwezKUsDKB/jzQuTY/PXQFbWtF8+L0DfR4e64ev65UCXJsob+64FVeXfCqfQFqdSYroSPMeQWO2D49fIlJiArh/QEpfDgwhcMncug5fB6v/rBR19aVKgGOLXRvcLD5w5CfC3NftztKietcvxI/DL2CXslxvP3zZq59R9fWlfI0LXQPOhWRAEk9YdUEyCl9J+JEhgTyap/GfDgwhYPHrbX113RtXSmP0UL3tGa3WWeOrvva7iS26Vy/Ej/+4wp6Nq7CWz9vpufweWzed8zuWEo5jha6p1VrAxVqw7IxdiexVWRIIK/1Teb9ASnsPXKSHm/PZdLSnRiHHquvlB0cW+hjrxvL2OvG2h3DmhKg2SDYuQj2rrM7je26JFVi+oOXk5xQjkcnr+LBCakcPZlrdyylHMGxhZ4QmUBCZILdMSyN+4N/ECz/2O4kXqFSRBnG3dmSYVfV4bvVv3PNW3NJ3XnY7lhK+TzHFvrENROZuGai3TEsoRWg/rWwcjzk6lmUAP5+wv2dajNpSCvy8g03jJjPe7O36CYYpYrBsYU+YukIRiwdYXeMPzUbVLBz9Bu7k3iVZtWimPb3y+mSVIkXpm9gyNhluglGqUvk2EL3OontIKpmqd85ejaRIYG8e3NT/tU9iZkb9tFz+Dx+23vU7lhK+Rwt9JJyeufojgWwb73dabyOiHBHu+p8emdLjpzIpdfweUxbXXrOsFXKHbTQS1LyTeAXCMt05+i5tKpRgW8fuJw6seHc++lyXpi2nlN5eiKSUq7QQi9JodFQv4fuHL2A2MgyTBjciltaVeW9OVsZMHoxB4/n2B1LKa/n2EKf3Gcyk/tMtjvGXzUbBCcPw7opdifxasEB/jzf6zJeubExS7cfopduV1fqghxb6NEh0USHRNsd46+qt4eoGrpz1EU3NItn4uBWZOXkcf2785m1cZ/dkZTyWo4t9DGpYxiTOsbuGH/1x87R+Y687qgnNKlanm/ub0t8VAi3j1nCmHnb9Hh1pc5CC90OjU/vHB1jdxKfEVeuLJPvbk3n+pV4Zuo6nvp6Dbm6s1SpMzi20L1aWAzU6warJkGenkTjqtDgAN67pRn3dKjJp4t2MHD0YjKz9O9PqdNcKnQR6SoiG0Vks4g8fpbXq4rILyKyQkRWiUg390d1mEb9IOsAbPnZ7iQ+xc9PeKxrPV69sTFL0w5x/Yh57DxY+uaaV+psLljoIuIPDAeuBpKA/iKSVGSxp4BJxpgmQD/gXXcHdZxaV0LZKFg5we4kPql3s3jG3tGC/UezuX7EfL0aklK4tobeAthsjNlqjMkBJgA9iyxjgIiCx5HAbvdFdKiAIGh4PWycBieP2J3GJ7WsUYEv7mlDkL8ffd5bwC96BIwq5QJcWCYO2Fno53SgZZFlngF+EJEHgFDgyrO9kYgMBgYDxMXFkZaWdsbrGRkZrmR2yYh21sRcRT+jJF1oPMExV1D51Acc+HU0x2r3KplQxeTO78gdAoG3eiTw+PQd3DFmCQ+1r0L3+uVd/vPeNh53cNqYnDYeKN6YEhMTz/maK4Xuiv7AGGPMqyLSGhgrIg2NMWcchmCMGQWMAkhJSTFnC3a+sL7ovOOpVg0W1iB690yiuwwtqUjF5m3fUSLwVa3q3Pvpcl6ZvZsc/xD+0aUOIuLan/ey8biD08bktPGAZ8bkyiaXXUDhK0XEFzxX2B3AJABjzAKgDGDrWT3vLnmXd5d4+aZ8EWjUF7b9CplF/0rVxQgLDuDDgSn0SYnnrZ8388jkVToHjCp1XCn0JUBtEakuIkFYOz2Lnre+A+gMICL1sQp9vzuDXqxJaycxae0kOyO4plEfwMBqH8jq5QL9/fhf70Y82Lk2k5elc/e4ZZzMzbM7llIl5oKFbow5BdwPzADWYx3NslZEnhWRawsWexi4S0RWAuOBQUZP5XNNVA2IbwErJ4L+lRWbiPCPLnV4tmcDZm7Yx4APF5N5Qo9VV6WDS9vQjTHTgGlFnnu60ON1QFv3RitFGveF7x6GPauhciO70zjCgNaJlA8J4qFJqfQbtZCPb29OxfAydsdSyqP0TFFv0OB6ayqAVV5yDVSH6NG4Ch8ObM72jOPcMGIB2zOO2x1JKY/SQvcGIVFQ+ypYPRnydZuvO7WvE8Nnd7Xi6Mlceo9YwNrdegKSci7HFvqsQbOYNWiW3TFc17gvHNsDW2fZncRxkhPK8fndrQn0F/qNWsjStIN2R1LKIxxb6D6n9t8gONKasEu5Xa2K4XxxTxtiwoK59cPF/PqbrQdhKeURji30V+a/wivzX7E7husCy0CDXrB+KuTotl5PqFKuLBOHtKZahRDuGLOUGWv32B1JKbdybKF/u+lbvt30rd0xLk7jfpB7HNb7WG4fEhMezMTBrUmqEsG9ny7nh02H7Y6klNs4ttB9UkIriKyqR7t4WGRIIOPubEmLxChe+HkX4xZutzuSUm6hhe5N/PysM0e3/gJHdXOAJ4UFB/DRbc1pVS2Mp75ew8jZW+yOpFSxaaF7m8b9wOTD6s/tTuJ4ZQL9ee6qqvRoXIUXp2/g9R836bVKlU9zbKGXDSxL2cCydse4eNG1Ia6ZNRWA8rgAf+GNvsnc2CyeN2f+xovfb9BSVz7LXdPnep3pN0+3O8Kla9QPpj8Ce9ZAbEO70ziev5/wv96NCA70473ZW8nOzeffPZJcnn5XKW/h2DV0n9awN/gFwCq9PF1J8fMTnuvZkDvaVWfM/DSe/GoN+fm6pq58i2ML/bnZz/Hc7OfsjnFpQitYJxqt+lynAihBIsJT19Tnvo41Gb94B8MmryRPS135EMcW+sxtM5m5babdMS6dTgVgCxHhkb/V4+Eudfhy+S4enLCCXL1QhvIRjt2G7vPqdIUykbByAtTqbHeaUueBzrUJDvTjv9M2kJuXz9v9mxIU4Nj1H+UQ+l+otwoItqbV3fAtZB+1O02pNLh9Tf7dI4kZa/dy76fLyD6lm7+Ud9NC92aN+0NuljW/i7LFbW2r81yvhvy0fh93j9VL2inv5thCrxBSgQohFeyOUTwJLaB8dWuzi7LNra2q8cL1lzFr037u+mSplrryWo4t9C/6fMEXfb6wO0bxiFhnjm6bA5npdqcp1fq3qMpLvRsxd/MBbh+zhKycU3ZHUuovHFvojtGoD2B0nnQvcGNKAq/1aczCrRkM+mgJx7O11JV3cWyhP/HTEzzx0xN2xyi+qBrWLIyrJoKekm6765rE80a/JizbfoiBoxdz9GSu3ZGU+oNjC31B+gIWpC+wO4Z7NO4L+zfA76l2J1HAtY2r8Hb/JqTuPMzA0Ys5oqWuvIRjC91RGlwH/kE6YZcX6XZZZd65qSmr0jMZ8KGWuvIOWui+oGx5qHu1NaVunhaHt+jaMJZ3b27K2t2Z3PrBIjKz9LtR9tJC9xWN+0PWAfjtB7uTqEKuahDLiJubsf73o9z84UIOZ+XYHUmVYo4t9PiIeOIj4u2O4T61ukB4ZVj2sd1JVBFXJlXivVubsWnvMW56fxGHjmupK3s4ttDHXT+OcdePszuG+/gHQPLNsPlHyNxldxpVRMd6FXl/QAqb9x+j//sLOailrmzg2EJ3pKa3WpenW+Ggf6gc5Io6MXw4MIVtB45z0/sLOXAs2+5IqpRxbKEP/X4oQ78fancM9yqfCDU6wIqxOk+6l7q8dgyjBzUnLeM4/UctZP9RLXVVchxb6Kl7Ukndk2p3DPdrOhAyd8KWX+xOos6hba1oRg9qTvqhE/R/fyH7jp60O5IqJVwqdBHpKiIbRWSziDx+jmX6iMg6EVkrIp+5N6b6Q71rIKQCLB9jdxJ1Hm1qRvPRbc3ZffgE/UYtZO8RLXXleRcsdBHxB4YDVwNJQH8RSSqyTG3gCaCtMaYBMNT9URVgzZPeuD9snA7H9tmdRp1HqxoVGHNbC/ZmnqTfqIXsydRSV57lyhp6C2CzMWarMSYHmAD0LLLMXcBwY8whAGOMNo0nNR0I+acg9VO7k6gLaFE9io9vb8H+o9n0HbWA3YdP2B1JOZgrhR4H7Cz0c3rBc4XVAeqIyDwRWSgiXd0V8FLVqVCHOhXq2B3DM2LqQNU2sPwTnbDLB6QkRvHJHS04eCyHvqMWkH4oy+5IyqHcdU3RAKA20AGIB+aIyGXGmMOFFxKRwcBggLi4ONLS0s54k4yMDDfFgScvexLgL59Rktw5nqJCq3Uj5ten2LPwc05WbuGxzynKk2OyQ0mNJwp4+ZqqDPs2jd7D5/L6tYlUiQjyyGfpd+T9ijOmxMTEc77mSqHvAhIK/Rxf8Fxh6cAiY0wusE1ENmEV/JLCCxljRgGjAFJSUszZgp0vrC/y2Hji7oQlLxO7awa07uOZzzgH/Y4u9XMgPq4KN3+wiIe/28n4u1qRGB3qoc9K9Mj72sVp4wHPjMmVTS5LgNoiUl1EgoB+wJQiy3yNtXaOiERjbYLZ6r6YF2/w1MEMnjrYzgieFVgWGvWF9VMg66DdaZSLGsZFMv6uVpzMzaPvqAVs3X/M7kjKQS5Y6MaYU8D9wAxgPTDJGLNWRJ4VkWsLFpsBZIjIOuAX4BFjjK2/J23K2MSmjE12RvC8pgMhL0evOepjkqpEMH5wK07lGfqOWsjmfUftjqQcwqXj0I0x04wxdYwxNY0x/1fw3NPGmCkFj40x5iFjTJIx5jJjjDZMSYhtCHHNYPnHunPUx9SLjWDC4FYYA/1GLWTjHi11VXyOPVO01Gg60Lqa0c5FdidRF6l2pXAmDmmFv5/Qb9QC1uzKtDuS8nFa6L6uYW8IjoRFI+1Ooi5BzZgwJg1pTUhQADe9v5CVOw/bHUn5MMcWenJsMsmxyXbH8LzgMGg2ANZNgcM7L7y88jrVKoQycUgryoUEcfMHi1iapju51aVxbKG/0fUN3uj6ht0xSkaLIdb94vfszaEuWXz5ECYOaUXF8GAGjF7Mgi3OO/ZaeZ5jC71UKZcASdfCsk8gWw+D81WVI8syYUgr4suXZdBHi5mzab/dkZSPcWyh3/LlLdzy5S12xyg5re6D7Eyd38XHVQwvw4TBrakZE8adHy/lp3V77Y6kfIhjCz39SDrpR9LtjlFyEppDfHNYOEIvfuHjokKDGH9XK+pXDufuccuYunK33ZGUj3BsoZdKre6FQ9tg0/d2J1HFFBkSyLg7W9K0WnkenLCCSUt1h7e6MC10J6l/LUQmwIJ37U6i3CC8TCAf39aCdrVjeHTyKsbM22Z3JOXltNCdxD8AWgyG7XPh95V2p1FuUDbIn/cHNONvDSrxzNR1DP9ls92RlBdzbKG3jm9N6/jWdscoeU0HQGCorqU7SHCAP8Nvakqv5Cq8PGMjL8/YgNGpHtRZuGs+dK/zwpUv2B3BHmXLQZNbYOlo6PIfCI+1O5FygwB/P17rk0zZoACG/7KF49l5PN09CT8/sTua8iKOXUMv1VrdbV2ibvH7didRbuTnJ/z3uobcdXl1xsxPY9jnK8nNy7c7lvIiji303pN603tSb7tj2COqBtTtZq2l5+o1LJ1ERHiyW32GXVWHL1fs4p5xyzmZq4epKotjCz0jK4OMrFJ8+nTre+HEQVg53u4kys1EhPs71ea5ng2YuWEvgz5azNGTuXbHUl7AsYVe6lVrC1Wawtw3IE//Z3eiW1sn8kbfZJamHeLmDxZx8HiO3ZGUzbTQnUoEOjwOh7dD6md2p1Ee0jM5jlEDmrFxz1FuHDmf3zN1E1tppoXuZLWvsq5oNOcVOKVrb07VqV4lPrm9BfuOZHPDiAVsP5RtdyRlE8cWeufqnelcvbPdMewlAh2ehMwdkDrO7jTKg1rWqMD4wa3IPpXPA99sY/mOQ3ZHUjZwbKH/64p/8a8r/mV3DPvV6mxN2jXnVTila25O1jAuki/vaUNYkD83vb+QXzbsszuSKmGOLXRVQAQ6PAFH0mHFWLvTKA+rWiGEd3pVp1bFMO78ZCmTl5WiGUeVcwv96k+v5upPr7Y7hneo2QkSWsKvr0HuSbvTKA+LCglgwuDWtK5RgWGfr2Tk7C06VUAp4dhCP5F7ghN6Uo3lj7X0XbD8E7vTqBIQFhzA6EHN6dG4Ci9O38Bz364nP19L3ekcW+iqiBodoGobmKtr6aVFUIAfb/ZN5ra2iYyet437x+tZpU6nhV5aiEDHJ+Do77BsjN1pVAnx8xOe7p7EU9fUZ9rqPXoCksNpoZcm1dtDtXYFa+m6Oaq0EBHuvLwG797clNW7Muk9Yj7bM47bHUt5gGMLvXud7nSv093uGN6n4xNwbK81cZcqVbpdVpnP7mzJoawcrn93Piv0WHXHcWyhD2szjGFthtkdw/sktrPW1Oe8Aif0f+jSJiUxii/vaUNocAD931/IjLV77I6k3Mixha7O46r/g5OHYdaLdidRNqgRE8aX97ahbmwEd49bxge/btXDGh3CsYXeYUwHOozpYHcM71S5ETQdaF0AY996u9MoG0SHBTPhrlb8LSmW579bz5NfrdGLZTiAYwtdXUCnf0FwGEx/DHTtrFQqG+TPuzc35d4ONRm/eAcDRy8mM0unWvZlLhW6iHQVkY0isllEHj/Pcr1FxIhIivsiKo8IrQAd/wnbZsOGb+1Oo2zi5yc82rUer97YmCVpB7nu3XlsO6BHwPiqCxa6iPgDw4GrgSSgv4gknWW5cOBBYJG7QyoPSbkDYurDjCf1MMZSrnezeD67qxWHT+TSa/g8FmwpxVf78mGurKG3ADYbY7YaY3KACUDPsyz3HPA/QE9D9BX+AXD1/+DwDpj/jt1plM2aJ0bx9b1tiQkP5tYPF/HZoh12R1IXyZVCjwN2Fvo5veC5P4hIUyDBGPOdG7MVS58GfejToI/dMbxfjSug/rXWyUaZu+xOo2xWtUIIX97bhra1onnyq9X886vV5JzSnaW+IqC4byAifsBrwCAXlh0MDAaIi4sjLS3tjNczMtz3a163mG4Af/mMkuTO8XhSQNLdVNk0g6yvH+bAFec/lNFXxuQqp40H3DOmpzvE8EFIPp8u2sHqHQd4pksCUSHFrotLot/RmRITE8/5mivf0C4godDP8QXPnRYONARmiQhALDBFRK41xiwt/EbGmFHAKICUlBRztmDnC3sxsnKzAAgJDHHL+10qd43HsxJh/4OEzXmJsA5/h2ptzr+0T4zJdU4bD7hnTC/UqE6rert47ItV3PfNdkbdmsJl8ZHFD3cJ9DtyjSubXJYAtUWkuogEAf2AKadfNMZkGmOijTGJxphEYCHwlzIvad0+7Ua3T7vZGcG3tPsHRMTD9Ech75TdaZSX6Jkcx+S72+Anwg0j5/P1Ct0s580uWOjGmFPA/cAMYD0wyRizVkSeFZFrPR1QlZCgEOj6X9izGua+bnca5UUaxkXyzf1taZxQjqETU3nu23V6EpKXcmmjmDFmGjCtyHNPn2PZDsWPpWyR1BMa9obZ/4M6f7POKFUK68zST+9syfPfruPDudtYlX6Yd25qSqWIMnZHU4XomaLqTN1egZAo+PoeOKXzZqs/Bfr78Z+eDXmzXzJrdh3hmrfmsnCr83ZY+jItdHWmkCjo8SbsXWOtqStVRM/kOL65vy0RZQO4+YNFes1SL+LYQh+UPIhByYPsjuGb6l4NyTdbx6anL7M7jfJCdSqFM+X+dnRtEMuL0zcwZOwyjpzUeWDspoWuzq7rCxBeBb6+W6cFUGcVFhzAOzc14V/dk/h5wz66vzWXVemH7Y5Vqjm20A9kHeBA1gG7Y/iuMpHQ8204sAl+ft7uNMpLiQh3tKvOxCGtOJWXT+8R83V+dRs5ttBvmHQDN0y6we4Yvq1mJ2sCrwXDYft8u9MoL9asWhTTHrycjnUr8vx367nj46V6MWobOLbQlZt0eRbKV7OOesk+anca5cXKhQTx3q3NeLZnA+b+doCr35yjR8GUMC10dX7BYdBrhDUj4zf36cUw1HmJCANaJ/LVfW0IDQrgpvcX8tqPmzilJyKVCC10dWHV2sCVz8C6b4hY+4ndaZQPaFAlkqkPtOO6JvG8NfM3bhi5QC+cUQK00JVr2vwdknpSftkbsHW23WmUDwgNDuDVPo15u38Tth04Trc3f+XTRdt1h6kHObbQ70m5h3tS7rE7hnOIQM/h5EYmwuTb4PDOC/4RpQB6NK7CjKHtSUkszz+/WsMdHy9l31G9Do4nOLbQ+zbsS9+Gfe2O4SzB4ezr+Lo1JcCkWyFX/6dUromNLMPHt7XgmR5JzNt8gK5v/MqMtXvsjuU4ji30nZk72Zmpa5HudioyEa4bCbtXwPRH7I6jfIifnzCobXW+faAdlSPLMGTsMv4xMZXDWXp4o7s4ttBv/epWbv3qVrtjOFP97nD5MFj+CSwbY3ca5WNqVwrnq3vb8vfOtZm6cjdXvjaH79fo2ro7OLbQlYd1fBJqdoZpj8CORXanUT4mKMCPh7rU4Zv721IxPJi7xy3jvs+Wk3Es2+5oPk0LXV0aP3/o/QFExsP4vrBvg92JlA9qUMW6eMawq+rww9o9dHl9DlNX7tYjYS6RFrq6dCFRcMuX4B8EY6/TI1/UJQn09+P+TrX57u+Xk1C+LA+MX8Fdnywl/VCW3dF8jha6Kp6o6lap5xy3Sv24nuqtLk2dSuF8cU8bnuxWj3mbM+jy2hxGzdnCqTxdW3eVYwv94dYP83Drh+2OUTrENoSbJkDmTvj0Bsg+Znci5aMC/P0Y3L4mPz7Unra1KvDfaRsY8uUWlu84ZHc0n+DYQu9Rtwc96vawO0bpUa0N3DgGfl8JE2/Ry9epYokvH8L7A1IYeUszMk/m0XvEfP751WoyT+hFNM7HsYW+8cBGNh7YaHeM0qXu1XDt27D1F/hqCOTn2Z1I+TARoWvDWD7pW4vb2lRn/OIddHplFhMW7yAvXzfDnI1jC33It0MY8u0Qu2OUPk1utqbcXfslfDtUS10VW0iQP0/3SGLK/e2oHh3K41+uptfweSzbftDuaF7HsYWubNT2QWj/iHXi0VdDIE9/TVbF1zAuks/vbs2b/ZLZfzSb3iMWMHTCCvZk6hQUpwXYHUA5VKenICgUfnrGOgLmho8gsIzdqZSPExF6JsdxZf1KjJi1hVG/buWHdXu5r2Mt7mhXnTKB/nZHtJWuoSvPafcP6PYKbJwGn/XRo1+U24QGBzDsb3X56R9X0K5WNC/P2EjHV2bx+dKdpXr7uha68qwWd0GvkZD2q3Wc+onDdidSDlK1QgijBqQwYXArKoYH88jkVVzz1q/M3rS/VJ5t6thCf6r9UzzV/im7YyiA5P5w48fWDI1jusOx/XYnUg7TqkYFvr6vLe/c1ISsnDwGjl7MrR8uZs2uTLujlSjHFvqVNa7kyhpX2h1DnZZ0rXXyUcZmGH0V7N9kdyLlMCJC90ZV+PGh9jzdPYm1uzPp/vZc7vt0Ob/tLR0XOHdsoafuSSV1T6rdMVRhta6EgVMg+yh80Bl++9HuRMqBggP8ub1ddWY90pEHOtVi1sZ9XPXGHIZOWOH465o6ttCHfj+Uod8PtTuGKiqhBdz1C5SvBp/eCPPehFK4rVN5XmTZQB6+qi6/PtaJwe1r8P3aPVz52mwe+XwlOw86c+Ivxxa68mLlEuD2GZDUE3582jpWXS9npzwkKjSIJ66uz6+PdmJg60S+Wbmbjq/M4tHJK9m631lHXrlU6CLSVUQ2ishmEXn8LK8/JCLrRGSViMwUkWruj6ocJSjUmvul41OwaiKM6QZHfrc7lXKwmPBgnu6RxJxHOnJzy6p8k7qbzq/N5r7PlrNu9xG747nFBQtdRPyB4cDVQBLQX0SSiiy2AkgxxjQCJgMvuTuociARuOIR6DvOukDGqA6wdZbdqZTDxUaW4T89GzL3sU7cfUVNZm/cT7e3fuX2MUt8fjoBV9bQWwCbjTFbjTE5wASgZ+EFjDG/GGNOb5RaCMS7N6ZytPo94M4foUwEfNILfviXztaoPC4mPJjHutZj3mOdeLhLHVbsOETvEQu4ceR8vl+zxydPUHLl1P84oPClaNKBludZ/g5genFCucN/O//X7gjqYlRqAINnw4wnYf5bsG2OdYm76Np2J1MOFxkSyAOda3PH5dUZv3gno+du4+5xy6gaFcLtbRO5MSWB0GDfmCXFrSlF5BYgBbjiHK8PBgYDxMXFkZaWdsbrGRnuu9pNFaoA/OUzSpI7x+MtPD6my4YSEtGICvOfQUa242CLxzhW+zpr84wH6Hfk/UpyPJ3jhSv6VGfutiN8viqDZ6au45UZG+iRVJ7rGlagYligWz6nOGNKTEw852uuFPouIKHQz/EFz51BRK4E/glcYYw566W7jTGjgFEAKSkp5mzBzhf2YszfOR+ANglt3PJ+l8pd4/EmHh9T4u3Q5Gr4agjR8/9D9KHlcM3rEBbjmY/T78jrlfR4atWAQZ1h2fZDjJ67jYkrf2fSqoNcWb8it7ZKpE3NCvj5FW8lwxNjcqXQlwC1RaQ6VpH3A24qvICINAHeA7oaY/a5PeUleHLmkwDMGjTL3iDq0kRUhlu/hgXvwMxnIW0udHkOmtzisbV1pYpqVq08zaqVZ+fBLD5bvIOJS3YyY+1eakSHcnOratzQNJ7IEPestbvDBXeKGmNOAfcDM4D1wCRjzFoReVZEri1Y7GUgDPhcRFJFZIrHEqvSw88P2v4d7pkHMfVhyv0w5hqdNkCVuISoEB7rWo8FT3Tijb7JlA8N4rlv19HyhZ945POVLNt+0CsmA3NpG7oxZhowrchzTxd6rJOmKM+JqQuDvoPUcdYRMCPawOUPQbuHdI51VaKCA/zp1SSOXk3iWLf7CGMXbmdK6i4+X5ZOjZhQ+qQkcH3TOCqG2/PfpZ4pqnyDnx80HQD3L4UG18Hs/8HItrD5J7uTqVIqqUoEL1x/GYv/eSUv3dCICqFBvDh9A61f+Jk7P17KD2v3kHMqv0Qz+caxOEqdFhYDvd+Hxv3gu4dhXG+o0RGueg5iL7M7nSqFQoMD6JOSQJ+UBLbsP8bnS9P5Ynk6P63fS7mQQK65rDK9msTRrGr5Yu9IvRDHFvobXd+wO4LypFqd4b5FsORDmPMSjLwcGve3Ln0XGWd3OlVK1YwJ4/Gr6zHsqjrM+W0/X6/YzRfL0/l00Q7iypWlV5Mq9EqOw1O7UR1b6MmxyXZHUJ4WEAyt77UuoPHra7BoJKz9ElrdC+2GQplIuxOqUirA349O9SrRqV4ljmef4od1e/hqxW5GzNrC8F+2cH+bWIbZdNiiT/ppq7VtVS9yUQqULW9tcml+J/z8HMx9DZaOhlb3QMsh1utK2SQ0OIDrmsRzXZN49h/N5ttVu6kd6pmpLRy7U/T5Oc/z/Jzn7Y6hSlL5atZ0AYNnQ7W2MOsFeP0y+Ok/cNxZZ08q3xQTHsxtbasTXy7YI+/v2EJXpViVZOj/Gdw919rWPvd1eKMhzPgnHN1rdzqlPEYLXTlX7GXQ52Nr52n9HrDwXavYv7obdqfanU4pt9NCV84XUxeuH2Udw950IKybAqOugNFdCUn7AfJO2Z1QKbfQQlelR4WacM0r8PB6+Nt/4chuKs56BN5sbB0lo5tjlI9z7FEu73V/z+4IyluViYTW90HLu9k7dyyVtn0BM/8DPz8Pdf4GTW6F2l3A33smXVLKFY4t9LrRde2OoLydnz8nqnaA9oOsCb9Sx0HqeNg4DcIqWWejJt8CMXXsTqqUSxy7yWXqxqlM3TjV7hjKV8TUgS7PwkProN94iEuB+e/A8ObWtU7nvw2Z6XanVOq8HLuG/uqCVwHoUbeHzUmUT/EPhHrdrNvRvbB6Eqz5An54yrpVbQ0Ne0NSL49dcEOpS+XYQleq2MIrQZsHrFvGFmtagdVfwLRhMP1R6+SluldD3W4QVd3utEppoSvlkgo1of0j1m3vOqvcN3xnXdR6xpNQMamg3K+BKk2s6X6VKmFa6EpdrEpJ1q3TU3BwK2ycDhumWWek/voqhFaEmh2hZmfrPqyi3YlVKaGFrlRxRNWwDoFsfR9kHYTffoDffrQuvLFqorVM7GUF5d4JElpAYFl7MyvHcmyhj71urN0RVGkTEmUd6ti4H+Tnw56VsHkmbPnFutj1vDfAPwjimkFiO2sbfEJLCAqxO7lyCMcWekJkgt0RVGnm52dtS6/SBNoPg+yjsH0+pM2F7fOsM1PnvAx+gRDX1Fpzj29h3YfH2p1e+SjHFvrENdavu30b9rU5iVJAcLh1Fmqdv1k/Zx+FHYsg7Ver6Be9Zx3rDhBZFRKaWwUf1wxiG+pmGuUSxxb6iKUjAC105aWCw6H2ldYN4FQ2/L4Kdi6C9MWwfYF1/DuA+EPF+ta0wJWToUpTa6eslrwqwrGFrpRPCQi21soTmv/5XOYu2L3Cuv2eah1Ns2Kc9Zr4QYVaUKkhVGpg7Xit1AAi9HqqpZkWulLeKjLOutXvbv1sDGTutOZy37Ma9q6FXcusY+JPKxNJbER1iLsMYuoV3OpCZDyIZ684r+ynha6UrxCBclWtW9K1fz5/MhP2rbdKfv8G2JFasDZf6EivwFDr5KgKtQrdFzzWa646hha6Ur6uTCRUbWXdgD1paSQmJsLxA7B/o1XyBzZZ0xfsXgHrvgaTX+jPl4PyidYtqvqfj8tVhYh4CAgq6RGpS+TYQp/cZ7LdEZSyV2i0dUtse+bzp3LgUBoc3AIZm63Hh9Jg7xprOoP83EILi3UYZWQ8RCZAuQTrPqJKwS0OQqJ1qgMv4dhCjw6JtjuCUt4pIMiaLvhs87zn58GR3XBoGxzeaW2zP7wTMndYa/frpxYpfKxj6SMqQ3gVa0KzsFjrH4HwWGte+fBYazqEkCjw8y+ZMZZSji30MaljABiUPMjWHEr5FD9/ay283DlOzMvPh+P74ehuq/iP7IYju/58vHeddWZs9pG//lnxs9bmQ2OsqYdDK1q/QYREFTwfDSEVrMchUdamIH/HVpRHOPZvSwtdKQ/w87PWwsMrWWfBnkvOcTi2F47use6P7Yfj++DYPusfhGP74OA2yMqAnGPnfp/gSOKCwiG84p8lX7ZcoftI63GZSCgTYd0HFzwuhZcQdGyhK6VsFBRqTVwWVePCy+aetIo964B1fzwDThyEE4cg6yDZB3YS6JdjTX6WscU6qufk4TN37J5NQFmr2IPDrVtQGASf/jnM+jkorOBx6J8/B4VAYEiRx6HWPDxefuinS4UuIl2BNwF/4ANjzItFXg8GPgGaARlAX2NMmnujKqUcKbDMn8fcn8WBtDTCEhPPfNIYa/qEk4etgj9x2NrMc/JIwX2mdcs+AtnHrGWzj8Lh7QXPHbV+i8jLcT2n+FvlHljWugWFWvcBZa0xFH5c+D4gGALKWPeB1s8B+RWBxAt94kW7YKGLiD8wHOgCpANLRGSKMWZdocXuAA4ZY2qJSD/gf4Cec6+U8gyRgk0sEcV7n1M51iafnOPWffYxyD0OOVmQm2U9n5tlvZZ7ouCWdeZ9Tpb1m0XuSTh14sz7vOyzfmyZ1k9Bo8uLl/0sXFlDbwFsNsZsBRCRCUBPoHCh9wSeKXg8GXhHRMQYY9yYVSml3CsgCAKirO3znmCMNU/PqZOFbtlkZZzwyMe5UuhxwM5CP6cDLc+1jDHmlIhkAhWAA+4IeSmm3TzNro9WSimLSMHmmDJnPJ1/PM0jH1eiO0VFZDAwGCAuLo60tLQzXs/IyCjJOB7ntPGA88bktPGA88bktPFA8caUWHR/QiGuFPouoPBBqfEFz51tmXQRCQAisXaOnsEYMwoYBZCSkmLOFux8YS/Gu0veBeDe5ve65f0ulbvG402cNianjQecNyanjQc8MyZXztddAtQWkeoiEgT0A6YUWWYKMLDg8Q3Az3ZvP5+0dhKT1k6yM4JSSpWoC66hF2wTvx+YgXXY4mhjzFoReRZYaoyZAnwIjBWRzcBBrNJXSilVglzahm6MmQZMK/Lc04UenwRudG80pZRSF0OnSFNKKYfQQldKKYcQu/Zdish+YHuRp6Ox8dh1D3DaeMB5Y3LaeMB5Y3LaeKB4YzpgjOl6thdsK/SzEZGlxpgUu3O4i9PGA84bk9PGA84bk9PGA54bk25yUUoph9BCV0oph/C2Qh9ldwA3c9p4wHljctp4wHljctp4wENj8qpt6EoppS6dt62hK6WUukRa6Eop5RC2FLqIdBWRjSKyWUQeP8vrwSIyseD1RSKSaENMl7kwnkEisl9EUgtud9qR01UiMlpE9onImnO8LiLyVsF4V4lI05LOeDFcGE8HEcks9P08fbblvImIJIjILyKyTkTWisiDZ1nGZ74nF8fjU9+TiJQRkcUisrJgTP85yzLu7TpjTInesCb42gLUAIKAlUBSkWXuBUYWPO4HTCzpnG4ezyDgHbuzXsSY2gNNgTXneL0bMB0QoBWwyO7MxRxPB+Bbu3Ne5JgqA00LHocDm87y353PfE8ujsenvqeCv/ewgseBwCKgVZFl3Np1dqyh/3FJO2NMDnD6knaF9QQ+Lng8Gegs4rWX23ZlPD7FGDMHa9bMc+kJfGIsC4FyIlK5ZNJdPBfG43OMMb8bY5YXPD4KrMe6clhhPvM9uTgen1Lw936s4MfAglvRo1Dc2nV2FPrZLmlX9Is745J2wOlL2nkjV8YD0Lvg197JIpJwltd9iatj9iWtC341ni4iDewOczEKfk1vgrUGWJhPfk/nGQ/42PckIv4ikgrsA340xpzzO3JH1+lO0ZIxFUg0xjQCfuTPf5GVd1gOVDPGNAbeBr62N47rRCQM+AIYaow5Ynee4rrAeHzuezLG5BljkrGu9NZCRBp68vPsKPSLuaQd57uknZe44HiMMRnGmOyCHz8AmpVQNk9x5Tv0GcaYI6d/NTbW3P+BIhJtc6wLEpFArPL71Bjz5VkW8anv6ULj8dXvCcAYcxj4BSg6qZZbu86OQvfJS9qdxwXHU2S75bVY2wd92RRgQMFRFK2ATGPM73aHulQiEnt6u6WItMD6/8JbVyAA6wgWrCuFrTfGvHaOxXzme3JlPL72PYlIjIiUK3hcFugCbCiymFu7zqUrFrmTcdgl7Vwcz99F5FrgFNZ4BtkW2AUiMh7riIJoEUkH/o21QwdjzEisq1d1AzYDWcBt9iR1jQvjuQG4R0ROASeAfl68AnFaW+BWYHXBNlqAJ4Gq4JPfkyvj8bXvqTLwsYj4Y/3jM8kY860nu05P/VdKKYfQnaJKKeUQWuhKKeUQWuhKKeUQWuhKKeUQWuhKKeUQWuhKKeUQWuhKKeUQWuhKKeUQWuhKKeUQWuhKKeUQWuhKKeUQWuhKKeUQWuhKKeUQWujK8URkqYistjuHUp6mha4creAqMA2B1GK8R5iI5ImIcfEW5bYBKHURSvwCF0qVsCQgmGIUOtb/JwOLPHcP0AYYBuwt9Hy2MeZgMT5LqUumF7hQjiYiA7Auyt3ZGPOzG993GdY/FuEFV2tXyna6yUU5XZOC+9TTT4hIORH5SkROisjgi33DgosZNwRWaZkrb6KbXJTTJQM7T28GEZFmwOeAAG2NMcsu4T0bAEHACneFVModdA1dOV0yBWvnInIPMA9YDzS7xDKHP9f6lxc3nFLupIWuHEtEEoFywG8i8hnwDvB/QPdi7rhsWnCva+jKq+gmF+Vkp9ekHwDyga7GmB+LLiQifkAm1gqOKbjfBNxnjJl3jvc9Beix7cqr6Bq6crLThf4+1qGLl51jubpAGJBgjAkDygMbgNeLLlhQ/o2B9caYk25PrFQx6Bq6crJk4IAx5j4RKQu8LCJpxpgviyzXDNh9ejOMMSZbRFYAdc7ynrWxyl+3nyuvo2voysma8Od27iHALGCciLQoslwzYBWAWBoBdwGvnOU9dfu58lpa6MqRRKQCEE9B8RpjcoHeQBowVUSqF1q8GdBRRA4Dx7GOihlvjPnsLG+tR7gor6WFrpzqdPH+sSZtjDkMXFPw4zQRKV+wTbwJcJMxppwxJgRoBTwqIlec430NsNJjyZW6RHrqvyrVRKQe1nHpCcaY9ELPHwAeMcZ8ZFs4pS6S7hRVpd3pHaLpACISCjyEdSboD3YGU+piaaGr0q4ZECsix7COVT8CzMeaFmCXrcmUuki6yUUppRxCd4oqpZRDaKErpZRDaKErpZRDaKErpZRDaKErpZRDaKErpZRDaKErpZRDaKErpZRD/D93ug1ObYaq0wAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "## plotting order parameters\n",
    "ax = plt.subplot(111)\n",
    "ax.plot(Tvals, avg_k_approx_vals/np.sum(Nelems), label = r'$\\langle k \\rangle$ (large $N$)')\n",
    "ax.plot(Tvals, avg_m_approx_vals/np.sum(Nelems), label = r'$\\langle m \\rangle$ (large $N$)')\n",
    "ax.set_xlabel(r'$k_BT$', fontsize = 18, labelpad = 10.5)\n",
    "ax.grid(alpha = 0.5)\n",
    "ax.axvline(x = kBTcrit_master(E0s, Dels, Evs, Nelems), color = 'g', linestyle = '--')\n",
    "\n",
    "# Hide the right and top spines\n",
    "ax.spines['right'].set_visible(False)\n",
    "ax.spines['top'].set_visible(False)\n",
    "\n",
    "# Only show ticks on the left and bottom spines\n",
    "ax.yaxis.set_ticks_position('left')\n",
    "ax.xaxis.set_ticks_position('bottom')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Metropolis Hastings simulation code "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Microstate transitions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "## dissociation operator\n",
    "def trans_dissoc(free_objs, bound_objs):\n",
    "        \n",
    "        # indices of non-empty\n",
    "        indxs = [i for i, x in enumerate(bound_objs) if x != \"-\"]\n",
    "        \n",
    "        # random choice for bound object\n",
    "        random_indx = random.choice(indxs)\n",
    "\n",
    "        ## new state vector \n",
    "        free_objs_new = copy.deepcopy(free_objs)\n",
    "        bound_objs_new = copy.deepcopy(bound_objs)\n",
    "\n",
    "        # putting empty slot\n",
    "        bound_objs_new[random_indx] = '-'\n",
    "        # appending previously bound object to free objects\n",
    "        free_objs_new.append(bound_objs[random_indx])\n",
    "\n",
    "        return free_objs_new, bound_objs_new\n",
    "\n",
    "\n",
    "## association operator\n",
    "def trans_assoc(free_objs, bound_objs):\n",
    "    \n",
    "        # random element to associate\n",
    "        elem = random.choice(free_objs)\n",
    "\n",
    "        # indices of empty spaces\n",
    "        indxs = [i for i, x in enumerate(bound_objs) if x == \"-\"]  \n",
    "\n",
    "        # random choice for empty space\n",
    "        random_indx = random.choice(indxs)\n",
    "\n",
    "        ## new state vector \n",
    "        free_objs_new = copy.deepcopy(free_objs)\n",
    "        bound_objs_new = copy.deepcopy(bound_objs)    \n",
    "\n",
    "        ## state\n",
    "        free_objs_new.remove(elem)\n",
    "        bound_objs_new[random_indx] = elem      \n",
    "\n",
    "        return free_objs_new, bound_objs_new\n",
    "\n",
    "## permutation operator\n",
    "def trans_perm(free_objs, bound_objs):\n",
    "    \n",
    "    Ncomp = len(bound_objs)\n",
    "    i1 = int(random.choice(range(Ncomp)))\n",
    "    i2 = int(random.choice(range(Ncomp)))\n",
    "\n",
    "    ## new omega vector \n",
    "    bound_objs_new = copy.deepcopy(bound_objs)\n",
    "    bound_objs_new[i2] = bound_objs[i1]\n",
    "    bound_objs_new[i1] = bound_objs[i2]\n",
    "    \n",
    "    return free_objs, bound_objs_new"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Logarithm of Botlzmann factor\n",
    "\n",
    "The logarithm of the Botlzmann factor for a microstate (i.e., the temperature normalized negative energy of the microstate) is defined as \n",
    "\n",
    "\\begin{equation}\n",
    "\\beta E(\\boldsymbol{k}, \\boldsymbol{m}) = \\sum_{i=1}^R(m_i \\ln \\delta_i + k_i \\ln \\gamma_i).\n",
    "\\label{eq:sim_en}\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def log_boltz(free_objs, bound_objs, mstr_vec, deltas, gammas, name_key):\n",
    "    \n",
    "    elem_set = list(set(mstr_vec))\n",
    "    count_dict = dict()\n",
    "    for elem in elem_set:\n",
    "        count_dict[elem] = bound_objs.count(elem)\n",
    "\n",
    "    bind_log_factor = 0\n",
    "    for elem in elem_set:\n",
    "        key = name_key[elem]\n",
    "        bind_log_factor += count_dict[elem]*np.log(gammas[key])\n",
    "        \n",
    "        \n",
    "    corr_log_factor = 0\n",
    "    for j in range(len(bound_objs)):\n",
    "        if bound_objs[j] == mstr_vec[j]:\n",
    "            elem = bound_objs[j]\n",
    "            key = name_key[elem]\n",
    "            corr_log_factor+=np.log(deltas[key])        \n",
    "    \n",
    "    return bind_log_factor+corr_log_factor    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Function to count the number of correctly bound particles"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "def m_calc(bound_objs, mstr_vec):\n",
    "    \n",
    "    num = 0\n",
    "    for k in range(len(mstr_vec)):\n",
    "        if mstr_vec[k] == bound_objs[k]:\n",
    "            num += 1\n",
    "    return num"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Checking logarithm of Boltzmann factor definition"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Energy for original ordering: -3.664820641553842\n",
      "Checking  energy value: -3.664820641553841\n",
      "Number of correctly placed elements: 4\n",
      "Number of bound elements: 4\n",
      "-----\n",
      "Energy after permutation and associaation values: -4.01912719922027\n",
      "Checking  energy value: -4.01912719922027\n",
      "Number of correctly placed elements: 3\n",
      "Number of bound elements: 5\n"
     ]
    }
   ],
   "source": [
    "# defining name key\n",
    "name_key0 = dict()\n",
    "key_list = ['A', 'B', 'C', 'D', 'E', 'F', 'G', 'H', 'I', 'J', ]\n",
    "for j in range(len(key_list)):\n",
    "    name_key0[key_list[j]] = j\n",
    "    \n",
    "# random energies    \n",
    "np.random.seed(2)\n",
    "q1 = np.random.rand(10)\n",
    "q2 = np.random.rand(10)\n",
    "\n",
    "# sample master list\n",
    "sample_master = ['A', 'B', 'C', 'D', 'E', 'F', 'G', 'H', 'I', 'J', ]\n",
    "bound_init_0 = ['-', '-', '-', '-', '-', '-', 'G', 'H', 'I', 'J', ]\n",
    "free_init_0 = ['A', 'B', 'C', 'D', 'E', 'F' ]\n",
    "\n",
    "print('Energy for original ordering:', log_boltz( free_init_0,bound_init_0, sample_master, deltas = np.exp(-q1), gammas = np.exp(-q2), name_key =  name_key0))\n",
    "e1 = -np.sum([q1[k] for k in range(len(sample_master)) if sample_master[k] == bound_init_0[k]])\n",
    "e2 = -np.sum([q2[k] for k in range(len(sample_master)) if sample_master[k] in bound_init_0])\n",
    "print('Checking  energy value:', e1+e2)\n",
    "print('Number of correctly placed elements:', m_calc(bound_init_0, sample_master))\n",
    "print('Number of bound elements:', np.sum([1 for elem in bound_init_0 if elem!='-']))\n",
    "print('-----')\n",
    "random.seed(1)\n",
    "free_init_0, perm_bound = trans_perm(free_init_0, bound_init_0)\n",
    "free_init_new, bound_init_new = trans_assoc(free_init_0, perm_bound)\n",
    "print('Energy after permutation and associaation values:', log_boltz(free_init_new, bound_init_new, sample_master, deltas = np.exp(-q1), gammas = np.exp(-q2), name_key =  name_key0))\n",
    "e1 = -np.sum([q1[k] for k in range(len(sample_master)) if sample_master[k] == bound_init_new[k]])\n",
    "e2 = -np.sum([q2[k] for k in range(len(sample_master)) if sample_master[k] in bound_init_new])\n",
    "print('Checking  energy value:', e1+e2)\n",
    "print('Number of correctly placed elements:', m_calc(bound_init_new, sample_master))\n",
    "print('Number of bound elements:', np.sum([1 for elem in bound_init_new if elem!='-']))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Metropolis Hastings algorithm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Metropolis Monte Carlo Algorithm \n",
    "\n",
    "## loads uniform random sampling \n",
    "runif = np.random.rand\n",
    "\n",
    "def met_assembly_grid(Niter, free_objs, bound_objs, mstr_vec, deltas, gammas, name_key):\n",
    "    '''\n",
    "    #################################################################\n",
    "    # function to sample using Metropolis \n",
    "    #  \n",
    "    # n_iter:  number of iterations\n",
    "    # initial_state: initial state for the start position for our chain\n",
    "    # gamma: energy cost for incorrect component\n",
    "    # temp: temperature \n",
    "    ##################################################################\n",
    "    '''\n",
    "    \n",
    "    # Initialize state values\n",
    "    free_objs_vals = [0]*(Niter+1)\n",
    "    bound_objs_vals = [0]*(Niter+1)\n",
    "    \n",
    "    # Set initial values\n",
    "    free_objs_vals[0] = free_objs[:]\n",
    "    bound_objs_vals[0] = bound_objs[:]\n",
    "        \n",
    "    # Initialize acceptance counts\n",
    "    # We can use this to tune our number of steps\n",
    "    accepted = 0\n",
    "    \n",
    "    # debugging code\n",
    "    debug_assoc, debug_dissoc, debug_perm = 0, 0, 0\n",
    "    \n",
    "    for i in range(Niter):\n",
    "    \n",
    "        # get current monomer and dimer states\n",
    "        current_free_objs = copy.deepcopy(free_objs_vals[i])\n",
    "        current_bound_objs = copy.deepcopy(bound_objs_vals[i])\n",
    "        \n",
    "        N_free = len(current_free_objs)\n",
    "        N_bound = len(current_bound_objs)-len(current_free_objs)\n",
    "        \n",
    "        u_trans = runif()\n",
    "        \n",
    "        if u_trans < 1/3: #first type of transition; monomer association \n",
    "            \n",
    "            if N_free < 1:\n",
    "                \n",
    "                log_alpha = np.log(1e-15) \n",
    "            \n",
    "            else:\n",
    "                # proposed new monomer and dimer states\n",
    "                new_free_objs, new_bound_objs = trans_assoc(current_free_objs, current_bound_objs)\n",
    "\n",
    "                # transition elements\n",
    "                log_init = log_boltz(current_free_objs, current_bound_objs, mstr_vec, deltas, gammas, name_key)\n",
    "                log_final = log_boltz(new_free_objs, new_bound_objs, mstr_vec, deltas, gammas, name_key)\n",
    "\n",
    "                # weight\n",
    "                num = N_free*N_free\n",
    "                den = N_bound+1\n",
    "\n",
    "                # Log-acceptance rate\n",
    "                log_alpha = log_final-log_init+np.log(num/den) \n",
    "            \n",
    "        elif 1/3 <= u_trans < 2/3: #second type of transition; bound monomer dissociation\n",
    "            \n",
    "            if N_bound <1:\n",
    "                \n",
    "                log_alpha = np.log(1e-15) \n",
    "                \n",
    "            else: \n",
    "                \n",
    "                # proposed new monomer and dimer states\n",
    "                new_free_objs, new_bound_objs = trans_dissoc(current_free_objs, current_bound_objs)\n",
    "\n",
    "                # transition elements\n",
    "                log_init = log_boltz(current_free_objs, current_bound_objs, mstr_vec, deltas, gammas, name_key)\n",
    "                log_final = log_boltz(new_free_objs, new_bound_objs, mstr_vec, deltas, gammas, name_key)\n",
    "\n",
    "                # weight\n",
    "                num = N_bound \n",
    "                den = (N_free+1)*(N_free+1)   \n",
    "\n",
    "                # Log-acceptance rate\n",
    "                log_alpha = log_final-log_init+np.log(num/den) \n",
    "            \n",
    "        elif 2/3 <= u_trans: #third type of transition; switching bounded elements\n",
    "            \n",
    "            if N_bound <2:\n",
    "                \n",
    "                log_alpha = np.log(1e-15) \n",
    "            \n",
    "            else:\n",
    "                \n",
    "                # proposed new monomer and dimer states\n",
    "                new_free_objs, new_bound_objs = trans_perm(current_free_objs, current_bound_objs)\n",
    "\n",
    "                # transition elements\n",
    "                log_init = log_boltz(current_free_objs, current_bound_objs, mstr_vec, deltas, gammas, name_key)\n",
    "                log_final = log_boltz(new_free_objs, new_bound_objs, mstr_vec, deltas, gammas, name_key)\n",
    "\n",
    "                # Log-acceptance rate\n",
    "                log_alpha = log_final-log_init\n",
    "            \n",
    "        # Sample a uniform random variate\n",
    "        u = runif()\n",
    "\n",
    "        # Test proposed value\n",
    "        if np.log(u) < log_alpha:\n",
    "            # Accept\n",
    "            free_objs_vals[i+1] = new_free_objs\n",
    "            bound_objs_vals[i+1] = new_bound_objs\n",
    "            \n",
    "            #log_current_prob = log_proposed_prob\n",
    "            accepted += 1\n",
    "        else:\n",
    "            # Stay put\n",
    "            free_objs_vals[i+1] = free_objs_vals[i]\n",
    "            bound_objs_vals[i+1] = bound_objs_vals[i]\n",
    "\n",
    "    # return our samples and the number of accepted steps\n",
    "    return free_objs_vals, bound_objs_vals, accepted"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Computing microstate averages from simiulations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "def avg_k(bound_objs_vals, Nmc):\n",
    "    \n",
    "    \"\"\"\n",
    "    Microstate average of number of bound objects\n",
    "    We only consider microstates near the end of theh chain to ensure\n",
    "    that the system has equilibrated\n",
    "    \"\"\"\n",
    "    \n",
    "    length = int(Nmc/50)\n",
    "    \n",
    "    ls = [0]*length\n",
    "    ls = np.array(ls)\n",
    "    for k in range(length):\n",
    "        ls[k] = len(bound_objs_vals[Nmc-length+k]) - bound_objs_vals[Nmc-length+k].count('-')\n",
    "    \n",
    "    return(np.mean(ls))\n",
    "\n",
    "# average number of correctly bound objects\n",
    "def avg_m(bound_objs_vals, mstr_vec, Nmc):\n",
    "    \n",
    "    \"\"\"\n",
    "    Microstate average of number of correctly bound objects\n",
    "    We only consider microstates near the end of theh chain to ensure\n",
    "    that the system has equilibrated\n",
    "    \"\"\" \n",
    "    \n",
    "    length = int(Nmc/50)\n",
    "    \n",
    "    ls = [0]*length\n",
    "    ls = np.array(ls)\n",
    "    for k in range(length):\n",
    "        ls[k] =  np.sum([1 for j in range(len(mstr_vec)) if bound_objs_vals[Nmc-length+k][j]==mstr_vec[j]])\n",
    "    \n",
    "    return(np.mean(ls))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Image grid for completely correct configuration"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAATwAAAEGCAYAAAD45CnNAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAATzUlEQVR4nO3deZhkVXnH8e8LyL64EiEaBxQN4oY7j8rySMQYlS3u4ckQQMHduCTiEx1DGARDXKKQAApxjVEIiGDigoCoiFs0iAIug1GUVZYBWYSTP95TTHGnqrtq6Oqe6vP9PE8/PXPr3lOn7vK75557u06UUpCkFqyz0BWQpPli4ElqhoEnqRkGnqRmGHiSmmHgSWrGegtdgblwzjnnzPmzNbvsskuv7Kkq9+KVc1fmIzZlzsuc5nK549C5K3Td5XNf5uTLjbktdP7ZwpPUDANPUjMMPEnNMPAkNcPAk9QMA09SMww8Sc0w8CQ1w8CT1AwDT1IzDDxJzTDwJDXDwJPUDANPUjMMPEnNMPAkNcPAk9QMA09SMww8Sc2IUuZ8OIiFsCg+hLSWc0yLhRQRz4uI404//fSFroqkKbAoWnjTOGrZNIzYNa2ji01DuRMZCQ0ctWwWU93Ck6RxGHiSmmHgSWqGgSepGQaepGYYeJKaYeBJaoaBJ6kZBp6kZhh4kpph4ElqhoEnqRkGnqRmGHiSmmHgSWqGgSepGQaepGYYeJKaYeBJasaiGNMCRy2T5oNjWiwkRy2TNI5F0cK7eOXct/CmaQSsSZU7TXWdj3LndCSwOrrYVNQVHLVMkqaNgSepGQaepGYYeJKaYeBJaoaBJ6kZBp6kZhh4kpph4ElqhoEnqRkGnqRmGHiSmmHgSWqGgSepGQaepGYYeJKaYeBJaoaBJ6kZBp6kZiyKMS1w1DJpPjimxUJy1DJJ41gULTxHLXPUsmkrd5rq2leuLTxJmhYGnqRmGHiSmmHgSWqGgSepGQaepGYYeJKaYeBJaoaBJ6kZBp6kZhh4kpph4ElqhoEnqRkGnqRmGHiSmmHgSWqGgSepGQaepGYYeJKasSjGtMBRy6T54JgWC8lRyySNY1G08By1bFW5vOPQuSv0ncuB6VsHkyo33rfPnJVZXnfKnJc5D+XawpOkaWHgSWqGgSepGQaepGYYeJKaYeBJaoaBJ6kZBp6kZhh4kpph4ElqhoEnqRkGnqRmGHiSmmHgSWqGgSepGQaepGYYeJKaYeBJasai+Ip3HMRHmg9+xftCchAfSeNYFC28aRzEZ04H24G7BtyZxCA+k6qrg/g4iM98m+oWniSNw8CT1AwDT1IzDDxJzTDwJDXDwJPUDANPUjMMPEnNMPAkNcPAk9QMA09SMww8Sc0w8CQ1w8CT1AwDT1IzDDxJzTDwJDXDwJPUDANPUjMWxZgWOGqZNB8c02IhOWqZpHEsihbeJEctm4rRxSZV7jTVta/caRgNrVfmB/ed20bTq07OQ2FC68AWniRNCwNPUjMMPEnNMPAkNcPAk9QMA09SMww8Sc0w8CQ1w8CT1AwDT1IzDDxJzTDwJDXDwJPUDANPUjMMPEnNMPAkNcPAk9QMA09SMww8Sc1YFGNa4Khl0nxwTIuF5KhlksaxKFp4Ex217IdzPLLWDsunp9xpqmtfuY5a5qhlw0x1C0+SxmHgSWqGgSepGQaepGYYeJKaYeBJaoaBJ6kZBp6kZhh4kpph4ElqhoEnqRkGnqRmGHiSmmHgSWqGgSepGQaepGYYeJKaYeBJaoaBJ6kZi2JMCxy1TJoPjmmxkBy1TNI4FkULz1HLmMiIXa6D5KhljlomSVPHwJPUDANPUjMMPEnNMPAkNcPAk9QMA09SMww8Sc0w8CQ1w8CT1AwDT1IzDDxJzTDwJDXDwJPUDANPUjMMPEnNMPAkNcPAk9SMRfEV7ziIjzQf/Ir3heQgPpLGsShaeHHs3LfwyiF3ld1sudNU12krd5rq2leuLTxJmhYGnqRmGHiSmmHgSWqGgSepGQaepGYYeJKaYeBJaoaBJ6kZBp6kZhh4kpph4ElqhoEnqRkGnqRmGHiSmmHgSWqGgSepGQaepGYYeJKasSjGtMBRy6T54JgWC6k3allEPJ/cGHP6Y7nTVddpK3ea6lrLfR5TbqoDr5Ryeinl5aWUSY3TOKkNPE3lTlNdp63caarrJMudN1MdePNgUkE6TeVOU12nrdxpqusky503i6UPT5JmZQtPUjMMPEnNMPAkNcPAk9QMA09SMww8Sc0w8CQ1w8CT1AwDT1Izmgy8iDg7IlasBfVYGhElInZd6LpoYUTEbhFxfkTcWPeFpRGxa+/fC12/ubQ2HHcjB15EbBwRr4+Ir0bEtRFxe0RcERFn1o203iQrujaJiL0iYtlC12M2EfHMiPiPiPi/iLi1HlTfiYjDI+JBC12/eyoilkTEsoh43Bose5+IeHtEfCsirouI2yLilxFxckTsExExgSqvVgfgFGAT4I3AfsC5k37fSaoZsXSh6zHMSCEVEQ8DzgAeDnwJOAK4GtgS2B04EXgk8JbJVHOtsxfwl8Cyha3GYBGxDvCvwIHAZcAngEuB9YEnAK8CDiK33zRbArwDWAH8z6gLRcSTgdPIz/9Z4OPADcAfAs8BTibX0TFzWdkBngTcGziglHJKX/1WABsBt0/4/Sfh9eT2OGnAa88iv2pqwcwaeBGxEfA5YFtg3/4NUx0ZEU8iN57WDsvIsPsksLSUclv/ixHxRjIo5kRE3AtYt5Ryy4DX1gU2KKXcPFfvd09ExAPJb/3YENillHJeZ5bDImIP4D7zUJ0H1t/X9k8spdwJrLYu59NM23RNdffDBVFKmfEHeA35jcLvmm3eznJ7AV8DbgJW1n/vOWC+FcDZwGPJ1uNK4ErgaDKQNwT+EfgVuROcC2zfKWNprePu5MF+GXAr8APgxQPe82xgxYDp2wEfBX4N3Fbr9m5gk86yZcDP0r55tgKOBX5Ry7kcOA7Ycki9d63/37v+/6Ah6/SHwE+o33IzZJ4tgd/Vum84xvZaUj/7FXXd/RRYDmzcmW9ZreMOwD8BvwTuAHbtbIe/q2Xc3ls35Nn9EOA7wM11W38F2G1Infat6/u6Ov/FwPvJlmrvvbo/Z8/yOY+u8x085v58IPDdum6vB74APH3AfIVs3ewEnEPu/9cAJwCbdvb71epfX9u1u0/V6fcDPlzLWwmcBezIgP25V48B9bvbPjfbNq2vv4hsCf+i7htXA6cCjxnwnoN+lsxy3O0MfLGu19/V9XzAsOMW2Jo8mf+27hf/DTx8lO04yiXtn9ffx40wLwAR8Urgg8CPgb+vk5cCp0bEK0op3bIeVD/wp4DPkE3fvwZ+T26EjYB3AfcH3lTL2b7kmbDfkWR/SO9SZH/gkxGxYSnlpFnq/ARyB7qOvBz8FRnCrwWeFhG7lFJuBw4n+z6fQfa59Hy9lvNHwDfIg/JD5EH/MPJA3y0inlhKuX5INU4HfgP8FXB8p35PJbsN3tY7Kob4M/Ik8ZEy4tk5Ih4CXABsQa67S8mD7q31sz+zlPL7zmIfJ3fOXoD8mgxNyBPUvepnuIEMKshAfQm5jU8ENgBeBnwxIvYppXy2r06HA4cCFwHvqeU/lAzBt5MnvuV1nuOAr9ZFr5jl4+5LnoT+bZb57hIRR5LdNRfU99sMeDnwlYjYs5RyZmeRx5FXRSeS3Qm7AgcAd9blIC/9/rT+fznwo1nqsAHZIHgcGagXAI+p064duuB4Bm1TgFeTIXscuX8+tNb7axHx+FLKpXW+/chtdTV5nPRcNewN67co/2ct92jgRuDFwAkRsW0p5W2dRTYht/355LbYBngdcFpEPKqUcseMn3CEM9s1wPVjnAnvQ559fgJs3jd9c/LgvxG494Az3Qs65XyH3EFOo69FQwZQAfYYcNa6DNiib/oWddq1wEbdM0Xn/b5PBvRmnem9Vld/C+4k6tl4wOc/jWyhPqgz/YlkgC+b5Wy7vE57ZGf54+vyW8+y/ns76z5jbLOP12We05n+7jr9gL5py+q0s4H1hrQeLmb1lmFvPb68M3094NvAz3vbGXhynfcsOq1U6teNlxlaQjN8zs3q/D8YY908ou6H5wHr903fmjw5riAv/fpbOXcCT+mUcwbZ2u1v5a22/Yd9LuCVddrbOvP2pnf35zVp4a22TevrmwyYtj3Z2jumM30FQ1rZdI47YF3y+Lyuf78mGwtfI1uZ23WWL8BbOuW+mU4mDPsZ5S7t5mRIjepPyBR+fynlht7E+u/3A5uSlzz9flVK+XRn2nnkzv3PpX6qqncm327Aex9b+lpP9d//QobwrsMqHBGPJs+WnwA2iIj7935qPW4iW50ziogtgOeSzf9bOuWsIE8Cs5VzPDVk+srdhLys+Hwp5fJZlt+8/r5hxrlWlb0O8Hzge2X1lsoR5MG794BF31tWb/X1HFtW77P7C3I/OrWzXu5NtmyXsGqbvqz+fmvptFJLNfsnG2isdVPtSe6HR5W+Pqi6HU4EHkJeVvb7Rinlm51pZ5HhvmScCvd5HhkA7+tMP4G8FJwLA7dpKeUmyEEtImLzut2uIk9sT7kH7/cE4I+AD/fv13U9H0VeSe3ZWeZOMkf6nVV/D8qEuxkl8G4gz4yj2qb+/uGA13rTtu1M//mAeX875LXe9PsNWGbQZcFFQ96z3/b19zvJDdn/cyUZ4H8ww/I9jyDX6QEDyrmqvj5jOaWUn5OXKfvVjmOAF5Lb4IQR6tA7mEfdZg8gT0Krba9SyrXkZc2gdXfJDGUOem37WqcrWH29LKvz9NbNdmTof3/W2o9n3HUDa7Y//2zAvNfU34P221HrcXkpZWX/xBoOg46fNTFwm0bEjhHxOfKEdT2rttujuWc3d9Zk3V7ePQkyxrodpQ/vQmDnej09aEPOhZmuu4e9Npe3t3tlHQ3815B5fjtk+qByPsbwPqLfjVDOccCnyZbXyWSA/oa8LJrNhfX3jmTfyKTMdNd10GtBHiQvnWG5C/v+3evwnjOllBsj4jLgjyNio1LKKNtiTcy0Py/oYxnMfMyvtt1qn/S55MniMLJVdxO5bd5Lnizn0z1at6ME3snkXZQDyU7C2fRCcQfgy53XHtmZZ65tT/ahjfuevU7XO0opXxrhfYYdiD+pr60/YjnD9PoBD4iIC4GnAUfOcAnZ7wzybvZ+EXF4KeXWWea/ijxz79B9oT4YuxVjPOM2g0vJ5zjP77ZSBriE7NB/LNk5P8yaBOIpwBvIDvZRbsT1788/7bw26f253wpg94jYtH/91auAbch+sH7XAvcdUM5MVzqD7E2G2vNLKV/pfyEi7kf24/UbZ5v0r9uuiazbUS5pTyBT/U0R0b2eBvIOZ70zC3m39SbgNRGxWd88m5GPuKys80zCIbUfrfeeWwAHkzvDOTMs9z2ydXFwRKy2Q0TEehHRv/OsrNPvtkOVUq4BzgT2qXdVu+VERDxgtg9R8m7wScAerHpe7kOzLVeXvZK82bCEvNO1/oB6bB4R76nz30n2oe0YEc/uzPq35D4yFy3Fj9Syjhj0YkT0X+p/ov5ePqT+vTN578AfdGAPcxQZ8kdFxE5D6vKsiHhx/e9nyYP4zX1dDETEVuRTAJeR+8+knU528r+uM/0g8uZc1yXAThGxcW9CPYHtP+b79lpUd2s9RcRBrHqOsN9KRt8e3yUfddm/Ph/ZK/terLoR0W3A3COztvBKKTdHxHPJlsOpEfEFMrCuIft/diMPzKPq/NdFxFvIx1K+GREn1aKWko9nvKIMfyzjnrq6vueJ9f/7k52iBw7oRL9LKaVExH5k5+cPIuLDZB/CxrXO+5CPaJxUFzmfvFV/TET07r59s/a/HULe6Dg3Ij5CHgzrkGfWPckDf9kIn+V4cqO/BDinrLr1P4plZMvsQODpEfHvZOtzffKxhheQj2a8oc5/KHmz6dSIOKbOuzN5o+RcxniEY5hSymfqdnl1RDyefGzjavKRpJ3I9bxtnfeC+ijI3wDfjYhPkZf025CPST2ZPIldRLZOXxkRN9dpV5ZSzmKIUspv6v58GnBeRJzKqku2rYFnA08ntyOllIsj4t3kYynn1rr0HkvZFHhZme1RiLlxAvAK4B/qXz71Hkt5Ibm9usfyB8iulbMi4qPkzaGDyIAeFFTDfJ681P1oRHyA7Np5GvkXKT8d8L7nk1cmh5F96ncCp/dufPQrpdwREa8mT6jfiojjyO35IuCpwPIx9/vZzXYbt+/W78bkAXJe/dC3kx3QZ5CXB+t25t+bfDbtpvrzdWCvAeWuYMBtbFbdKl/Smb6kTl824Fb77uSNh94Dkv8LvHS22+N90x9C3tVdQQbCNeTjMUcAD+6bbx3yWbPeA5rdRwjuT7ayLiEvL6+rdXkffY+bMOSxhL7Xv1xf32/U7dRZfneyL/CX9fPcWD/PYcBWnXm3IZ+Tu7LO+zNmfvB4yYD3m/Hz1Hn2I++031DXzQryMvNFA+Z9Cfl4wo11H/ox2W/U/3jIc8iWwi2M8OBx33L3JVvP3yY74m+r6+kz5OVbd/6DyJPXLbXuXwSeMWC+cR4HGbi+GP7g8QPIk+61dX2cRZ7Avg1cNOA938yqh/B/RD7fOageQ7dpfX1n8ri/kdyXzwAexeDHu7Yku8GuJcPurnIHzV+n71LXZ2+f+B4zPHg8YPoSOpkw7GdRjEsb+cfKJ5JP7J+9sLWZOxFxJtn62bpMroNdU6z+6d7V5BVGt0tCHU1+PdQ0qJctewAfM+wEd/1de9fB5OXqpPrFF5VmvtJpWkTEU8i7za8lL7OOXtgaaS1yfERsSHYP3Uq2/l9K9uGN/KefLbOFt/Y5hPwD8c3JDvEVC1sdrUW+ADyY/GKG95J9fSeQX2Iwzl9DNWtR9OFJ0ihs4UlqhoEnqRkGnqRmGHiSmmHgSWqGgSepGf8PTXzYZW6IINgAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# defining master_list\n",
    "master_list =compute_master_list()\n",
    "# testing plot\n",
    "imshow_list(master_list, title = 'Completely Correct Configuration');\n",
    "# defining Nelems\n",
    "Nelems = np.zeros(8)\n",
    "key_list = list(rgb_map.keys())[:-1]\n",
    "name_key_ = dict()\n",
    "for j in range(len(key_list)):\n",
    "    name_key_[key_list[j]] = j\n",
    "    Nelems[j] = master_list.count(key_list[j])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 9.,  9., 10.,  5.,  7.,  6.,  3., 51.])"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# displaying copy-number counts of the various elements\n",
    "Nelems"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Simulating system"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Temperature Run: 1 ; Current Time: 62.85 secs\n",
      "Temperature Run: 2 ; Current Time: 125.33 secs\n",
      "Temperature Run: 3 ; Current Time: 184.85 secs\n",
      "Temperature Run: 4 ; Current Time: 236.61 secs\n",
      "Temperature Run: 5 ; Current Time: 283.27 secs\n",
      "Temperature Run: 6 ; Current Time: 341.44 secs\n",
      "Temperature Run: 7 ; Current Time: 406.55 secs\n",
      "Temperature Run: 8 ; Current Time: 471.81 secs\n",
      "Temperature Run: 9 ; Current Time: 529.68 secs\n",
      "Temperature Run: 10 ; Current Time: 583.34 secs\n",
      "Temperature Run: 11 ; Current Time: 633.34 secs\n",
      "Temperature Run: 12 ; Current Time: 682.38 secs\n",
      "Temperature Run: 13 ; Current Time: 730.92 secs\n",
      "Temperature Run: 14 ; Current Time: 780.55 secs\n",
      "Temperature Run: 15 ; Current Time: 831.61 secs\n",
      "Total Simulation Run Time: 831.6078369617462 secs\n"
     ]
    }
   ],
   "source": [
    "## Generate lf for each temperature from .03 to 2.0 in npoints steps\n",
    "t0 = time.time()\n",
    "\n",
    "# number of steps for MC algortihm\n",
    "Nmc = 30000\n",
    "    \n",
    "# parameter definitions    \n",
    "R = 8\n",
    "np.random.seed(24)\n",
    "Del_bar, sigma_D = 3.0, 1.0\n",
    "Dels = np.random.randn(R)*sigma_D+Del_bar\n",
    "E0_bar, sigma_E = 14.0, 2.0\n",
    "E0s = np.random.randn(R)*sigma_E+E0_bar\n",
    "Evs = np.ones(R)*0.001    \n",
    "    \n",
    "# initial monomer and dimer states; \n",
    "# system in microstate of all correct dimers\n",
    "random.seed(0)\n",
    "free_objs_0 = []\n",
    "bound_objs_0 = random.sample(master_list, len(master_list))\n",
    "mstr_vec = copy.deepcopy(master_list)\n",
    "\n",
    "# temperature limits\n",
    "Tmin = .05\n",
    "Tmax = 3.0\n",
    "\n",
    "npoints = 15 #number of temperature values\n",
    "navg = 5 # number of times we run simulation at each temperature; 50 in paper\n",
    "temp_vals = np.linspace(Tmin, Tmax, npoints).tolist()\n",
    "\n",
    "# list of dimer values \n",
    "sim_k_vals = [0]*npoints\n",
    "# list of correct dimer values\n",
    "sim_m_vals = [0]*npoints\n",
    "# accepted list \n",
    "accepted_list = [0]*npoints\n",
    "# saved list for plotting\n",
    "saved_list = dict()\n",
    "\n",
    "for k in range(npoints):\n",
    "    \n",
    "    fin_k_vals = [0]*navg\n",
    "    fin_m_vals = [0]*navg\n",
    "    fin_accepted = [0]*navg\n",
    "    \n",
    "    for j in range(navg): \n",
    "        \n",
    "        # make copy of initial monomer and dimer states \n",
    "        free_objs_copy = copy.deepcopy(free_objs_0)\n",
    "        bound_objs_copy = copy.deepcopy(bound_objs_0)\n",
    "        \n",
    "        # defining helper functions\n",
    "        gammas_ = gamma_func(E0s, Evs, temp_vals[k])\n",
    "        deltas_ = delta_func(Dels, temp_vals[k])  \n",
    "\n",
    "        # metroplois generator\n",
    "        free_list, bound_list, accepted = met_assembly_grid(Nmc, \n",
    "                                                            free_objs_copy,\n",
    "                                                            bound_objs_copy,\n",
    "                                                            mstr_vec,\n",
    "                                                            deltas_,\n",
    "                                                            gammas_,\n",
    "                                                            name_key_) \n",
    "        \n",
    "        \n",
    "\n",
    "        \n",
    "        # averaging final states to compute observables\n",
    "        fin_k_vals[j] = avg_k(bound_list, Nmc)\n",
    "        fin_m_vals[j] = avg_m(bound_list, mstr_vec, Nmc)\n",
    "        fin_accepted[j] = accepted\n",
    "    \n",
    "    # saving every 5 temperatures\n",
    "    if (k+1)%5 ==0 or k ==0:\n",
    "        saved_list[k] = ['white' if x=='-' else x for x in bound_list[-1]]   \n",
    "    \n",
    "    # averaging over computed equilibrium averages\n",
    "    sim_k_vals[k] = np.mean(np.array(fin_k_vals))\n",
    "    sim_m_vals[k] = np.mean(np.array(fin_m_vals))\n",
    "    accepted_list[k] = np.mean(np.array(fin_accepted))\n",
    "    \n",
    "    t_prelim = time.time()\n",
    "    print(\"Temperature Run:\",str(k+1),\"; Current Time:\", round(t_prelim-t0,2),\"secs\")\n",
    "    \n",
    "\n",
    "t1 = time.time()\n",
    "print(\"Total Simulation Run Time:\",t1-t0,\"secs\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Simulated image grid at various temperatures"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAgQAAAIRCAYAAAA1GaW+AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAjrklEQVR4nO3de5Cld1kn8O+TGUKAKRcFNaIgeBt1oAymVNDSRkvB1RrBxF3BVcAqRWPUaHnbnUWRQmdLLZRApaKDtYoXFKkE3LNecZWWAhUVUkCE8YrcFjaAomMMJJPf/vG+PXPS0+dM9znv6e53+vOpOtWn38tzfuf2zLd/73vOVGstAMDBdtleDwAA2HsCAQAgEAAAAgEAEIEAAIhAAABkRiCoquNVdWoymbQkLi4u++MyKvqIi8u+vMxUF/kegrkrgV1Vez2ABekjsH/M7COHt7P3+vr6IKNYW1vb9zWn650+s3S5JMnRI+evH6Sa0/Vy9sTyBZPk0Mnx1SRJUjdeM0iddsOt566voo8MMc7pMe7X9+dYao6xj6zisdwNziEAAAQCAEAgAAAiEAAAEQgAgAgEAEAEAgAgAgEAEIEAAIhAAABEIAAAIhAAABEIAIAIBABABAIAIAIBABCBAABIUq21eevnrgR2Ve31ABakj8D+MbOPbDlDUFXHq+rUZDJZ3ZCAS5o+AuNihgDGwwwBsKyZfeTwtnY/e2KYYRw6ee7q+vr6ICXX1tbOXT99Zvl6R4+cvz5EvYNcc9VjzO0DvS6PnX9druK1TmcV7/m68ZpBarYbbj13fQyv/YNUc+V9ZBXv+ZH2EScVAgACAQAgEAAAEQgAgAgEAEAEAgAgAgEAEIEAAIhAAABEIAAAIhAAABEIAIAIBABABAIAIAIBABCBAABIUq21eevnrgR2Ve31ABakj8D+MbOPbDlDUFXHq+rUZDJZ3ZCAS5o+AuNihgDGwwwBsKyZfeTwdvY+fWaYURw9sv9rjmGMY6k5XS9nTyxfMEkOnTx3dRX3exXjpLNfX6erqDmGMY6l5hjGuLnmWPuIkwoBAIEAABAIAIAIBABABAIAIAIBABCBAACIQAAARCAAACIQAAARCACACAQAQAQCACACAQAQgQAAiEAAAEQgAACSVGtt3vq5K4FdVXs9gAXpI7B/zOwjW84QVNXxqjo1mUxWNyTgkqaPwLiYIYDxMEMALGtmHzm8nb1PnxlmFEeP7P+aYxjjWGqOYYy7UZPOSh7b208MU/TYyXNXvfb3V837PN9nB3q+Dw37fCf3HWfdPEzNdt0wdbbLSYUAgEAAAAgEAEAEAgAgAgEAEIEAAIhAAABEIAAAIhAAABEIAIAIBABABAIAIAIBABCBAACIQAAARCAAACIQAABJqrU2b/3clcCuqr0ewIL0Edg/ZvaRLWcIqup4VZ2aTCarGxJwSdNHYFzMEMB4mCEAljWzjxzezt6nzwwziqNH9n/N6Xp5zonlCybJc0+eu7pf7/cqak7XqxuvWb5gknbDraOrSWe/vk5XUVMfGa7mqseYswM9P4fOPz9j7SNOKgQABAIAQCAAACIQAAARCACACAQAQAQCACACAQAQgQAAiEAAAEQgAAAiEAAAEQgAgAgEAEAEAgAgAgEAEIEAAEhSrbV56+euBHZV7fUAFqSPwP4xs49sOUNQVcer6tRkMlndkIBLmj4C42KGAMbDDAGwrJl95PB29j59ZphRHD2y2pp5zonlCz735LD1dqHmfn1+puvVjdcsXzBJu+HW0dWks19fp5tr6iPL0UeS3D7Q83Ps5MW3GZCTCgEAgQAAEAgAgAgEAEAEAgAgAgEAEIEAAIhAAABEIAAAIhAAABEIAIAIBABABAIAIAIBABCBAACIQAAAJKnW2rz1c1cCu6r2egAL0kdg/5jZR7acIaiq41V1ajKZrG5IwCVNH4FxMUMA42GGAFjWzD5yeDt7nz4zzCiOHpn65Tknhin63JPD1hy63khrDvGcTz/fq3gN3XTtMP8+Xn/L+X+vVvJaJ4k+chBr6iPL2e0+4qRCAEAgAAAEAgAgAgEAEIEAAIhAAABEIAAAIhAAABEIAIAIBABABAIAIAIBABCBAACIQAAARCAAACIQAAARCACAJNVam7d+7kpgV9VeD2BB+gjsHzP7yJYzBFV1vKpOTSaT1Q0JuKTpIzAuZghgPMwQAMua2UcOb2fv02eGGcXRI1O/3H5imKLHTg5bc+h6u1BzfX19kJJra2vnrg/xnE8/36t4Dd107TD/Pl5/y/l/r1byWifJePrIEO+n6ffSWPrIKmqOoY/s1/ud7H4fcVIhACAQAAACAQAQgQAAiEAAAEQgAAAiEAAAEQgAgAgEAEAEAgAgAgEAEIEAAIhAAABEIAAAIhAAABEIAIAIBABAkmqtzVs/dyWwq2qvB7AgfQT2j5l9ZMsZgqo6XlWnJpPJ6oYEXNL0ERgXMwQwHmYIgGXN7COHt7P36TPDjOLokalfbj8xTNFjJ4etOVXvQN3vTTWHuO/T93sVj+VN1w7z7+P1t5z/92olzzlJVvPYrq+vD1JzbW3t/C8j6CP79n4n+kjG20ecVAgACAQAgEAAAEQgAAAiEAAAEQgAgAgEAEAEAgAgAgEAEIEAAIhAAABEIAAAIhAAABEIAIAIBABABAIAIAIBAJCkWmvz1s9dCeyq2usBLEgfgf1jZh/Zcoagqo5X1anJZLK6IQGXNH0ExsUMAYyHGQJgWTP7yOFt7X3zMKNo1+3/mmMY41hqjmGMu1GTzlieL6/9/VVzDGPcXPOma4fJ7tffsrtZ2kmFAIBAAAAIBABABAIAIAIBABCBAACIQAAARCAAACIQAAARCACACAQAQAQCACACAQAQgQAAiEAAAEQgAACSVGtt3vq5K4FdVXs9gAXpI7B/zOwjW84QVNXxqjpVVc/rdx7sUlVfs99rjmGMY6k5hjGOpWZVHc+I6CP7f4xjqTmGMY6l5rw+smUgaK1NWmvPSnLlrB2XsIqmNnTNMYxxLDXHMMax1BxVINBHRjHGsdQcwxjHUnNngWDKZOCBjKXmGMY4lppjGONYaq5ijLthDI/tKmqOYYxjqTmGMY6l5sx6FzuHAAA4AHzKAAAQCAAAgQAAiEAAAEQgAAAiEAAAEQgAgAgEAEAEAgAgAgEAEIEAAIhAAABEIAAAIhAAABEIAIAIBABABAIAIAIBABCBIFX16Kq6p6q+YmrZE6qqVdUzp5Y9uao+UlWfPtDt/kVVvXmIWsDu20Hv0E8YhQMfCJL8dJLXttZeNW+j1tpvJnlzkp9Y9gar6nCSRye5bck6R6rqbN9stnP5mGXHPmcsl1XV91bV26rqrqp6Z1U9v6oetIoac+7jmWHvGcy0rd6xFf3komNZqp/09+VEVb25qv61qt5fVa+rqmdWVW2x/cdX1c/2t/ORqnpHVd1YVQ8e/M7tY4f3egB7qaoen+Qrkjxl06o/TvKAJHdvWn5jkpdU1bHW2u1L3PRnJ7l/lnwDp3v+nrFp2XVJvjDJ9yd539TyD7fWPrjk7c3zM0m+O8krkjw/yWf1vz+2qr68tXbvCmq8JsmpTcs2P2cwuB32Dv1k5xbuJ1V1WZLfSTfulyR5UZIHJnlakl/oa/3Q1PYfl+TPkjwsyc8leUu6gHVdki+pqi9qrd059B3cl1prB/aS5JeT3JHkftvc/kiSf0vyoiVv9+lJWpIvW8F9+ssk/57k8C4+jseS3Jvklk3Lv6u/n98wdI1+2S/u9WvI5WBedto7ZtTQT7a+zaX6SZLH99v9zKbllyf5+yT/vGn5C/rtn7Zp+dP65c/e69fbbl0O7CGDfprtKUn+oLV296Z1FxzzS5LW2pl0f5V+3ZI3/9j+521Tt/ngqnpFPz32rEWKVtX90iXbN7XW7llyjDvxtCSV7o017cVJ7kzyjauqUVWXV9WRHYwVlrLT3qGf7Niy/eSj+p/vmV7YWvtIkvenC2HTvjRd6Pn1TctfluSuJN+8nUFfCg7yIYOr0yX01+9wvz9J8qSq+szW2tsWvO2rkryz9VNuVXV1kpenexN8UWvtLxeseyxdCn7jxTbsp9V2cgzwg232NN3npUv093ksW2t3VdVt/fqLWaTG16VrDoeq6o50b+Bnt9Y+tI3bg0Ut2ju2op9caNl+8vok/5zkB6vq7ekOBzww3eGQq5N8+6bt75/krtZPC0zd3r1V9e9JPqWqHtpae//F7tTYHeRA8Nn9z7/b4X4b2x9Lsswb+DVJUlXXpTte9n+SfFNb7rjcxl8Kb9jGto9I8g87qP2oJG+fse5hSd7fWvvwFuveneQLq+ryPqHPstMar0/X9P423V8EX5XkO5OsVdUX9n99wSos2ju2op9caKl+0lr7p6r6miQ/n+Q3plb9a5JrW2uv3LTL7UmOVtVVrbXbNhZW1VVJPrr/9RHpZhcuaQc5EHxs/3Onb5gP9D8/bpEbrapHJnlwkr+pqpcm+fokP5rkxzYn1AV8bv/zook+yXvTnRS1Xe+ds+6BSbZ68ybdlNvGNvMCwY5qtNa+YNM2v1RVb0ry40lu6H/CKizaO7ain1xoiH5yJt3Jgf8ryevSzV5cn+SlVfXkdt9Phrwg3SGg36iq7+n3O9YvvzvJ/frbu+Qd5ECw8Wa54CMoF7Gx/aJvto3U/V3ppsW+sm3xsaV+Cu5D6T4a2vqff53k+tbaa+fUvifdx5nmaq3dleQPdjz6rd2Z2Q3tiqltVl3jp5I8J8lXRyBgdRbtHVvRTy60VC+oqsekCwHf21r72anlv5buH/sXV9WnttbOJklr7TVV9dQkL0zyW/3mZ9PNMNye5GuT/Mvid2c8DuxJhenOEE52dtxrevs75m4128Yb+MXpjl09ZsZ2R9Mdp3x4a+1Iuqmrt6WbDrxA/4b/nCRv7d+cc1XVoaq6cgeXQ3PKvSfJQ6vq/lus+8R003/z0vwgNfoTvN6T5KEXuS1YxqK9Yyv6yYWW7QXfmy44vHx6Yes+OvhbST45ySM3rXt5kk9K93h+SZKHtda+vV92T7pDk5e8gzxD8Jb+506/KezTNu2/U1ele0FfX1UPSPJTVfX21tqtm7a7Osl7No4BttY+XFVvTPIZM+p+ero3/HaO9yXJwzPcMb8/T/LEJJ+f/lhmklTVFenu7x9vo/7SNfptPynJn27j9mBRi/aOregnF1q2F3xi/3Or0HF4089z+hmD26Zu78p0AWG9HZDvITjIgeCN6aaBHrfD/R6X5H2ttdML3u5jc/6Y3LelS6u/UlVPaK1Nn1V7dZI3JUn/zVqPSfKtSX5kRt2dHO9Lhj3m97IkJ5J8T6bewOnG+8Akv7qxoP8o06cmubO19o4FazyktfaBXOh56V7Tk4veG1jcor1jK/rJhZbtJ3+VLlA8M8lPTm374CRPTvJPuchf/P0MyQvThYqDc/hxr78IYS8v6b616p+S3H/T8iekO872zE3Ll/oikSQP6ev+xNSyB6d7Ab8vyaOmlv9xuhNo/jnd8bJ7kzxvTu2f7Gt/8R49li/qb//WJN+S7tvF7k7y6iSXTW33yH67Vy9R42fSfVzrZLqPEH1/kj/s9/3TJA/Y69eWy6V92Unv0E8Wum8L95N0oegD/X385b5HnEg3g9GSfMcWz8NfpfuH/1uSfF+Sv+i3PbHXr7Vdfdz3egB7eue7KamW7qMo08tnvYGf0S9/9IK39+X9/k/dtPxR/Rv4remO7V2W7iMy12wa64eTrM2o/ar+DfBRe/RYHurfSKf7cb473Xe9H9m03QVv4AVqPDnJ7/Xr7+qb6m39m/6KvX5duVz6l530Dv1kofu2VD9JN2vwkiTvShck/iVdKLpmi9u6PMmvpQsMd6X79MjvJXnSXr/OdvtS/QNyYFXV7yZ5UGvti7ex7RuSvL21ds2Kx/SZ6d7MD2+tvWtq+fuT/EBr7RdWefvAxe2kd8zYXz9hXznInzLY8H1JHl9VT5y3UVU9Jd3XeP7QvO0GsnEC0Lv6235QVf1wuiT7+7tw+8DFbat3bEU/YT86yCcVJkla97+MXfRxaN23W12+8gF1rk5yZXX/le+96aa7Xpfua0jfvUtjAObYbu+Yse8ro5+wzxz4QwYAgEMGAEAEAgAgAgEAEIEAAIhAAABEIAAAIhAAABEIAIAIBABABAIAIAIBABCBAACIQAAARCAAACIQAAARCACACAQAQAQCACACAQAQgQAAyIxAUFXHq+rUZDJpSVxcXPbHZVT0EReXfXmZqVqbu37uSmBX1V4PYEH6COwfM/vI4d0cxaqdPrN8jaNHhq13kGtO17vp2mH+Lbv+lql/W86eGKRmDp0cpg6XhDH0kZVYwfupbl6+XLtu+Rpsj3MIAACBAAAQCACACAQAQAQCACACAQAQgQAAiEAAAEQgAAAiEAAAEQgAgAgEAEAEAgAgAgEAEIEAAIhAAABEIAAAklRrbd76uSuBXVV7PYAF6SOwf8zsI1vOEFTV8ao6NZlMVjck4JKmj8C4mCGA8TBDACxrZh85vJujWLW6efka7brla8xz+swwdY4e2f81p+sN8dwk931+1tfXB6m5trY2SB3Yyn59f26uCU4qBAAEAgBAIAAAIhAAABEIAIAIBABABAIAIAIBABCBAACIQAAARCAAACIQAAARCACACAQAQAQCACACAQCQpFpr89bPXQnsqtrrASxIH4H9Y2Yf2XKGoKqOV9WpyWSyuiEBlzR9BMbFDAGMhxkCYFkz+8jh7ex907XD9KHrbznfF06fGaRkjh45f32ImkPXO8g1p+utwirud55zYpiizz05TB1Gb7++P8dSc9V9hPOcVAgACAQAgEAAAEQgAAAiEAAAEQgAgAgEAEAEAgAgAgEAEIEAAIhAAABEIAAAIhAAABEIAIAIBABABAIAIAIBAJCkWmvz1s9dCeyq2usBLEgfgf1jZh/Zcoagqo5X1anJZLK6IQGXNH0ExsUMAYyHGQJgWTP7yOHt7H36zDCjOHpk6pezJ4YpeujkuatDjHMMY0zuO879WnO63vr6+vIFk6ytrZ27vl/v9+aadA7qY3tQ73cyzPt+1e95znNSIQAgEAAAAgEAEIEAAIhAAABEIAAAIhAAABEIAIAIBABABAIAIAIBABCBAACIQAAARCAAACIQAAARCACACAQAQJJqrc1bP3clsKtqrwewIH0E9o+ZfWTLGYKqOl5VpyaTyeqGBFzS9BEYFzMEMB5mCIBlzewjh7ez9+kzw4zi6JHz1+vGawap2W649XzNmweod93UL2dPLF8wSQ6dPHd1LI/lEOOcHuMq7vdYnh8Y2lj6SG4f4D167OTFt2EQTioEAAQCAEAgAAAiEAAAEQgAgAgEAEAEAgAgAgEAEIEAAIhAAABEIAAAIhAAABEIAIAIBABABAIAIAIBABCBAABIUq21eevnrgR2Ve31ABakj8D+MbOPbDlDUFXHq+rUZDJZ3ZCAS5o+AuNihgDGwwwBsKyZfeTwdvY+fWaYURw9stqa6+vrS9dbW1s7d71uvGbpeknSbrh1kDqz7NfnZ9XP9ypqjuU5pzPEez657/t+6Nf+Kozl/cS4OKkQABAIAACBAACIQAAARCAAACIQAAARCACACAQAQAQCACACAQAQgQAAiEAAAEQgAAAiEAAAEQgAgAgEAECSaq3NWz93JbCraq8HsCB9BPaPmX1kyxmCqjpeVacmk8nqhgRc0vQRGBczBDAeZgiAZc3sI4e3tfvtJ4YZxrGT566ePjNMyaNHzl+vm5ev164btt7mmuvr64PUXFtbO//LcwZ6fp578uLbLGgVz/cqjGWcY3RgH9uR9M8xWEVP5jwnFQIAAgEAIBAAABEIAIAIBABABAIAIAIBABCBAACIQAAARCAAACIQAAARCACACAQAQAQCACACAQAQgQAAiEAAACSp1tq89XNXAruq9noAC9JHYP+Y2Ue2nCGoquNVdWoymaxuSMAlTR+BcTFDAONhhgBY1sw+cng7e58+M8wojh6Z+uX2E8MUPXby3NUhxnmfMY7E+vr6IHXW1tbOXa+bl6/Xrjt/fRWvobHUBPavm64dJmdff8v4c6+TCgEAgQAAEAgAgAgEAEAEAgAgAgEAEIEAAIhAAABEIAAAIhAAABEIAIAIBABABAIAIAIBABCBAACIQAAARCAAAJJUa23e+rkrgV1Vez2ABekjsH/M7CNbzhBU1fGqOjWZTFY3JOCSpo/AuJghgPEwQwAsa2YfObybo5h2+swwdY4emfrl9hPLFzx28tzVVYxx397v5D73/aA+lit5fhiXswO89g+dvPg27NhY3p9jGedmTioEAAQCAEAgAAAiEAAAEQgAgAgEAEAEAgAgAgEAEIEAAIhAAABEIAAAIhAAABEIAIAIBABABAIAIAIBABCBAABIUq21eevnrgR2Ve31ABakj8D+MbOPbDlDUFXHq+rUZDJZ3ZCAS5o+AuNihgDGwwwBsKyZfeTwtnY/e2KYYRw6ef767QPVPHa+5ukzy5c7emTqlxWMsW4epmS77vz1/Vpzut4Qz02y6flZgZuuHebf3Otv8W8gq1M3XjNInXbDreeur6+vD1JzbW1tkDq7xXv+PCcVAgACAQAgEAAAEQgAgAgEAEAEAgAgAgEAEIEAAIhAAABEIAAAIhAAABEIAIAIBABABAIAIAIBABCBAABIUq21eevnrgR2Ve31ABakj8D+MbOPbDlDUFXHq+pUVT2v33mwS1V9zX6vOYYxjqXmGMY4lppVdTwjoo/s/zGOpeYYxjiWmvP6yJaBoLU2aa09K8mVs3Zcwiqa2tA1xzDGsdQcwxjHUnNUgUAfGcUYx1JzDGMcS82dBYIpk4EHMpaaYxjjWGqOYYxjqbmKMe6GMTy2q6g5hjGOpeYYxjiWmjPrXewcAgDgAPApAwBAIAAABAIAIAIBABCBAACIQAAARCAAACIQAAARCACACAQAQAQCACACAQAQgQAAiEAAAEQgAAAiEAAAEQgAgAgEAEAOWCCoqkdX1T1V9RVTy55QVa2qnrnAsidX1Ueq6tMHGt9fVNWbh6gF7J7t9pYd1Bu0t/Q19RfmOlCBIMlPJ3lta+1VQxRrrf1mkjcn+Ylla1XV4SSPTnLbknWOVNXZvhFt5/Ixy459zlj+W1W9vKr+vr+tty9Q4+Or6mer6p19g3xHVd1YVQ/eYttZ9/HMEPcH5ti3vSUZpr/sl95SVUer6ler6q1V9aGqurOq3lZVP11Vn7CDOpdV1ff2+97V95jnV9WDVjHuMTi81wPYLVX1+CRfkeQpm1b9cZIHJLl7gWVJcmOSl1TVsdba7UsM8bOT3D9LBoJ0z+kzNi27LskXJvn+JO+bWv7h1toHl7y9eU4m+WCSNyR58E53rqqPS/JnSR6W5OeSvCVdU7suyZdU1Re11u7ctNtrkpzatGzzcwaD2WFv2YmheksyTH/ZL73lk5J8QpJXJHlXknuSPCbJs5I8taquaq39v23U+Zkk393XeX6Sz+p/f2xVfXlr7d5VDH5fa60diEuSX05yR5L7DVz3SJJ/S/KiJes8PUlL8mUruO9/meTfkxze5cf8U6auvyXJ23e4/wv6x+Rpm5Y/rV/+7E3LW5Jf3O3XlsvBvuz33tLXWkl/2aveMmMs/6m/jz+4jW2PJbk3yS2bln9XX+Mb9vr+7MXlQBwy6KfLnpLkD1prd29at/A5BEnSWjuT7q/Sr1tymI/tf942dZsPrqpX9NNZz1qkaFXdL91f1W9qrd2z5Bh3pLX290uW+NJ0zebXNy1/WZK7knzzVjtV1eVVdWTJ24aL2klvqar/2P/+3TNq/UlV3dG/Z4fsLckK+ste9pYZ/rH/+dHb2PZpSSrdHx3TXpzkziTfONywxuOgHDK4Ol3afv2K6v9JkidV1We21t62YI2rkryz9dNsVXV1kpene9F+UWvtLxeseyzJ5UneeLENq+qyJDs57vfBttpptfsnuav10X1Da+3eqvr3JJ9SVQ9trb1/avXXpXszH6qqO9KFh2e31j60wnFycO2kt/x+kvem+2v9hdMr+pMHH5fkhZuCxRC9JVlNf9l2b+lvc9D+UlVXpHvsr0h3SGTjfIvf3kbtz0s3Q3Cf5621dldV3davP3AOxAxBuhdLkvzdiupv1D22RI2r0qf3qrouyWuTvDXJ1UuEgeT8XwZv2Ma2j0g39bndyyOWGNd23J7ko6vqqumF/e8bfwVMj+H1SX40XSh4RpI/TPKdSV5jxoAV2XZvaa2dTfIrSa6uqs/etPrp/c+XbFo+RG9JVtNfdtJbkuH7y7f0270zye+lO0/pG1trr9nGWB6W5P2ttQ9vse7dSR5aVZdvo84l5aDMEHxs/3NVJ7l8oP/5cYvsXFWPTPdi/puqemmSr0/3D9uPbf7reAGf2//cTop/b7qTo7brvTsfzo68IN107G9U1fekOw/hWL/87iT3S/LAjY1ba1+waf9fqqo3JfnxJDf0P2FIO+0tL0l3At7Tk/zXJKmqSjer9ZbW2uZ/XJfqLX39R2Y1/WUnvSUZvr+8Msnb0s0SPDbJ1yR56DZrPzDJVmEg6Q5HbmzzkW3WuyQclECw8aKvFdXfqLvom2sjaX9Xummsr2xbfHypn3L7ULqZndb//Osk17fWXjun9j3pPsI0V2vtriR/sOPRr0hr7TVV9dR006u/1S8+m+Tn080efG2Sf7lImZ9K8pwkXx2BgOHtqLe01t5SVW9I8l+q6kQ/Jf4lSR6Z5Ae32GXZ3pKsrr9su7ckw/eX1tq70n3KIEleWVW3JPnzqnpga+1/XGT3OzM7ZF0xtc2BclAOGdzR/1zVZ+436t4xd6vZNt6wL0533PwxM7Y7mi4NP7y1diTdtPnb0n185gL9G/xzkry1fzPOVVWHqurKHVwO7fB+7lhr7eXpPmb02HSN82GttW/vl92T5G8vsv/dSd6T7f/lADuxSG/5pXSv3y/rf396uqD7K1tsu2xvSVbQX3baW/p9VtpfWmtvSjdb8R3b2Pw96Q4L3H+LdZ+Y7nDCgZodSA7ODMFb+p+DfevXJp+26XZ26qp0L8Drq+oBSX6qqt7eWrt103ZXJ3nPxolBrbUPV9Ubk3zGjLqfnu4Nvt1jfA9P8g87GPejkrx9B9svpD/2etvG71V1Zbomt94u/B6C++hPPPqkJH+6yjFyYC3SW16abubq6VX12nTnvLyqtfZ/t9h22d6SrKa/7LS3JLvTXx6Q7YWzP0/yxCSfn+6THEnO9Yur0n2HxIFzUALBG9NNLT9uRfUfl+R9rbXTC+7/2Jw/DvdtST45ya9U1RNaa9NnwV6d5E3JueOOj0nyrUl+ZEbdvT7Gty39x5c+NcmdrbV3XGTby9IdQjiUqUMAVfWQ1toHttjleele55Mhxgqb7Li3tNbuqKrfSXJNun94PioXnky4Ydnekqymv+y0tyQD9ZequrK1dsG6qvrSdB+DfPWm5Vv1l5clOZHkezIVCNLd3wcm+dUdjPOScSACQWvtbFXdmuQpVXX/GWeWLqQ/e/2Lk/zPBfd/SLq/YF/aj/Xuqro2yeuSTKrqca21jVR9dZLPr6p/TvdxnyuS/Hhr7aUzyu/oLOChj/FV1Telaz5Jd/LV5VX17P73f2yt/XJ//RPTnfG8nuQJU/tvfJzrFen+svgP6T4/fHWS/95a+6Opm3t2VT0uyR8leUe6v16+Kt13GfxZkhcNdb9gwxK95SXpToJ7frrj9q/cvMGyvaWvsar+stNPGAzZX26u7iuK/zDddw9ckW7sT03yr0m+b9P2F/SX1tqbq+qmJN/ZP3+/nfPfVLie/vE6cPb6m5F265JuaqgluXbT8if0y5+502X98mf0yx+94Li+vN//qZuWPyrdV4G+Nd2xvMvSvdiv2XSfPpxkbUbtV6U7ieij9ugxf3V/37a6vHpqu0duXtYvvzzJr6ULA3elO5P795I8aYvbenK/7t39tv+W7jDDiSRX7PXrz+XSveykt0ytuzzdJwhakhfPqLtUb+lrrKS/7GVvSfKfk/zvdB83vCvdl5e9LV3of8QW28/qL4fShYfT/f18d7r/k+LIXr+m9upS/QNzIFTV7yZ5UGvtiwes+YZ0X8l7zVA1Z9zOZ6Z78z68dWfXbix/f5IfaK39wipvH5htzL2lvy39hQPzKYMN35fk8VX1xCGKVdVT0h2z+qEh6l3Exgk/7+pv+0FV9cPp/tL4/V24fWC2MfeWRH8hB+Qcgg2t+x/DBrvPrbVXpnvD7Iark1xZ3X/le2+6E5lel+5rR9+9S2MAtjDy3pLoLyQH65ABALC1g3bIAADYgkAAAAgEAIBAAABEIAAAIhAAABEIAIAIBABAkv8PTdCik0sgFRwAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 648x648 with 4 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# figure parameters\n",
    "rows, cols, idx = 2, 2, 0\n",
    "fig, axes = plt.subplots(nrows=rows, ncols=cols, figsize=(9,9))\n",
    "\n",
    "# list of keys for saved snapshots of image\n",
    "img_key_list = list(saved_list.keys())\n",
    "\n",
    "for i in range(rows):\n",
    "    for j in range(cols):\n",
    "        if idx < 4:\n",
    "            axes[i, j].imshow(color_to_rgb(list_to_matrix(saved_list[img_key_list[idx]])))\n",
    "            ax = plt.gca()\n",
    "            \n",
    "            # making ticks invisible\n",
    "            axes[i, j].set_xticks([])\n",
    "            axes[i, j].set_yticks([])\n",
    "\n",
    "            # Minor ticks\n",
    "            axes[i, j].set_xticks(np.arange(-0.5, 11, 1), minor=True)\n",
    "            axes[i, j].set_yticks(np.arange(-0.5, 10, 1), minor=True)\n",
    "            axes[i, j].tick_params(axis='y', colors='red')\n",
    "            \n",
    "            # labeling images\n",
    "            itimes = 'i'*(1+idx) if idx<3 else 'iv'\n",
    "            \n",
    "            # Gridlines based on minor ticks\n",
    "            axes[i, j].grid(which='minor', color='w', linestyle='-', linewidth=3)\n",
    "            axes[i, j].set_title(fr'({itimes})  $k_BT = {round(temp_vals[img_key_list[idx]],2)}$', fontsize = 18, y = -.2)\n",
    "\n",
    "            # making spines invisible\n",
    "            axes[i, j].spines['right'].set_visible(False)\n",
    "            axes[i, j].spines['top'].set_visible(False)\n",
    "            axes[i, j].spines['left'].set_visible(False)\n",
    "            axes[i, j].spines['bottom'].set_visible(False)  \n",
    "            \n",
    "            idx +=1\n",
    "    \n",
    "# plt.savefig('general_grid_assembly_grid_plots.png', bbox_inches='tight', format = 'png')      \n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Comparing analytical and simulation results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "<ipython-input-2-a582446d2d66>:4: RuntimeWarning: divide by zero encountered in true_divide\n",
      "  phi_func = lambda x, z, gamma, delta: x*(1+ 1/(z*gamma))/(1-delta)\n",
      "<ipython-input-2-a582446d2d66>:4: RuntimeWarning: invalid value encountered in multiply\n",
      "  phi_func = lambda x, z, gamma, delta: x*(1+ 1/(z*gamma))/(1-delta)\n",
      "/Users/mobolajiwilliams/opt/miniconda3/lib/python3.8/site-packages/scipy/optimize/minpack.py:175: RuntimeWarning: The iteration is not making good progress, as measured by the \n",
      "  improvement from the last ten iterations.\n",
      "  warnings.warn(msg, RuntimeWarning)\n",
      "<ipython-input-2-a582446d2d66>:4: RuntimeWarning: overflow encountered in true_divide\n",
      "  phi_func = lambda x, z, gamma, delta: x*(1+ 1/(z*gamma))/(1-delta)\n",
      "/Users/mobolajiwilliams/opt/miniconda3/lib/python3.8/site-packages/scipy/optimize/minpack.py:175: RuntimeWarning: The iteration is not making good progress, as measured by the \n",
      "  improvement from the last five Jacobian evaluations.\n",
      "  warnings.warn(msg, RuntimeWarning)\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAa8AAAFLCAYAAABhvgOVAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAABln0lEQVR4nO3dd3hUxfrA8e/sphfSEyCEkgCKQKgiUgSlI1XlUlWkiYKiV+HaCz9FvKBXbCgIggWsIEhR7KDYEAOEIhogpAAhnfRNdn5/JFlII4GUzW7ez/PsAzs7e847e5J9M+fMmVFaa4QQQghbYrB2AEIIIcSlkuQlhBDC5kjyEkIIYXMkeQkhhLA5kryEEELYHEleQgghbE69SF5Dhw7V8fHxGmhQj9JtVkpppZTV46rLNjeEh7S5YTwaWpvrqL0VcrjYi3UlMTGRgoICa4dR50q3eebMmVaKpO7IcW4YpM32z9rtrRc9r6pKSUkhKCiIqKgoTpw4gVKKPXv2AHDgwAGCg4PJzMy0cpSXb8WKFaxYscLaYVyyix2X0s+PHDli88dJCGF9NpW8Fi1axPDhwwkLCyMkJIRTp07RuXNnADp27EjPnj158cUXrRtkA3Sx41L6+ZVXXinHSQhRbTaTvLKysnjrrbeYPn06AEajkcaNG+PgcP7M5x133MHy5cvJz8+3VpjV8scff/DHH39YO4xLUtlxscfjJISwPptJXtu2bUMpRe/evYGyp6MABg8eTHJyMt9//72Voqye7t270717d2uHcUkqOy72eJyEENZnM8lr165ddOvWDaVUhXWcnJzo3LkzP/zwQx1G1rBV5biUJsdJCFFdNpO8oqOjadq0aaX1mjZtyokTJ2o/IAFU/biUJsdJCFEdNpO8srOzcXFxqbSeq6sr2dnZdRCRgKofl9LkOAkhqsNmkpe/vz8pKSmV1ktOTiYgIKAOIhJQ9eNSmhwnIUR12Ezy6tKlC4cOHaq0XmRkJF27dq2DiARU/biUJsdJCFEdNpO8hgwZwuHDh0lKSqqwzokTJ4iLi2Pw4MF1GFnDVpXjUpocJyFEddlM8urYsSM9evTggw8+qLDO+vXrGTx4MC1atKjDyBq2qhyX0uQ4CSGqq17MbVhVTz75JPPmzWP27Nm0bNkSrc/P25ibm8vy5ctZv369FSNsmC52XOQ4CSFqg830vACGDh3KnDlziI2NLfNadHQ0jz76qOVmWVF3LnZcSouLi5PjJISoNpvqeQHce++95Za3bduWtm3b1nE0olhFx6W00NBQ+vXrV8vRCCHsnU31vIQQQgiwwZ6XPWsI63kJIURNkORlReboaGIWPYfB1ZUmzy2yybW8hBDCGiR5WYE5M5Ozr79O1pq1ULQaqTblEfzSSyij0crRCSFE/Vfvr3ntjkqk9+Jv2R2VaNP7ANj9z1nunb2Uw4OHkbxqNZjNeI0di6FRI8599TU7Zs8usXTIZe2jrtpSR/sRQojyVCl5KaXmKqX2KKVylVJrKql7v1LqtFIqXSm1WinlfLnB7Y5KZPqaPcSlZjN9zZ5a+aKsi30A/LJzL9HTZ3DX96swJJ0lv82VuL7yMk2fW0TI66+hnJ0Z+tZbXH311Ze9j7pqS13tRwghKlLVnlc88Ayw+mKVlFJDgIeAAUALIBR4+nICK/6CzDYVnlbLNhXU+BdlTe9Dm0yYEhLI+esvMn/+mbStW0l+733+fPBR3GbfRuczRznn6MrLnW5mQqeZ7HMPBsCte3eCX3yBq1xcuMrZmeR166zfFq0xmU1kmbIoMBeU2s+v5DlEYXT/G5PzfmZ8spL//rSWj49+zP6z+zEVmC5rn+WRHp4QojzqwtkPKq2s1DNAM6311ApeXwec0Fo/UvR8APC+1rrxxbbbvXt3vXHjRkJCQgDIWxyGU07FX1aZ+f780XEuX0d/dfF4i5rW2rs1Y9qMAa1Ba7TWbDz4FT+e/AVl1ijAYAaD1hg0OBVAkIsjbiiMBZrW7i0I92qHzsvFnJlJQWYmMaePkpWWhFNuAU65ZpxzC3DJMV80nh3Nr2Z1+xtJc/bAOWgzDm4naeXvjodz4aXHLr8mMerjWLSCj6e04HAnbwDuDL+T65tfX2JbC35YwIn0EwBk5OZzMjkTs+VYnv/XoGBB90e5rWvJ99+y+RZSc1MxazMFuqDwYS78N9+cj8l8PgFtGr2JUO/Q8wmyIAvPK56ssJ3ORmfa+7WnU2AnugR04Zom1+Dm6AZATEyM5ThX5sKE7OpoZNXU7vQK86/Se+uTS2mzvZA22786am+Fq9zW9ICN9sCmC57vA4KUUn5a6yrN3Lo7KpFeF0lcALnHswj4ZAUTqxzWcWIpmeh6FD2q5jTJ/FqixKfocSGzgnOukO5W/FCku0KqDuXnRsP5y/f8XH4Gp0QMLrFEZwAZhWUHW0PmdQYm7jQzdl00hwtiONjCQGpuapmIotKiOJpy1PJcOUNFQz0Wf7mP1l4dS3zxJ2QlkJJbtaVMTGZTqZ7dxX9scgty2Zuwl70Je3mbt/ni5i8syQvgVMYpPj/2Od7O3vi4+ODt7E2AawCBboGWehX1JG01gQkhalZNJy8PIO2C58X/9wRKJC+l1CxgFkBwcLBlVvL7PzhcKk2U5ehWwJ/Nm6AdT1dYx9KfVODj5EObRm2Ld8xvsRk0cdhHtLMRswKtsPxbYACTEfIdIN9Y+P9Qnyvo2bQvODqiXN1Qbq5sPPsFP2bsJccZsp0KH1kuoFXZPxRyE9qSl1S1SWg39lJ4ZyqG/aGZ/4mZpyYrkpOTiYmJKVHPZKr6qbnc/ALu/2Avn9zWrpwPqGJGZcRBORB3Ko6HP4u2JBIwkp/VErQRzE5osyNoR1D5OLnHoh3OH2pfJ18KkguISSmMPykpiWNJx3hl3yvl7tPDwQN3gw+nUlwx+zfCMTcIU3JfoDCB3fH2b/z3xlZ0beZR5fZb26XMuG8vpM32ry7ae7GeXU0nrwyg0QXPi/9/rnRFrfUKYAUUnjb08/MjJCSE/01whXcvvhPP4Fx8//0Gwf65VQqqkXMjWjQ6nzwSoxIJXdeWwEqGpe8fvoHwYG98XX0J9ggu8dq4c30ZlJN68ffHpfLMlsOYcsp+0eacGYk6m4Ozg4HHR1xFeDNvOgZ0BGDw6QjMTy/D7dufeX6TB/5jryaw1EF80eNFcgty2R+XysLPD5Obby5KSBcmz8JLms7aj//d3pWQkPM9lg9HfgiA0WDEoAwYlRGjwViYsAwOOBocMajzl0T/51KqJxQ9u0ybXB2NvHFjd9o2hX0J+/gz4U+cHZxp3rx5iXoxeTFl3lssIz+DDDIwuIMTUJAdYkleALn5moe//5JuHf6ic5PWNG/UnDbebWjj0wZ3R/cKt2ttDel0UjFps/2zZntrOnkdBDoBHxU97wScqeopQ6DKp4SGtrv8eQx7hflDQQGNCwouWq9j574VvhbiGUKI58UPXMcAaOXZrvBLn5L70nkBODkoVt3ao0ybOwR1wvzSG8TceSdZP/9C+uz78Vn3Po6Nz186bO3TGoD2/tDC44oSieVCFV0rauLR5KKxl9YrzJ9VU7tXeT8DWgxgQIsB5W4r1DuU6R2mk5qbSnJOMik5KZzNPktCVkKJa20AZlOjMu/Pd4wmIvlHIpJ/LFEe7BFMG582tPFuQ1uftrT1aUuLRi0wGuTeOSHsTZWSl1LKoaiuETAqpVyAfK11fqmq7wBrlFLvUzhC8TFgTc2Fa3t6bbiWw8aEii9IbQiE+X+XKTY4OdHslVc5OX0aOfv2c3LadFq89y4Ovr5l91FBYqnpQQ41tZ/2fu1p79e+TLlZm9lxJIr7P9lJnkrG4JCOOd+rTD1H5+RytxuXEUdcRhzfx3xvKRsVNopn+zx70Xh2RyUy/+P9LBkXLtfThLARVR0q/xiQTeEw+ClF/39MKdVcKZWhlGoOoLX+Avgv8B1wEogGKh6W1hBkJlz260YPd5q/+SbObduSd+wYJ2fMoOBcmTOwwPnE4upYmCVra3Rebe7HoAwMbdeGVZNG45TbEVNqTwoy2pWo4+po5On+9/Df6/7LnM5zGN5qOK29W+Ogyv87rJ1vuzJlS35fwoKdC/jor4/YGBnB9DW/yz1rQtiYKvW8tNZPAU9V8HKJCzpa6xeBF6sVlbAwenvTfNVbnJg8hdxDh4mZfRfN31qJwdW1TN3ixFLbvYja3s/l9PBMBSaOpR3j79S/+Tvlb46mHOVw0mHCA8LL1P3m5DfEZcSx/fh2AAwtvHDJDMOU2Zrp76ax6tYB0gMTop6rn3MbugdevMfiHlh3sdQDDgEBNF+9mugpU8j+4w9i77mXZq+/hsHJqUzdXmH+/PTQDbUeU23vp3QCq6yH52h05ArfK7jC9wpLWXn3MCZmJxKXEVeizOCYhsF7L47eewGY9e0KBkT1YWLHYXQL6oajwbEGWyaEqAn1M3mVcw2ooXNqFkzz1auInnIrmT/+SPyD8wl+8QWUQ/08hDWhuj08Vc5tC97O3jzcaTnPfL0d7fI3RrdjKGPJUavKKYFvT23g21Mb+OLmL8qMNBVCWF+9n5i31lTWe6uHvTvn0NDCU4YeHpzbsYNTTzyJNl98Vg9bV9zDq6nTeA4GB177Io+sxF5kx95OxtEnyDx+N7kJg8nPDEObz/8xYDA1LZO4YtJj+PDIhyTnlD9oRAhRN+z3z/bK2GjvzuWqqwh58w1OTp9B2oYNGFycCXr88XJ7GaJ8S8aFX3A9zYg5pzl5Oc0h6QZQJoxuJ3BpdITxXbuUee+X0V+ybO8ynvvtOXo17cWNoTdyfcj1JWYQEULUvobb87Jhbt260ezVV1FOTqSsW8+ZZxeVe31HlK/0iMkStCNOeVeycsQzPNb3zjIvfxfzHQAFuoBdcbt4aNdD9P+oP//Z+R92xu4sc5/ahWSSYSFqjiSv2lZLpyc9+vSm2SsvoxwdSXnvPRIWL5YEdgkqSmAXGxiitWZU6Ci6BJbskWXnZ7Pt+DbmfDOHgR8P5MU9L3Iy/WSJOrKMjBA165Jmla8tpWeVbyhKz8pcfOrvUo7Jue++I/beeWAy4TttGoHzH6zXpxDr28zblztzfey5WLYf387WY1uJSosqt86aoWvoFtSNTT8f5qFt0bV6A3l9U9+Oc11oaG229qzy0vOqR7p27UrXrl0v6T2e119Ps5f+Bw4OJK9ezdkXX5Qe2CUo7oEFe7teUjJp5tmMmeEz2Th6Ix+P/Jip7acS6Ha+Fx3oGkingE7sjkpkwdbjZabUqo316YRoSKTnZUU1+ZdL+ldfEXff/VBQgN/sOwmYN69e9sDs+a/TAnMBP8X/xMd/fcxV/lfRyWNciRutHTwO4uD9B6ak6yjIbgEou+2B2fNxrkhDa7O1e14Nd7ShnWk0aBC88AJxDzxA0htvogxGAu69x9phNShGg5Hrml3Hdc2uQ2tNn+e/K9HjcvT5BQePv3H0PERBVnPykvuSfa498z/eXyc3lgthT+S0oT1Y0gae8qLRL/8i+JqzoDSJr7/O2VuC4SmvwtdFnVJKsWRcuGVAiHJIx+j+j+V1o9tJXJu9j0fYC9zY5xjZ+dnWClUImyTJqx5RSl3eqb4LptJq1DyHpj1TCxNYZCMSD3pUPjmwqBXF19OcHRQ6vxFZx+4jL7U72nx+hKNySmJ91DIGfzKY5fuWcy6v/ImXhRAlSfKyQ14tsml6TSoozdkDRQlMWEWvMH/+e2MrXB2NmPOCyD11C5lR/8GcfANuDp6Weqm5qbwe8TpDPx3Kiv0rZNCNEJWQ5FWPaK1r7EvLq2U2TXukAkUJbMXKGtmuuHRdm3mUuKfMRfmwctTjfPevr3mox0MlpqBKz0vncNLhejnYRoj6RJKXHfNqlU2Ta1IBzdkXXyTprbesHVKDVd6QfDdHNya3m8yWsVt4ts+zlpW5Z3eaXeb90hMToiQZbWjnvFtlg4ZTv/uSsPQFUAb8pk+zdlgNUkXLyDgYHBgVNorhrYbzx5k/SizrApBlymLytslMvHIiN7W5CQeD/NoKIT2veqRbt25069atxrfrHZpNk2f+D4CEJUtIentNje9DVJ+DwYFrmlxTpnz9kfX8k/oP//fL/3Hz5pvZGbtTemKiwZPkVY/s3buXvXv31sq2vW++mcb/txCAhOefJ3nt2lrZj6hZWmu+PPGl5fmxtGPM+WYOM3fM5HDS4XLfY40JgFNSUggKCiIqKooTJ07QvHlz9uzZA8CJEydQSlX4/MCBAwQHB5OZmVln8QrbJ8nLHlRx8l+fceNo/PTTAJx5bjHJ77xb25GJalJK8c6wd5jXdR7uju6W8l9P/8r4LeN59MdHOZ152lJurQmAFy1axPDhwwkLCyMkJIQ9e/bQuXNnAEJCQjh16lSFzzt27EjPnj158cUX6yRWYR9keigrqomJeS9HygcfcPqpwiQW9Nhj+E6ZXKv7u1BDm0IHaq7NSdlJLN+3nE+OfkKBPj9zh4vRhRkdZ3Cl20hmv7O/zicAzsrKIjg4mM8//5w+ffoAl97mLVu2MGvWLE6ePImDja4O3tB+tq09PZT0vBognwkTCHricQDOPPMMye+/b+WIRFX4ufrxWM/H2DB6A/1D+lvKcwpyeDXiVeZ8O41sU8n1xOpiAuBt27ahlKJ3794Al3zaEGDw4MEkJyfz/fff11qcwr5I8mqgfCdNIuixxwA483/PkLJ+vZUjElUV6hXKKze8wuohq7nS90pLeW5aOOX9Std2Atu1axfdunWr1r1pTk5OdO7cmR9++KEGIxP2TJJXA+Y7ZTJBjzwCwOmnF5LywYdWjkhciqsbX80HN36AS9rNFGQ1x5Tcu5xahaegs00FzP94f63EER0dTdOmTau9naZNm3LixInqByQaBEleDZzvbbcS9PBDAJx+6ilSPvrIyhGJS2E0GFl241x0/FxK37ZpcDqDW6tlGF2P4+poZMm48FqJITs7GxcXl2pvx9XVlexsmaBYVI0kL4Hv7bcT+NB/ADj9xJOkfPyxlSMSl6Jw9o6rLdNPFdI4N/4Mo8tp3Fq+ycC+P9GpuWut7N/f35+UlJRqbyc5OZmAgIAaiEg0BJK8BAB+U6cSuGABAKcff4KUD6UHZkuKp5+yLMHilIjRJd7y+venPmPsprHsjttd4/vu0qULhw4dqvZ2IiMjL3klcdFwSfISFn7T7iBw/nwATj/5pAzisDEXJjCdF0DByQfp6HOt5fVTmae48+s7eezHx0jLTaux/Q4ZMoTDhw+TlJR02ds4ceIEcXFxDB48uMbiEvZNkpcowW/6tPOnEJ9eSPJ7MozelpSYAPjWgbw/8k0W912Mt7O3pc6mqE2M2TSGb6K/qZF9duzYkR49evDBBx9c9jbWr1/P4MGDadGiRY3EJOyf3KRsRda6Sbkqkt99jzPPPgtA0MMP4Xv77TWy3YZ2IyfUjzYnZSex+LfFfHHiixLlw1oO49Gej+Ll7FWt7X/xxRfMmzePQ4cOYTQaL6nNubm5tGnThvXr11vuFbNF9eE41yW5SVlY1OR6XtXle+uU8zcyP7eYpNVvWzkiUR1+rn4s6beEZdcvI8D1/KCI7Se2cyip+terhg4dypw5c4iNjb3k90ZHR/Poo4/adOISdc8252ERdcJ30iSU0YHTTz5Jwn//C+YC/GbMsHZYohpuaH4D3Rt3Z+nvS9n4z0YmXTmJa5teW/kbq+Dee++9rPe1bduWtm3b1kgMouGQ5CUuymf8v1BGA6cef4KEpS+g8wvwn32ntcMS1dDIqRELey9kUItBdG/cvczruQW5OBudrRCZEFUnpw3rkdpaz6u6vG+5hSbPPgtKcfall0j430v15vSmuHx9m/XF1aHkvV+pOamM2jiKtw68RYG5oIJ3CmF90vOqR2prLa+a4H3TWNSX/yF+p4GkN9/E/M2LBHVLo8R0du6BMP9vq8UoqkdrzcJfFhKfGc+yvcv4Me5HFvddTGP3xtYOTYgypOdVj+zZs6fETNv1jVeTMzTrnYIyaFL+cSf+F2+0+YIKmQlWi01UX6YpkzOZZyzP/zjzB+M+H8fO2J1WjEqI8knyqkfq62nDC3k2yyGkXxLKwUx6tBuxP/kgZ5fsg4eTB2uHreXuTndjUIVfDam5qcz5Zg5Lf1+KqeD8civWWK1ZiAtJ8hKXzD0ojxbXJ2FwMpMR50rMD34UmC5/OQxRfzgYHLir812sHrKaQLfzK3SvPbSW27+4ndhzsVZbrVmIC0nyqkdmzZrFrFmzrB1Glbj6mWhxQyJGlwKyEpw5+b0fBbmSwOxFt6BufDLyE65rdp2l7EDiAW7adAszPl5tWa25Lha7FKI8krzqkZUrV7Jy5Uprh1FlLt75tByQiKN7PjlJTkR/609+onyJ2QsfFx9eveFVHuz+IA6qcGxXdkEmxibvYnSLstSTBCasoUrJSynlq5TaqJTKVEpFK6UmVVDPWSn1hlLqjFIqWSn1uVIquGZDFvWJk2cBLQYk4tTIRG6aIzGz78KclWXtsEQNUUpxe/vbeTD8ZbTJBwBTensKskJL1JMEJupaVXterwF5QBAwGViulGpfTr15wLVAONAUSAFeqYE4RT3m6GamxQ1JhT2wyEjiHngQXSCjOOzJ61+YyDh2L3nJ15Jz6mbKm3KuNldrFqK0SpOXUsoduBl4XGudobX+EdgM3FpO9VbAl1rrM1rrHOBDoLwkJ2yRe2CFLzm4mAkZCkYvLzK++44zzz4rNzLbkSXjwnE1epB7ZjSY3Uq9mo+D1x5cHVWtrdYsRGlVuUm5LZCvtT56Qdk+oF85dVcBy5RSTYFUCntp26sbpKgnKrkB2RlodtMfnLxjGinr1uMY3Ay/6dPqJjZRq4qXWpm+Zo9lsEYx56BtOPnuppPvKTqG9LVShKKhqUry8gDSS5WlAZ7l1P0biAHigALgADC3vI0qpWYBswCCg4OrtZCdraqozTExMXUcSQ0KDMRp/oPkPruIhCVLSHN2wrHf+b9z5DjbrhAnWDy8BQu2Hic3v7BXbXT/CyffwtWZ9yX/xC2f3cJjHR/DM7e8rwf7Zi/Huarqor0XW3KlKskrA2hUqqwRcK6cuq9R+Ae4H5AJLKCw53VN6Ypa6xXACihcz8vPz69BrYVTrLw22/zncOutJOXlkbBkKXnP/5fGV16JW/fzE8DafPsug720OSQEAgIDLD0wp7wrGND0Fr6J/wSA+Ox4Htj7APddeR+TQsod12XX7OU4V5U121uVARtHAQelVJsLyjoBB8up2xlYo7VO1lrnUjhYo4dSyr/akQqb4jttGj6TJqJNJmLmzCX32HFrhyRqSInVmqdew0uDnmTJdUssk/xm52fzXORzvPLnK5hLzB8mRM2pNHlprTOBDcBCpZS7Uqo3MBp4t5zqvwO3KaW8lFKOwN1AvNZaxs82MEopgh55BI/rr8eclkbMrFlyD5gd6RXmz08P3UCvsMK/S4e2Gsq64eto0aiFpc6K/Su477v7yDRlWitMYceqOlT+bsAVSADWA3dprQ8qpfoqpTIuqPcgkEPhta+zwHBgbA3GK2yIcnAg+IWluHTogCk2lujbp2I+e9baYYla0tqnNe8Pf5/eTc+viPxdzHdM2TaFmHQbvo4r6qUqJa+i04BjtNbuWuvmWut1ReW7tNYeF9RL0lpP1loHaq29tdZ9tNa/1Vbwov4zuLkR8sZynNu0IS8qiuz77pdTiHbMy9mLVwe8yk0hN1nK/kn9hyd/ftKKUQl7JNNDiVrn4O9Pi3ffwbVzZ3RCAtGTJ5MdWd4lU2EPHAwOzGgzg2f7PIujwRFvZ28W9lpo7bCEnZHFKOuR+ryWV3UZvb1pvnoV/9x5JwW/7+Hk7bfT7PXXcb+mh7VDE7VkVNgoWjZqSV5BHs08m1k7HGFnpOdVj9jCel7VYXBzw+Xpp2k0fDjmzExiZs7k3NdfWzssUYvCA8Lp3rh7mfLdcbtJzyt9+6gQVSfJS9Qp5ehI06VLCofR5+URe+88Uj/dYO2wRB3ae2Yvc7+dy63bbiXmnAzkEJdHklc9YkvreVWHMhgIevxx/O++G8xmTj36KEmrVslciA3Aubxz3PfdfZjMJo6lHWPy1sn8mfCntcMSNkiSVz1ia+t5VYdSioB77yHokUcASFiylDPPLpLZ6O2cp5MnD/V4CCeDEwApuSlM/3I6W49ttXJkwtbIgI165M0337R2CHXO97ZbcfD3I/4/D5Hy3nvknzlN0yVLMLi4WDs0UUuGhw6nqUdT5n03j+ScZExmEw/teojo9Gju6nQXSsmK3KJy0vOqRxrKacPSGg0fTsiqtzA0asS5r77m5NQ7yE9JsXZYohZ1DuzM+8PfJ9Tr/KKWy/ct5+EfHyavIM+KkQlbIclL1AvuPXrQ8v33cGjShOyICKInTCTPlmfXF5Vq5tmMd4e/y7VNrrWUbT22ldlfz7aMRNwdlUjvxd/KCs2iDEle9ciKFStYsWKFtcOwGuc2bWj5wQc4X3kledHRnJgwkewDB6wdlqhFjZwa8drA17il7S2Wst9P/87UL6ay8+9TTF+zh7jUbKav2SMJTJQg17zqkTvvvBOgQZ46LOYYFEiLqyOJS8kn80wS0RNuIbhXCp7BuecruQdWujCmsB2OBkee6PkEwR7BLNu7DIBuPsO48519loUvs00FTF+zh1VTu1smAxYNm/S8RL1jNCUQ0i8Zr5ZZ6AIDsT/6knzUHctI+swEq8Ynap5SihkdZ/Bc3+cY2mwS7+0IKbNic3ECkx6YAEleop5SBmhyTSr+7c+BVpzZ68XpP7yQ5aHsm6/uyZbvu5RJXKh8QBKYOE+Sl6i3lIKAjudo2jMFZdCk/uNOzA++FOTJUGp7Nf/j/eUkrjzcWryJk/9XgCbbVMD8j/dbJT5Rf0jyEvWeV8tsmt+QiNG5gMwzLpz42p+8kyetHZaoBUvGhePqaLygpADX4PUYXWNwDvgG58YbcXUsrCcaNklewia4+ZtoOSgRZy8TeemOnPjXeLJ+/93aYYka1ivMn1VTu59PYKqg8FHEyec3uvX4nK4tPCrYgmgoJHkJm+HkUUCLgYm4N8mhIDWV6GnTSd2w0dphiRpWIoFpJ7JjbseU2sXy+r7kn5j11SzSctOsGKWwNklewqYYHTUhfZPxvf02MJk49cgjnFn8PDo/39qhiRpUsgdmRCVOYHDweMvrfyb8ye3bb+d05mnrBSmsSpKXqH/cAy/6svIMJOjhh2n81FPg4EDymjXEzJolU0rZmeIEFuztyqqpPXhh4GPM7z7f8npUWhRTtk0hKjXKilEKa5GblEX9U8UbkH0mjMe5dRix8+4jc/fPnBj3L5q99houV7St5QBFXekV5s9PD91geX5b+9vwd/Xn0Z8eJd+cz5msM9y2/TZeHfAqXQK7XGRLwt5Iz0vYNLfu3Wn1yce4tG+PKTaWExMnkv7lDmuHJWrR8NDhvDbgNdwc3ABIz0tn2d5lsh5cAyPJS9g8xyZNaPH+ezQaNRKdlUXcvHkkLFuGNssdzfaqV9NerB66Gl8XX1p5teJ//f8nS6k0MHLasB5piOt51RSDiwtNn38elyvbkbB0KUnL3yD38BGaLvkvRk9Pa4cnakF7v/a8N+w9HAwO+Lj4WDscUcek51WPNNT1vGqKUgq/aXcQsnIFBi8vMr7/nuO33ELOX39ZOzRRS0IahdDEo0mZ8v1n98tpRDsnyUvYHY/evWn18Uc4X3klpuiTnBg/gbTNm60dlqgjX5z4ginbprDwl4UUmAsqf4OwSZK86pGGvp5XTXJq3pyW69fhNWYMOieH+AX/4dTTT2POk1V67VlkYiQP73oYjeaTo58wf+d8WZnZTknyqkfuvPNOy5peovoMrq40eW4RjRc+jXJ0JHX9B0RPuRVTfLy1QxO15ArfKxjacqjl+VfRX3H3N3eTacq0YlSiNkjyqkdmzpzJzJkzrR2GXVFK4fOvf9Fi3TocmzYlZ/9+jt90Mxk//mTt0EQtcDQ48myfZ5nSboql7NdTvzJzx0yZTsrOSPKqR+S0Ye1x7diBlp9+gnufPhSkphIzcyZnX30NXSDXROyNQRlYcPUC7ulyj6XsQOIBpn4xlYQsWcjUXkjyEg2Gg48PIW++gf/cuQAkvvoqJ6fPIP/sWStHJmqaUopZ4bN47JrHUBTe//VP6j/cvv12Ys7FWDk6URMkedUjf/zxB3/88Ye1w7BrymgkYO4cQt5aidHPj6xffuHYmLFk7t5t7dBELRh/5Xie6/scDqrwltbYjFhu3347f6dUbQoyUX9J8qpHunfvTvfu3a0dRoPg0bs3rTZuwO2aayhISuLk9BkkvPSSzE5vh24MvZFlNyzD2egMQGJ2IsfSjlk5KlFdkrxEg+UYGEjz1avwv2cuKEXSG29ycuodmM6csXZoooZd1+w6lg9cjoejB49f+zhDWg6xdkiimiR5iQZNGY0EzJlD87ffxhjgT9aePRwfM5aMH36wdmiihl3d+Gq2jN3CuLbjrB2KqAGSvIQA3K/pQehnn+HeuzcFKSnE3Dmb04sWYc7NtXZoogb5ufqVKUvLTeP13z6j9+Jv2R2VaIWoxOWQ5CVEEQc/P0JWriDggX+DgwMp77zLiX+NJ/dvubhvrzJNmUzZMoPlhx8nQX3L9DV7JIHZCJlVXjRMS9pAZtl7fhTgD7gPDyRuXyi5f/3F8VvGEfTQf/CeMEGW3bAzC759lhMZRwBwabyJ3IRcpq+BVVO70yvM37rBiYuSnpdomMpJXBdydUsgdMMGvG66CZ2by+mnFxI7Zy75KSl1FKCobbujEvl2d3cKskMsZc6BX1DgvZXpa36XHlg9J8lLiAoY3N1puuhZgl98AYOnJxnffsvxUaPJ+EmmlrJ1u6MSmb5mD9m5zmSdnEF+ZqjlNWf/7zH7bmD6mt8kgdVjkryEqESj4cMJ/Wwjrt26kX/2LDHTZxQO5sjJsXZo4jLN/3g/2aaiqcHMzmTH3EH+uXaW1518f0H7f8CDH/9ppQhFZaqUvJRSvkqpjUqpTKVUtFJq0kXqdlVK7VRKZSilziil5tVcuEJYh2NwMC3WriFg3r2WwRzHb76F7MiD1g5NXIYl48JxdTSeL9COZMdOwZTWyVLk6P0nYe03yJIq9VRVe16vAXlAEDAZWK6Ual+6klLKH/gCeBPwA1oDO2omVCGsSzk44H/XXbRcvx6n0FDyoqI4MWECicuXy8wcNqZXmD+rpnYvmcAwkhM/nryUHpaSiOQfmfvNXLJMWXUfpLioSpOXUsoduBl4XGudobX+EdgM3FpO9X8DX2qt39da52qtz2mtD9dsyEJYl2vHDrTa8Ck+t94K+fmcXfYy0ZOnkHfihLVDE5eg/ARmwJB0C4OCx1tKMvNlLbD6qCo9r7ZAvtb66AVl+4AyPS+gJ5CslNqtlEpQSn2ulGpeE4E2BLKel+0wuLjQ+NFHaL56FQ6NG5O9bx/Hxt5Eyvr1aK2tHZ6ootIJzNXRyKqpV/PCgEeZ23kuV/hcwesDXsfN0c3KkYrSVGW/aEqpvsDHWuvGF5TNBCZrrfuXqnsUCAQGAQeA/wLdtNa9y9nuLGAWQHBwcLctW7bg51f27nd7lpSUJG22kqbr+mHMTqrw9QJXP+InVW2KKH3uHLmvvkr+t98BYOzaFecH/o0hMLBG92NL6stxrqq9sRks+jaGR24IoWszD0t5bkGuZULfytham6urLtobEhJS4Y2VVUleXYCftNZuF5Q9APTXWo8sVXcfsFdrfUfRcz8gEfDWWle4jGn37t31xo0bCQkJqaiKXYqJiZE225H07ds5/fRCClJTMbi7E/jQf/COnEGl9zU/ZX8r/Nrzcd4ctZkugV0I8SzZPntuc3nqqL0V/vZU5bThUcBBKdXmgrJOQHnDrPYDF2ZDOX9yCWQ9L9vWaNgwQrd8jueggZgzMzn9+BPE/OCLKUvuSLEXm6M28+iPj3L79tuJSo2ydjgNWqW/VVrrTGADsFAp5a6U6g2MBt4tp/rbwFilVGellCPwOPDjxXpd4jxZz8v2Ofj7E/zyyzRdsgSDlxeZp104tj2Q1GOuyKUw25aSk8KzvzwLwNnss0z9YiqHkg5ZOaqGq6p/Et4NuAIJwHrgLq31QaVUX6VURnElrfW3wCPA1qK6rYEK7wkTJXXt2pWuXbtaOwxRTUopvEaOIGzL53gEZ2M2GTj1mw8xO6UXZst8XHx4dcCruDkUXkFJzU1l+pfT2Xtmr5Uja5iq9JuktU7WWo/RWrtrrZtrrdcVle/SWnuUqrtcax2stfbRWo/UWsfURuD2SE4b2heHgACa9Umhac8UDE5mMk+5cGxbICn/uEkvzEZd3fhq3hr8Fo2cGgGQYcrgzq/uZHfcbitH1vDIn4FC1CKlwKtlNqHDEgp7YfkGTu/x5uR3fuSdM1a+AVHvdAzoyNtD38bPpXCkXU5BDnO/nctPCTLnZV2S5CVEHXB0NdOsTwrBvZIxOheQleDMsS8CSTriLrNz2KC2Pm1ZO2wtTdybAGAym3gu8jk2R222cmQNhySvekQpJetF2Rv3QMt/lYJGzXMIHZ6AV8ssdIEiIcKLExMmkvPXX1YMUlyOFo1asHboWlo0agGAGTOP/vgoHx750MqRNQyyGKUQtWn+32Xuh3EAmgKNdu7k1JNPkRMZyfGbb8Fv5gz8Z8/G4Fy1m2KF9TXxaMKaoWu486s7OZpyFKMy4u8qi1jWBel5CWElHtddR+jnn+MzaRLk55O0/A2OjxpN5i+/Wjs0cQn8Xf1ZPWQ17bza8X+9/48BLQaUW293VCK9F38ra4TVEEleQliR0cOdxk88Tov338OpdRh50dGcnDqV+IcfkVWbbYiXsxf/7fpfRoaNLPf14sUv41Kzmb5mjySwGiDJS4h6wK1bN0I3bCBg3r0oJyfSNm7k2PAbSdu8WSb6tRFGVXb0aGJ2Ig98/QzT1/xqWfwy21QgCawGSPISop5QTk7433UXrTZ9hluPHhSkpBC/4D/ETJ9B3smT1g5PXKK03DRu3TqDHXEfQtAaUCbLa5LAqk8GbAhRzzi3akXztWtI2/gZCc8/T+bu3RwbOQr/2XfiO306Biensm9a0gYyEyreqHsgzP+79oIWZSz7dT2xmYXzHzp4HsE1ZDXZsbeD2QU4n8BWTe1OrzAZ5HGppOclRD2klML7prGEbt+G1+hR6Nxczi57meOjx5D5yy9l33CxxFWV10WN+2J3G3ITr7c8d3A/jlvzlSjj+cUts00FzP94vzXCs3mSvISoxxx8fWn6/PM0X7MGp9BQ8o4f5+TUO4h7cD75Z89aOzxxEUvHdcKYOpycM8MsZUbXOFxbvIlyKJyr3NXRyJJx4dYK0aZJ8hLCBrj3vIbQzzYScN99KGdn0rdsIWr4jSSvW4cuKLB2eKIcxas0O5y7gZxTN6F14QQERucE3Fosx9UtWU4ZVoMkLyFshHJywn/2nYRu3YJ7v+swnzvHmYX/x4nxE8hOcrR2eKIclgSWeS05cRPQuvAr1+CUik/rlfh4S+/5cknyEsLGODVrRsgbbxD88jIcGjcmJzKSE1/5c+o3L/Jz5Ve6vilOYI45XcmOuR1tLvxD45wphWlfTuOvZJka7HLIT3o9Iut5iapSStFo8GDCtm7Bb+YMUJB6zJ1jW4uWXDFbO0JxoeIE1tixMws6v4CnkycA7fza0dKrpXWDs1EyVL4ekbW8xKUyuLsT+MADeMUu4vReL7LOOHN6jzepx9xo3C0NVz9T5RsRdaJXmD8/PXQDANe0eJuX/3yZ5/s+j7NR5rK8HJK8hLADzk19ad4ogXMxLpz504ucZCdOfOWPd2gWAZ3O4eArgwLqkyt8r+C1Aa9ZOwybJslLCHsw/28U0AjwyMwkcflyktasJfWYO+mJTQi491588vNRDvIrX5+tjlxNdn42d3e6W5ZHqoRc86pHZD0vURMM7u4EPvggoZs+w71XL8zp6Zx55hmO33Qzmb/9Zu3wRAU+Pfop//vjf7yx7w2e+eUZCsxyC8TFSPISwk45h4URsuotmr36Co7BweQePcrJ224n7t//xnTqlLXDExfQWvPNyW8szz86+hHzd84ntyDXilHVb5K86hGttcwgLmqUUgrPgQMJ3boF/3vvQbm4kL5tO1HDbyTxjTcw58qXY32glGLZ9csY1ur8bBxfRX/FXV/fRUZehhUjq78keQnRABhcXAi4+27Ctm7Bc+hQdHY2Z19axrEbR5D+1VfyR1M94Gh0ZHHfxUxpN8VS9vvp35n25TQSs2X2+dIkeQnRgDgGB9Pspf/RfM0anNu0wRQbS9w993Lyjmnk/HXU2uE1eAZlYMHVC5jXdZ6l7HDyYW7bfhsx52KsGFn9I8mrHunWrRvdunWzdhiiAXDveQ2tNm6g8ZNPYPTyIuuXXzg+diynFy6UFZytTCnFjI4zeLrX0xhU4Vd0zLkYbt12K0eSj1g5uvpDklc9snfvXvbu3WvtMEQDoRwc8Jk4kbAvv8BnyhRQipR164kaOozkd95Fm+QGZ2u6qc1NvNT/JctNzEk5Scz5eo4M4igiyUuIBs7o7U3jxx4l9LONhUPr09I4s2gRx8aMJWPXLmuH16Bd3/x63hz0Jp6OnjgoB/6v9//JjBxF5I5FIQQAzm3aELLqLTK++44zzz9PXlQUMTNn4d63L0H/WYDzxmGyWrMVdAvqxtpha4lKi6JXcC9rh1NvSM9LCGGhlMLzhhsI/fxzAhcswODpSeauXRwbPYbTO3PIz73ITfSyWnOtaePThqEth5YpP515usHezCzJSwhRhsHJCb9pdxD25Rd4T5wAWpPytwdRW4JI/ssd3TC/L+uV05mnmbJtSpmbmXdHJdJ78bfsjrLv4fWSvIQQFXLw9aXJk0/S6rONuAflYDYZOPOnF8e+CORcrAtye5h1ZJoyufubuzmTdYavor9i1o5ZpOWmsTsqkelr9hCXms30NXvsOoFJ8hJCVMqlbVtC+ifT7LoknDxN5J1zIPZHX05+50dOilw6r2suRhd6NO5heb43YS/jNk9m+ntfkW0q7BZnmwrsOoFJ8hJCVIlS4Nk0l9BhZwnqmobRyUxWgjPHvwwg/ldvTNnydVJXjAYj/7n6PzzQ7QFL2amsaAzNXsHgEmsps+cEJj9tQohLogzg2zaTsBFn8L0iAwyQdtyNqC2BnH3tNczZ2dYOsUFQSjG1w1Smt30UrY0AGBzO4dbiTYwehyz17DWBSfISQlwWo5MmqEs6YcMS8GyWjS4wkPjKq0QNHYZpxw602WztEBuEj38IJPvkdHSBKwDKYMK12bs4+vxsqZNtKmD+x/trLYaUlBSCgoKIiorixIkTKKXYs2dPld6bm5tL8+bNq1y/mCQvIUTVuAeWW+zkWUCzPik0H56Py1VXkX/mDLlLlnL85lvI/OWXOg6y4VkyLhwnUxuyTtyFOc8XAKU0Lo034Ry4BTDj6mhkybjwWoth0aJFDB8+nLCwMEJCQjh16hSdO3eu0nudnZ2ZP38+//nPfy5pn3KlVQhRNZXcgOwOtDSbSf/8c04tXUru4cOcnHoHHv37Ezj/QZzDwuomzgamV5g/q6Z2Z/qaPWSduBvXkLUYXQsn8VXGLFwdHVg1tTu9wvxrZf9ZWVm89dZbfP755wAYjUYaN258SduYPHkyDz74IAcPHqR9+/ZVeo/0vOoRWc9L2DplMOA1ejRuq1cTcN99GNzcyPj+e46NGs2pp58mPynJ2iHapeIE5mLwIit6Jqb09uRnhqESx7Fq6tW1lrgAtm3bhlKK3r17A5Q4bWg2mwkJCeGVV14p8Z6jR4+ilLLM5err60vv3r1Zv359lfcryUsIUeOUiwv+s+8kbMeXeI8fD1qTuv4DogYPIfHNFZhzcqwdot0pTmCuDq7kxE2G03ewamrPWk1cALt27aJbt24oVXb2FYPBwMSJE3n//fdLlL///vu0a9eOrl27Wsp69OjBDz/8UOX9SvISQtQaB39/mjz9FKGbN+He7zrMmZmc/d//iBo6jNTPPpNBHTWsOIEFe7uz6vbeJRKX1pqFPy9kZ+zOGt1ndHQ0TZs2rfD1KVOm8OuvvxIVFWUpW7duHVOmTClRr2nTppw4caLK+61S8lJK+SqlNiqlMpVS0UqpSZXUd1JKHVZKxV6snihJ1vMS9sq5dWuav/kmzd9ejXO7duSfPs2phx4uHNTx88+Vb0BUWa8wf3566IYyPa4V+1fw8dGPuefbe3jn4Ds1dokiOzsbFxeXCl8PDw+nY8eOlt5XcSKbPHlyiXqurq5kX8JtFlXteb0G5AFBwGRguVLqYlfV5gNnqxyFAGQ9L2H/3K+9llaffkKTxc/h0Lhx4aCOO6ZxctYsco7KSs61JcuUxcZ/NgJg1maW7FnCwl8WYjJXf802f39/UipZwHTKlCmW5PX+++/Tp08fWrRoUaJOcnIyAQEBVd5vpclLKeUO3Aw8rrXO0Fr/CGwGbq2gfitgCvBclaMQAOzZs+eS73UQwtYogwHvMWMI+2I7Af/+NwZ3dzJ37uL4mLGcevxxTGdkdvqa5uboxvvD36dTQCdL2SdHP+Gur+4iLTetWtvu0qULhw4dumidSZMm8c8///DLL7/w4YcfljllCBAZGVniGlhlqtLzagvka60v/LNoH1BRz+sV4BFAbrO/RHLaUDQkBhcX/GfNJOyrHfhMmgQGA6kff0LU0KGcffllCjIyrR2iXfFz9WPVkFXcGHqjpezX078yedtkjqcdv+ztDhkyhMOHD5N0kZGkzZo1o1+/fsyePZu0tDTGjRtXps6uXbsYOrTssi8Vqcp9Xh5AeqmyNMCzdEWl1FjAqLXeqJTqf7GNKqVmAbMAgoODL9pweyVtbhikzVVwx1TcBg0kd9VqCn78kcTXl5O0bj1Ot92Kw7BhKIf6f0uqrRznu1vcja/25d3j7wIQnR7NxC0TWdB+AVf7XV3l7RS319vbm86dO7N8+XJuv/124uPjATh9+jQxMTGW+sOHD2f+/PkMHTqUjIwMMjIyLK/98ccfpKSk0LNnzxLvCQkJqXD/qrKLdkqpLsBPWmu3C8oeAPprrUdeUOYORADDtdZ/FyWv97TWzSr7ELp37643btxYbqDp6ekkJCRgMlX/3Gx9k5+fj8MFv5TFPwx+fn7WCqnWlW5zRRwdHQkMDKRRo0Z1EFXtiomJuegvoT2qTpuz9u4l4b9LyI6IAMCpVSsCH/g3HgMGlDscu76wteP85YkveezHx8gpKLxtQaG4r9t93NH+jip9zhe294svvmDevHkcOnQIo9F4ybGMGzeOLl268Mgjj5R+qcJAqvLnzFHAQSnVRmtdfIt9J+BgqXptgJbArqKGOwFeSqnTQE+t9Ykq7KuE9PR0zpw5Q3BwMK6urvX6B/dy5OXl4eTkZHlefL2rXbt21gqp1pVuc3m01mRnZxMXFwdgFwlMVJ3bN+NpcUUC59xdSNjXiLzjx4mdew+u/rkEdk7HrYVPpbN9iMoNaTmEEM8Q5n03j9OZp9FoPjzyIePajsPTqcyJtYsaOnQoc+bMITY2tsxAjMrk5uYSHh7O/ffff0nvq/Sal9Y6E9gALFRKuSulegOjgXdLVY0EQoDORY8ZwJmi/8dwGRISEggODsbNzc3uEpeomFIKNzc3goODSUiQi/cNTmYCSkGjkBzChicQ1DUVo3MB2YnORH8dQOyXJnKPX/41GnHeVX5X8cGNH9A1sCuuDq68fMPLl5y4it17772XnLigcG7Dxx9/HFdX10t6X1VPJN8NrAYSgCTgLq31QaVUX2C71tpDa50PnC5+g1IqGTBrrU+Xu8UqMJlMl9wgYT9cXV3t8nSxqLrC5Vey8GqVTdJhD5L/cudcrCvnRo7C51/j8L/7bhz8a3cGCXvn5+rHW4Pf4q+Uv7jC9wprh1NlVbrPS2udrLUeo7V211o311qvKyrfpbX2qOA931fleldlpMfVcMmxF8WMjprA8HOE3ZiAd2gmmM2krFtP1OAhnH39dcxZWdYO0aY5Gh3p4N+hTPm3J7/lwyMf1ss5V2V6KCGEzXB0M9OkRxqhmz7Do39/zFlZJL78Cv8MGULKBx+i8/OtHaLd+CflHx7e9TDP/PoMj/30GDn59Ws+SkleQgib49ymDSFvLKf5O2tx6diRgrOJnH7qKY6NHEX6V1/Vy56CrVm+bzlZ+YU92s1Rm7l1+63EnLus4Qu1QpJXLWrfvj3ff/99ne0vJCSEP//8s872J4S1uffoQcuPPiT4pf/h2KI5ecePE3fPvURPnETWH39YOzyb9myfZxnTeozl+ZHkI4zfMp6dsTvZHZXILe8cZndUotXis+vktTsqkd6Lv63VD/jHH3+kV69eeHl5Wdak+f333wE4ePAg/fv3r7V9L1u2jPvuuw8oXIb71KlT5Q6zb9WqFbGxZedIjo+Pp1mzwsuSLVu2JDAwkMzM87MavPXWW7UavxA1QSlFo6FDCduyhaAnHsfo50d2RATRk6cQc/cccv/5x9oh2iQXBxcW9lrIk9c+iaPBEYBzeeeY880cZm5+ljPncpm+Zo/VEpjdJq/dUYlMX7OHuNTsWvuA09PTGTFiBPfccw/JycnExcXx5JNP4uzsXOP7Ks+mTZsYM2YMAAcOHKB169blzu48cuRINm/eXKZ827ZtJaZjKSgoYNmyZbUWrxBV4h54Wa8rR0d8J00i7Msv8Z8zB+XmRsa333Js1GjiH3sM0+nLHvjcYCmluKXtLawdupbG7udXRzb4fo1ryBqyC85ZLYHZZfIqTlzZpgIAsk0FtfIBHy2aBXvixIkYjUZcXV0ZPHgw4eHhQGFv5uuvv7b8f8mSJYSHh+Pu7s706dM5c+YMw4YNw9PTk4EDB5KeXnoWroqlpKRw6NAh+vbtC8D+/fvp0KFwtFBWVhaTJk3ipptuIiMjg1GjRvHZZ5+V2ca2bdsYPny45fn8+fNZunQpqampl/NxCFEz5v8NT6VV/KjkBmWjhzsB98yl9Zdf4D1xAhgMpH3yKVFDhpKwdCkFadWbiLYh6hjQkQc7LMec1dpSZnA+jVLmWvt+rYzdJa/SiatYbXzAbdu2xWg0cvvtt7N9+/ZKlwX49NNP+eqrrzh69Ciff/45o0aNYtGiRZw9exaz2cyHH35Y5X1v27aNIUOGWKZiOXDgAB07duT48eP07t2bK664gk8//RQPDw/69etHREQEaRf80ppMJnbu3MmgQYMsZd27d6d///4sXbr0Ej8JIeofh4AAmjz5JGFbPsdz2FB0bi5Jb63in0GDSVq1SlZzvgS7oxK57/1/yIyeRm5if7TZSHbsFHRB4Z1S1khgdpW8KkpcxWr6A27UqBE//vgjSilmzpxJQEAAo0aN4syZM+XWv+eeewgKCiI4OJi+ffty9dVX06VLF1xcXBg7dix//fVXlfe9adMmRo8ebXm+f/9+Tp06xfXXX89TTz3Fk08+ablPytHRkQEDBrB9+3ZL/Z07d9KpUyc8PUveTb9w4UJeeeUVzp6V5diEfXBq2ZJm//sfLT/+CLdrrsGcnk7CkqWFqzl/+qkMr6+C+R/vL/peNZB3diiZUfMx5zQvUSfbVMD8j/fXWUx2lbzOf8AVq+kPuF27dqxZs4bY2FgiIyOJj4+3DKIoLSgoyPJ/V1fXMs+VUlWe1zAyMpKOHTsChXMBRkZGsnHjRu66664SSa1YeHg4Bw4csDwvfcqwWIcOHRgxYgSLFy+uUhxC2ArXjh1pvuZtQlauPL+a86OPcWz0GM59/bUMr7+IJePCcXU8P+GuzvcuU8fV0ciSceF1FpNdJa/SH3B5avMDvvLKK5k6dSqRkZGX9X6j0Yi7u3uV6g4aNIgdO3YAcLxonrevv/6aF154odwFLXfs2MHgwYMtzytKXgBPP/00K1eutEyMK4S9UErh0bcPrT79hKZL/otjs2bkRUURO/ceoidMJPO336wdYr3UK8yfVVO7V/j96upoZNXU7vQKq7upuuwqedX1B3zkyBFeeOEFyzD0mJgY1q9fT8+ePWtk+xczatQoNm3aBBSeMgwPD6djx46sWLGCsWPHcurUKUvdlJQUDh8+TJ8+fYDCZJebm1thL69169aMHz+el19+udbbIYQ1KIMBr5EjCdu2laDHHiscXr9vHydvu52Ts2aRc+SItUOsdyr6frVG4gI7S15Qtx+wp6cnv/76K9dccw3u7u707NmTDh068MILL1zW9nJycjhx4gQAw4YNY9GiRRXW7devH/v27SMtLY0DBw5YRjiOGTOGWbNmMWbMGHKKLkhv3bq1xOCOrVu3VtjrKvbEE0+UuOdLCHuknJzwnTK5cHj9PXMxuLuTuXMXx8feRNyD88mLqT8zStQHpb9frZW4oAqLUdaFihajPHz48GWvbXXh4A1rfsC1adKkSYwcOZKJEydetN64ceOYNGkSY8eOBQpXNJ07d26lCaw2VGU9rwtV52egvrC1RQprgk20eUkbyCy55E5+joGkQx6k/OOONitwcMDnX+Pwmz0bx8CL339mE22uIbujErn/g738b0LX2v5erXB2brvreRUr/gsh2NvVLhMXwLRp00hOTq60npubW4nrXf379+f666+vzdCEqP8yy64V5+BiJqhrOqHDE/BqmVVi9vqEF16Ue8SK9Arz55Pb2ln1e9VukxcUfsA/PXSDXSYugIEDBzJnzpxK661du7bEQJAFCxbIOmlCXISTRwFNe6YWzl4/cAA6J4eklSv5Z9BgElesxJydbe0QGzy7Tl5CCFEdzm3aEPLqq7T8YL3lHrGzL77IP4MHk7xuHTovz9ohNliSvIQQohKunTsX3iO26i1c2ren4GwiZxb+H1HDhpP62WfogovfXypqniQvIYSoAqUUHr170/KTjwl+6SWcQkMxxcVx6qGHOTZqNPm7dsmNznXIwdoBCCGELSlcgmUInoMGkrb5cxJffZW8qChY+H+c+HQDAffdh3uf3oXTs5UzorEE98BKJxoW5ZPkVQ1nzpwhISHBMk1TVZw7d67SOQy9vLxo06ZNdcMTQtQiZTTiPXYMXjcOJ+WTT0h49TVyDh4kZuZMXLt3I3DePNwulrjg4olNXJQkr2rIzs6+5FF7rq6uXHnllZbnsbGxaK1L3B/i6OhYYzEKISrgHlh5r6gKlJMTvpMmkXH11bjv3EnSipVk7/mD6Ftvwz3Il4Dwc7j6mWooaFFMklc1ZGdn4+XldUnvcXBwwMPDw/LcZDLh5eVVokwIUQdq+HSdcnHBb/p0vMePJ3ntWpLfXkPmGcj8ygWPpjkEdEzHxUdmsK8pMmDjMmmty/S8EhIS2Lt3b5VuHAYwm83k5uaWu/qxEMI2GT08CJgzh9Zff4XfVedQDmYy4l04/mUgsT/5kJsmfYaaIMnrMuXm5mI2m3Fzc0NrzcmTJ4mPj6dt27b4+voC0L59e77//vsKt1E896DcMCyE/TF6exMYfo7WIxLwvSIDZdCci3Hl2PYA4nZ7k5suSaw6JHldpuzsbAwGAz///DPdunWjffv23HDDDQwZMoTff/8dgIMHD9K/f/+LbgMuP3ktW7aswrXDKtOqVSvLbPgXio+Pp1mzZgC0bNmSwMDAEhP0vvXWWxdtkxCiJAcXM0Fd0gkbcQaf1plggPSTboVJ7GdvcouWNBKXxv5Sfx0NTc3KyiI3N5dRo0bx5JNPMnfuXMxmM7t27cLZ2ZnExERiY2NxdHTEbDbj7OxMWFiYZWZ3KExejo6OODhc3mHYtGkT8+fPJzY2lmbNmmEymUhPT8fPz6/S944cOZLNmzdz9913lyjftm0bQ4cOtTwvKChg2bJlPPLII5cVoxCikKObmcbd0/Brl0HiIQ9Sj7mRHu1G+o0j8Bo5Ev+778KpRQtrh2kz7K/nVUdDU7Ozszl27Bhaa+bMmYOzszOurq4MHjyY8PBwsrOzGTVqFKdOnaJjx45cf/31LFq0iPDwcNzd3Zk+fToxMTHcc889eHp6MnDgQFJSUqq8/5SUFA4dOsTgwYNp1qwZq1atYsiQIdx33334+PjQtm1bDh06xLJly2jevDn+/v5s2LDB8v5Ro0bx2Wefldlu6UUq58+fz9KlS0lNTa3OxyVEw1TOiEVH9wKaXJ1G2I0JeLc1g1KkbdpE1PAbiX/kUfJOnrRCoLbH/pJXHcnOzqZnz54YjUYmTJjA1q1bSySfrKyswpsUKbw+BrBlyxa++uorjh49yueff860adN46KGHOHv2LGazucTijyaTiaioKA4dOkRkZCTp6ekAHDt2jOjoaFauXEnfvn2Jjo4mPT2dPXv2sHfvXq677jp++OEHrrrqKoYNGwZAVFQUjz/+OM8884xl+/369SMiIoK0C2bJNplM7Ny5k0GDBlnKunfvTv/+/Vm6dGktfIpC2Ln5f8NTaeU+nJYm02TzX4R9sR2vm24CIG3DBqKGDSf+oYfJK1rbT5RPktdlKCgoIDc3l8DAQHbu3InZbGbGjBkEBAQwatQozpw5Q3Z2Nvn5+Rw/fpyDBw9iNBq59957CQoKIjg4mD59+tC+fXuuvvpqXFxcGDt2LH/++SdQOJLx6NGj+Pr6ctVVV3HVVVfh5uYGFCZFo9HInj17mDhxItnZ2bi5uXHo0CFmz57NpEmTLKsqd+zYkXnz5uHo6EiHDh3Izz8/TNfR0ZEBAwawfft2S9nOnTvp1KkTnp6eJdq7cOFCXnnlFc6ePVsHn64QDYtTSAhNFz1L2LateBWtuZf22WeFPbH//EeuiVVAktdluHCgRXh4OOvWrWPLli38+OOPxMfHc++992I0GnF0dKRVq1aEhoZSUFBAUFCQZRtOTk74+vpahsm7urqSkZEBQGpqKs7Ozvj4+ABgMBhwcHDAbDZjMplo2rQpkZGRtG/fnoKCAhwcHNi/fz/XXnutZfDHoUOHGDFihGV/hw4dKnFzNEB4eDgHDhywPC99yrBYhw4dGDFiBIsXL66Jj08IUQ6nFi1o+twiwrZvw+vmm8BgIG3TZo7dOIK4+QvIPXas4jcvaQNPeVX8WGJ/M/ZI8roMpQdaeHh40KxZM5ycnJg0aRKRkZEl7t1ydXUtM2FnQdEs1OWNNMzOzi73puWcnBzc3d0xGAwMGjSIbdu24eLiQnR0NHl5ebRq1QqDofCQRkRE0LlzZ8t79+/fX+I5wI4dO0osUllR8gJ4+umnWblyJXFxcRf5ZIQQ1eXUvDlNn32WsC+24z1uHBgMpH/+OcduHEHs/feTU970cg1wGipJXpchICCATp06ceTIEV544QViY2MJCgoiICCATz/9lM6dO5dIXomJiZakUszZ2ZmmTZuWO9LQwcHB0ruDwmtRUHjKsDjZjRo1is2bN+Pm5sa+ffto3749zs7OAKSnpxMdHU14eLhlG/v27aNTp06W5ykpKRw+fJg+ffoAcPz4cXJzc2nXrl25bW7dujXjx48vcV1OCFF7nJo1o8n/LaT1l1/gPX48ODhwbvsXHB89hpi755C9f7+1Q7QqSV7V4Onpya+//so111yDu7s7PXv2pEOHDixYsIDk5GRMJhPR0dGYTKZLGg7v7+/P1KlTmTdvHgcPHuTcuXMAlutbUDjgIjIyEpPJxL59++jSpQtms5mDBw/y888/07p1a0tds9lMZGRkiZ7X1q1bGTJkiGXo/tatWyvsdRV74oknStzzJYSofY7BwTR5+ilaf7UDn1tvRTk7k/Htt5z413hOTptOVtF9pQ2Nqg/rz3Tv3l1v3LixxOS0AIcPH66wJ1ChBrQEwaRJkxg5ciQTJ0685PeOGzeOSZMmMbboAvHw4cOZO3dupQmsuvLy8nBycqpy/cv6GahnYmJiyvxs2ztpc+3JT0wkee1aUt5fhzkrCwBX/1z8r8rAvUkuRYOcy3oqrYIXLk8dtbei1tjhTcp2kpiqYtq0aZUur1IRNze3Ete7+vfvz/XXX19ToQkhaomDvz+BDzyA3/TpJL/3Psnvvkt2YhoxO51x9jbh1+4cjUJyUHZ+Xs3+klcDMnDgQAYOHHhZ7127dm2J5wsWLKiJkIQQdcTo7U3A3Dn4Tp1K6vS2JP3lTm6qI/E/+3J2fz5+7TLwapWFwVj5tmyRJC8hhLBhRg93/Npl4NM2g7TjbiQd8cCU4cDpPd6cjfTE94pMfMIysbccJslLCCFsnXsghswEfFpn4R2axblYFxIPeZKb6sjZfY1IOuSJt/cSfG+7FcfGja0dbY2Q5CWEELbugmv9CmgEeGpN5o8/kbRyJVm//Uby6tUkv/MOXjcOx3faNFyuuMJq4dYEO7+kJ4QQDZNSCo++fWjxzlpafvwxjYYPA7OZtE2bOT56DCenzyBz9+4yEyjYiiolL6WUr1Jqo1IqUykVrZSaVEG9+UqpSKXUOaXUcaXU/JoNVwghxKVy7diB4BdfJGzHl4X3irm5kfnTT5ycNp3jY28idcNGzHl51g7zklT1tOFrQB4QBHQGtiql9mmtD5aqp4DbgP1AGLBDKRWjtf6ghuK1K6XvedqzZw9QOJO7EELUNKdmzWj86CMEzLmblA8+JPm998g9coRTjzxCwgsv4DNxIj4TJ+BQ3pqApe6hLXOHVx3fQ1tpz0sp5Q7cDDyutc7QWv8IbAZuLV1Xa/1frfVerXW+1vovYBPQu6aDFkIIcfmM3t74z76T1t9+Q5PnnsP5yispSEoi8dVX+af/9cQ/8mjZORTr2fyJVel5tQXytdZHLyjbB/S72JtU4WJWfYE3K3h9FjALIDg4mKSkpDJ18vPzybOxruyluHCJkgs1xDZfrH5MTEwtRVM3yvvZtnfSZhvSvRsO3bqi9u/HtGEDBT//QtqGDaRt2ICxc2ccx4zG2LMnzauwqZr+Xb3YDB5VSV4eQHqpsjTAs5y6F3qKwp7d2+W9qLVeAayAwumh/Pz8ygSakZFxSVMJ2aLy2tcQ21wRBwcHu5hmyB7acKmkzTameXMYMYK86GiS33uftE8/pSAigoKICByaNiHR3wPvsCwcnM0VbqIu21+VARsZFI68vFAj4FxFb1BKzaXw2teNWuvcyw9PWFtISIhlkUwhhP1zatGCxo8+QusfvifokYdxatGC/PhTnN3fiH82BRH/izfZyY7WDrNKyeso4KCUunA1s05A6cEaACilpgEPAQO01rHVD7F+a9myJV9//bW1wyghJSUFpRTXXnttifLZs2dz//33V/i+ZcuWcd9995XYzqlTp8qdGLdVq1bExpY8vPHx8TRr1gwo/FwCAwNLzEL/1ltvMWjQoMtpkhCijhk9PfG97TZCt28jZOVKPJrmoM2QdsKNEzsCOP6VP+finK0WX6XJS2udCWwAFiql3JVSvYHRwLul6yqlJgOLgEFa64ss+ykudKnXgSoTERFB48aNOXToEKdPn7aU//nnn2UWpLzQpk2bGDNmjOX5gQMHaN26dYm1yYqNHDmSzZs3lyjbtm0bQ4cOtTwvKChg2bJll98QIYTVKYMBj759CLkumbAbE/C9IgODo5mcJCfy0q03z0VVb1K+G3AFEoD1wF1a64NKqb5KqYwL6j0D+AG/K6Uyih5v1GzItmHx4sWEhYXh6enJVVddxcaNG0u83rJlS5YuXUp4eDju7u789ttvTJ48mX79+jFu3DjGjx/PY489ZqkfHx/PzTffTEBAAK1atbroopARERF0796dQYMGsWnTJqAwkRw4cIAuXbqU+56UlBQOHTpE3759LWX79++nQ4cOQOFCmJMmTeKmm24iIyODUaNG8dlnn5XYRumVmOfPn8/SpUtJTU2t0mcmhKjfnDwLCOqSTpvRZ2h8dSpeoVlWi6VKyUtrnay1HqO1dtdaN9darysq36W19rigXiuttaPW2uOCx+zaCr4+CwsLY9euXaSlpfHkk08yZcoUTp06VaLOhx9+yNatW0lISGDcuHGMGDGCb775hokTJ5ZIdmazmZEjR9KpUyfi4uL45ptveOmll/jyyy/L3XdxD2vMmDGWBHPkyBHMZnOFa2Nt27atxOKUUNjz6tixI8ePH6d3795cccUVfPrpp3h4eNCvXz8iIiJISytcI8hkMrFz584SpwW7d+9O//79Wbp06WV9hkKIesQ90PJfg4PGJywLB2dd7ut1webmNnw94nWW71tepbo3t7mZp3o9VaLsqd1P8enfn1b4nrs63cXdne+uTohA4WKPxcaPH89zzz3Hb7/9xujRoy3lc+bMISQkhJ07d5Kfn8/cuXNRStGzZ0969Ohhqff7779z9uxZnnjiCQBCQ0OZOXMmH3zwAUOGDCmz74iICEaPHs0NN9zA7NmzOXfuHBEREbRv3x5Hx/IvtG7atIkJEyaUKNu/fz9KKa6//nqWLVtWInZHR0cGDBjA9u3bmTBhAjt37qRTp054epYchLpw4UJ69+7NvHnzLuHTE0LUO6VuQLb2gqM2l7xsxTvvvMOLL77IiRMngMJh/4mJiSXqFA9uiI+PJzg4mFatWlleu/CHIjo6mvj4eLy9vS1lBQUFJU7xFcvNzeXw4cN07twZHx8fevTowfbt2yu93hUZGUnHjh0tz7XWREZGcuzYMf7973+XSFzFwsPDOXDgABMmTChzyrBYhw4dGDFiBIsXL7b5FZGFEPWHTMxbC6Kjo5k5cyavvvoqSUlJpKam0qFDhzITYKqi9bqbNGlCXFxcidcvvNkvJCSEVq1akZqaanmcO3eObdu2ldl3ZGQkbm5uhIaGAlhOHf75558VXu8CGDRoEDt27LA8P378OABff/01L7zwgmXqqgvt2LHDshpzRckL4Omnn2blypXExcVVuH8hhLgUNtfzurvz3dU6rfdUr6fKnEqsLpPJRE5OjuV5WloaSikCAgIAePvtt4mMjKzw/ddeey1Go5EXXniBmTNn8v333/Pbb7/Rv39/AHr06IGnpyfPP/889957L05OThw+fJjs7GyuvvrqEtv6888/CQ8PtyTGUaNG8cQTT2AwGHj66acrjGHUqFE8//zzzJkzByg8ZRgeHk7Hjh1ZsWIFY8eO5bfffqNJkyZA4QCPw4cP06dPH44fP05ubm6FPavWrVszfvx4Xn75ZcsAECGEqA7pedWA4cOH4+rqanl89NFHPPDAA1x77bUEBQVx4MABeveueIpHJycnNmzYwIoVK2jSpAnvvfceI0aMwNm58B4Ko9HIli1biIiIoFWrVvj7+zNjxgzLYIkLRURElDg92LJlS1q2bElqaiqdOnWqMIZ+/fqxb98+yzYPHDhAeHg4UNh7mzVrFmPGjLEk6a1bt1oGeGzdurXCXlexJ554osQ9X0IIUS1aa6s/unXrpk+ePKlLO3ToUJkye5Kbm1vi+cGDB/XBgwe11lr36NFDr169uk7jmThxol63bl2V6t5yyy16w4YNWmuthw0bprdu3Vql95Vuc2Xs4WegvJ9teydttn911N4K84b0vOqJH374AV9fX9q2bcvatWvZv39/iRt+68K0adNITk6uUl03NzfL9a7+/ftz/fXX12ZoQghRgs1d87JXf/31F//617/IzMwkNDSUTz75xHJ9qa4MHDiQgQMHVqnu2rVrLf9fsGBBbYUkhBDlkuRVT8yaNYtZs2ZZOwwhhLAJctqwHtmzZ0+5Q9KFEEKUJMlLCCGEzZHkJYQQwuZI8hJCCGFzJHkJIYSwOZK8hBBC2BxJXkIIIWyOJC8hhBA2R5KXEEIImyPJq5patmzJ119/be0wSkhJSUEpxbXXXluifPbs2dx///0Vvm/ZsmXcd999l7XPVq1aERsbW6Y8Pj7esuhmy5YtadasWYnZ5d966y3L0i9CCFFVkrzqgfz8/BrdXkREBI0bN+bQoUOcPn3aUl7ZasqbNm1izJgxl7XPkSNHsnnz5jLl27ZtKzHBcEFBAcuWLbusfQghRDFJXrVk8eLFhIWF4enpyVVXXcXGjRtLvN6yZUuWLl1KeHg47u7u/Pbbb0yePJl+/foxbtw4xo8fz2OPPWapHx8fz80330xAQACtWrXi5ZdfrnDfERERdO/enUGDBrFp0yagMGkcOHCgwtWUU1JSOHToEH379rWUrVq1iiFDhnDXXXfh4+ND27ZtOXToEMuWLaN58+b4+/uzYcMGoHAxy88++6zMdkuvsPzvf/+bpUuXkpqaWulnKIQQFZHkVUvCwsLYtWsXaWlpPPnkk0yZMoVTp06VqPPhhx+ydetWEhISGDduHDfddBO7d+9m4sSJJZKd2Wxm5MiRdOrUibi4OL755hteeuklvvzyy3L3XdzDGjNmjCWhHDlyBLPZXOFqx9u2bbMsLlls37597Nmzh1tuuYXExEQ6duzIsGHDAIiKiuLxxx/nmWeeAQoXs4yIiCixQKbJZGLnzp0MGjTIUta1a1f69+/P0qVLL+HTFEKIkmxqVvnDV5b/xVvT2h05XO1tjBs3zvL/8ePH89xzz/Hbb78xevRoS/mcOXMICQlh586d5Ofns2jRIpRSdOrUiR49eljq/f7775w9e5YnnngCgNDQUGbOnMkHH3zAkCFDyuw7IiKC0aNHc8MNNzB79mzOnTtHREQE7du3x9HRsdx4N23axIQJE0qU7du3j4cffpgBAwYAcNVVV5Gbm8u8efMA6NChg+WUp6OjIwMGDGD79u2W7ezcuZNOnTrh6elZYrsLFy6kd+/elu0IIcSlkp5XLXnnnXfo3Lkz3t7eeHt7ExkZSWJiYok6xQMZ4uPjCQ4ORilleS0kJMTy/+joaOLj4y3b8vb2ZtGiRZw5c6bMfnNzczl8+DCdO3fGx8eHHj16sH379kqvd0VGRtKxY8cSZfv372fEiBGW54cOHSrz/Morr7Q8Dw8P58CBA5bnpU8ZFuvQoQMjRoxg8eLFFcYjhBAXY1M9r5roEdWF6OhoZs6cyTfffMO1116L0Wikc+fOaK1L1CtOVk2aNCEuLg6ttaUsJiaGsLAwoDCRtWrVir///rvSfUdGRuLm5kZoaCiA5dThmTNnGDt2bIXvGzRoEDt27KBNmzaWNuTl5dG2bVtLnYiICObPn295vn///hIJcceOHTz11FOW59u2bbNcEyvt6aefpmvXrjzwwAOVtkkIIUqTnlcNMJlM5OTkWB5paWkopQgICADg7bffJjIyssL3Fye4BQsW8Msvv7Bp0yZ+++03y+s9evTA09OT559/nuzsbAoKCoiMjOT3338vs60///yT8PBwSxIcNWoU27Ztq7TnNWrUKMvgDig8ZdixY0cMhsIfkfT0dKKjowkPDy9Rp1OnTkDhgI/Dhw/Tp08fAI4fP05ubm6F19hat27N+PHjLzrwRAghKiLJqwYMHz4cV1dXy+Ojjz7igQce4NprryUoKIgDBw7Qu3fvCt/v5OTEhg0b2LRpEzfccAPvvfceI0aMwNnZGQCj0ciWLVuIiIigVatW+Pv7M2PGjBKDI4pFRESUSFItW7akZcuWpKamWhJNefr168e+ffss29y3b1+J7ezfv5/WrVvj5uYGFA4iiYyMtNTZunVriQEfW7duLfeU4YWeeOKJEvd8CSFElWmtrf7o1q2bPnnypC7t0KFDZcrsSW5uboWv9ejRQ69evboOo9F64sSJet26dZf13ltuuUVv2LDB8nzYsGF669atZepdrM3lsYefgfJ+tu2dtNn+1VF7K8wb0vOqJ3744QdOnz5Nfn4+a9euZf/+/SVu7q0L06ZNIzk5+bLe6+bmxuDBgy3P+/fvz/XXX19ToQkhRAk2NWDDnv3111/861//IjMzk9DQUD755BOaNGlSpzEMHDiQgQMHXtZ7165dW+L5ggULaiIkIYQolySvemLWrFmWwQ5XXXWVlaMRQoj6TZJXPZKVlWXtEIQQwibINS8hhBA2R5KXEEIIm1Pvk5cuNSuFaDjk2AshKlKvk5ejoyPZ2dnWDkNYSXZ2doUTCQshGrZ6nbwCAwOJi4sjKytL/gpvQLTWZGVlERcXR2BgoLXDEULUQ/V6tGGjRo2AwlnXTSaTlaOpefn5+Tg4nD8ExbPOHz5sGxMQX47Sba6Io6MjQUFBlp8BIYS4UL1OXlCYwOz1CywmJqbE0ifF93fZcy+zdJuFEOJyVOm0oVLKVym1USmVqZSKVkpNqqCeUko9r5RKKno8ry5cpEoIIYSoAVXteb0G5AFBQGdgq1Jqn9b6YKl6s4AxQCdAA18Bx4E3aiJYIYQQAqrQ81JKuQM3A49rrTO01j8Cm4Fby6l+O/CC1jpWax0HvABMrcF4hRBCiCqdNmwL5Gutj15Qtg9oX07d9kWvVVZPCCGEuGxVOW3oAaSXKksDPCuom1aqnodSSulSoxCUUrMoPM0IkNG8efMkILFKUdsPf8pps51fJiy3zXZO2twwNLQ210V7v9Bal7s2VFWSVwZQerhfI+BcFeo2AjJKJy4ArfUKYEXxc6XUHq119yrEYzekzQ2DtLlhaGhttnZ7q3La8CjgoJRqc0FZJ6D0YA2KyjpVoZ4QQghx2SpNXlrrTGADsFAp5a6U6g2MBt4tp/o7wL+VUsFKqabAA8CaGoxXCCGEqPL0UHcDrkACsB64S2t9UCnVVymVcUG9N4HPgQNAJLC1qKwqVlRexe5ImxsGaXPD0NDabNX2KnuezUEIIYR9qtcT8wohhBDlkeQlhBDC5tRp8mqIcyReQpufUkqZlFIZFzxC6zre6lJKzVVK7VFK5Sql1lRS936l1GmlVLpSarVSyrmOwqxRVW2zUmqqUqqg1DHuX2eB1hCllLNSalXRz/M5pVSEUmrYRerb/HG+lDbby3EGUEq9p5Q6VXTsjiqlZlykbp0e57rueV04R+JkYLlSqrwZOC6cIzEcGAncWUcx1rSqthngQ621xwWPY3UWZc2JB54BVl+sklJqCPAQMABoAYQCT9d6dLWjSm0u8nOpY/x97YZWKxyAGKAf4AU8BnyklGpZuqIdHecqt7mIPRxngOeAllrrRsAo4BmlVLfSlaxxnOsseTXEORIvsc12QWu9QWv9GZBUSdXbgVVa64Na6xTg/7DBYwyX1Ga7oLXO1Fo/pbU+obU2a623UDgBd5kvNezkOF9im+1G0XHLLX5a9Agrp2qdH+e67Hk1xDkSL6XNACOVUslKqYNKqbtqPzyrKu8YByml/KwUT13popRKLDoF87hSqt6vqVcZpVQQhT/r5U1IYJfHuZI2gx0dZ6XU60qpLOAIcArYVk61Oj/OdZm8amSOxFqKrbZcSps/AtoBAcBM4Aml1MTaDc+qyjvGUP5nYy92Ah2AQAp75BOB+VaNqJqUUo7A+8BarfWRcqrY3XGuQpvt6jhrre+m8Hj1pXDCitxyqtX5ca7L5FUrcyTWc1Vus9b6kNY6XmtdoLXeDSwDbqmDGK2lvGMM5f882AWt9TGt9fGi004HgIXY8DFWShkonGknD5hbQTW7Os5VabO9HWeAou+lH4FmQHlnher8ONdl8mqIcyReSptL04Ct9TQvRXnH+IzWukFcNypis8e46CzIKgoHIt2stTZVUNVujvMltLk0mz3O5XCg/GtedX6c6yx5NcQ5Ei+lzUqp0UopH1WoB3AvsKluI64+pZSDUsoFMAJGpZRLBef73wGmK6WuUkp5Uzh6a03dRVpzqtpmpdSwomslKKWuBB7HBo9xkeUUnuYeqbXOvkg9uznOVLHN9nKclVKBSqkJSikPpZSxaEThROCbcqrX/XHWWtfZA/AFPgMygZPApKLyvhSeFiyup4D/AslFj/9SNJWVrT0uoc3rKRytlkHhhdF7rR37Zbb3Kc6PSip+PAU0L2pb8wvq/hs4Q+F1wbcBZ2vHX5ttBpYWtTcTOEbh6SRHa8d/Ge1tUdTGnKL2FT8m2+txvpQ229FxDgB+AFKLjt0BYGbRa1Y/zjK3oRBCCJsj00MJIYSwOZK8hBBC2BxJXkIIIWyOJC8hhBA2R5KXEEIImyPJSwghhM2R5CWEEMLmSPISQghhcyR5CVGLilZYPmDtOISwN5K8hKglRfMbdgAiqrENj6Il5XUVH7411gAh6jGbXSBNCBtwFeBMNZIXhb+jt5cquwvoBTxI4VxyxXK11snV2JcQNkPmNhSiliilbgPWAgO01t/W4Hb/oDAxemqt82tqu0LYEjltKETt6VL0b0RxgVLKWym1USmVo5SadakbLFrFtwOwXxKXaMjktKEQtaczEFN8Kk8p1Q34mMIlf3prrf+4jG22B5yAP2sqSCFskfS8hKg9nSnqdSml7gJ+Ag4D3S4zccH53tze6gYnhC2T5CVELVBKtQS8gb+VUuuAV4FngRHVHFTRtehf6XmJBk1OGwpRO4p7SPcAZmCo1vqr0pWUUgYgjcI/JHXRv0eBOVrrnyrYbj6Fq9oK0WBJz0uI2lGcvFZSOFy+YwX1rgA8gBCttQfgAxwB/le6YlGi6wQc1lrn1HjEQtgQ6XkJUTs6A4la6zlKKVdgiVLqhNZ6Q6l63YD44lOJWutcpdSfQNtyttmGwkQn17tEgyc9LyFqRxfOX5e6E/geeE8p1aNUvW7AfgBVKByYCSwtZ5tyvUuIIpK8hKhhSik/oBlFSUZrbQJuBk4AnyulWl1QvRtwvVIqFcikcHTieq31unI2LSMNhSgiyUuImlecZCw9JK11KnBj0dNtSimfomtYXYBJWmtvrbUb0BNYoJTqV8F2NbCv1iIXwkbI9FBCWIlS6koK7/sK0VrHXlCeCMzXWr9tteCEqOdkwIYQ1lM8WCMWQCnlDvybwhk0dlgzMCHqO0leQlhPN6CxUiqDwnvB0oHdFE4dFWfVyISo5+S0oRBCCJsjAzaEEELYHEleQgghbI4kLyGEEDZHkpcQQgibI8lLCCGEzZHkJYQQwuZI8hJCCGFzJHkJIYSwOf8PAAR7m+tIIYUAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 504x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize = (7,5))\n",
    "ax = plt.subplot(111)\n",
    "\n",
    "# simulation results\n",
    "plt.plot(temp_vals,np.array(sim_k_vals)/np.sum(Nelems), \n",
    "        label = r'Sim. $\\langle k \\rangle$/N',\n",
    "        markersize = 7.5,\n",
    "       marker = 'D',\n",
    "       linestyle = '')\n",
    "plt.plot(temp_vals,np.array(sim_m_vals)/np.sum(Nelems), \n",
    "        label = r'Sim. $\\langle m \\rangle$/N',\n",
    "        markersize = 7.5,\n",
    "       marker = 's',\n",
    "       linestyle = '')\n",
    "\n",
    "# large N analytical results\n",
    "k_avg_approx_vals = [k_avg(T, E0s, Dels, Evs, Nelems)/np.sum(Nelems) for T in Tvals]\n",
    "m_avg_approx_vals = [m_avg(T, E0s, Dels, Evs, Nelems)/np.sum(Nelems) for T in Tvals]\n",
    "plt.plot(Tvals, k_avg_approx_vals, label = r'Large $N$ $\\langle k \\rangle$/N', linestyle= '--', linewidth = 3.0)\n",
    "plt.plot(Tvals, m_avg_approx_vals, label = r'Large $N$ $\\langle m \\rangle$/N', linewidth = 2.0 )\n",
    "ax.axvline(x = kBTcrit_master(E0s, Dels, Evs, Nelems), color = 'k', linestyle = '-.', linewidth = 2)\n",
    "\n",
    "\n",
    "plt.legend(loc = 'best', fontsize = 12)\n",
    "# plot formatting\n",
    "ax.set_xlabel(r'$k_B T$', fontsize = 18)\n",
    "plt.xlim([-0.01,3.2])\n",
    "plt.ylim([0,1.1])\n",
    "plt.grid(alpha = 0.45)\n",
    "\n",
    "\n",
    "# Hide the right and top spines\n",
    "ax.spines['right'].set_visible(False)\n",
    "ax.spines['top'].set_visible(False)\n",
    "\n",
    "# Only show ticks on the left and bottom spines\n",
    "ax.yaxis.set_ticks_position('left')\n",
    "ax.xaxis.set_ticks_position('bottom')\n",
    "\n",
    "# increase label size\n",
    "ax.tick_params(axis='both', which='major', labelsize=12)\n",
    "ax.tick_params(axis='both', which='minor', labelsize=12)\n",
    "ax.text(kBTcrit_master(E0s, Dels, Evs, Nelems)-.2, 0.25, r'$k_BT_{crit}$', color='black', fontsize = 14.5,\n",
    "        bbox=dict(facecolor='white', edgecolor='none', pad=5.0))\n",
    "\n",
    "\n",
    "for i in range(4):\n",
    "    ax.text(temp_vals[img_key_list[i]], sim_k_vals[img_key_list[i]]/np.sum(Nelems)+.05,'('+'i'*(1+i)+')' if i<3 else '(iv)', fontsize = 14 ) \n",
    "    \n",
    "# plt.savefig(f'general_grid_assembly.png', bbox_inches='tight', format = 'png')         "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Total Notebook Runtime: 13.918 mins\n"
     ]
    }
   ],
   "source": [
    "print('Total Notebook Runtime: %.3f mins' % ((time.time()-nb_start)/60))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}