{
 "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 the Derangement Model \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. For the derangement model, particles can only exist \\textit{on} thhe grid. Each particle type has a collection of \"correct\" locations on the grid. A particle of type $j$ attaches to this correct location with a Boltzmann factor of $\\delta_j$. Here we want to use simulations of this system to affirm analytical calculations of the average number of correctly bound particles as functions of temperature. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Numerical representations of analytical work"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### System Partition Function \n",
    "\n",
    "The partition function for the system is \n",
    "\\begin{equation}\n",
    "Z_{\\boldsymbol{n}}(\\{\\beta\\Delta_k\\})  = \\int^{\\infty}_{0} dx\\, e^{-x} \\, \\prod_{k=1}^{R}\\left(e^{\\beta \\Delta_k}-1\\right)^{n_k}\\, L_{n_k} \\left( \\frac{x}{1-e^{\\beta\\Delta_k}} \\right).\n",
    "\\label{eq:znr}\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# delta function definition\n",
    "delta_func = lambda Del, T: np.exp(Del/T)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# parameters for the model going forward\n",
    "np.random.seed(42)\n",
    "Dels0 = np.random.randn(10)+2 # randomly generated delta from normal distribution\n",
    "Nelems0 = np.random.randint(1,10,10) # randomly generated integers "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.03200812850398603"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "## integrand of partition function\n",
    "integrand = lambda x, deltas, ms: np.exp(-x)*np.prod([eval_laguerre(m, x/(1-delta))*(1-delta**(-1))**m for m, delta in zip(ms, deltas)])\n",
    "\n",
    "# test function \n",
    "integrand(.40, np.exp(Dels0), Nelems0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "189563.4375182251"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "## partition function for \n",
    "ZN_betalambda = lambda deltas, ms: quad(integrand, args =  (deltas, ms), a= 0,b =np.inf)[0]\n",
    "error_ZN  = lambda deltas, ms: quad(integrand, args =  (deltas, ms), a= 0,b =np.inf)[1] \n",
    "\n",
    "# test function\n",
    "ZN_betalambda(np.exp(Dels0), Nelems0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Calculating $\\langle m \\rangle$ exactly\n",
    "\n",
    "From this partitition function, we can calculate the average number of elements in their correct placement through the equation\n",
    "\n",
    "\\begin{equation}\n",
    "\\langle m \\rangle  = \\sum_{j=1}^{R} n_{j} e^{\\beta \\Delta_j} \\frac{Z_{\\boldsymbol{n}_{j}}(\\{\\beta\\Delta_k\\}) }{Z_{\\boldsymbol{n}}(\\{\\beta\\Delta_k\\})},\n",
    "\\label{eq:ord_zn}\n",
    "\\end{equation}\n",
    "\n",
    "where $\\boldsymbol{n}_{j}$ is $\\boldsymbol{n}$ with 1 subtracted from the $j$th component: $\\boldsymbol{n}_{j}= (n_1, \\ldots, n_{j}-1, \\ldots, n_R)$. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "original n vector: [5 1 6 9 1 3 7 4 9 3]\n",
      "(n vector)_j value for j=3: [5 1 6 8 1 3 7 4 9 3]\n"
     ]
    }
   ],
   "source": [
    "# function to get n_j\n",
    "mj_vals = lambda j, ms: np.array([ms[k] if k != j else ms[k]-1 for k in range(len(ms))])\n",
    "\n",
    "num = 3\n",
    "# test function\n",
    "print('original n vector:', Nelems0)\n",
    "print(f'(n vector)_j value for j={num}:', mj_vals(num, Nelems0))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# average of ell\n",
    "avg_m_exact = lambda T, Dels, ms: np.sum([ms[j]*ZN_betalambda(delta_func(Dels, T), mj_vals(j, ms))/ZN_betalambda(delta_func(Dels, T), ms) for j in range(len(ms))])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "32.24972879913478"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# testing evaluation\n",
    "avg_m_exact(1.0, Dels0, Nelems0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Thermal limits of disorder\n",
    "\n",
    "At what temperature, will all the particles in the system settle into their correct location? In other words, with $\\langle m \\rangle$ representing the number of elements in their correct location, at what temperature does the system settle into the microstate  $\\langle m \\rangle = N$?\"\n",
    "\n",
    "In the accompanying paper, we found that this temperature $\\beta_c = 1/k_BT_c$ is defined implicitly by \n",
    "\n",
    "$$ 1 = \\sum_{j=1}^R n_j e^{-\\beta_c\\Delta_j}$$ \n",
    "\n",
    "\n",
    "We also found the necessary (but not sufficient) condition for the system to settle into the correct macrostate is \n",
    "\n",
    "$$ k_BT_c \\leq \\frac{1}{\\ln N} \\sum_{j=1}^r \\frac{n_j}{N} \\Delta_j$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.561392354979486"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# implicit thermal constraint\n",
    "\n",
    "def therm_constr(beta,Dels, Ns):\n",
    "    \n",
    "    F =  1-np.sum(Ns*np.exp(-Dels*beta))\n",
    "    \n",
    "    return F\n",
    "\n",
    "# critical temperature\n",
    "kBTcrit = lambda Dels, Ns: 1.0/fsolve(therm_constr, x0 = 0.5, args = (Dels, Ns))[0]\n",
    "\n",
    "\n",
    "# Test function \n",
    "kBTcrit(Dels0, Nelems0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.6809298387017404"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# necessary condition for temperature\n",
    "kBTlimit = lambda Dels, Ns : 1.0/(np.sum(Ns)*np.log(np.sum(Ns)))*np.sum(Dels*Ns)\n",
    "\n",
    "# Test function \n",
    "kBTlimit(Dels0, Nelems0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# checking inequality\n",
    "kBTcrit(Dels0, Nelems0) < kBTlimit(Dels0, Nelems0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Large Temperature Limit\n",
    "\n",
    "We found that in the limit of large temperature the total number of correctly ordered components has the behavior\n",
    "\n",
    "$$ \\lim_{T \\to \\infty} \\langle m\\rangle  = \\dfrac{ \\sum_{j=1}^R n_j^2}{ \\sum_{j=1}^R n_j}.$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "## ell limit\n",
    "m_limit = lambda Ns : np.sum(Ns*Ns)/np.sum(Ns)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Equations of Large $N$ approximation \n",
    "\n",
    "The order parameter for this system can be approximated as \n",
    "\\begin{align}\n",
    "\\langle m\\rangle & =   \\sum_{j=1}^R \\frac{n_j }{1- \\delta_j^{-1}} \\frac{\\displaystyle L_{n_j-1} \\left( \\bar{\\sigma}_{j}\\right)}{\\displaystyle  L_{n_j} \\left( \\bar{\\sigma}_{j} \\right)},\n",
    "\\label{eq:jtox}\n",
    "\\end{align}\n",
    "where $\\bar{x}$ is defined from \n",
    "\\begin{equation}\n",
    "N- \\bar{x} = \\sum_{j=1}^{R} n_j \\frac{\\displaystyle L_{n_j-1} \\left( \\bar{\\sigma}_{j}\\right)}{\\displaystyle  L_{n_j} \\left( \\bar{\\sigma}_{j} \\right)}\n",
    "\\label{eq:barx}\n",
    "\\end{equation}\n",
    "with $\\bar{\\sigma}_{j} \\equiv{\\bar{x}}/(1-\\delta_j); \\delta_j = e^{\\Delta_j/T}$ for notational simplicity  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "# function that determines xbar\n",
    "def xbar_func(x, T, Dels, Ns):\n",
    "    \n",
    "    R = len(Ns)\n",
    "    deltas_ = delta_func(Dels, T)\n",
    "    sigmas_ = x/(1-deltas_)\n",
    "    \n",
    "    RHS = np.sum([Ns[j]*eval_laguerre(Ns[j]-1, sigmas_[j])/eval_laguerre(Ns[j], sigmas_[j]) for j in range(R)])\n",
    "    \n",
    "    LHS = np.sum(Ns) -x\n",
    "    \n",
    "    return LHS-RHS\n",
    "\n",
    "# function that determs ellavge for large N\n",
    "def avg_m_approx(T, Dels, Ns):\n",
    "    \n",
    "    R = len(Ns)\n",
    "    xbar = fsolve(xbar_func, x0 = np.sum(Ns), args = (T, Dels, Ns))\n",
    "    deltas_ = delta_func(Dels, T)\n",
    "    sigma_bars = xbar/(1-deltas_)\n",
    "    \n",
    "    return np.sum([Ns[j]*deltas_[j]/(deltas_[j]-1)*eval_laguerre(Ns[j]-1, sigma_bars[j])/eval_laguerre(Ns[j], sigma_bars[j]) for j in range(R)] )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Example analytical plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "# temperature vals\n",
    "Tvals = np.linspace(0.1, 3.0, 50)\n",
    "\n",
    "# parameters for the integral\n",
    "np.random.seed(42)\n",
    "Dels0 = np.random.randn(10)+2\n",
    "n_vals = np.random.randint(1,10,10) # randomly generated integers \n",
    "\n",
    "# computed averages for each temperature\n",
    "avg_m_exact_vals = np.array([avg_m_exact(T, Dels0, Nelems0) for T in Tvals])\n",
    "avg_m_approx_vals = np.array([avg_m_approx(T, Dels0, Nelems0) for T in Tvals])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYsAAAEYCAYAAACtEtpmAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAzoUlEQVR4nO3dd3hUdfbH8feZSSMQWuhJIIB0pEYEXKSqgAKrIKI/C64KiqxiF9eyYlvXsqwKrFgRVwHRdQFRLIi6UgQ1ICjd0KQlEHpIO78/ZtQYEyZtcjM35/U8eTJz7zd3PidXc7hdVBVjjDHmVDxOBzDGGFPxWbMwxhgTkDULY4wxAVmzMMYYE5A1C2OMMQFZszDGGBOQK5vFwIEDFfjN144dO343LZS/CqunT58+2qdPH8fzlWVNofrltnrcWJPb6imDmgrlymaRmpr6u2k5OTkOJAket9UD7qvJbfWA+2pyWz0QvJrCgrJU45h7773X6QjGGBeyZuEyAwYMcDqCMcaFXLkbqjJLTk4mOTnZ6RjGGJexLQuXmTBhAgBLlixxNIcxxl1sy8IYY0xA1iyMMcYEZM3CGGNMQHbMIpTk5sLetZDyBfW+XwT7kyGqJjQ5C5r09H0ZY0wQONosRORl4AJgn6q2L2C+AP8EBgPHgdGq+k35pnTQySOwZy3sToaU/8G2L+HEQQDCYxpDm6GQkQ6bFsHqNwB4tH0NaHA6bPkUmvd1LrsxxlWc3rJ4FXgOeK2Q+YOAFv6vM4Fp/u/ukp0Jh3ZA2hbYswb2fOf7fmDrr2NqNoFW50PTXpDYi10Hs0hMTPTNU4XUjbDtS3puWwpbP4NZl8GNX0HNBEdKMsa4i6PNQlU/F5HEUwwZBrymvme/LheRmiLSUFV3l0/CUso+CccPwIkDv/1+dB8cTIH0bXBwGxzexW9uy1KzCTTsAB0v820lNOwA1Rv9dtkHU359LQJ1W0HdVizNbA3VBtNz5ThYdA9cMrMcCjXGuJ3TWxaBxAE78rzf6Z8WlGaxevEcvCumFjrfQy6CIuQiqnj8r8M0i3DNJEIzCdcswvWk7z1ZhS7rYFgdDoY3Ij2yHYcbnMvRKo04XjWB3PrtqFmrLvWqR1K3WiR1YyKJCvcWuYZ77rkHgCWTboPFD8PmT+C0/kX/JRhjTAEqerMoMhEZA4wBiIuLIyUl5Tfz09LSAi7jYFoq9bJPFDr/11bhQfGSK7732cRwUiLIlHCyJIKT+F6foAqHJYbDxHBIYjhMNdKJ4SDVOZ4bRmaWkpmhZB3IJTNHycxRYLv/61fVI720qVeFLvFV6RJXleaxURw8cKDAjBkZGQCkxA0jLuY1mHcLu4bNBW9EwPqdVpR1FErcVg+4rya31QOlq+mXXdsFqOjNYheQd6d7vH/a76jqdGA6QFJSkhZU9Kl+Eb7544BxJQpaFrJzcjlwLJN9R06y/8hJ9h3JYP+Rk+xKP8FXPx5g2rK9ANSKDqdDgyqc06EOg09vSO2qvzaCqKgoABKbt4Sh/4B/DyfxpwXQ61ZHaiquQOso1LitHnBfTW6rB4JTU0VvFvOA8SIyC9+B7UMhc7yiBMK8HupVj6Je9agC5+85lMHSLal8uTmNzzfs4bN31/L3D9Zzx8DWXNatMV6P/PYHWgyA1hfA509Ah5FQI74cqjDGuJHTp86+CfQB6ojITuABIBxAVf8FLMR32uxmfKfOXu1M0oqhQY0oLuoSz0Vd4vnxx+qciIzl4fe+57531/LWqh08NOx3Zx/DeY/ClDN9B7tHFnbSmTHGnJrTZ0NdGmC+AjeWU5yQIiK0bVSdf197JvPX7ObhBd/zx6lf0m/ojYzukfjrwFpNoNdt8OnDdu2FMabE7HYfIU5EGNqxEZ/c1ps/ndWUJalVufnjdN5bk2dvXc8/Q+1msPAO3zUdxhhTTNYsXCImKpz7LmjLPR2ziNq3jj+/+Q1LNuzzzQyPgkF/h7RNsHyKs0GNMSHJmoXLzJz2D7yr/0OrBtUZ/8a3bNhzxDejxTm+K8A/e8J3YaAxxhSDNQsX8nqEl65KIjrCyzUzVpJ69KRvRr97IesYrHrJ2YDGmJBjzcKlGtWswotXJZF69CRjXltFRlYO1G8LzfvDVy/4bkVijDFFZM3CxTrE1+QfIzvxzfZ07py7BlWFnuPh6F74bq7T8YwxIcSahcsNOr0hd5zXinmrf+Kfn2yCZn2hXjtYNsV3t1pjjCmCin4Ftymm559//nfTxvVpztb9x5j88Saa1qnKsB43wn/HwZbFdpNBY0yR2JaFy7Rq1YpWrVr9ZpqI8OhF7enWtDZ3zl3DroTzoVp9WPacQymNMaHGmoXLzJ8/n/nz5/9uemSYl6dHdkSBpz5JgW5jfFsWe9eVe0ZjTOixZuEyTz31FE899VSB8+JrRXP1WYn8J3kXP8SNgLAqvmMXxhgTgDWLSmZcn9OoWSWchz/dg3b6P1gzB47scTqWMaaCs2ZRydSoEs6f+7Xgy81prGgwCnKzfdddGGPMKVizqIQu796EJrHRPPDFCbTVYN8V3ZnHnI5ljKnArFlUQhFhHu48rzUb9h5hce2RcOIgJL/hdCxjTAVm11m4zMyZM4s0bvDpDejcuCYTVx5nWcPOeJdPhaQ/gccb5ITGmFBkWxYuk5CQQEJCQsBxIsJfBrdh39FMPqxxMRzYChsXlUNCY0wosmbhMrNnz2b27NlFGpuUWJvz2tXn7u+bkFOtAay0A93GmII53ixEZKCIbBCRzSJydwHzm4jIJyKyRkSWiEi8EzlDxbRp05g2bVqRx981sDXHsj18Wu0C30V6qZuCmM4YE6ocbRYi4gWmAIOAtsClItI237AngddUtQMwCXisfFO6W7O61bjszMbcs70r6gmHlS86HckYUwE5vWXRDdisqltVNROYBQzLN6YtsNj/+tMC5ptS+nO/FqR7apFcvY/vrKiTR5yOZIypYJxuFnHAjjzvd/qn5bUauMj/+kIgRkRiyyFbpVE3JpIRXeN5NPVsOHkYVs9yOpIxpoIJhVNnbweeE5HRwOfALiAn/yARGQOMAYiLiyMlJeU389PS0oKds1wVVk9GRgbA7+oPZHCzCK78qhk7q7Wg3pdT+anuABApbcxiqSzrKJS5rSa31QOlqykxMbHQeU43i11A3vM84/3TfqGqP+HfshCRasBwVU3PvyBVnQ5MB0hKStKCij7VLyIUFVTPggULAKhTp07xlgUMWnuMaZvO4ZGMqSTqdmjau/Qhi6kyrKNQ57aa3FYPBKcmp3dDrQRaiEhTEYkARgHz8g4QkToi8nPOicDL5ZwxpNSpU6fYjeJnY3s3Y25GN06E14SvppdtMGNMSHO0WahqNjAeWAT8AMxR1XUiMklEhvqH9QE2iMhGoD7wiCNhQ8Srr77Kq6++WqKf7RBfk67NGzI7py+6YSGk7wj8Q8aYSsHpLQtUdaGqtlTV5qr6iH/a/ao6z/96rqq28I+5VlVPOpu4YitNswC4vndzXjjex/d47lW2EWeM8XG8WZiKpVeLOtRs1JylYWeg38yArAynIxljKgBrFuY3RISxvZsz5fgA5HgarHvH6UjGmArAmoX5ncHtG7CzRld2eBPQFc/j2ydljKnMrFmY3wnzerju7OY8nzEA2Z0MO1c5HckY4zBrFi6zcOFCFi5cWOrlXNw1gSWR/Tgu0bB8ahkkM8aEMmsWLhMdHU10dHSpl1MlwsvFPdswM6sv+v1/IX17GaQzxoQqaxYuM3XqVKZOLZstgSt7NGGWDCZXgeX/KpNlGmNCkzULl5kzZw5z5swpk2XVqhpB325dWJDTndxvZkDGoTJZrjEm9FizMKd0ba+mvJx7Pp7Mo/D1DKfjGGMcYs3CnFKjmlVo3vEslms7cpdPg5wspyMZYxxgzcIEdH3v5jyfNRjPkZ9g3X+cjmOMcYA1CxNQy/oxhLU8h63Ekfvls3aRnjGVkDULl1myZAlLliwp8+WO7dOC6VmD8OxdAylflPnyjTEVmzULUyRJibXZHj+Eg1Qnd+mzTscxxpQzaxYu8+STT/Lkk08GZdnX9GnDq1nn4Nn0IezfEJTPMMZUTNYsXGbBggW/PFq1rPVtVY9ltS/kJBHosilB+QxjTMVkzcIUmccjjOrbmbnZvchNfhOO7nc6kjGmnDjeLERkoIhsEJHNInJ3AfMbi8inIvKtiKwRkcFO5DQ+Qzo2YkH0hXhzM2HlC07HMcaUE0ebhYh4gSnAIKAtcKmItM037F58z+buDIwC7BaoDgr3eji39x9YlJNE9rJpdgsQYyoJp7csugGbVXWrqmYCs4Bh+cYoUN3/ugbwUznmCzlVqlShSpUqQf2MS85I4JWwkYRlHoYVzwf1s4wxFYPTzSIO2JHn/U7/tLz+ClwuIjuBhcCfyydaaHr//fd5//33g/oZ0RFh9Dq7Px/ldCX7y+ds68KYSiDM6QBFcCnwqqo+JSI9gJki0l5Vc/MOEpExwBiAuLg4UlJSfrOQtLS0copbPpyup2+c8KB3BOdkTuTgosc51HFMqZfpdE1lzW31gPtqcls9ULqaEhMTC53ndLPYBSTkeR/vn5bXNcBAAFVdJiJRQB1gX95BqjodmA6QlJSkBRV9ql9EKCqonoceegiA++67L+if37fveXz08Vv0Xfc6tc67G6KqB/6hACrDOgp1bqvJbfVAcGpyejfUSqCFiDQVkQh8B7Dn5RuzHegPICJtgCjAztksxCeffMInn3xSLp91ZY8mvBp+CWGZh+ArO3ZhjJs52ixUNRsYDywCfsB31tM6EZkkIkP9w24DrhOR1cCbwGhVu5NdRRAdEUbvPufwcU5n37GLk0ecjmSMCRKntyxQ1YWq2lJVm6vqI/5p96vqPP/r71X1LFXtqKqdVPVDZxObvK7onsiMiFGEnUyHr6Y7HccYEySONwsT2qpEeOnd51wW53Qi64tnbOvCGJeyZuEysbGxxMbGlutnXt69CTMiRhGemQ5f2VXdxriRNQuXefvtt3n77bfL9TOjwr306TeQT3M6kvW/Z+Dk0XL9fGNM8FmzMGXi0m6NmRk5ivCTB1HbujDGdaxZuMzEiROZOHFiuX9uVLiXPv0HsySnI9lfTIYT6eWewRgTPNYsXGbZsmUsW7bMkc++5IwEXo66Em/mIfSLpx3JYIwJDmsWpsxEhnk5/9xzeSenF7nLp0H6dqcjGWPKiDULU6ZGdE3g3VpXk50LOR9PcjqOMaaMWLMwZcrrEcYO6cUL2YPwrn0Ldn3jdCRjTBmwZuEy8fHxxMfHO5qhV4u6rG16NQeoTtYH94LdncWYkGfNwmVef/11Xn/9dadjcOsFSUzOGk74ji9h4wdOxzHGlJI1CxMULevHkNv1KrZqQzI/uBdysp2OZIwpBWsWLjNhwgQmTJjgdAwAbj6nLU/r/xFxcDN8M8PpOMaYUrBm4TLJyckkJyc7HQOAujGRtOkzihW5rcn85BG7yaAxIcyahQmqa3o1Y3rk1URkpKH/m+x0HGNMCVmzMEEVFe5lyOAhzMvpQc7S5+xCPWNClDULE3RDOzbi3TpjyMpRshfe6XQcY0wJWLNwmZYtW9KyZUunY/yGxyOMG9aHyVkXEbbxfdjwvtORjDHF5HizEJGBIrJBRDaLyN0FzP+HiCT7vzaKSLoDMUPG9OnTmT694j3eNCmxNkc6Xsem3Diy5t8OmcedjmSMKQZHm4WIeIEpwCCgLXCpiLTNO0ZVb/E/e7sT8CzwTrkHNWXijvNP52/eMYQf3Yl+9oTTcYwxxeD0lkU3YLOqblXVTGAWMOwU4y8F3iyXZCFqzJgxjBkzxukYBapVNYKB5w9nbs7Z5C59BvZvcDqSMaaInG4WccCOPO93+qf9jog0AZoCi8shV8jauHEjGzdudDpGoUZ0jeeDhuM4qpFkzbvF7htlTIgIczpAMYwC5qpqTkEzRWQMMAYgLi6OlJSU38xPS0sLdr5yVVg9GRkZAL+rvyK5okcznnh7FA/veIn9i6dwrPkFQOVZR6HMbTW5rR4oXU2JiYmFzitRsxARD3A24AWWFPYHvAh2AQl53sf7pxVkFHBjYQtS1enAdICkpCQtqOhT/SJCUUH1REVFFTqvokhMhK/2Xsu3Sz+j3cqnqdvz/6BKLf+8REezlTW31QPuq8lt9UBwairybigRiRaRi0TkNWAfvt1BHwH7ReR1ERkhItWK+fkrgRYi0lREIvA1hHkFfHZroBbgzPNCTZkb368Vz0bfgDfjoD0kyZgQcMpmISL1ReRaEVkApAJzge7AK/i2LM4CXgSSgDn4GsdCERkjIg0CfbiqZgPjgUXAD8AcVV0nIpNEZGieoaOAWaq2gzuQTp060alTJ6djBFQlwssVfxzKjOxz8Xz9CuxY6XQkY8wpBNoN9ZP/+yrgIeC/qvp9vjHLgDv9//r/IzAUmAZMFZEqqpp1qg9Q1YXAwnzT7s/3/q8Bchq/yZMnOx2hyPq2rsf8ljeye8tK6r49FrngDacjGWMKEWg31DggXlXPVNXHCmgUv1DV9ar6N1Xtie+MphuAkh7LMJXEXX/sxgNyA+HpW6jx9bNOxzHGFOKUzUJVn1fV3cVdqKruUdUXVDW35NFMSVx++eVcfvnlTscosvrVoxg09FJmZg+gxg//hm1LnY5kjClAsa+zEJFYEekrIr1EpF4wQpmS27lzJzt37nQ6RrFc2DmOlS0msEPrkvn29ZB5zOlIxph8itUsRGQYkAJ8DCwBdovIDhH5r4g8ICJDRKTAi+qMKYyIcP/wbvxVbiDs8HZyPrw/8A8ZY8pVca+zeBzIAP4CHAbaAF3wnRU1xD9G8V1/YUyR1akWSZ+zB/DK4q+4ZtWL0HYINOvjdCxjjF9xm0UCcJ+qPpN/hv92HF2BTmWQy1RCZzerztS9t7B1fTLx74wjYvxyiKrudCxjDMU/ZvEDIAXNUNVtqvpO/tNeTfnq0aMHPXr0cDpGif1lWBceCb8J79HdZH9wj9NxjDF+xW0W/wQu999a3FRAjz32GI899pjTMUqsRnQ4V468mOnZ5xOWPBM2fuh0JGMMxWwWqjoT+B/wjojEBieSqex6t6zL7s63sD43gay3x8LhYp+9bYwpY8U9GyoOaIzvYPYuEflARP4qIkPtLKiKYfjw4QwfPtzpGKV21wUdebTqXWSfPEbWW3+CXLu+0xgnFXc31Cv4GkUy8CXQGbgf+A+wXUT2iMjCwn/cBFtaWporbrtcNTKMOy8fygPZ1xC+Yym65HGnIxlTqRW3WZwFvKCqXVS1v6rWx7elcSHwML57SHUq24imsmofV4PTL7ieuTlnw+d/h62fOR3JmEqruM3iIPB13gmqulNV56nqA6p6gao2Krt4prK7/MzGLG11N1u0IZlvXQNH9zkdyZhKqbjN4h2gVzCCGFMQEeHBEd14JPpOck8cIuut6yDXbjlmTHkrbrN4FugoIpcEI4wpvf79+9O/f3+nY5SpmKhwbrv8Ih7OvorwbUvI/d8/nI5kTKVT3Cu4NwDHgTf8DWMOsFxVU8o6mCmZ++67z+kIQdE+rgatzh/PvIVruWDxI9CkJzQJ3YsPjQk1xd2yeALfWVBp+B509AawRUTSROQjEXlcREaWcUZjALi8exM+a/kXtmtdMmddaddfGFOOintR3l2qep6q1sN3n6ihwIPA50AL4A7gzeIsU0QGisgGEdksIncXMmakiHwvIutExB6ndgqDBg1i0KBBTscIChHhrxd358HoiWSdOEzmG5dCVobTsYypFIr9PIufqeouVV2gqpNU9UJVTQRigXOKugz/bUOmAIOAtsClItI235gWwETgLFVtB0woaebK4MSJE5w4ccLpGEETExXOxNHDuTt3PBF7viX7vzeBPZrdmKArcbMoiKoeVNXFxfiRbsBmVd2qqpnALGBYvjHXAVNU9aD/M+zcyUquZf0Yho0aw9PZIwhbOxtdao9jNSbYTtksRGRcSW7jISINRWSsiARqRnHAjjzvd/qn5dUSaCkiX4rIchEZWNw8xn0GtK1PVP+7WJjTDf3oAdj8sdORjHG1QGdD/RN4VkS+Bd4F5qnqmoIGikh7fFsFQ4EkIBd42f+9tBlbAH2AeOBzETldVdPzff4YYAxAXFwcKSkpv1mIG26BkVdh9WRk+Pbh568/FBR3HQ1KDOfJhDtpuvNWmr95FfuG/pvsGonBCVcCbvtvDtxXk9vqgdLVlJiYWOi8QM2iPr57QQ0D7gYeFJHtwH+BeUCOf94wIBHfU/Q+Aq4F5qtqVoDl78J3oPxn8f5pee0EVviX9aOIbMTXPFbmHaSq04HpAElJSVpQ0af6RYSiguoZMWJEofNCQXFzP3V1Y/489a88fuBm6n56G5HXL4aoGsEJVwKhuh5OxW01ua0eCE5Np9xNpKoHVHWGql4E1MF3uuxiYBS+53AvBi4HPgMuAuqo6h9V9RVVTS3C568EWohIUxGJ8C93Xr4x7+LbqkBE6uDbLbW1SNVVQrfffju333670zHKTVS4l4dGn8/EsDvwpv9I5pw/QU6207GMcZ0iH+BW1QxVna+q1wANgT8AvYEGqvonVf2vqhbrNBxVzQbGA4vwPYVvjqquE5FJIjLUP2wRkCYi3wOfAneoqvu2HU2JNagRxQ2jr+LBnKuJ2Pox2fMn2BlSxpSx4l7BDYCqKrC0LAKo6kJgYb5p9+d5rcCt/i8TQJ8+fQBYsmSJoznKW6eEmmwbfivPzd3P+OSZ5NaIw9N3otOxjHGNMj111hgnDesUR+Q5D/BW9tl4PvsbuupVpyMZ4xrWLIyrXNe7OVu6P8qnOR3RBbfAhvedjmSMK1izMK5z5+D2fND2b3yXm0j27NGwY2XAnzHGnJo1C+M6Ho/w8MjuTE/4Gzuza5A5cwSkbnI6ljEhzZqFy4wcOZKRI+3Gv+FeD3+/sj+PxT7C4ZPKyVeHwaH8l/AYY4rKmoXLjBs3jnHjxjkdo0KoGhnGo9cM5S/R95N59AAnXz7fbmtuTAlZs3CZ48ePc/z4cadjVBix1SK597rLuC38frIP7ebkyxfYc7yNKQFrFi4zePBgBg8e7HSMCiWhdjR/uf4qbg27l9z0HZx86QI4VpQbDBhjfmbNwlQKTWKrcvf1f+IW7z3owR99WxjHDzgdy5iQYc3CVBpN61Tljuuv5Tbv3ZC2mYyXh8CJg07HMiYkWLMwlUrzutW4ZewYbvfchSd1PRmv/BFOpDsdy5gKz5qFqXROqxfDTWOv507P7Xj3reXkS4Ph6H6nYxlToVmzcJnRo0czevRop2NUeC3qx3DDmPHc6rmL3NRNnHzhXEjfEfgHjamkrFm4jDWLomvVIIZbx93IrREPcDJ9Dyenn2NXehtTCGsWLpOamkpqqp0WWlRN61Tl/huv4c5qj3L02DEyXzgXdq92OpYxFY41C5cZMWLEL49WNUXTsEYVHhv3f9xb+0lSM4SslwbDtjJ5XIsxrmHNwhigVtUInrx+OI81nMz2zBiyZ1wIGxc5HcuYCsOahTF+VSPDePLa83m+2XN8n92Q3DdGkbv8eadjGVMhON4sRGSgiGwQkc0icncB80eLyH4RSfZ/XetETlM5RIZ5eezK/rzTYTqf5HTG88GdZC+4A3JznI5mjKMcbRYi4gWmAIOAtsClItK2gKGzVbWT/+vFcg1pKh2vR3hg+BlsP+d5XsweRNiq6Zx8/RI4edTpaMY4xukti27AZlXdqqqZwCxgmMOZQtoNN9zADTfc4HSMkCciXHN2CxIv+ycP5l5D2NZPODH9XHsmhqm0whz+/Dgg75VQO4EzCxg3XETOBjYCt6jq766eEpExwBiAuLg4UlJSfjM/LS2tjCJXDIXVc+aZvl9f/vpDQUVcR6dFgw69jlsX1uPR1H9wfMrZpJ/3HJmxbQL+bEWsp7TcVpPb6oHS1ZSYmFjoPKebRVHMB95U1ZMiMhaYAfTLP0hVpwPTAZKSkrSgok/1iwhFBdWzY4evjyYkJJRzmrJREddRYiJ0bNmUiS/HcdfBB6j/3mi8w/6JdBxVhJ9NDHq+8ua2mtxWDwSnJqd3Q+0C8v5Vi/dP+4WqpqnqSf/bF4Gu5ZQtJF1xxRVcccUVTsdwnXrVo3h83GU80+x5VmU3Rf4zlsx5t0J2ptPRjCkXTjeLlUALEWkqIhHAKGBe3gEi0jDP26HAD+WYz5hfVInw8rcrB7Cm7wxezB5MxDcvceLFQfaoVlMpONosVDUbGA8swtcE5qjqOhGZJCJD/cNuEpF1IrIauAkY7UxaY3wHvsf0bUXbq59loucWdPd3ZEz5A6R86XQ0Y4LK6S0LVHWhqrZU1eaq+oh/2v2qOs//eqKqtlPVjqraV1XXO5vYGOjZvA4333w3E2Mn89OJMHJfHULOl8+BqtPRjAkKx5uFMaGqQY0onhg3itmdX+OjnM54P/oLGa+NsOd7G1eyZuEyt912G7fddpvTMSqNiDAPEy88k4yLZvCoXo38+BkZz5wJWxY7Hc2YMmXNwmWGDBnCkCFDnI5R6QzrHM/lNz3KXbUns+NEJMy8kGornrSzpYxrWLNwmQ0bNrBhwwanY1RKjWOjeeLGy3i/x5v8O6c/dX6YyfF/9YO0LU5HM6bUrFm4zNixYxk7dqzTMSqtcK+HmwZ1pMWfXuROz+1k7t9K1pSzyF3xAuTmOh3PmBKzZmFMEHRrWpsrLrmMp5q/wrKs0/C8fzvHXzwfDqY4Hc2YErFmYUyQxER6mXTFuRwaPptJcj25u74l89nu5Cx/3rYyTMixZmFMEIkIQzrFceNtk3i82Sssy2qB94M7OfbCIDiw1el4xhSZNQtjykFstUgeumoQJ0bOYZJnHLk/rSbruR5kfzEZcrKcjmdMQKFw11lTDPfee6/TEcwpDDy9Id2bP8jT7/Sl54bHOOeTBzj29RtUvegZaNzd6XjGFMq2LFxmwIABDBgwwOkY5hRqRkfwwOXnEnHFHP4SOZH0g6nw8nkcf+t6OOa+5ysYd7Bm4TLJyckkJyc7HcMUQe+Wdbnv9jt4t+e7vJA7hPC1c8j4R2eyV82wA+CmwrFm4TITJkxgwoQJTscwRRQV7uXG8zpw3s3TmRT3L1ZnNiRswU0cndoHtq9wOp4xv7BmYUwF0Dg2mknXXcyhke/yUPhNHN2/HV4+l6OvXw4HtzkdzxhrFsZUFCLCue0bcsedf+W/veYzVYfj3bSI7GeSyHj/Psg47HREU4lZszCmgokK9zJ2wOlcfPs0JreZxbzsbkSteIYTT3Uke8WLdqqtcYQ1C2MqqLoxkUwc1Z92N87igfrP8d3JuoS9fxtHn+pM7uo5dhDclCvHm4WIDBSRDSKyWUTuPsW44SKiIpJUnvlCzaOPPsqjjz7qdAxThlo1iOHBG67g+GXz+Wu1B9h+1IPnP9dxePKZ5P7wnj2dz5QLR5uFiHiBKcAgoC1wqYi0LWBcDHAzYKeHBNCzZ0969uzpdAwTBH1a1+f+W28hZcT7PBR1O6nph/HMvoxDz/VGtyyxpmGCyukti27AZlXdqqqZwCxgWAHjHgIeBzLKM1woWrp0KUuXLnU6hgkSj0cY3CGOe+68l+Qhi3g8YhzHUncgM4dxaEo/dOOH1jRMUDjdLOKAHXne7/RP+4WIdAESVPW98gwWqu655x7uuecep2OYIPN6hIvOSOSWOx9hybmLeCLsOo7u34a8cTGHnjmL3O/n2zENU6Yq9L2hRMQDPA2MLsLYMcAYgLi4OFJSUn4zPy3NXbdRKKyejAzfxlf++kNBZVlHZa1nXARZl45j3vqRHPzmHS5Le4cacy7nQHRTMruOJaPpOeApm//VbR1VfKWpKTExsdB5TjeLXUBCnvfx/mk/iwHaA0tEBKABME9EhqrqqrwLUtXpwHSApKQkLajoU/0iQlFB9URFRRU6LxSEau7ClGc9LZo3I3tgEgtX38Csj15l+NHZtPjibo6seobIP4wnIukqiKxW6s+xdVTxBaMmp3dDrQRaiEhTEYkARgHzfp6pqodUtY6qJqpqIrAc+F2jMMb4hHk9DO3ShDvvuJ+tIz/mser3sv5YNSI+uoeMv7fm2Hv3wuHdTsc0IcjRZqGq2cB4YBHwAzBHVdeJyCQRGepkNmNCmccjnNe+EXffcjvyp0U8Hvcsn2a2Ieqr58h+uj3pb1wDu1c7HdOEEKd3Q6GqC4GF+abdX8jYPuWRKZRNnjzZ6QimAhERkhJrk3TdlWxLG85zi/9HnbUvceGG+bBxLgdqd6b62TcQ1v5CCItwOq6pwBxvFqZsderUyekIpoJqEluVmy8+j0Pn92PW0rUcXTGDoakLqf3uGI69dzfa5Sqq9bwOasQFXpipdKxZuMzHH38MYA9AMoWqER3OnwZ0JqdfJ5asv4u3l7xDp91v0Xf5ZHJW/JP0hAHU+sO1eFoMAI/X6bimgrBm4TIPP/wwYM3CBOb1CP3bNqR/2xvZljaaf322jOg1M7hg+xI8b37IkYj60PkyYrpfDbWaOB3XOMyahTGGJrFVGXfRADKG9OXD77az5Yu5dE6dz9nLJ5O7YjIH6vWg5h+uQaLbOx3VOMSahTHmF1HhXoZ2aQpd7mBb2jimf7kSb/K/Gbz3E8LeuYZ6nmjSThtK7Z5XIo17gMfps+9NebFmYYwpUJPYqlw/tA/Z55/N5xv28PaX7xG/47+ct+E/yMZZHIpsBB0vocaZV0Bsc6fjmiCzZmGMOaUwr4d+bRvRr+11fLehHwvTstizYi4dD3zAWSsmw1f/IK16GyI7XUy1LhdDzcZORzZBYM3CZZ5//nmnIxgXi4n0cnHP5tDzXnYevJXXVqwmM3k2Z6Z/RqfPJ8Hnk9hfsyPRXS6maqcRUL2h05FNGbFm4TKtWrVyOoKpJOJrRXP1wB4wsAcb9x7hpRUryV37Dmcd+Jy2i+8ld/F9pNbqRHSHP1Kt0x+hVqLTkU0pWLNwmfnz5wMwZMgQh5OYyqRl/RhaDu2HDunL97sP8+KKZXi+f5fuaUtp+9kD8NkDpMa0JrzdUGp0GQ51W4Hv5qAmRFizcJmnnnoKsGZhnCEitGtUg3YXDkT/eB7r9xzhlVWryPl+Pp0PfUHX5X+H5X/nYFQC2c3PJbbLMDyJPcEb7nR0E4A1C2NMUIgIbRpWp82QfjCkH9vTjvPGt2s4tnoeLQ79jx5rZ+BZ9xLHPdU4FNeb2l2GEtnqXIiu7XR0UwBrFsaYctE4NprLBnSHAd05cCyTRd//yJ5vPqDuT4v5w/YviNzxHrl42F+jPWEtz6F2x8FIoy52LUcFYc3CGFPualeNYOgZreCMVmTl/JmVP6Yyf9USwrZ+TIeDq+jw1dPIyqc46q3B4Ua9qNlhENGt+tvZVQ6yZmGMcVS410PP0+rR87SRwEh2HjzOf9Zu4uB3H1Bv7//osf1zoncsgPcgtUpTMpucTd0O5xHerBdEVXc6fqVhzcJlZs6c6XQEY0olvlY08b06Qq+OZOXcwTcpaby/Zjm65VOaHl7JGT+8Sfj6GeTgYX/19tC0F3Xb98fbpDtEVHU6vmtZs3CZhISEwIOMCRHhXg9nNq/Lmc2HAEM4kpHF0k172LFmCRHbP6dN+jecnjwN7+opZBPG/urt/M2jH2GNzyyTZ44bH2sWLjN79mwALrnkEoeTGFP2YqLC6X96Apx+BXAF+45k8NHGHexd9xmRO5fSKn01HZKnErb6Od+WR7VW5MT3ILZtb6KanQXV6jpdQshyvFmIyEDgn4AXeFFV/5Zv/vXAjUAOcBQYo6rfl3vQEDFt2jTAmoWpHOrFRDGoawvo2gK4lv1HTvLxxu3s+/5zwnetoNnh1XT6YSaR618GIDWyMScaJBHT4ixqtjwLNNLZAkKIo81CRLzAFOAcYCewUkTm5WsGb6jqv/zjhwJPAwPLPawxpsKrGxPJwF+axzUczshi+dY97Pp+ObptGQ0PfUunlA+pue0d+BjqSVV21uqAp3E36rQ5i4iEM+w6j0I4vWXRDdisqlsBRGQWMAz4pVmo6uE846sCWq4JjTEhq3pUOL3bJkDbBOBiMrNz+eGnQ3y6fjXHty6j6t6vaZO6gZZpz+BN/icAaRHxHK3TkeimZ1C7ZU+8jTpAeBVnC6kAnG4WccCOPO93AmfmHyQiNwK3AhFAv4IWJCJjgDEAcXFxpKSk/GZ+WlpamQSuKAqrJyMjA+B39YeCyrKOQpkbaqoBdG2ZCC0TSUs7l+NR1Xnrp1SO7lhDZOpa6h9bz+m7llH3p/fgS8jGy77IJhyp0RpPw/ZEx51Odu1WaFjF3IVVmnWUmJhY6Dynm0WRqOoUYIqIXAbcC1xVwJjpwHSApKQkLajoU/0iQlFB9URFRRU6LxSEau7CuK0ecF9NiYmJdG3XAugBQG6usjX1KO9t3MihLSsI35tM/aPraZfxBbH7FsJqXwNJrdKMjDrtqNK4M7HNuxLWqGOFue4jGOvI6WaxC8h7rme8f1phZgHTgpooxM2dO9fpCMaENI9HOK1eDKfV6wp/6ApAVk4um/YcYdnWDRzZ+hXhe9dQ79h62mz/lLo73oUvfT+bFt6II7XaEtbodGo160zVhI5Qo7ErblnidLNYCbQQkab4msQo4LK8A0Skhapu8r89H9iEKVSdOnWcjmCM64R7PbSNq0HbuG7QqxsA2Tm5pKQd4+utWzi89WvY+x21Dv1Aiz3fEb/vY0j2/ewJiSat6mlk1WlDdHwHajftSHjDdiF3IN3RZqGq2SIyHliE79TZl1V1nYhMAlap6jxgvIgMALKAgxSwC8r86tVXXwVg9OjRjuYwxu3CvB7/Fkgn6N4JAFVl7+GTfLFjN/u3JpP103dUObCeBoe30PrIfKqnzIb/+X4+3RvLoZjTyK3Tmuj404lNbE9Y/dZQpZZjNZ2K01sWqOpCYGG+affneX1zuYcKYdYsjHGOiNCgRhQNajSF9k2BCwHIzM5ly74jLEvZxJHta8jdt56YwxuJO5BCi4PfUmXzjF+Wke6N5XC1pmTXbkmVuHbUbtyeyIatoVp9Rx8Y5XizMMYYt4sI89CmUQ3aNEqCnkm/TD+emc2mPYf4KWUjR3etg33rqXZkCw0ObuO09Heo9uMbv4w9IdEcqNKEjJrNCavbkurxbagR3wZPbHOIiA56DdYsjDHGIdERYXRoHEuHxj34+WwsgIysHFJSj7Jz22aO7PyBnP0biErfSu1j20g8tpRGPy2A1b8u50BYPQ5HNyG7VjNOxPcGF54NZYwxJp+ocC+tG9agdcOuQNdfpv98TGTZT3s5sGM9GXs24DmwhapHU6iXvpOmh+bzVVYspwchkzULY4wJEb8eE2kCbZoA5/0yLzsnl50HjhO791RXH5ScNQuXWbhwYeBBxhjXCfN6SKxbDY4F58pyaxYuEx0d/ANdxpjKJ/QvKzS/MXXqVKZOnep0DGOMy1izcJk5c+YwZ84cp2MYY1zGmoUxxpiArFkYY4wJyJqFMcaYgKxZGGOMCUhU3feUUhHZD2zLN7kOkOpAnGBxWz3gvprcVg+4rya31QOlqylVVQcWNMOVzaIgIrJKVZMCjwwNbqsH3FeT2+oB99XktnogeDXZbihjjDEBWbMwxhgTUGVqFtOdDlDG3FYPuK8mt9UD7qvJbfVAkGqqNMcsjDHGlFxl2rIwxhhTQtYsjDHGBOS6ZiEiA0Vkg4hsFpG7C5gfKSKz/fNXiEiiAzGLrAj1jBaR/SKS7P+61omcRSUiL4vIPhFZW8h8EZFn/PWuEZEu5Z2xOIpQTx8ROZRn/dxf3hmLS0QSRORTEfleRNaJyM0FjAmZ9VTEekJqPYlIlIh8JSKr/TU9WMCYsv1bp6qu+QK8wBagGRCB7ym1bfONGQf8y/96FDDb6dylrGc08JzTWYtR09lAF2BtIfMHA+8DAnQHVjiduZT19AEWOJ2zmDU1BLr4X8cAGwv47y5k1lMR6wmp9eT/vVfzvw4HVgDd840p0791btuy6AZsVtWtqpoJzAKG5RszDJjhfz0X6C8iUo4Zi6Mo9YQUVf0cOHCKIcOA19RnOVBTRBqWT7riK0I9IUdVd6vqN/7XR4AfgLh8w0JmPRWxnpDi/70f9b8N93/lP1upTP/Wua1ZxAE78rzfye//o/hljKpmA4eA2HJJV3xFqQdguH9XwFwRSSifaEFT1JpDSQ//7oL3RaSd02GKw7/rojO+f7nmFZLr6RT1QIitJxHxikgysA/4SFULXUdl8bfObc2iMpoPJKpqB+Ajfv2XhKkYvgGaqGpH4FngXWfjFJ2IVAPeBiao6mGn85RWgHpCbj2pao6qdgLigW4i0j6Yn+e2ZrELyPsv63j/tALHiEgYUANIK5d0xRewHlVNU9WT/rcvAl3LKVuwFGUdhgxVPfzz7gJVXQiEi0gdh2MFJCLh+P6w/ltV3ylgSEitp0D1hOp6AlDVdOBTIP8NAMv0b53bmsVKoIWINBWRCHwHdeblGzMPuMr/egSwWP1HgCqggPXk2088FN/+2FA2D7jSf7ZNd+CQqu52OlRJiUiDn/cTi0g3fP/PVdR/nAC+M52Al4AfVPXpQoaFzHoqSj2htp5EpK6I1PS/rgKcA6zPN6xM/9aFlfQHKyJVzRaR8cAifGcSvayq60RkErBKVefh+49mpohsxndgcpRziU+tiPXcJCJDgWx89Yx2LHARiMib+M48qSMiO4EH8B2cQ1X/BSzEd6bNZuA4cLUzSYumCPWMAG4QkWzgBDCqAv/j5GdnAVcA3/n3iQPcAzSGkFxPRakn1NZTQ2CGiHjxNbY5qrogmH/r7HYfxhhjAnLbbihjjDFBYM3CGGNMQNYsjDHGBGTNwhhjTEDWLIwxxgRkzcIYY0xA1iyMMcYEZM3CGGNMQNYsjDHGBGTNwhhjTEDWLIwxxgRkzcIYY0xA1iyMMcYEZM3CmFIQkVUi8p3TOYwJNmsWxpSQ/+lj7YHkUiyjmojkiIgW8at2mRVgTDG46uFHxpSztkAkpWgW+P4fvCrftBuAnsDtwN4800+q6oFSfJYxJWYPPzKmhETkSmAG0F9VF5fhcr/G14hiVDW7rJZrTGnYbihjSq6z/3vyzxNEpKaI/EdEMkRkTHEXKCLh+HZtrbFGYSoS2w1lTMl1Anb8vGtIRLoCbwECnKWqX5dgme2ACODbsgppTFmwLQtjSq4T/q0KEbkB+BL4AehawkYBv26tfFPacMaUJWsWxpSAiCQCNYFNIvIG8BzwCHBBKQ9Cd/F/ty0LU6HYbihjSubnLYA/A7nAQFX9KP8gEfEAh/D9w0z93zcCN6rql4UsNxuwazdMhWJbFsaUzM/N4gV8p8+eXsi4VkA1IEFVqwG1gPXAP/IP9DeWjsAPqppR5omNKQXbsjCmZDoBqap6o4hUAZ4QkRRVfSffuK7ATz/vmlLVkyLyLdCygGW2wNdY7HiFqXBsy8KYkunMr8cVxgJLgNdFpFu+cV2BNQDi0wG4DniygGXa8QpTYVmzMKaYRCQWiMf/R11Vs4DhQAowX0Sa5hneFegrIunAMXxnT72pqm8UsGg7E8pUWNYsjCm+n/+o/7IFoKrpwPn+twtFpJb/GERn4DJVramq0UB34E4R6V3IchVYHbTkxpSQ3e7DmCARkdb4rrtIUNWdeaanAneo6iuOhTOmmOwAtzHB8/PB7Z0AIlIVuBXfFdofOhnMmOKyZmFM8HQFGojIUXzXYhwGluK7FcguR5MZU0y2G8oYY0xAdoDbGGNMQNYsjDHGBGTNwhhjTEDWLIwxxgRkzcIYY0xA1iyMMcYEZM3CGGNMQNYsjDHGBPT/9MxY+BqI8IUAAAAASUVORK5CYII=\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_m_exact_vals/np.sum(Nelems0), label = 'Exact')\n",
    "ax.plot(Tvals, avg_m_approx_vals/np.sum(Nelems0), label = 'Large N')\n",
    "ax.set_ylabel(r'$\\langle m \\rangle$', fontsize = 18)\n",
    "ax.set_xlabel(r'$k_BT$', fontsize = 18, labelpad = 10.5)\n",
    "ax.grid(alpha = 0.5)\n",
    "ax.axvline(x = kBTcrit(Dels0, Nelems0), color = 'k', linestyle = '--')\n",
    "# ax.set_ylim([0,50])\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": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "## permutation operator\n",
    "def perm(init_ls):\n",
    "    Ncomp = len(init_ls)\n",
    "    i1 = int(random.choice(range(Ncomp)))\n",
    "    i2 = int(random.choice(range(Ncomp)))\n",
    "\n",
    "    ## new omega vector \n",
    "    fin_ls = copy.deepcopy(init_ls)\n",
    "    fin_ls[i2] = init_ls[i1]\n",
    "    fin_ls[i1] = init_ls[i2]\n",
    "    \n",
    "    return(fin_ls)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Function to count the number of correct elements in a bound vector"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "## number of correctly placed elements\n",
    "def m_calc(vec, mstr_vec):\n",
    "    \n",
    "    num = 0\n",
    "    \n",
    "    for k in range(len(mstr_vec)):\n",
    "        if mstr_vec[k] == vec[k]:\n",
    "            num += 1\n",
    "    \n",
    "    return num"
   ]
  },
  {
   "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{m}) = \\sum_{i=1}^R m_i \\ln \\delta_i.\n",
    "\\label{eq:sim_en}\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "# log of boltzmann factor\n",
    "def log_boltz(bound_objs, mstr_vec, deltas, name_key):\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 corr_log_factor    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Checking logarithm of Boltzmann factor definition"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Energy for original ordering: -3.58801017672414\n",
      "Checking energy value: -3.58801017672414\n",
      "Number of correctly placed elements: 10\n",
      "-----\n",
      "Energy after switching first and last values: -3.58801017672414\n",
      "Checking energy value: -3.58801017672414\n",
      "Number of correctly placed elements: 10\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",
    "\n",
    "# sample master list\n",
    "sample_master = ['A', 'B', 'C', 'D', 'E', 'F', 'G', 'H', 'I', 'J', ]\n",
    "\n",
    "print('Energy for original ordering:', log_boltz(sample_master, sample_master, np.exp(-q1), name_key0))\n",
    "print('Checking energy value:', -np.sum(q1))\n",
    "print('Number of correctly placed elements:', m_calc(sample_master, sample_master))\n",
    "print('-----')\n",
    "changed_list = perm(sample_master)\n",
    "print('Energy after switching first and last values:', log_boltz(changed_list, sample_master, np.exp(-q1), name_key0))\n",
    "print('Checking energy value:', -np.sum(q1*(np.array(changed_list)==np.array(sample_master))))\n",
    "print('Number of correctly placed elements:', m_calc(changed_list, sample_master))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Metropolis Hastings algorithm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Metropolis Monte Carlo Algorithm \n",
    "\n",
    "## loads uniform random sampling \n",
    "runif = np.random.rand\n",
    "\n",
    "def met_gen_derang(Niter, bound_objs, mstr_vec, deltas, 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",
    "    # Initialize trace for state values\n",
    "    state_vals = [0]*(Niter+1)\n",
    "    \n",
    "    # Set initial values\n",
    "    state_vals[0] = bound_objs\n",
    "        \n",
    "    # Initialize acceptance counts\n",
    "    # We can use this to tune our step size\n",
    "    accepted = 0\n",
    "    \n",
    "    for i in range(Niter):\n",
    "    \n",
    "        # get current permutation\n",
    "        current_bound_objs = state_vals[i]\n",
    "        \n",
    "        # proposed new permutation; generated from random integer sampling\n",
    "        new_bound_objs = perm(current_bound_objs)\n",
    "        \n",
    "        #sets current log-probability\n",
    "        log_current_prob = log_boltz(current_bound_objs, mstr_vec, deltas, name_key)\n",
    "                \n",
    "        # Calculate posterior log-probability with proposed value\n",
    "        log_proposed_prob = log_boltz(new_bound_objs, mstr_vec, deltas, name_key)\n",
    "\n",
    "        # Log-acceptance rate\n",
    "        log_alpha = log_proposed_prob- log_current_prob\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",
    "            state_vals[i+1] = new_bound_objs\n",
    "            log_current_prob = log_proposed_prob\n",
    "            accepted += 1\n",
    "        else:\n",
    "            # Stay put\n",
    "            state_vals[i+1] = state_vals[i]\n",
    "\n",
    "    # return our samples and the number of accepted steps\n",
    "    return state_vals, accepted"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Computing microstate average from simiulations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "def avg_m(state_vals, mstr_vec, n_mc):\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(n_mc/50)\n",
    "    \n",
    "    ls = [0]*length\n",
    "    ls = np.array(ls)\n",
    "    for k in range(length):\n",
    "        ls[k] = m_calc(state_vals[n_mc-length+k], mstr_vec)\n",
    "    \n",
    "    \n",
    "    return(np.mean(ls))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Image grid for completely correct configuration"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "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": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 9.,  9., 10.,  5.,  7.,  6.,  3., 51.])"
      ]
     },
     "execution_count": 23,
     "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": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Temperature Run: 1 ; Current Time: 49.38 secs\n",
      "Temperature Run: 2 ; Current Time: 96.82 secs\n",
      "Temperature Run: 3 ; Current Time: 144.21 secs\n",
      "Temperature Run: 4 ; Current Time: 187.47 secs\n",
      "Temperature Run: 5 ; Current Time: 226.98 secs\n",
      "Temperature Run: 6 ; Current Time: 263.06 secs\n",
      "Temperature Run: 7 ; Current Time: 299.26 secs\n",
      "Temperature Run: 8 ; Current Time: 335.83 secs\n",
      "Temperature Run: 9 ; Current Time: 372.2 secs\n",
      "Temperature Run: 10 ; Current Time: 405.93 secs\n",
      "Temperature Run: 11 ; Current Time: 441.39 secs\n",
      "Temperature Run: 12 ; Current Time: 473.89 secs\n",
      "Temperature Run: 13 ; Current Time: 503.44 secs\n",
      "Temperature Run: 14 ; Current Time: 533.46 secs\n",
      "Temperature Run: 15 ; Current Time: 560.87 secs\n",
      "Total simulation run time: 560.8682470321655ecs\n",
      "------\n",
      "\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 = 20000\n",
    "\n",
    "# binding energy parameters\n",
    "R = 8\n",
    "Del_bar, sigma_D = 4.0, 2.0\n",
    "Dels = np.random.randn(R)*sigma_D+Del_bar\n",
    "    \n",
    "random.seed(0)\n",
    "\n",
    "# initial monomer and dimer states; \n",
    "# system in microstate of all correct dimers\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.05\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_m_vals = [0]*npoints\n",
    "\n",
    "# accepted list \n",
    "accepted_list = [0]*npoints\n",
    "\n",
    "# saved list for plotting\n",
    "saved_list = dict()\n",
    "\n",
    "for k in range(npoints):\n",
    "\n",
    "    final_m_vals = [0]*navg\n",
    "\n",
    "    for j in range(navg):\n",
    "        \n",
    "        # make copy of initial monomer and dimer states \n",
    "        bound_objs_copy = copy.deepcopy(bound_objs_0)        \n",
    "        \n",
    "        # defining helper functions\n",
    "        deltas_ = delta_func(Dels, temp_vals[k])        \n",
    "\n",
    "        # make copy of list \n",
    "        master_list_copy = copy.deepcopy(master_list)\n",
    "\n",
    "        # metroplois generator\n",
    "        bound_list, accepted = met_gen_derang(Nmc, bound_objs_copy, mstr_vec, deltas_, name_key_)\n",
    "\n",
    "        # calculate final j value\n",
    "        final_m_vals[j] = avg_m(bound_list, master_list, Nmc)\n",
    "\n",
    "    # saving averaged values        \n",
    "    sim_m_vals[k] = np.mean(np.array(final_m_vals))\n",
    "\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",
    "    time_prelim = time.time()\n",
    "    print(\"Temperature Run:\",str(k+1),\"; Current Time:\", round(time_prelim-t0,2),\"secs\")\n",
    "\n",
    "t1 = time.time()\n",
    "print(\"Total simulation run time: %secs\" % (t1-t0))\n",
    "print(\"------\")\n",
    "print()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Simulated image grid at various temperatures"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAgQAAAIRCAYAAAA1GaW+AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAkgklEQVR4nO3de7Ckd1kn8O/DDCHAlKKiRhTE66gDZXBKDVo6aHlbrRFM3BVcBaxa0eyo0fK2jmhMqbPrWijRSkUHS8ULXqgE3F6vWMpIgYIKFJCF8YLIbXGDKDrGQDL89o/3PZnOyemec7rfPqffcz6fqq7T5708/Xv7dD/z7V+/3VOttQAAB9uD9noAAMDeEwgAAIEAABAIAIAIBABABAIAIDMCQVWdrKqzk8mkJXFxcVmPy6joIy4ua3mZqS7zPQRzVwK7qvZ6AAvSR2B9zOwjh7ez97lz5wYZxYkTJ9a+5nS98xeWLpckOXrk0vWDVHO6Xi6eXr5gkhw6M76aJFnNc/6W64bJSKduu5RZ9JH1qqmP7B7nEAAAAgEAIBAAABEIAIAIBABABAIAIAIBABCBAACIQAAARCAAACIQAAARCACACAQAQAQCACACAQAQgQAAiEAAACSp1tq89XNXAruq9noAC9JHYH3M7CNbzhBU1cmqOjuZTFY3JGBf00dgXMwQwHiYIQCWNbOPHN7O3ufOnRtkFCdOnFhpzfMXlq939Mil60PUO8g1Vz3GXDw9TNFDZ1ZbkyT6iJrrUW9zTX3kEicVAgACAQAgEAAAEQgAgAgEAEAEAgAgAgEAEIEAAIhAAABEIAAAIhAAABEIAIAIBABABAIAIAIBABCBAABIUq21eevnrgR2Ve31ABakj8D6mNlHtpwhqKqTVXV2MpmsbkjAvqaPwLiYIYDxMEMALGtmHzm8nb3PXxhmFEePrH/NMYxxLDWn6+Xi6eULJsmhM/ddXcVxr2KcdNb1cbqKmmMY41hq6iO7x0mFAIBAAAAIBABABAIAIAIBABCBAACIQAAARCAAACIQAAARCACACAQAQAQCACACAQAQgQAAiEAAAEQgAAAiEAAASaq1Nm/93JXArqq9HsCC9BFYHzP7yJYzBFV1sqrOTiaT1Q0J2Nf0ERgXMwQwHmYIgGXN7COHt7P3+QvDjOLokfWvOYYxjqXmGMa4GzXpjOXv5bG/XjXHMMbdqLkbnFQIAAgEAIBAAABEIAAAIhAAABEIAIAIBABABAIAIAIBABCBAACIQAAARCAAACIQAAARCACACAQAQAQCACACAQCQpFpr89bPXQnsqtrrASxIH4H1MbOPbDlDUFUnq+rsZDJZ3ZCAfU0fgXExQwDjYYYAWNbMPnJ4O3ufvzDMKI4eWf+a0/Vy4+nlCybJTWfuu7qux72KmtP16uZrly+YpN1w++hq0lnXx+kqauojw9XUR3aPkwoBAIEAABAIAIAIBABABAIAIAIBABCBAACIQAAARCAAACIQAAARCACACAQAQAQCACACAQAQgQAAiEAAAEQgAACSVGtt3vq5K4FdVXs9gAXpI7A+ZvaRLWcIqupkVZ2dTCarGxKwr+kjMC5mCGA8zBAAy5rZRw5vZ+/zF4YZxdEjq62ZG08vX/CmM8PW24Wa6/r3ma5XN1+7fMEk7YbbR1eTzro+TjfX1EeWo4+Mt484qRAAEAgAAIEAAIhAAABEIAAAIhAAABEIAIAIBABABAIAIAIBABCBAACIQAAARCAAACIQAAARCACACAQAQJJqrc1bP3clsKtqrwewIH0E1sfMPrLlDEFVnayqs5PJZHVDAvY1fQTGxQwBjIcZAmBZM/vI4e3sff7CMKM4emTqlxtPD1P0pjPD1hy63khrDvE3n/57r+IxdMt1w/z7eOq2S/9ereSxThJ95CDW1EeWs9t9xEmFAIBAAAAIBABABAIAIAIBABCBAACIQAAARCAAACIQAAARCACACAQAQAQCACACAQAQgQAAiEAAAEQgAAAiEAAASaq1Nm/93JXArqq9HsCC9BFYHzP7yJYzBFV1sqrOTiaT1Q0J2Nf0ERgXMwQwHmYIgGXN7COHt7P3+QvDjOLokalf7jg9TNFjZ4atOXS9kdYc4m8+/fdexWPoluuG+ffx1G2X/r1ayWOdJPrIQaw5hj5y7ty5QWqeOHHivutj7SNOKgQABAIAQCAAACIQAAARCACACAQAQAQCACACAQAQgQAAiEAAAEQgAAAiEAAAEQgAgAgEAEAEAgAgAgEAEIEAAEhSrbV56+euBHZV7fUAFqSPwPqY2Ue2nCGoqpNVdXYymaxuSMC+po/AuJghgPEwQwAsa2YfObydvc9fGGYUR49M/XLH6WGKHjszbM2pegfquDfVHOLYp497FfflLdcN8+/jqdsu/Xu1kr85SQ7Y80kfSaKPLGu3+4iTCgEAgQAAEAgAgAgEAEAEAgAgAgEAEIEAAIhAAABEIAAAIhAAABEIAIAIBABABAIAIAIBABCBAACIQAAARCAAAJJUa23e+rkrgV1Vez2ABekjsD5m9pEtZwiq6mRVnZ1MJqsbErCv6SMwLmYIYDzMEADLmtlHDm9r71uHGUW7fv1rjmGMY6k5hjHuRk06Y/l7eeyvV80xjHE3au4GJxUCAAIBACAQAAARCACACAQAQAQCACACAQAQgQAAiEAAAEQgAAAiEAAAEQgAgAgEAEAEAgAgAgEAEIEAAEhSrbV56+euBHZV7fUAFqSPwPqY2Ue2nCGoqpNVdbaqfrjfebBLVX3lutccwxjHUnMMYxxLzao6mRHRR9Z/jGOpOYYxjqXmvD6yZSBorU1aa89KctWsHZewiqY2dM0xjHEsNccwxrHUHFUg0EdGMcax1BzDGMdSc2eBYMpk4IGMpeYYxjiWmmMY41hqrmKMu2EM9+0qao5hjGOpOYYxjqXmzHqXO4cAADgAfMoAABAIAACBAACIQAAARCAAACIQAAARCACACAQAQAQCACACAQAQgQAAiEAAAEQgAAAiEAAAEQgAgAgEAEAEAgAgAgEAEIEgVfW4qrq3qr54atmTqqpV1TOnlj25qt5fVZ800O3+RVW9fohawO7bQe/QTxiFAx8IkvxEkpe31l4yb6PW2m8leX2SH1v2BqvqcJLHJXntknWOVNXFvtls5/Khy459zlgeVFXfUVVvqqq7q+ptVfWcqnr4KmrMOcYLwx4ZzLSt3rEV/eSyYxmin3xkVf1Mv+/7q+qtVXVzVT1ii22/r6peWFVv7o/tLUMez1gc3usB7KWqemKSL07ylE2r/iTJQ5Pcs2n5zUmeX1XHWmt3LHHTn5bkIVnyCZzu7/eMTcuuT/I5Sb4ryT9MLX9fa+09S97ePD+Z5NuSvCjJc5J8av/7E6rqi1prH1hBjZclObtp2ea/GQxuh71DP9m5pfpJVX1EklcmeVSSn03yhnSh6fokn19Vn9tau2tqlzNJ3pPk1UkeMeyhjEhr7cBekvxykjuTPHib2x9J8m9JfnrJ2316kpbkC1dwTH+Z5N+THN7F+/FYkg8kuW3T8m/tj/Nrh67RL/vFvX4MuRzMy057x4wa+snWtzlEP3luv+3TNi1/Wr/82ZuWf/zU9TckecteP8b24nJg3zLop9mekuQPW2v3bFr3gPf8kqS1diHdq9KvXvLmn9D/fO3UbT6iql7UT489a5GiVfXgdCn4da21e5cc4048LUmlexJOe16Su5J83apqVNUVVXVkB2OFpey0d+gnOzZEP/mCdEHm1zct/40kdyf5humFrbU3LzLQ/eYgv2VwPF1Cf9UO9/vTJF9aVZ/SWnvTgrd9dZK3tX7KraqOJ3lhuifB57bW/nLBuseSXJHkNZfbsKoelGQn7wG+p82epvvMdIn+fvdla+3uqnptv/5yFqnx1emaw6GqujPdk/3ZrbX3buP2YFGL9o6t6CcPNEQ/eUiSu1v/kn+qxgeq6t+TfHxVPbK19u4djHnfO7AzBOned0uSv93hfhvbH1vitq9On+ar6vokL0/yxiTHl3jyJpdeKbx6G9s+Jt2U53Yvj5lT61FJ3t1ae98W696R5JFVdcVlxrPTGq9K8kPpQsEzkvxRkm9J8jIzBqzYor1jK/rJAw3RT+5I8iFVdfX0wv73D5kaM1MO8gzBh/c/d3pizD/2Pz9ikRutqsemO2nlr6vqBUm+Jt0/bD+yOc0u4DP6n5dN9Enele6kqO1615x1D0uy1ZM36abnNrZ5/1A1WmufvWmbX6qq1yX50SQ39D9hFRbtHVvRTx5oiH7y3HRv6/xmVX17uvMCjvXL70ny4L4GUw5yINh4stQO99vYftEn20bq/tZ002Jf1rb42FI/BffedLM4rf/5V0lOtdZePqf2vek+zjRXa+3uJH+449Fv7a7MbmhXTm2z6ho/nuTGJF8RgYDVWbR3bEU/eaCle0Fr7WVV9dQkP5Xkt/vFF5P8XLrZg69K8i/LD3V/OchvGdzZ/9zpZ2k3tr9z7lazbTyBn5fufa7Hz9juaLr3KR/dWjuSbprrTek+jvMA/RP+05O8sX9yzlVVh6rqqh1cDs0p985003gP2WLdR6eb/puX5gep0Z/g9c4kj7zMbcEyFu0dW9FPHmiIfpLW2guTfEy6++jzkzyqtfbN/bJ7k/zN5WocNAd5huAN/c+dflPYJ27af6euTveAPlVVD03y41X1ltba7Zu2O57knRsnCrXW3ldVr0nyyTPqflK6J/x23u9Lkkcn+bsdjPvjkrxlxro/T/IlST4r3VnTSZKqujLd8f7JNuovXaPf9mOS/Nk2bg8WtWjv2Ip+8kBD9JMkSWvtYu7/6Yur0gWEc+3+30NADnYgeE26KaNrdrjfNUn+obV2fsHbfUIuvSf3TUk+NsmvVNWTWmvTZ9UeT/K6JKmqSpf8vzHJD86ou5P3+5Jh3/P7jSSnk3x7pp7A6cb7sCS/urGg/yjTJyS5q7X21gVrfFhr7R/zQD+c7jE9uezRwOIW7R1b0U8eaIh+8gD9rMdPJTkUbylu6cAGgtbaxaq6PclTquohM85ovZ/+7PXPS/Lzi9xmVX1YulewL+jHcE9VXZfkFUkmVXVNa20jZR9P8llV9c/pPvpzZZIfba29YEb5nZwRPOh7fq2111fVLUm+pb9PfyeXvlnsXPrj7X10ujOgzyV50oI1nl1V1yT54yRvTfdK5svTffb4lUl+eojjgq0s0ju2op/MrLV0P+nv21el+6bDv0vywem+3+B4ku9vrf3x9G1W1denC1NJd9LoFVX17P73v2+t/fIQx7b29vqbkfbykm5KqiW5btPyJ/XLn7lp+TP65Y9b8Pa+qN//qZuWf1y6rwV9Y7r39h6U5F+TXLtprO9LcmJG7ZekO6nog/bovjyU5DuTnO/H+Y503/V+ZNN2j+3vg5cuUePJSX6/X393um97e226VxVX7vXjymX/X3bSO/SThY5tqX6SLvT8WrowcHe6T4T8fpIvnXF7L+3rbHV56SqOcR0v1d8ZB1ZV/V6Sh7fWPm8b27463VdaXrviMX1Kuifzo1trb59a/u4k391a+4VV3j5weTvpHTP2109YKwf5UwYbvjPJE6vqS+ZtVFVPSfc1nt+7C2PaOAHo7f1tP7yqfiBd6v2DXbh94PK21Tu2op+wjg7sOQQbWve/jF32fmitvTjdE2g3HE9yVXX/le8H0p3A9Ip0X0P6jl0aAzDHdnvHjH1fHP2ENXPg3zIAALxlAABEIAAAIhAAABEIAIAIBABABAIAIAIBABCBAACIQAAARCAAACIQAAARCACACAQAQAQCACACAQAQgQAAiEAAAEQgAAAiEAAAEQgAgMwIBFV1sqrOTiaTlsTFxWU9LqOij7i4rOVlpmpt7vq5K4FdVXs9gAXpI7A+ZvaRw9va+9ZhRtGuv3T93Llzg9Q8ceLEoDWn652/sHS5JMnRI5euH6Sa0/Vy8fTyBZPk0Jn7rq7iuFfxWKeziuf8uvYmfaRTN1+7dL12w+2XfllBH8mNA9W8abW9aTc4hwAAEAgAAIEAAIhAAABEIAAAIhAAABEIAIAIBABABAIAIAIBABCBAACIQAAARCAAACIQAAARCACACAQAQAQCACBJtdbmrZ+7EthVtdcDWJA+AutjZh/Zcoagqk5W1dnJZLK6IQH7mj4C42KGAMbDDAGwrJl95PB29j5/YZhRHD1y6XrdOkzNdv2l6+fOnVu63okTJ+67fst1w/TfU7dd6oeruC9z8fQwRQ+due/qEOOcHuNYjnsVj0s6Qzw/k/s/R1fxuBrDY/8g1Vz1GOvmawep2W64/dIvK+hNu8FJhQCAQAAACAQAQAQCACACAQAQgQAAiEAAAEQgAAAiEAAAEQgAgAgEAEAEAgAgAgEAEIEAAIhAAABEIAAAklRrbd76uSuBXVV7PYAF6SOwPmb2kS1nCKrqZFWdnUwmqxsSsK/pIzAuZghgPMwQAMua2UcOb2fv8xeGGcXRI+tfcwxjHEvN6Xq5eHr5gkly6Mx9V1dx3OfOnRuk5okTJwaps5+s6+N0FTWn69Wty9dLknb9pevretyba95y3fIZ9tRtU3lyBX1kJc/5Gwca501nLr/NgJxUCAAIBACAQAAARCAAACIQAAARCACACAQAQAQCACACAQAQgQAAiEAAAEQgAAAiEAAAEQgAgAgEAEAEAgAgAgEAkKRaa/PWz10J7Kra6wEsSB+B9TGzj2w5Q1BVJ6vq7GQyWd2QgH1NH4FxMUMA42GGAFjWzD5yeDt7n78wzCiOHlltzbr52qXrtRtuv+/6WI57XWuOYYy7UZPOWP5eQz/2h+hLyep70yrGmTtOL1/w2Jn7rq7r33tzzVwc4LiT5NCZy28zICcVAgACAQAgEAAAEQgAgAgEAEAEAgAgAgEAEIEAAIhAAABEIAAAIhAAABEIAIAIBABABAIAIAIBABCBAACIQAAAJKnW2rz1c1cCu6r2egAL0kdgfczsI1vOEFTVyao6O5lMVjckYF/TR2BczBDAeJghAJY1s48c3s7e5y8MM4qjR9a/5nS93Hh6+YJJctOZ+66u63GvouZ0vbr52uULJmk33D66mnTW9XG6ipqr7iO3XDdMNjx126WsNor78o6B7stjl+7LunWYku36qV8uDjTOQ2cuv82AnFQIAAgEAIBAAABEIAAAIhAAABEIAIAIBABABAIAIAIBABCBAACIQAAARCAAACIQAAARCACACAQAQAQCACACAQCQpFpr89bPXQnsqtrrASxIH4H1MbOPbDlDUFUnq+rsZDJZ3ZCAfU0fgXExQwDjYYYAWNbMPnJ4O3ufvzDMKI4eWW3N3Hh6+YI3nRm23i7UXMl9eXGAcR66NMa6+drl6yVpN9x+6Zc7Brovj634viTJePrIEDWn691y3TA57tRtU7lqJH1k6Ptybf/dSO7f54fon8n9euhucFIhACAQAAACAQAQgQAAiEAAAEQgAAAiEAAAEQgAgAgEAEAEAgAgAgEAEIEAAIhAAABEIAAAIhAAABEIAIAk1Vqbt37uSmBX1V4PYEH6CKyPmX1kyxmCqjpZVWcnk8nqhgTsa/oIjIsZAhgPMwTAsmb2kcPb2fv8hWFGcfTI1C83nh6m6E1nhq05Ve9AHfemmkMc+/Rxr+K+vOW6Yf59PHXbpX+v6tZBSqZdP0yd/WQVj4G1fY6OsY/cMVDNY+PqI+vaP5NN49wFTioEAAQCAEAgAAAiEAAAEQgAgAgEAEAEAgAgAgEAEIEAAIhAAABEIAAAIhAAABEIAIAIBABABAIAIAIBABCBAABIUq21eevnrgR2Ve31ABakj8D6mNlHtpwhqKqTVXV2MpmsbkjAvqaPwLiYIYDxMEMALGtmHzm8nb3PnTs3yChOnDhx6Zc7Tg9SM8fO3Hf1luuW75enbpvqXSsY41hqnr+wfLmjRy5dH6Le5pq5ONBxHxr2uJNN4yTJau7bunWYmu36S9eHfuyPpX+uouYQf5+h/zbJpsfQzdcOUrPdcPt918faR5xUCAAIBACAQAAARCAAACIQAAARCACACAQAQAQCACACAQAQgQAAiEAAAEQgAAAiEAAAEQgAgAgEAEAEAgAgAgEAkKRaa/PWz10J7Kra6wEsSB+B9TGzj2w5Q1BVJ6vq7GQyWd2QgH1NH4FxMUMA42GGAFjWzD5yeDt7n78wzCiOHpn65Y7TwxQ9dmbYmlP1DtRxb6o5xLFPH/dK7ssbBzrum1b8NyfJau7bunWYmu36qV/0keFqXhyg5qHx3Zdj7SNOKgQABAIAQCAAACIQAAARCACACAQAQAQCACACAQAQgQAAiEAAAEQgAAAiEAAAEQgAgAgEAEAEAgAgAgEAEIEAAEhSrbV56+euBHZV7fUAFqSPwPqY2Ue2nCGoqpNVdXYymaxuSMC+po/AuJghgPEwQwAsa2YfObyt3S+eHmYYh87cd7VuHaZku/7S9SFqDl1vc83zF4apefTIamsOfl/efO3yBZO0G26/7/q6Hndy/2OnM5a/1xDjnB5jbhyof960/v1zc82h78uVHPdIetNucFIhACAQAAACAQAQgQAAiEAAAEQgAAAiEAAAEQgAgAgEAEAEAgAgAgEAEIEAAIhAAABEIAAAIhAAABEIAIAk1Vqbt37uSmBX1V4PYEH6CKyPmX1kyxmCqjpZVWer6of7nQe7VNVXrnvNMYxxLDXHMMax1KyqkxkRfWT9xziWmmMY41hqzusjWwaC1tqktfasJFfN2nEJq2hqQ9ccwxjHUnMMYxxLzVEFAn1kFGMcS80xjHEsNXcWCKZMBh7IWGqOYYxjqTmGMY6l5irGuBvGcN+uouYYxjiWmmMY41hqzqx3uXMIAIADwKcMAACBAAAQCACACAQAQAQCACACAQAQgQAAiEAAAEQgAAAiEAAAEQgAgAgEAEAEAgAgAgEAEIEAAIhAAABEIAAAIhAAADlAgaCqHldV91bVF08te1JVtap65gLLnlxV76+qTxpofH9RVa8fohawGtvtIzuoN2gf6WvqJSzkwASCJD+R5OWttZcMUay19ltJXp/kx5atVVWHkzwuyWuXrHOkqi72zWk7lw9dduxzxvJ9VfXCqnpzf1tvWaDGR1bVz1TV2/qm+daqurmqHrHFtkeq6nRVvb6q/rWq3l1Vr6iqZ1ZVDXFMkDXuI8kwvWRd+khVHa2qX62qN1bVe6vqrqp6U1X9RFV91A7qPKiqvqPf9+6+nzynqh4+Y/tZx3lhuKNbT4f3egC7oaqemOSLkzxl06o/SfLQJPcssCxJbk7y/Ko61lq7Y4khflqSh2TJQJDu7/mMTcuuT/I5Sb4ryT9MLX9fa+09S97ePGeSvCfJq5M8Yqc7V9VHJHllkkcl+dkkb0jX6K5P8vlV9bmttbv6bR+U5HfTHefzk/x0kocleVqSX0jyqUm+d7nD4aDbYR/ZiaH6SDJML1mXPvIxST4qyYuSvD3JvUken+RZSZ5aVVe31v7fNur8ZJJv6+s8J10/+LYkT6iqL2qtfWCLfV6W5OymZYv+fcejtbbvL0l+OcmdSR48cN0jSf4tyU8vWefpSVqSL1zBsf9lkn9PcniX7/OPn7r+hiRv2eH+z+3vk6dtWv60fvmzp5Y9sV/2k5u2vSLJm5P8824eu8v+vKx7H+lrraSX7FUfmTGW/9gf4/dsY9tjST6Q5LZNy7+1r/G1W+zTkvziXh/nXlz2/VsG/RTaU5L8YWvtnk3rFj6HIElaaxfSJcmvXnKYT+h/vnbqNh9RVS/qp7ietUjRqnpwulfVr2ut3bvkGHektfbmJUt8QboG9Oublv9GkruTfMPUsg/qf75z0xjen+Td6ZotLGwnfaSq/kP/+7fNqPWnVXVn//wcso8kK+gle9lHZvj7/ueHbGPbpyWpdC8wpj0vyV1Jvm7WjlV1RVUdWWSAY3UQ3jI4ni6Bv2pF9f80yZdW1ae01t60YI2rk7yt9VNvVXU8yQvTPZA/t7X2lwvWPZbuVfJrLrdhP+2+k/cC39O2nmobykOS3N36yL6htfaBqvr3JB9fVY9srb073d/2n5N8T3+uwivTvWXwjHR//29e4Tg5GHbSR/4gybvSvVr/qekV/cmD1yT5qU3BYog+kqyml2y7j/S3OWgvqaor0933V6Z7S2TjfIvf2Ubtz0w3Q3C/v1tr7e6qem2/fitfnS4sHKqqO9O9EHl2a+2927jN0dr3MwTpHkBJ8rcrqr9R99gSNa5On+ir6vokL0/yxiTHlwgDyaVXC6/exraPSTcdut3LY5YY13bckeRDqurq6YX97xuvDB6TJK21f0rylenOWfjNdK8g3pjkVJLrWmvPW/FY2f+23UdaaxeT/EqS41X1aZtWP73/+fxNy4foI8lqeslO+kgyfC/5L/12b0vy++nOSfq61trLtjGWRyV5d2vtfVuse0eSR1bVFZuWvyrJD6ULBc9I8kdJviXJy/b7jMFBmCH48P7nqk58+cf+50cssnNVPTbdA/yvq+oFSb4m3YPxRza/Ol7AZ/Q/t5Ps35XuhKntetfOh7Mjz003RfubVfXt6c5DONYvvyfJg9PNAmy40G/zv5K8It0rlFNJXlBVT24DnRXOgbXTPvL8dCfgPT3Jf0uS/tMuX5fkDa21zf+4LtVH+vqPzWp6yU76SDJ8L3lxkjelmyV4Qrrw/8ht1n5Ykq3CQNK99bixzfs3FrbWPnvTdr9UVa9L8qNJbuh/7ksHIRBsPBFW9dGzjbqLPuE20ve3ppva+rKt/vHqp+Hem25Wp/U//yrJqdbay+fUvjfdx5rmaq3dneQPdzz6FWmtvayqnppuyvW3+8UXk/xcutmDr0ryL0lSVY9PFwK+o7X2Mxs1qurX0oWE51XVJ/Sv3GARO+ojrbU3VNWrk/znqjrdT4l/fpLHJvmeLXZZto8kq+sl2+4jyfC9pLX29nSfMkiSF1fVbUn+vKoe1lr775fZ/a7MDllXTm1zOT+e5MYkX5F9HAgOwlsGd/Y/V/WZ+426d87daraNJ/Hz0r1v/vgZ2x1Nl5Af3Vo7km7a/E3pPlLzAP2T/tOTvLF/gs5VVYeq6qodXA7t8Dh3rLX2wnQfPXpCumb6qNbaN/fL7k3yN/2m35Huyf3CTfvflS5MfGy6RgyLWqSP/FK6x+oX9r8/PV2o/ZUttl22jyQr6CU77SP9PivtJa2116Wbrfiv29j8neneFnjIFus+Ot3bCe/fYt3m27xno9ZOxjo2B2GG4A39z8G+CWyTT9x0Ozt1dboH5amqemiSH6+qt7TWbt+03fEk79w4Wai19r6qek2ST55R95PSPem3+77fo5P83Q7G/XFJ3rKD7RfSv6p/7cbvVXVVusZ3rv8HP+me2EmyVWM5vOknLGKRPvKCdK8sn15VL0/3nvRLWmv/d4ttl+0jyWp6yU77SLI7veSh2V44+/MkX5Lks9J9kiPJfScqXp3uOyQuq9/+Y5L82Q7HOSoHoUm+Jt3U8jUrqn9Nkn9orZ1fcP8n5NJ7c9+U7tXsr1TVk1pr02fGHk/yuuS+9yIfn+Qbk/zgjLp7/b7ftvQfafqEJHe11t56mW0flO4thEO5/7Td/0n3pH9mkv85tf0jkjw5yT/l0mwCLGLHfaS1dmdV/W6Sa9P9w/NBeeDJhBuW7SPJanrJTvtIMlAvqaqrWmsPWFdVX5DuY5Av3bR8q17yG0lOJ/n2TAWCdMf7sCS/uqnGh7XW/jEP9MPp/r2cXP5wxmvfB4LW2sWquj3JU6rqITPONl1If8bp5yX5+QX3/7B0qfMF/Vjvqarr0r0fPqmqa1prG0n7eJLPqqp/TvcRoCuT/Ghr7QUzyu/ozOCh3/erqq9P15CS7oSsK6rq2f3vf99a++X++kenOwv6XJInTe2/8RGvF6V7tfHB6T5TfDzJ97fW/njq5p6bbjr2f/TnE7w83auHb0z3TWennD/AMpboI89PdxLcc9K9b//izRss20f6GqvqJTv9hMGQveTW6r6i+I/SfXLoynRjf2qSf03ynZu2f0Avaa29vqpuSfIt/d/vd3LpmwrPpb+/pjy7qq5J8sdJ3ppuduTL030vyivTfQvq/rWX34q0W5d000Ut3UfQppc/qV/+zJ0u65c/o1/+uAXH9UX9/k/dtPzj0n096BvTvb/3oHRPgGs3HdP7kpyYUfsl6U4s+qA9us9f2h/bVpeXTm332M3L+uVXJPm1dGHg7nRnd/9+ki+dcXufkK75vj3dpxD+Jd2rsmtXdYwuB+uykz4yte6KdJ8gaEmeN6PuUn2kr7GSXrKXfSTJf0ryv9N93PDudF9U9qZ0/yg/ZovtZ/WSQ+nCw/n+ON+R7v+kOLJFjSf3feYd/W3+W7q3LE8nuXKvH4OrvlR/J+x7VfV7SR7eWvu8AWu+Ot1X8l47VM0Zt/Mp6Z7Qj27dGbcby9+d5Ltba7+wytsHOmPuI/1t6SXMdBA+ZbDhO5M8saq+ZIhiVfWUdO9j7cZ/mrNxEtDb+9t+eFX9QLpXH3+wC7cPdMbcRxK9hDn2/TkEG1r3v4gNdryttRenexLthuNJrqruv9/8QLrp8Fek+yrSd+zSGODAG3kfSfQS5jgwbxkAALMdpLcMAIAZBAIAQCAAAAQCACACAQAQgQAAiEAAAEQgAACS/H/J4LIXx7MVbQAAAABJRU5ErkJggg==\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('derang_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": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAdYAAAFLCAYAAABr4q3HAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAABPBElEQVR4nO3dd3hUZfbA8e9JJSQBQoIgIbRIbwEiUkSxISCIilhRUQQVFXUV17qWZdEVsKy6KD9RrKBgRRB1VRCMjZJA6FJC6KSSninv749JhiQkkDKTmSTn8zzzwNx5773nzQTOfd9777lijEEppZRSruHj6QCUUkqp+kQTq1JKKeVCmliVUkopF9LEqpRSSrmQJlallFLKhTSxKqWUUi5U7xPriBEjDGAOHjxogAb3qiv9HjZsmBk2bFiD7HtD/t6179p3D/bbbfzcuXFvkJKSAoDNZvNwJJ5RV/o9aNAgl2+zrvTdHbTvDVND7bu39bvej1grKz09nZYtW7Jr1y727t2LiLB27VoANm3aRGRkJDk5OR6Osv567rnneO655zwdhlJK1Zgm1iIzZ85k1KhRREdHExUVxaFDh4iJiQGgV69eDBw4kBdffNGzQSqnUx0IFb9PSEgA9MBIKVW7NLECubm5vPXWW0yaNAkAX19fWrVqhZ/fiZnyW2+9lblz52K1Wj0VZr02btw4xo0bV+n2pzoQKn7fo0cPQA+MlFK1SxMrsHz5ckSEIUOGAJw0AgIYPnw4aWlprFy50kNR1m+DBg2q9HnW0x0I6YGRUsqTNLECq1evpn///ohIhW0CAgKIiYlh1apVtRhZw/HQQw/x0EMPVart6Q6Eyk4Fgx4YKaVqjyZWICkpidatW5+2XevWrdm7d6/7A1KnVJkDobL0wEgpVVs0sQJ5eXk0atTotO2CgoLIy8urhYganmHDhjFs2LBKta3sgVBZemCklKoNmliBiIgI0tPTT9suLS2NFi1a1EJE6lQqeyBUlh4YKaVqgyZWoG/fvmzZsuW07RITE+nXr18tRKROpbIHQmXpgZFSqjZoYgUuvfRStm7dSmpqaoVt9u7dy4EDBxg+fHgtRqbKU9kDobL0wEgpVRs0seK4z3HAgAEsWrSowjYLFy5k+PDhtGvXrhYjU+WpzIFQWXpgpJSqLZpYizz11FP85z//wWaz0b59e4wxxMbGAlBQUMDcuXN58sknPRylgsodCJWlB0ZKqdqiibXIiBEjuPvuu9m/f/9JnyUlJfH4448775tUnneqA6Hi93369AH0wEgpVbvq/dNtqmLatGnlLu/cuTOdO3eu5WjUqZQ8EDrdKFQPjJRStUkTq6qzKjoQKksPjJRStUkTq/IKo0eP9nQISinlEppYlVeobJ1gpZTydnrxklJKKeVCDTqxxu1KYcjzPxK3K6Ve7cub2ew2MvIzSDqexO7M3VjsFqBqtYKVUsqbeXwqWETuASYCvYCFxpiJp2j7APB3oDGwBLjLGFNQnf3G7Uph0oK15FlsTFqwlvkTYxkcHVGdTXnVvor3N33xRmaN701UgOu2a4yhwFZAjiWHXEsuOdYcAnwD6Ni0Y6l2q5JX8VPyT2RbsjlecJzMwkwyCzI5XnicrMKsUm17R/Tmw8s+ZOLEiaX2U/bJNSX75M6fnVJK1ZTHEytwEJgBXAoEVdRIRC4FHgEuLFrnc+CZomVVUjLRAW5NeO7YlzEGu7FjNVZsdhuN/RuX3t+7v1AgR5i0cCd3DIpggG8yhbZCCmwFFNgKyLfmO/8e4h/CdV2vK7X9Xw78wlub3iLPmkeeNY9ca64zmdqMrVTbnmHnsG/LDaUS3ta0rXy689NK9eWssLMASiXWNxLeYPme5fRr2Y8hrYfgU9CZez/YVmsHJkopVRMeT6zGmM8ARCQWaHOKprcA840xm4va/xP4kCom1rKJrljZ/7QnvzuArKJpytN54aJXadv2XOd7S2EuN344GJuBQquN6CgwRZ+Zotc/fwS/VT4gwkdXfU1o0xNdP3A4nptW3IINgx2wAXbAjsECWEsM5oKM8MfEjaX6dkbAWtLafgHA/P2OV0XaSuBJiTXjcDxrj6ytVN/3H9zBgYy8Uj+7UHxPuU6ofyhNAptgtVvp2rwrACkpjinyiIgINqduZu/xvew9vpfPdn6GMT7QOoqA7M4U5HRh0gI78ycO0OSqlPJKHk+sVdAD+LLE+wSgpYiEG2MqVTS2oqRarGRy3WbPJcOncg/Szi84Xuq9MXa2+hTt45Q5xg6AzVZYen1bIcfEXql9W4q2UbJvjf0tpFVqbSiw5p+0LCg/q5yWDv7GEGy3E2w3NDZ2WlrTSab0z+7s7AyeSEkjxG6nid1OU58AmgY0pWmjMEKCIvALjIBgx8tIGABXX301ACtXrmRH+o5S+xSx49c4Cb/GSQTyPXZrMJOXd2HawPHc1v9S/H38K9lbpZRyv7qUWEOAzBLvi/8eCpRKrCIyBZgCEBkZSXJyMqmpqTyw7GiFSbVYnsXGA4vWE3CmASqXWFNTU0hulOx8b7HkVmq9YgcOHyQr50RySE05/QVOPsbgZ8AfwxdxW/j78r0UWB3jYh9bAF0KCvHHEGAM/gYCjKGRMQQWv+yOP5v4hZCcnFxq220Lgnjr0BGCjCGoKIEG2x0JtWwK+93elW+K/p5nsXHrO3+wpEcW12Zll2iVDznHgWTKyo8cQkrjHhQUOE6VJycnM/+Ma0j46zvezcxnd1A2BY1SEDHOdXz8cqDJel7bsp6U9GPc3O3icn9GVSnSX99o3xumhtr36vQ7KirKDZE41KXEmg00KfG++O8nDa+MMfOAeQCxsbGm+Af40nVtTjliBQjy9+Wl6/rRzPIv7JWcCo5ufxFBjZs739ttVmYkP8hbq/dgsZYceZ5I1GKEAF8/7hzWkS6d+uPnf+L08pktm/O/0JfwFcFHfPAVH8ef+ODn44uf+OIjJy7oHvLBYWdSBdhb0Jsue8IoAAqLJqF9MLQIDeDF62MAA8YOxoBfIJT9BQu9irNaRILdBnar42UrBGsBSccy+GrdHnzshfhjY585o9SqBVbD4q25nNW0LUH2PMjLgFP8HINadSYqKorAwECg6Jd90w6i9n3LaIBMOCr+fBF0JiuDQtjRuIACP8f2jC2QZevP5PHhJ+K32C0czD5IuybtTmyvgdK+N0wNte/e1O+6lFg3A32AT4re9wGOVHYaGGBwdATzJ8ZWmFyD/H1LXBhzRbUD9fH1Y+z5E2nRpuKp5yB/X+bfXP5FOH4BwbQ8q/xRWHlmjS+9nxyCWG9Kl/AL8vdl/rWx0KES5yWbRTle5bjh+R85UNj3lKu/V3gBP1hG8csjFzqSd0EW5KZAbhrkpkJOiuPPrMPQbtDJG0jfW+rtGcbClNx9TMkFeypsDQhgeXAwRwlg0jmlJ71/O/gbU3+YSr8z+jHijBGMjxyPn09d+jVXStV1Hv8fR0T8iuLwBXxFpBFgNcZYyzR9D1ggIh/iuCr4CWBBVfdXUXItnVRdo7b2VbUDhpqZNb53pUb9s8b3drwRgUZNHK/mHStcp5SBd0G7QRz7ay05e9bRXg45P/IBehQW0qOwEEiHkLxSq369+2sA1h9dz/qj63lnzztc1/U6xnUaR9PAplXpqlJKVYs3FIh4AsjDcXXvhKK/PyEibUUkW0TaAhhjVgAvAD8B+4Ak4Knq7LA4EQX5O64sckdSre19ld1PsUA/cUsSL7ufYi7pX/QFcO4DtJj4IQdvXkOs7R2uLXiSf1pu5HPbEHab1ifalhnZC4LfiVlxDuUc4qV1L3Hx4ot59tdn2ZWxq/pxKaVUJYgx5vSt6rDY2Fizdu1akpOTT5qDr82iA7W1r5JXBwf5+/L8qHaMHdTNrfspVpOkWlx1aeXKlafcl3MfLSyw71foeVXpxke2cGTeuXzcJIQloSGk+558ADC49WBu73U7sS1jTypEUV+U9/veUGjfG17fq9lvt/3j94YRq8cMjo7gl0curJX7IWtrX8UjyshmQcyfGEu/NiFu3U9tjvqL+zQ4OgKanHlyUgU4FE9Lm51p6Zl8n3yAZ4+l0qWg9O1McQfjuO3b20g6nuTyWJVSqkEn1vrKU0ncnfurdJ9iboAHEuGif+DTpANXZuew+OBh3j50hItzcvEpmqG5gGDap+51W7xKqYbL4xcvqbqtOOHVVMmShjXWtA0MfZDD7cYT5XsMiV/I2YlLOPtoCsl+vsxv2pSrsw7D1qXQcZhzta92fUWATwCXtLsEX59TV49SSqmKaGJVXsGlibWYCET2d7wunQk7vyNqwwc8vWMFYOCcO51N86x5zFk7h7T8NDo27chdMXcxvN3wUvcLK6VUZWhiVV6hZK1gt/ALgG6jHa+03bB7JUR0cn786Y5PSct33BO7O3M301dNZ17YPO7uczcXtr2w3l7kpJRyPT0cV17h6quvdtYLdrvmHSH2tlKLLut4GZOb9CTYfqJS1s70ndy/8n6u/fpaViWvor5fQa+Ucg0dsSqv8OCDD3p0/2GNwphmCeDm5IO82zSUD5uEkufjOO7cmraVe368h14Rvbg75m4Gtx6sI1ilVIV0xKq8wpgxYxgzZoxHY4jrP4e7LM/QM6UjK5IPMjHjOI1KjGA3pWzizv/dyRc7P/dglEopb6eJVXmF7du3s337do/tv7gIxa+WaCZZpnNT3gy6pnRmWfIhJmQeJ8DumAYOs9kYvuZNyDzgsViVUt5Np4KVV7jjjjuA8isvuVt5VaQ2m/bcbbmfaOsBplq+YmnGr7wdFkoHi4Vg8YXgFs62yceTybZk0y3c9RWulFJ1j45YVYM3ffHGCh8qsMtE8qDlLq7Jm03PjLO5MbsAznvQcZVxkZfWv8Q1X1/D31b+TWsRK6U0sSo1a3zvCh8qUCzV70zOuGEuTFsPfW5wLt+aupXvk74H4Puk77nyyyt5ZPUjWi5RqQZME6tq8Kr0xJ5mbUuNVhv7N+biiBPPpzUYlu1extjPL+epuKc4mH3Q7fErpbyLJlalqDi5nu7hAu2atOOl8MEsOnSMobknng1rw85nOz/jss9GMeO3GRzNPerW+JVS3kMTq1JFqv3EnoF30uPOP/hv5CjeP3yMc/LynR9ZjY2Pt3/MqE9H8Mn2T05aNW5XCkOe/5G4XSku7YtSynM0sSpVQrWf2NOsLVz+H2Km/MZbZw5n/uFj9M0/kWAL7BY6/P42HN3qXFZ8NfKBjDwmLViryVWpekITq1Jl1Oixe2HtYezrDJgcx7tnXMzcI8foUVDAOXn5nL37V7BbgeKk+jt5tiwA8iw2Ta5K1ROaWJVyh+YdkSvncu6kX1gYPow5x1Kh2xho1cs5UrUG/0HIWS8QEPE/8MnX5KpUPaEFIpRX8HStYLcJj0bGzaPp+dPBx/dEMQprAcERPyK++QS2+B9BzVeTlzqMvPTBTFqw1u0PjldKuY8mVuUVPF0n2O2KHlE3fd6P5FlsiH8GYj9xBbLdt4DAM74lqPkq8tKG8dASQ9zfR3oqWqVUDehUsPIKnq4VXFuKi1EYSwSD9p7HjKOptLFYnJ/b/fIJPGMFRD7LO4nvkGvJ9WC0Sqnq0MSqvMIdd9zhrBdcn5W8ped7+zm8nvYP7kpqw7PHUom0WJ3tsmzZvLjuRUYuuYSPt33swYiVUlWlU8HKK8ycOdPTIdSa4uQ6acFadliimGZ5gC6p+5iW9Rn5oZt5M6wph/wc/zTTCo9z5M83oWk3OLO3hyNXSlWGjliVVxg8eDCDBw/2dBi1pmwxin1+HQi7eSHjbvqeZSH9eTIljVZWK8F2O7fs2QB7fi61vtVuLW+zSikvoIlVeYW4uDji4uI8HUatKrcYRcse+F/7AddM+I5lIf1569BRmjZqDrG3Otfbd3wfly65lHc3v6vnYJXyQjoVrLzCY489BnjmeayeVFyM4iQtexBw7Qf0PLIZMvZBQLDzo3kb53E07yiz187m7YR53NZnCuM7j6exf+NT7ituVwrTF29k1vjeeiuPUm6kI1alvFnLHtDlxG03BbYC1h5Z63yfZjnO7LWzGbn4YuZvmk+OJafczWj5RKVqjyZWpeqQQN9AvrrgvzyZmk4r64nzrGmWLF5e/zKXfnIh8xLeJMd6IsE6i1IUPcxdKzwp5V6aWJWqYwKad+Sam39iWbNzeTI1nTNLJNhMay6vxr/GrWsmMDf+v6zeebhUUi2myVUp99HEqlRdFNGJgHH/xzW3rGJZ+DCeSc0odR9str2AleveZOF7805KqsU0uSrlHppYlarLwqPxv2IuV926hqUthzMjNYO2RZWc7k45zLV8X2YFe6l3eRYb0xdvrKVg3S89PZ1+/fqxa9cu9u7di4iwdq3jnPTp3m/atInIyEhycso/T61UZWliVao+CGuH/5hXGDvpN748czQvH8tgaF4+bzKuRCM7jTu8SmDLLxG/DMDxMPdZ4+tP4YmZM2dywQUXEB0dTVRUFIcOHSImJgbgtO979erFwIEDefHFFz0TvKo3NLEqVZ80jcTvstl0HfMVMmo290yc4CxC4Re6Cd9Ghwho/iuhZ/2bJpGLmXltq3pz601ubi5vvfUW1157LQC+vr60atUKv6IqVqd7D3Drrbcyd+5crFYtwKGqTxOrUvWQvXEEDJhcqsKTX8iJhxwYMZgm63h63c1M//F+tqVt82C0rrF8+XJEhLPPPhs4/dRv2fcAw4cPJy0trcHdT61cSxOr8gozZ85sUPWCa1NxcpWU63j4oB+xefnOz+wYViT/wPil47lrxSTWHVnnwUhrZvXq1fTv3x8RqfY2AgICiImJYdWqVS6MTDU0mliVV2hotYJrmyO5ns07zOTJqDtZkGljSG5eqTZrjvzBxBUTuWXptezO2O2hSKsvKSmJ1q1b13g7rVu3Zu/evTUPSDVYmliVV2iItYJr2+DoCH58dCQdL3uQ/lPjeeOcf/Bxth/Ds3MQY5ztNqdspun3z0DWYQ9GW3V5eXk0atSoxtsJCgoiLy/v9A2VqoAmVuUVHnvsMWe9YFUL/AKg3810v2stcy54ia8Km3FVVjZ+xnBldjbhu34qVZ84Iz/D6wv+R0REkJ6eXuPtpKWl0aJFCxdEpBoqTazKK7z55pu8+eabng6j4fHxhe5jaT95Nc+MeodvpC13ZGTCgCkQGOps9uqGV7l0yXD+u+E10vNrnrzcoW/fvmzZsqXG20lMTKRfv34uiEg1VB5PrCLSXEQ+F5EcEUkSkRsqaBcoIm+IyBERSRORpSISWdvxKvfo0qULXbp08XQYDZcIdBxGq1uWE3HbD3DOXc6PUvJS+OKvL8gozGTuxjcZ/smFzIx7lv1Z+z0Y8MkuvfRStm7dWqNR6969ezlw4ADDhw93YWSqofF4YgVeBwqBlsCNwFwR6VFOu/uAQUBvoDWQDrxaW0Eq91q6dClLly71dBgKILIfBIc73x7IPkCLwGbO9/nGysKdi7nss1H87X93s/GYd1Ru6tWrFwMGDOCrr76q9jYWLlzI8OHDadeunQsjUw2NmBIXLdT6zkWCcSTInsaYHUXL3gcOGGMeKdN2LpBljHm46P1lwIvGmFMOc2JjY83atWtJTk4mKirKLf3wZnWl38OGDQNc+zzWutJ3d3B1363blvH9imm809iPrYEBJ33eN6wrt/S5k2FRw/D18XXZfqtqxYoV3H333ezYsQNf36rFUVBQQKdOnVi4cCFDhgxxU4Tu1VB/56vZ7+rfl3Uanh6xdgasxUm1SAJQ3oh1PjBERFqLSGMco9tvaiFGpRo8v66XMfKuBD7udT/zjsOgMlfNbkjfxv0r7+fRFbeBBw/WR4wYwc0338z+/VWfpk5KSuLxxx+vs0lVeQ+/0zdxqxDgeJllmUBoOW13AsnAAcAGbALuKW+jIjIFmAIQGRlJcnIyqamproq5Tqkr/S4oKAAcR56uUlf67g5u63ub0bRpPYKZe77j8JZ3WGg/wvKQYKxFRRkuSfyWrOypZAzy3BXeY8eOxcfHp8q/S0FBQYwaNcqlv4O1raH+zlen3+4c2Xs6sWYDTcosawJkldP2dSAQCAdygIdxjFjPKdvQGDMPmAeOqeDiH2BDnCKButHvwMBAwPWx1oW+u4tb+97uDiLOn8K/9v3GfXEvsfDo78QFBXJhbh6+Z19PaNG+bXYbz//xPKOjR9OnRR/3xVOGfu8Njzf129OJdQfgJyKdjDE7i5b1ATaX0zYGeNwYkwYgIq8Cz4pIhDFGHyipVG0TgXaDOKPdIO5L3cW03+YifonQfihxu1KYvngjN1yQyaLti1i0fRG9m0ZzY+8pXNL+Evx9/D0dvVJu49FzrMaYHOAzHAkyWESGAGOB98tp/idws4g0FRF/YCpwUJOqUl4gPBq5bDbc+g1xu1OZtGAtBzLymLthgbPJxsxd/H313xmxaBjz4t8gNa9hTluq+s/TFy+BI0EGAUeBhcBdxpjNIjJURLJLtHsIyMdxrvUYMAq4sraDVUpVrDip5llsAOQeuoxhWeBf4oKmo5bjvJrwOpd8ciGP/fQ3Nh3b5KlwlXILjydWY0yaMeYKY0ywMaatMeajouWrjTEhJdqlGmNuNMacYYxpZow51xjzh+ciV0qVFLcrpVRSBbAXRJJ0YDIPJbXhrowswq0nPrNgZ+m+77lh+Q3c8Pnl7EzbUd5mK7XfIc//SNwunbxS3sHjiVUpVfeVl1SLbTTRPJZ3L4tSnuPtppczM7OQ3vkFpdpsz9hNxOLbIOWvau33QEYekxas1eSqvIImVuUVtFZw3TZ98cZyk2pJ+yxNuOmv4YyZupEPB81goTWcy7Oy8TeGkTk5hKXtgyZnOtvvydzDz/t/xmYvf7tlk3mexcakBWtZvz+73PZK1RZNrMoraK3gum3W+N4E+Z+60lGQvy+zxvcGv0DofQ09J63kX1d9zv+aDuKerHzoO6HUE3Xe2/Ied/9wN5d9PIy3ylzsVNEIOc9i4+Fle3TkqjzK07fbKAXgrBM8ZswYD0eiqsPxIPXYCqeDg/x9mT8xlsHREaU/aB1D8yv/D3LToMTINLswm2W7lwFwoDCDVxJe5/WEuVwSeS7dwq7khS8s5Fns5cZSYDVMWrC2/P0pVQt0xKq8wpw5c5gzZ46nw1A1UJxcy45cK0yqJTVuDiEnnoFqtVu5psNomtpOJE8rdr458DMvJj5ASNtnadT8J8S3/GnfPIuN6Yu94+EAquHRxKq8wpIlS1iyZImnw1A1VDa5ViqplqNZo2Y8FHMPP7S/jpnHrfQpc7FTfkA2/i2/pUmnGTSPnA9YS33unHZWygM0sSqvEBERQUSETtvVB8XJNbJZUM2mYxs3J/DCJxkzdSMfnDeHJbTh2uNZhNhPjGLtAr39N/Ge/yx8cCwP9BOdBlYepYlVeYUFCxawYMECT4ehXGRwdAS/PHKha5Kbrz90G0OXW77hiQkr+aHN1fzzuJWYolHsuKxs8gnAjg9B/r68cFkHMuQPVuxdQaGtsOb7V6qK9OIl5RWKk+rEiRM9GofycmHtaXzJM1xx4ZNcsfM71q15na45+7nXdoFz2jnSP4cpf7zCwZyDNPUJZFSHUVzZ7Xq6hXfzdPSqgdDEqpSqe3z9oOso+ncdxdpNifz19SHmX9OXwdERLE34lYM5BwHItBewcNfnLNz1OV0an8nYbjcyKno04UHhHu6Aqs90KlgpVafF9urJz49e4px2jgqOYmrrC4m0lL6gaXvuIV5YN5uLPrmAe5bdzLd7vqXAVlDeJpWqEU2sSql6JTwwnLuGPsvyPtN5yxrGZdk5BJa44MmGYVXKBh76+SFuWzwKjh/yYLSqPtLEqpSqfxo1xefs2zhn0s88f803/Nh6DE9lWemXn1+q2bBDO2DxLaWW5VpyK9ysFvxXlaGJVSlVv7XoQpNLn+fqqYm8e9EbLG8cw12ZOURZLIzOzoHe15RqPunbSVz3xZW8n7iAlLwTCVQL/qvK0ouXlFINg48vnHURUWddxNSCLO7a/CWSuBh6XOVsknQ8icTURAA2r5vD7HVzGBjWnc4RY3hrRQh5Fn/gRMF/vV9WlUcTq1Kq4QkMRfpNgH4TSi3enLIZP/HFahx1i+1AXPoW4tK3ENhBCM1pT1rmQKzZ3cizBGhyVeXSqWCllCoyquMoVvZ/kqey7fTPK30+1uZjyAvdQ1CbhYR1fpqWLRc6R646LaxK0sSqvILWClbeommPcVx91yYWjHiHBw9Ec2dqDl0LSldwsvrYucp3NfP852jBf3USnQpWXkHrBCuv4uMD7c+l87VduXPBb8SmbOCeRqvJD93JDyGB7Anw59KcXFaZdqUK/j+86mH8MFzSYRSDIwcT6Bvo4Y4oT9DEqryCljRU3mhwdARvTBzIpAV+/Jjbj8DcQs4/toHLg1bT2XqYJ2SI8xxrZkEm3yd9j9VYWbp3BcHix7AzB3JJ56sZEjmERn6NPN0dVUs0sSqvoIlVeavSD3EP4Dv7OXyXcw7/9bfw2sRznRcu/XLgF6zmRLWnHGNl2cE1LDu4hiDx5dwWfbmo8zjOizqf0IBQT3VH1QI9x6q8wsqVK1m5cqWnw1CqXOU9Z7ZkUgUY2WEki3rex23Hc2ljsZRaP8/Y+P7oWh5Z8yjnLRzC9OUTwVa65GJZWoyi7tLEqpRSlXC658yKCD36384Dd2xk+cCZfNKoG7cfz6V9Yekka8Uge36GN4eWWl6ybrEWo6jbdCpYeYXZs2cD8NBDD3k4EqUqVvyc2VMKCEZ6XkW3nlfRrTCXabt+ZPfmT/jh8K/8L8CHrYEBXJyTC93OKbXaU3FPseVIPF1DYlm6rg15lkjAR4tR1EGaWJVX+PrrrwFNrKqeCWiMdBtNdLfRRNssTEn6hQOJn9A85UfoNtrZzGK38PP+n8kqzGJPzgH82kALqz9kdyI9qx95OZ00udYhmliVUqo2+PpDx2FEdhwGJZ62A7Ancw8F1tL3yub7WaDZFoKabcHXLjTOa8kdn8fyr5E3cHmPXrUYuKoqPceqlFK1zcfH8SrSOawz7ZPv5Z5DAVx5PJvmNlup5jYfQ1bwYWjxNU/+eQM5f/wfWPJqO+o6LT09nZYtW7Jr1y727t2LiBgRia3MuiISKCL7KtteR6xKKeUFHrz2MiYtaElgRgbDfDbQJfg3soOT+a2xPzsCA5ztzs7PJ/j7f0Dfm5zLtqVtY3PKZjr6dCSKKE+E7/VmzpzJqFGjiI6OxuY4cDkTqNRVYcaYAhGZBfwbuOh07TWxKqWUFyh5v+wXlvMg6zz8s6yc7bON2wJ+p1HIVjYFWzgvNw86DgP/EwUnvtr1Fe9veR+AzutacG7UBZzbcQQxLWLw9/X3UI+8R25uLm+99RZLly4FwNfXF2PM4Spu5kNgtoj0MMZsPlVDTaxKKeUlShejsGHBjzh7TzbY+jD/yljuapoBO7+F8E6l1ltzYI3z7zvyj7Fj5ye8vfMTGosvA8K6MajDCIa0HUbb0LaISC33yvOWL1+OiDBkyBAA9u7dS4cOHQxwNrAeSAJeMMa8WryOiHQGtgP9jTHrjTFpIvILcD3wxKn2p+dYlVLKi5RXjMJ5NXDEWTDobug83NnebuyMO+sqzik0+BlTalu5xsbKtESeWzeb0Z+PZuRHQ1i7dQmUaVffrV69mv79+5d7UGGMsQMLgRvLfHQjsNUYs77Esj+A80+3P02sSinlZU5XjKIkH/Hhlh638NbYxSxrcRX/4Qyuyco5qfoTwAFrFq0+nQKJnzqXGWN4Z93/GPz8t/W2EEVSUhKtW7c+VZMPgHNEJLrEshuKlpd0EGh/uv2ddipYRPYC63AMl9cD64wxR0+3nlJKqeqrVDGKYiJwRjdsvSdzQdSzXFCQDXvXsG/HUuIOxPGL/Th/BDUiwmajjdUGHc5zrrpk05+8mPgAvi2Fv393Bhd1GcH4PiPpGtYVXx9fN/WuduXl5dGyZcsKPzfGbBSRTThGqc+KyDlANI7zqqU2BQSdbn+VOcdaAFwBXAkYABE5xIlkuw5Yb4w5WIltKVUurROslAsFhkCXEbTtMoK2wHXpSVh2/Y/Df30H7c+CkDMAR+nEp777At8Wjlt6Mhof4dPkd/k0+V1CxZ/YsC4MaH8JA9qcy1nNzsJH6uYkZ0REBOnp6adr9gEwCXgWR4JdY4xJKtOmOXDsdBs6bWI1xnQRkVCgL9AP6F/0ugwYw4lke5QTI9p/nG67SimlaklYO/xjJxEVO8m5qLgesTXYl7DCRhwPyC+1Spax8FNaIj+lJcL6l2juE8g1bS7k7nOfLXVFcl3Qt29f5xO0TuEj4DkRGQhcCzxZTpueOPLcKVXq8MMYk2WM+dkY87Ix5iZjTHegCXAe8ADwPo77gUYAj1dmm0qVNHv2bGe9YKWUexUn1TyLDUvGQJrvvpmb9nRjylFhZHYOEVbbSeuk2QuwJyyCT24utXzD0Q3sTN+J3dhPWsdbXHrppWzdupXU1NQK2xhj9gOrgDeApsDicpoNBVacbn/Vvt3GGJMLrBGR40A4jsuWhaIRrFJV8euvv3o6BKUajOmLN5JnOZE8N5mObMrvCPnQNDWbAT5bGB0YT0DobvY0yuXPRoFk+PoyID8f2pZ+eMDzfzzPltQtNPMJoH+zLvRvdyH9IgfRJawLfj7ecUdnr169GDBgAIsWLeLuu+8+VdMPgPnA58aYUnPHIjIIR8Jdcrr9VavXInI2MK7o1bFo8e/AdODTitZTqiKffqq/NkrVllnjeztHrGVlEsL39gGssQ5i/thY7m1psO9Zxc5dK2h/fAO0O9fZNqswi21p2wDIsBfyQ9omfkjbBBteIVj86BPann5thtIvaii9InrRyK8RcbtSmL54I7PG967VBwo89dRT3Hfffdx55520b98eY8xJ994YY94G3q5gE38DZhljTltLstKJVUTOxZFIrwLaAHZgNfAyjuyuFy8ppVQdULYQRVml7p0FfHpdTZdeV5/ULrswm4tan8vafT+R7lv6CuIcYyXu+F/EbfkLtryDH8IbLW7mvl87cswSWOtP6xkxYgR33303+/fvp127dlVaV0QCgY3AS5VpX5nbbeYCY4GWgBX4Afgn8IUxpn7e9KRq3aOPPgrAc8895+FIlGoYKkquZZPqqZwZciYvnjcLs/M7du36lnVH1rKuMJX1jQI54lc6vYRZLQz4459geR0IdDxn9sPlDI1NYHiHWPqeOYAOTTu49crjadOmVWs9Y0wBjrxXKZUZsd6BI6G+BzxjjNlbrcgqICLNccxpD8dxAdSjxpiPKmjbD8cIuR+QA8w0xrziyniUZ+g5VqVqX9nkWpWk6hTQGOlxBWf1uIKzgGvzMjDJf3Bw9/9Yf+g31uUeYH2AP10KC0m2t+AYYc5VrQE7+fXYd/x67DsAQsWf3qHt6H3mOfSJGkqvFr1oEtDExb12v8pOBfsBNwM3i0gSJ4pFrMdxD2tNCka8DhTiGBHHAMtEJKFskWMRicBxNdYDOE4eB+CYklZKKVVNxcnVZec9g5ohnYcT2Xk4kUDEjgPsfH8x4fbNvEvpBwL4Nt5b6n2WsfDL8b/45fhfsN1RmyE6IIze4d0Z3uUazm1XyYIZHlaZxNoUxwix+B7WfjgKRlxF6YIRJRPtV5XZuYgE4zhv29MYk43jKuOvgJuAR8o0/xvwrTGmuBJGAbC1MvtRSilVsSpVeaqihz7bzgFLJ36h00mfFaYN5aLCA/g1TiK+USBpvidXetpVmM6uQ7/QYusKzu18HYx83vnZtrRtNAtshvGy2seVKRCRhePenlXFy4oSYgwnEm1/YCQwGkeyrWwdrM6A1Rizo8SyBMovcjwQ2CQiccBZOK5CvtsYs6+S+1JKKVXLTnUFsj0/kmX5D3Kmfw7zL7QTWriehEN/sDFnPwn+PmwPCMBWVDi/V34+NG5eav2n4p5iS+oWbmk2kIfa/l+t9KcyqnW7jTEmB/il6AWAiAThSLZ9q7CpEOB4mWWZQGg5bdvgSOKXAJuAF3A8kWBI2YYiMgWYAhAZGUlycvIpbwyuz+pKvwsKCgBITk522TbrSt/dQfveMHlj36MC4PlR7Xh42R4KrCePLAP9hEdG9SK0TQgwhD7doY/dxm3pO7EfWc+eY3+wJWsnfQoKOBrQloKi/yMKbAVsT9sOQGRB+f93rN+fzcwfk3nswij6tQkpHVeU+x4IX5mrgpOBL4peK40xJx92AEX39vxa9KqsbBwVnEpqAmSV0zYPx209fxbF9QyQIiJNjTGZZWKZB8wDiI2NNcU/QHf+IL1ZXeh3YGAg4PpY60Lf3UX73jB5Y9+joqDFGS2qdgVyu/bAJbSmaPSUlw7+jcHP8X/F4ZzD9G8RQ+LhP+nZ7pKT+h23K4VHlm8mz2LjkeVJtXprT2Wua/4SxznV74GjIvK+iFwpIo1dsP8dgJ+IlJx87wOU93T2jZSu6uRdk+pKKaUqdMrnzFZGUJgzqQK0Cm7F/IvnEhf7DM1axZZqWrJkI+C4tWfB2lp7LN5pE6sx5h5jTBSOc5zzgFgc1ZWOiciXIjJRRMKrs/OiKeXPcDymJ1hEhuC4Z/b9cpq/A1wpIjEi4o+jQPKasqNVpZRS3qkqz5mtFP8gfHtdDT4nrjYum1SL1WZyrfSduMaYP4wxjxpjugHdgRlAKxz3oB4WkZ9EZJqItK1iDFNxPN/uKI5zpncZYzaLyFARyS6x/x+Bx4BlRW3PwvEgWlUPhIeHEx5ereMzpVQdUnwFsjumZStKqsVqK7lWq8SFMWabMeY5Y8w5QFsc95bagNnAHhFZLyIjKrmtNGPMFcaYYGNM2+LiEMaY1caYkDJt5xpjIo0xYcaYMcYY113pojzq008/1XrBSqkaKftwgfLkWWxMX7zRrXHUuHaUMeaAMeY1Y8zFOIo83ArsxfHcOqWUUqpWzBrf23kOtyJB/r7MGt/brXG4tCijMSbdGPOeMeYqY4w+XFNV2qOPPuqsF6yUUtVR9gKpsqpVsrEaqpVYRcRHRIaJyEUiUtliEEpVKDU11SvvwVNK1S0VJdfaSqpQtcfGNQZG4Lj1ZhRQXAIjQ0SW47jPdUVRaUKlqmTevHmeDkEpVU+45OECNXDKEauItBSR20XkaxxPnlmC47abd4DzcNy3+xaOW3A+wXELznIRmSIirdwbulJKKVU+l9/aUwWnG7EWP7x8LY5n0X1pjNlSps2vwMMi0hXHaPZyYC7wXxEJMsZYXBivqqemTJkC6MhVKeU67ny4wKmcLrFOBb4yxhw63YaMMduA54Hni0arY3DcgqPUae3YseP0jZRSqg44ZWI1xrxZnY0aYw4D3vOoAaWUUqqWVPnpNkXlC3sDVmB7DR9yrpRSStUrVUqsIjIW+ABoXGLZQUo85BzHg84PuDJIpZRSqq6o6oj130A+8DiO56h2w/GM1CE4zqlC1R50rpRSStUrVU2sUcCTxpj/lP1ARNoB/XE87FwppZRqkKqaWLcCUt4HxpgkIAnHY+CUUkqpBqmqJQ1fASZoGUOllFKqfFUasRpj3heRAcBnInKbMUaLuyqX6Ny5s6dDUEopl6jqVcGROJ6/OgY4ICIrgd9wXA28Tq8GVtWlFZeUUvVFVc+xvgNcDMQD6UBfYDiOK4ERkWM4brcZ5cIYlVJKqTqjqol1CPB/xpg7iheISBsct9z0L3r1c114qqHQWsFKqfqiqok1HVhXcoExZj+wH/jKVUGphic8PNzTISillEtUNbF+BgwFdFihXOq5557zdAhKKeUSVb3d5lWgj4hc645glFJKqbquqol1O9AR+EhEPhOR60SkvevDUg3NuHHjGDdunKfDUEqpGqvqVPAsHCUL++J4qPkVgBGRDE4U4V9njPnEZRGqBiE1VW+JVkrVD1UtEPH34r8X3dPaF8dVwH2LXhfhuPVGE6tSSqkGqcrPYy1WVAziAPB18TIRCcORYJVSSqkGqdqJtTzGmHTgR1duUymllKpLTnnxkohMLZryrRIROVNE7hCRql4cpZRSStVpp0t8rwD7RGStiDwhIr0raigiPUXkcRH5HUfBiNfQB54rpZRqYE43FdwSR8H9scAjwDMisg/4EkelJVvRZ2OB9kA+8D1wO7DUGGNxT9hKKaWUdzplYjXGpAHvAu+KSCPgEhy32FwHTMNxBXAqjguY/gZ8Z4zJc2fASimllDer9MVLxph8YCmwVEQEGIRjKjnOGGN3U3yqgRg0aJCnQ1BKKZeo1lXBxhgDxLk4FtWAaa1gpVR9oVftKqWUUi6kiVV5Ba0VrJSqL1xaIEKp6tJzrEqp+kITq/IKDz30kKdDUEopl9CpYKWUUsqFNLEqrzBs2DCGDRvm6TCUUqrGNLEqpZRSLuTxxCoizUXkcxHJEZEkEbnhNO0DRGSriOyvrRiVUkqpyvKGi5deBwpx1CWOAZaJSIIxZnMF7acDx4DQ2glPKaWUqjyPjlhFJBgYBzxpjMk2xqzBUdz/pgradwAmAFqmRymllFfy9FRwZ8BqjNlRYlkC0KOC9q8CjwFa6F8ppZRX8vRUcAhwvMyyTMqZ5hWRKwFfY8znIjLsVBsVkSnAFIDIyEiSk5NJTU11ScB1TV3pd0FBAQDJycku22Zd6bs7aN8bpoba9+r0Oyoqyg2ROHg6sWYDTcosawJklVxQNGX8AjCqMhs1xswD5gHExsaa4h+gO3+Q3qwu9DswMBBwfax1oe/uon1vmBpq372p355OrDsAPxHpZIzZWbSsD1D2wqVOOB6kvtrxxDoCgKYichgYaIzZWzvhKqWUUqfm0cRqjMkRkc+AZ0XkdhxXBY8FBpdpmgiUPBwZDLwG9MNxhbCq40aPHu3pEJRSyiU8PWIFmAq8DRwFUoG7jDGbRWQo8I0xJsQYYwUOF68gImmA3RhzuNwtqjpHawUrpeoLjydWY0wacEU5y1fjuLipvHVWAm3cGphSSilVDZ6+3UYpQGsFK6XqD4+PWJUCmDhxoqdDUEopl9DEqryCJlalVH2hU8HKK6SkpJCSkuLpMJRSqsZ0xKq8wtVXXw3AypUrPRuIUkrVkI5YlVJKKRfSxKqUUkq5kCZWpZRSyoU0sSqllFIupIlVKaWUciFNrEoppZQLaWJVSimlXEjvY1VewR33r3rTg4+VUg2HjliVUkopF9LEqpRSSrmQJlallFLKhfQcq1INiMViYf/+/eTn53s6FLexWq1kZ2d7Ooxq8/X1pVmzZkRERODjo2OfukgTq1INyP79+wkNDaV9+/aIiKfDcYvCwkICAgI8HUa1GGOwWCwcOXKE/fv307ZtW0+HpKpBD4eUakDy8/MJDw+vt0m1rhMRAgICiIyMJCcnx9PhqGrSxKqqbNGiRQQGBmKxWKq03h9//IGInPLl4+NDVlaWmyJXgCbVOkCngOs2nQpWVZaQkED37t3x9/ev0npnnXUWv/76q/P9f/7zH77++mu+++475zJ/f39CQ0NdFqtSStU2PSxSVZaQkEDfvn2rvF7z5s0ZOHCg85WZmUnPnj1LLevfv78bIlauELcrhSHP/0jcrpRa3/eHH37I8OHD3bLtiRMn8sQTT1R7/ZCQEHbv3u3CiFRdp4lVVVl8fDwxMTHO9z/++CPh4eHcf//92Gy2Sm9n06ZN9OzZ0w0RKleL25XCpAVrOZCRx6QFa92WXNesWcPgwYNp2rQpzZs3Z8iQIfz555/ceOONpWY2PGXYsGG89dZbpZZlZ2fTsWNHD0WkvJEmVlUlx44d49ChQ87E+uqrrzJ69Gief/55Xn75ZXx9fSu1nczMTJKTk+nVq5cbo1WuUJxU8yyOg6Y8i80tyfX48eOMHj2ae++9l7S0NA4cOMBTTz1FYGCgS/ejlLtpYlVVkpCQgIjQvXt3Jk+ezNNPP83y5cuZPHlylbaTmJgIoCNWL1c2qRZzR3LdsWMHANdffz2+vr4EBQUxfPhwevfuzYIFCzj33HOdbUWE//73v3Tq1InQ0FCefPJJdu3axeDBg4mIiOCaa66hsLAQ4KR1i9f/66+/ToohPT2d0aNH06JFC8LCwhg9ejT79+8H4PHHH2f16tXcc889hISEcM8995y0rczMTG6++WZatGhBu3btmDFjBna7vVQcDz30EGFhYXTo0IFvvvnGZT8/5T00saoqiY+Pp3nz5lx11VWsWbOG33//nWHDhjk//+KLLwgPDycmJobu3bvTp08ftmzZctJ2ihNrZUasixcv5tZbb3VZH1TlVJRUi7k6uXbu3BlfX19uueUWvvnmG9LT00/Z/ttvv2XdunX89ttvvPDCC0yZMoUPPviAXbt2kZiYyMKFC6scg91u59ZbbyUpKYl9+/YRFBTkTKD/+te/GDp0KK+99hrZ2dm89tprJ61/7733kpmZye7du1m1ahXvvfce77zzjvPz33//nS5dupCSksLDDz/MpEmTMMZUOU7l3TSxqipJSEjAGMPq1at57bXXOOuss0p9Hh8fz7333kt8fDxbtmzh7LPPLvc/uE2bNtGyZUsiIiJOu89169bRr1+/SsdYlfO8qmLTF2+sMKkWy7PYmL54o0v216RJE9asWYOIMHnyZFq0aMHll1/OkSNHym3/8MMP06RJE3r06EHPnj0ZPnw4HTt2pGnTpowcOZINGzZUOYbw8HDGjRtH48aNCQ0N5fHHH2fVqlWVWtdms7Fo0SKee+45ZxGOBx98kPfff9/Zpl27dkyePNl5AHHo0KEK+6fqLk2sqkoSEhKYNm0aEyZM4NZbb+XYsWOlPo+Pj3de2ZuRkUFiYiJDhw49aTuJiYkVjlb379/PZZddRq9evZg8eTJ//vmnM7EePnyYG264gYEDB9K9e3c+/fRTAO68807uvfdezjvvPMaPH893333HoEGD6Nu3Lz169GD16tXO7d9+++089NBDXHzxxbRt25YZM2Y4P9u3bx8jR46kd+/eTJ48mXPPPZfNmzfX7IdWR80a35sg/1OfMw/y92XW+N4u22e3bt1YsGAB+/fvJzExkYMHD3L//feX27Zly5Yn4ggKOul9dcoa5ubmcscdd9CuXTuaNGnCeeedR0ZGRqUO1lJSUrBYLLRr1865rF27dhw4cMD5vlWrVs6/N27cGKBOl19U5dPEqiqtoKCAbdu20bt3b+bNm0dERATjx4/HarU628THx/PII4/Qu3dvWrVqxbhx48q9TSIxMbHc86s2m43LL7+ce++9l02bNnH++eezatUqYmJisNlsTJgwgenTp/Pbb7+xcuVK7rvvPud+c3Jy+Omnn/jss8/o378/cXFxbNiwgVdeeYVZs2Y595GQkIDVauW7777jzz//5PXXXwcc5eSuv/567rvvPjZu3Mh5553Hhg0b6Nq1q6t/lHXC4OgI5k+MrTC5Bvn7Mn9iLIOjTz/rUB1du3Zl4sSJztMG1RUcHExubq7z/eHDhytsO2fOHLZv387vv//O8ePH+fnnnwGc07WnKq4RERGBv78/SUlJzmX79u0jMjKyRvGrukcTq6q0LVu2YLFY6NWrF0FBQXz++eds3ryZBx54AHCMULOysti8eTMbN25k06ZNzJgxg+PHj5fazqFDh0hNTS13xPrNN9/QunVrRowYATjOwXbq1Ing4GCWL19OQkICt956KzExMQwfPpzGjRtjt9vZtm0bc+bMcV6VvGTJEoYMGUKfPn2YMmUKQUFBgCNx//XXX8yYMQMfHx9sNhvh4eEALF++nLCwMOe+e/ToQe/evSt9pXN9VFFydUdSLf4Oiy8WSk5OZuHChQwcOLBG2+3Tpw+bN28mPj6e/Px8nn766QrbZmVlERQURLNmzUhLS+OZZ54p9XnLli0rvGfV19eXa665hscff5ysrCySkpJ48cUXmTBhQo3iV3WPJlZVaQkJCQQHBxMdHQ04prk+/vhj5s6dy4IFC4iPj6d79+7O9p06dUJEyMzMLLWdTZs2AeVfERwfH1/qfGrJaeCNGzfy8MMPEx8f73zt2LGDHTt20LlzZ8LCwgBHUv3iiy9YtmwZCQkJXHXVVfTp0weA7du307lzZ0JCQpz7K/5sw4YNpQpfbNy4sUrnduurssnVXSPV0NBQfv/9d8455xyCg4MZOHAgPXv2ZM6cOTXabufOnfnHP/7BxRdfTKdOnU66Qrik+++/n7y8PCIiIhg4cKDzIKvYfffdx5IlSwgLC2PatGknrf/qq68SHBxMx44dOffcc7nhhhu47bbbahS/qoOMMfX61b9/f2OMMfv27TMNUW32+6WXXjJTpkxxvv/4449Njx49qrSNuXPnmrFjxxpjjDlw4IBp3769mT17tjHGmLfffttceOGFprCw0BhjzMGDB83hw4fNwoULS+33iSeeME899ZQxxph169aZkJAQs3z5cmOMMR999JG54447nG1nzJhhZs2aZYwx5o033jBXX321McaYQ4cOmbPOOsu89dZbVYrfW1T0vW/ZsqXa2/zlr2Nm8HM/mF/+OlbtbdSGgoICT4fgEtX5rvT/uSpxW97RWsHKZeLj41mxYgW///47/v7+REZG8tVXX1VpGzfeeCOLFi2ie/fuREZG0qRJE+eo8cYbb2TlypV069aNkJAQwsPD+eijj4iPjy810rz55psZO3YsS5cuZciQIYSHhzs/T0hIKFU1asOGDdx5550A3HDDDSxcuJBu3boRExNDkyZNOP/882v4U6k/BkdH8MsjF3o6DKW8nph6fg9VbGysWbt2LcnJyURFRXk6nFrXUPtdHdnZ2c4p4mXLlvHKK694RRm96qjoe9+6dSvdunXzQES1py4/j7Wk6nxXDfXfezX77bbHPOmIVakir732Gu+99x6BgYFER0fz7rvvejokpVQdpIlVqSKPPPIIjzzyiKfDUErVcXpVsFJKKeVCmliVUkopF9LEqpRSSrmQJlallFLKhTyeWEWkuYh8LiI5IpIkIjdU0G66iCSKSJaI7BGR6bUdq1JKKXU6Hk+swOtAIdASuBGYKyI9ymknwM1AGDACuEdErqu1KJVSbtW+fXuCgoIICQlxvoqfhepqK1eupE2bNm7ZtlIevd1GRIKBcUBPY0w2sEZEvgJuAkrd92CMeaHE2+0i8iUwBFhUW/Eqpdxr6dKlXHzxxZ4OQ6ka8fR9rJ0BqzFmR4llCcAp68iJ49lNQ4E3K/h8CjAFIDIykuTkZFJTU10TcR1TV/p9zTXXAPDJJ5+4bJupqanOJ9c0NBV971arlcLCwlqOpvIsFstJ8d1zzz0cO3aMjz/+GIDHHnuMdevWsWLFCjIyMrj11lv5888/sVqtDBo0iFdeecX5TNS0tDT+/ve/8/3335OXl8fQoUN59913GTlyJAUFBc5KW4mJibRu3bp2O3saVquV5OTkKq1TV/69u1p1+u3OClWeTqwhwPEyyzKB0NOs9zSOaex3yvvQGDMPmAeOkobFP8CGWOoL6ka/AwMDAdfHWhf67i7l9T07O7t0ub+fnoNVz1dug/1ugcv/U3rZV9Ng/SkqVJ3/CFzwaOW2D/j7+59UjvDll18mJiaGjz76iOjoaOeTlAIDA/Hz82PSpEksWbIEm83GbbfdxoMPPuisUT1p0iRCQkLYvHkzISEhxMXFERYWxjfffMOECROcj6jzRn5+ftX6/W2ov/Pe1G9PJ9ZsoEmZZU2ArIpWEJF7cJxrHWqMKXBjbKoWzZw509MhKC9wxRVX4Od34r+lWbNmMXnyZN5//31GjhxJaGgor776qvP8aHh4OOPGjXO2f/zxx7ngggsAx3N/v/nmG1JTU52PFNSHKqja4OnEugPwE5FOxpidRcv6AJvLaywit+E493qeMcZ7DzVVlQ0ePNjTISgv8MUXX5R7jvWcc86hY8eOHD161HnaACA3N5cHHniAFStWkJ6eDjgeVm6z2UhOTqZ58+bOpKpUbfFoYjXG5IjIZ8CzInI7EAOMBU76X1ZEbgRmAhcYY3bXaqDK7eLi4gBNsLXugkerNFV7ksv/c/L0sBu8/vrrFBQU0Lp1a1544QUefdQR85w5c9i+fTu///47rVq1cj5C0BhDVFQUaWlpZGRk0KxZs1Lbc1ymoZR7eMPtNlOBIOAosBC4yxizWUSGikh2iXYzgHDgTxHJLnq94YF4lRs89thjPPbYY54OQ3mhHTt28MQTT/DBBx/w/vvv88ILLxAfHw84RqdBQUE0a9aMtLQ0nnnmGed6Z555JiNHjmTq1Kmkp6djsVj4+eefAWjZsiWpqalkZmZ6okuqnvN4YjXGpBljrjDGBBtj2hpjPipavtoYE1KiXQdjjL8xJqTE607PRa5c6c033+TNN8u9yFs1IGPGjCl1H+uVV17JhAkT+Pvf/06fPn3o1KkTM2fO5KabbqKgoID777+fvLw8IiIiGDhwICNGjCi1vffffx9/f3+6du3KGWecwcsvvwxA165duf766+nYsSPNmjXj4MGDHuitqq/0Qef1XEPtN2jf9UHndZs+6LzyvO1B5x4fsSoFjsIAS5cu9XQYSilVY56+KlgpwHERCjimApVSqi7TEatSSinlQppYlVJKKRfSxKqUUkq5kCZWpZRSyoU0sSqllFIupIlVKaWUciFNrEoppZQLaWJVSimlXEgTq/IKWitYAbRv357//e9/ng6jlPT0dESEQYMGlVp+55138sADD3goKuXNNLEqr9ClSxe6dOni6TBUPWC1Wl26vfj4eFq1asWWLVs4fPiwc/mGDRuIiYlx6b5U/aCJVXkFrRWsTuX5558nOjqa0NBQunfvzueff17q8/bt2/Pvf/+b3r17ExYWxh9//EHfvn0JDQ1l/PjxXHvttTzxxBPO9gcPHmTcuHG0aNGCDh068J//VPxM2fj4eGJjY7nkkkv48ssvAbDZbGzatIm+ffu6p8OqTtNawcoraK1gz/hv/H+ZmzC3Um3HdRrH04OfLrXs6bin+XTnpxWuc1efu5gaM7UmIQIQHR3N6tWradWqFYsXL2bChAn89ddfnHnmmc42CxcuZNmyZTRq1IjY2Fj+9re/MXXqVJYuXcp1113Hww8/DIDdbmfMmDGMHTuWhQsXsn//fi6++GK6dOnCpZdeetK+i0emXbp04cMPP+SOO+5g27Zt2O32ev+kIFU9OmJVXmHJkiUsWbLE02EoLzV+/Hhat26Nj48P1157LZ06deKPP/4o1WbatGlERUWxceNGrFYr06ZNw9/fn6uuuooBAwY42/35558cO3aMf/zjHwQEBNCxY0cmT57MokWLyt13fHw8MTExXHbZZaxevZqsrCzi4+Pp0aMH/v7+bu23qpt0xKq8QkREhKdDUF7svffe48UXX2Tv3r0AZGdnk5KSUqpN8fM4Dx06RGRkJCJy0mcASUlJHDx4kGbNmjmX2Ww2hg4detJ+CwoK2Lp1KzExMYSFhTFgwAC++eYbPb+qTkkTq/IKCxYsAGDixIkejaOhmRoztUZTtU8Pfvqk6WFXS0pKYvLkyfzwww8MGjQIX19fYmJiMMaUalecSFu1asWBAwcwxjiXJScnEx0dDTiSbIcOHdi5c+dp952YmEjjxo3p2LEjAFdccQVffPEFR44c4corr3RlN1U9olPByissWLDAmVxVw2axWMjPz3e+MjMzERFatGgBwDvvvENiYmKF6w8cOBBfX19ee+01rFYrX375Zalp4wEDBhAaGsq///1v8vLysNlsJCYm8ueff560rQ0bNtC7d29ngr788stZvny5jljVKWliVUp5lVGjRhEUFOR8ffLJJzz44IMMGjSIli1bsmnTJoYMGVLh+gEBAXz22WfMnz+fZs2a8cEHHzB69GgCAwMB8PX15euvvyY+Pp4OHToQERHB7bffTmZm5knbKj6/Wqx9+/a0b9+ejIwM+vTp4/K+q/pBp4KVUl6j+Bxqef71r39Ver3Y2Fji4+Od788555xSV5y3bt2ahQsXnjae11577aRlJberVHl0xKqUqndWrVrF4cOHsVqtvPvuu2zcuJERI0Z4OizVQOiIVSlV72zfvp1rrrmGnJwcOnbsyJIlS0rd86qUO2liVUrVO1OmTGHKlCmeDkM1UDoVrJRSSrmQJlallFLKhTSxKqWUUi6k51iVV9A6wbWnZEUi5Z3sdrunQ1A1oIlVeQWtFVw7GjVqRGpqKuHh4ZpcvZAxBovFwpEjRwgODvZ0OKqaNLEqr6C1gmtHmzZt2L9/P8eOHfN0KG5jtVrx86u7/7X5+fnRtGlTPdisw+rub5+qVzSx1g5/f386dOjg6TDcKjk5udTTbJSqbZpYlVdYuXKlp0NQSimX0KuClVJKKRfSxKq8wuzZs5k9e7anw1BKqRrTxKq8wtdff83XX3/t6TCUUqrGNLEqpZRSLqSJVSmllHIhTaxKKaWUC2liVUoppVzI44lVRJqLyOcikiMiSSJyQwXtRET+LSKpRa9/i9ZkU0op5WW8oUDE60Ah0BKIAZaJSIIxZnOZdlOAK4A+gAG+B/YAb9RapEoppdRpeHTEKiLBwDjgSWNMtjFmDfAVcFM5zW8B5hhj9htjDgBzgIm1FqxSSilVCZ6eCu4MWI0xO0osSwB6lNO2R9Fnp2unlFJKeYynp4JDgONllmUCoRW0zSzTLkRExBhjSjYUkSk4po4BskVkOxABpLgk6rqlTvXbxafN61TfXUz73jA11L5Xp98rjDEj3BGMpxNrNtCkzLImQFYl2jYBsssmVQBjzDxgXsllIrLWGBNbs3Drnobab9C+a98bnobad2/rt6engncAfiLSqcSyPkDZC5coWtanEu2UUkopj/FoYjXG5ACfAc+KSLCIDAHGAu+X0/w94G8iEikirYEHgQW1FqxSSilVCZ4esQJMBYKAo8BC4C5jzGYRGSoi2SXavQksBTYBicCyomWVNe/0Teqlhtpv0L43VNr3hser+i3lnKJUSimlVDV5w4hVKaWUqjc0sSqllFIuVG8Sa0OuOVyFvj8tIhYRyS7x6ljb8bqKiNwjImtFpEBEFpym7QMiclhEjovI2yISWEthukVl+y4iE0XEVuY7H1ZrgbqYiASKyPyi3/MsEYkXkZGnaF9vvveq9L2+fe8AIvKBiBwq+i53iMjtp2jr0e+93iRWStccvhGYKyLlVWYqWXO4NzAGuKOWYnSXyvYd4GNjTEiJ1+5ai9L1DgIzgLdP1UhELgUeAS4C2gEdgWfcHp17VarvRX4t852vdG9obuUHJAPnA02BJ4BPRKR92Yb18HuvdN+L1KfvHeA5oL0xpglwOTBDRPqXbeQN33u9SKwNueZwFfterxhjPjPGfAGknqbpLcB8Y8xmY0w68E/q8HcOVep7vWKMyTHGPG2M2WuMsRtjvsbxMI6T/oOlnn3vVex7vVP0PRYUvy16RZfT1OPfe71IrDTsmsNV6TvAGBFJE5HNInKX+8PzCuV95y1FJNxD8dS2viKSUjR99qSIeLrimsuISEsc/wbKKxZTr7/30/Qd6uH3LiL/FZFcYBtwCFheTjOPf+/1JbG6pOawm2Jzt6r0/ROgG9ACmAz8Q0Sud294XqG87xzK/xnVNz8DPYEzcMxsXA9M92hELiIi/sCHwLvGmG3lNKm333sl+l4vv3djzFQc399QHMWFCspp5vHvvb4kVrfUHK4jKt13Y8wWY8xBY4zNGBMHvAJcXQsxelp53zmU//tRrxhjdhtj9hRNHW4CnqUefOci4oOjQlshcE8Fzerl916ZvtfX7x2g6P+vNUAboLxZN49/7/UlsTbkmsNV6XtZBqirI/WqKO87P2KMaVDnJ4vU+e+8aHZpPo6L9cYZYywVNK1333sV+l5Wnf/ey+FH+edYPf6914vE2pBrDlel7yIyVkTCxGEAMA34snYjdh0R8RORRoAv4CsijSo4j/QeMElEuotIMxxXUy6ovUhdr7J9F5GRRefiEJGuwJPU4e+8yFwcpzTGGGPyTtGu3n3vVLLv9e17F5EzROQ6EQkREd+iK3+vB34op7nnv3djTL14Ac2BL4AcYB9wQ9HyoTimeovbCfACkFb0eoGi0o519VWFvi/EcRVpNo6T/9M8HXsN+/00J64OLH49DbQt6mPbEm3/BhzBcT76HSDQ0/HXRt+B2UX9zgF245gS9Pd0/DXod7uivuYX9bP4dWN9/96r0vd6+L23AFYBGUXf5SZgctFnXve9a61gpZRSyoXqxVSwUkop5S00sSqllFIupIlVKaWUciFNrEoppZQLaWJVSimlXEgTq1JKKeVCmliVUkopF9LEqpRSSrmQJlal6igRWSsimzwdh1KqNE2sStVBRXWBewLxNdhGiIjYRMRU8tXcZR1Qqh6r8w++VaqB6g4EUoPEiuPf/y1llt0FDAYewlFrtViBMSatBvtSqsHQWsFK1UEicjPwLnCRMeZHF253HY6kHWqMsbpqu0o1JDoVrFTd1Lfoz/jiBSLSTEQ+F5F8EZlS1Q2KiD+O6eWNmlSVqj6dClaqbooBkounZ0WkP7AYx2MRhxhj1lVjmz2AAGCDq4JUqiHSEatSdVMMRaNVEbkL+AXYCvSvZlKFE6Pg9TUNTqmGTBOrUnWMiLQHmgE7ReQj4DXgX8DoGl5g1K/oTx2xKlUDOhWsVN1TPLK8F7ADI4wx35dtJCI+QCaOA2hT9OcO4G5jzC8VbNcK6L2xStWAjliVqnuKE+v/4bjlplcF7boAIUCUMSYECAO2AS+VbViUhPsAW40x+S6PWKkGREesStU9MUCKMeZuEQkCZonIXmPMZ2Xa9QcOFk8PG2MKRGQD0LmcbXbCkYT1/KpSNaQjVqXqnr6cOA96B7AS+EBEBpRp1x/YCCAOvYHJwOxytqnnV5VyEU2sStUhIhIOtKEoARpjLMA4YC+wVEQ6lGjeH7hARDKAHBxXES80xnxUzqb1imClXEQTq1J1S3ECdI4sjTEZwGVFb5eLSFjROdO+wA3GmGbGmMbAQOBhETm/gu0aIMFtkSvVQGhJQ6XqIRHpiuO+1ihjzP4Sy1OA6caYdzwWnFL1nF68pFT9VHzh0n4AEQkG/oajstJ3ngxMqfpOE6tS9VN/oJWIZOO41/U4EIej3OEBj0amVD2nU8FKKaWUC+nFS0oppZQLaWJVSimlXEgTq1JKKeVCmliVUkopF9LEqpRSSrmQJlallFLKhTSxKqWUUi6kiVUppZRyof8HcQfAd51/JRQAAAAASUVORK5CYII=\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",
    "# plot of simulation results\n",
    "plt.plot(temp_vals,sim_m_vals/np.sum(Nelems), \n",
    "        label = f'Simulation',\n",
    "        markersize = 7.5,\n",
    "       marker = 'D',\n",
    "       linestyle = '')\n",
    "\n",
    "# plot of exact analytical results \n",
    "m_exact_vals = [avg_m_exact(T, Dels, Nelems) for T in Tvals]\n",
    "plt.plot(Tvals, m_exact_vals/np.sum(Nelems),linestyle = '--', linewidth = 3.0, label = f'Exact')     \n",
    "\n",
    "# plot of large N analytical results\n",
    "m_largeN_vals = [avg_m_approx(T, Dels, Nelems) for T in Tvals]\n",
    "plt.plot(Tvals, m_largeN_vals/np.sum(Nelems),linestyle = '--', linewidth = 3.0, label = f'Large $N$')     \n",
    "\n",
    "# plot of critical temperature \n",
    "plt.axvline(x = kBTcrit(Dels, Nelems), color = 'k',  linestyle = '-.')   \n",
    "\n",
    "plt.legend(loc = (.5,0.05), fontsize = 12)\n",
    "# plot formatting\n",
    "ax.set_xlabel(r'$k_B T$', fontsize = 18)\n",
    "plt.xlim([-0.01,3.1])\n",
    "plt.ylim([0,1.1])\n",
    "plt.ylabel(r'$\\langle m \\rangle/N$', fontsize = 18)\n",
    "# plt.yaxis.set_label_coords(-0.1,.5)\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",
    "\n",
    "ax.text(kBTcrit(Dels, Nelems)-.2, 0.25, r'$k_BT_{derang}$', color='black', fontsize = 15,\n",
    "        bbox=dict(facecolor='white', edgecolor='none', pad=5.0))\n",
    "\n",
    "for i in range(4):\n",
    "    ax.text(temp_vals[img_key_list[i]], sim_m_vals[img_key_list[i]]/np.sum(Nelems)+.05,'('+'i'*(1+i)+')' if i<3 else '(iv)', fontsize = 14 )\n",
    "\n",
    "# plt.savefig(f'derang_model.png', bbox_inches='tight', format = 'png')    \n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Total Notebook Runtime: 9.734 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
}