kmc.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Two-Transition Path Kinetic Monte Carlo Simulation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np \n",
    "from matplotlib import pyplot as plt \n",
    "import scipy as sp"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We first define the Hamiltonian in order to properly simulate the system.\n",
    "\n",
    "We have a system with two basins $A$ and $B$, each with 3 microstates labeled as below:\n",
    "\n",
    "$$ \\{0,1,2\\} \\in A; \\{3,4,5\\} \\in B $$\n",
    "\n",
    "We arbitrarily define our reaction landscape such that Basin $A$ has a lower energy than Basin $B$.\n",
    "\n",
    "Let the energy of basin $A$ defined as $E_A$, and the energy of basin $B$ defined as $E_B$.\n",
    "\n",
    "The free energy of each basin $i$ can be written as:\n",
    "\n",
    "$$ F_i = E_i - \\beta^{-1} \\log \\Omega $$\n",
    "\n",
    "where $\\Omega$ is the number of microstates in each basin (here it is simply 3) and $\\beta^{-1} = k_B T$.\n",
    "\n",
    "Based on transition state theory, the reaction rate for a reaction $X \\to Y$ is given by the Arrhenius rate equation:\n",
    "\n",
    "$$ k_{X \\to Y} \\propto e^{-\\beta \\Delta F_{X \\to Y} }$$\n",
    "\n",
    "where $\\Delta F_{X \\to Y}$ is the height of the activation energy barrier (difference in free energy) in the direction $X \\to Y$. \n",
    "\n",
    "Note that we also have transitions between microstates within a basin, albeit with a relatively small activation energy energy barrier. \n",
    "\n",
    "To specify our system entirely, we must thus define $F_I, F_{II}, F_A, F_B$ where $F_I$ and $F_{II}$ are the free energies of the saddle points in transition paths $I$ and $II$, and $F_A$ and $F_B$ correspond to the free energies of basins $A$ and $B$. In addition, we must also include the reaction rates within each basin: for simplicity, we make these activation energy heights the same across all basins:, on the order of $k_B T = \\beta^{-1}$. \n",
    "\n",
    "To tune the height of the activation energy barrier for our two transition paths $I$ and $II$, we define two constants $\\alpha_1$ and $\\alpha_2$, which scale the activation energy barrier for transition paths $I$ and $II$ respectively. \n",
    "\n",
    "To simulate the system realistically, even though the true microstates all have the same energy, we add a noise term $\\epsilon \\sim \\mathcal{N}(0,(.01\\cdot\\beta^{-1})^2)$ to the underlying true energy ($E_A$ or $E_B$) of our system.\n",
    "\n",
    "Thus, we can now generate our 6-state MSM and the corresponding KMC transition matrix based on the energies we impose on our system.\n",
    "\n",
    "A simple implementation of a 6-state MSM KMC simulation is given below, where we have the connections:\n",
    "\n",
    "$$ 1 \\to 3, 2 \\to 4 $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Free energy of microstate 0: -0.008948107850882564\n",
      "Free energy of microstate 1: -0.010613674409235347\n",
      "Free energy of microstate 2: -0.008703739012656626\n",
      "Free energy of microstate 3: 0.29040840056337686\n",
      "Free energy of microstate 4: 0.2903698349330016\n",
      "Free energy of microstate 5: 0.28978714378995396\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([[0.00000000e+00, 3.67879441e-01, 3.66981556e-01, 0.00000000e+00,\n",
       "        0.00000000e+00, 0.00000000e+00],\n",
       "       [3.61802909e-01, 0.00000000e+00, 3.60919855e-01, 2.23734481e-06,\n",
       "        0.00000000e+00, 0.00000000e+00],\n",
       "       [3.67879441e-01, 3.67879441e-01, 0.00000000e+00, 0.00000000e+00,\n",
       "        3.38584881e-04, 0.00000000e+00],\n",
       "       [0.00000000e+00, 4.53999298e-05, 0.00000000e+00, 0.00000000e+00,\n",
       "        3.67879441e-01, 3.67879441e-01],\n",
       "       [0.00000000e+00, 0.00000000e+00, 6.73794700e-03, 3.67737593e-01,\n",
       "        0.00000000e+00, 3.67879441e-01],\n",
       "       [0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 3.65601050e-01,\n",
       "        3.65742073e-01, 0.00000000e+00]])"
      ]
     },
     "execution_count": 42,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# generate 6-state MSM and corresponding transition matrix based on rates\n",
    "from itertools import permutations \n",
    "\n",
    "beta = 10\n",
    "\n",
    "e_a = 1/beta\n",
    "e_b = 4/beta\n",
    "\n",
    "alpha1 = 10\n",
    "alpha2 = 5\n",
    "\n",
    "def gen_MSM(alpha1,alpha2):\n",
    "    # define free energies \n",
    "    f_A = e_a - 1/beta * np.log(3)\n",
    "    f_B = e_b - 1/beta * np.log(3)\n",
    "\n",
    "    # define free energies of states in basin A and B\n",
    "    f_basinA = np.random.normal(loc=0,scale=.01/beta,size=3) + np.ones(3)*f_A\n",
    "    f_basinB = np.random.normal(loc=0,scale=.01/beta,size=3) + np.ones(3)*f_B\n",
    "    \n",
    "    f_system = np.array([f_basinA, f_basinB]).flatten()\n",
    "    \n",
    "    for i in range(len(f_system)):\n",
    "        print(f\"Free energy of microstate {i}: {f_system[i]}\")\n",
    "    \n",
    "    # we now calculate reaction rates based on these free energies assigned to our microstates:\n",
    "    # our connections are as follows: \n",
    "    # 0<->1, 1<->2, 2<->0, 3<->4, 4<->5, 5<->3, 1<->3 (transition path 1), 2<->4 (transition path 2)\n",
    "    \n",
    "    deltaf_saddle1 = alpha1 / beta\n",
    "    deltaf_saddle2 = alpha2 / beta\n",
    "    deltaf_basin = 1 / beta\n",
    "    \n",
    "#     print(deltaf_saddle1, deltaf_saddle2)\n",
    "    \n",
    "    # generate probability transition matrix \n",
    "    P = np.zeros((6,6))\n",
    "    for (i,j) in list(permutations(list(range(6)),r=2)):\n",
    "        \n",
    "        if i == j: # both states are the same\n",
    "            rate_ij = 0\n",
    "        elif i in list(range(3)) and j in list(range(3)): # both states are in basin A\n",
    "            f_activation = max(f_system[i]+deltaf_basin,f_system[j]+deltaf_basin)\n",
    "            rate_ij = np.exp(-beta*(f_activation-f_system[i]))\n",
    "#             print(f'{i},{j},{f_activation}') \n",
    "        elif (i not in list(range(3)) and j in list(range(3))) or (i in list(range(3)) and j not in list(range(3))): # one is in A, the other is in B\n",
    "            if (i == 1 and j == 3) or (i == 3 and j == 1):\n",
    "                f_activation = max(f_system[i]+deltaf_saddle1,f_system[j]+deltaf_saddle1)\n",
    "                rate_ij = np.exp(-beta*(f_activation-f_system[i]))\n",
    "            elif (i == 2 and j == 4) or (i == 4 and j == 2):\n",
    "                f_activation = max(f_system[i]+deltaf_saddle2,f_system[j]+deltaf_saddle2)\n",
    "                rate_ij = np.exp(-beta*(f_activation-f_system[i]))\n",
    "            else:\n",
    "                rate_ij = 0\n",
    "#             print(f'{i},{j},{f_activation}') \n",
    "        else: # both states are in basin B\n",
    "            f_activation = max(f_system[i]+deltaf_basin,f_system[j]+deltaf_basin)\n",
    "            rate_ij = np.exp(-beta*(f_activation-f_system[i]))\n",
    "#             print(f'{i},{j},{f_activation}') \n",
    "        \n",
    "        P[i,j] = rate_ij\n",
    "        \n",
    "    return P\n",
    "        \n",
    "gen_MSM(alpha_1,alpha_2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Free energy of microstate 0: -0.011256228866484388\n",
      "Free energy of microstate 1: -0.00886734782603811\n",
      "Free energy of microstate 2: -0.009429632157401929\n",
      "Free energy of microstate 3: 0.2884027322764716\n",
      "Free energy of microstate 4: 0.2899496686670319\n",
      "Free energy of microstate 5: 0.28966178922016966\n",
      "[[0.00000000e+00 3.59195378e-01 3.61220766e-01 0.00000000e+00\n",
      "  0.00000000e+00 0.00000000e+00]\n",
      " [3.67879441e-01 0.00000000e+00 3.67879441e-01 2.32288456e-06\n",
      "  0.00000000e+00 0.00000000e+00]\n",
      " [3.67879441e-01 3.65816717e-01 0.00000000e+00 0.00000000e+00\n",
      "  3.37551317e-04 0.00000000e+00]\n",
      " [0.00000000e+00 4.53999298e-05 0.00000000e+00 0.00000000e+00\n",
      "  3.62232371e-01 3.63276666e-01]\n",
      " [0.00000000e+00 0.00000000e+00 6.73794700e-03 3.67879441e-01\n",
      "  0.00000000e+00 3.67879441e-01]\n",
      " [0.00000000e+00 0.00000000e+00 0.00000000e+00 3.67879441e-01\n",
      "  3.66821915e-01 0.00000000e+00]]\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(6, 6)"
      ]
     },
     "execution_count": 56,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# define transition matrix\n",
    "P = gen_MSM(alpha1,alpha2)\n",
    "print(P)\n",
    "P.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [],
   "source": [
    "# implement KMC without rejection and define helper functions\n",
    "def kmc(P: np.array, epochs: int,init_state: int):\n",
    "    \n",
    "    init_vec_t0 = np.array([int(j == init_state) for j in range(P.shape[0])]) # one-hot vector \n",
    "    \n",
    "    R_ki = np.cumsum(P,axis=1) # returns matrix of cumulative functions \n",
    "    \n",
    "    Q_k = R_ki[:,-1] # total rates\n",
    "    \n",
    "    curr_vec = init_vec_t0\n",
    "    \n",
    "    states_list = [init_state]            # to keep track of state values\n",
    "    states_matrix = [init_vec_t0]           # matrix to keep track of state vectors\n",
    "    times = [0]                           # to keep track of timesteps\n",
    "    \n",
    "    prob_vector = init_vec_t0\n",
    "    \n",
    "    for _ in range(epochs):\n",
    "        \n",
    "        k = np.argwhere(curr_vec == 1)[0][0]\n",
    "        \n",
    "        u = np.random.uniform() # generate random number from uniform distribution\n",
    "        \n",
    "        # we now want to find i such that R_{k,i-1}< u * Q_{k} <= R_{ki} \n",
    "        \n",
    "        i = np.searchsorted(R_ki[k,:], u*Q_k[k])\n",
    "        states_list.append(i)\n",
    "        \n",
    "        curr_vec = np.array([int(j == i) for j in range(P.shape[0])]) # transition state\n",
    "        states_matrix.append(curr_vec)\n",
    "        \n",
    "        v = np.random.uniform() # draw uniform number for timestep\n",
    "        times.append(times[-1] + 1/Q_k[k] * np.log(1/v))\n",
    "        \n",
    "    return np.array(times), np.array(states_list), np.array(states_matrix).T\n",
    "\n",
    "def probabilities(state_matrix):\n",
    "    return np.sum(state_matrix,axis=1) / state_matrix.shape[1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 68,
   "metadata": {},
   "outputs": [],
   "source": [
    "times, states, state_matrix = kmc(P,epochs=100000,init_state=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 73,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAEWCAYAAABmE+CbAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO29d5wkWXXn+zsRGVlp2lRWd483PQ4rhGtGDHaAwYqVVhJ6kpZlMRIjIQPsSuKhhx4MEtLKffah1fIwkhAgEB5pJQTCSFgBM9MDjJ+BmR7X49pUVXdXZWZFZMTdP+LeiBs3w6WrzMo638+nPpUZGRlx80bEPffYS0IIMAzDMIw17QYwDMMwswELBIZhGAYACwSGYRhGwgKBYRiGAcACgWEYhpGwQGAYhmEAsEBgmNIQ0RoRXZjz+c1EdPkmNmkiEJEgooun3Q5m82GBwAwNEV1FRB8eYP/LiejwCOcb6fujIoTYIYQ4JNvyASJ6h/H5Y4UQX51K4xhmDLBAYOYKIqpMuw3TZLv/fmY0WCAwhRDR/01E9xPRKSK6nYieR0QvAvD/APg5aUq5Xu77aiK6Ve57iIh+WW5vAvg8gLPk/mtEdBYRWUT0ZiK6k4iOE9EniGgppQ1Z37+KiD5FRB8mopMAXkVElxLRt4lolYgeJKL/RURV7ViCiH6FiH5IRCtE9C4iIvnZxUT0NSI6QUTHiOjjxvcuJqIrAbwcwJtkO/5Jfn43EV0hXy8Q0TuJ6AH5904iWpCfXU5Eh4noN4noiGzjq3P6/wIi+rrs0y/L9n5YfrZftusXieheAP8mt3+SiB6Sv+PrRPRY7XgfIKL3ENGX5DG/RkTnG6e9Iq1/mDlHCMF//Jf5B+CRAO4DcJZ8vx/ARfL1VQA+bOz/4wAuAkAAng2gDeBJ8rPLARw29n8jgO8AOAfAAoD3AvhoRlvSvn8VAA/Af0Q4wakDeDKApwKoyPbeCuCN2ncEgM8CWARwHoCjAF4kP/sogLfIY9UAPMP43sXy9QcAvMNoy90ArpCvf0/+rtMA7APwLQC/r/2OntzHAfAS2U+tjN/9bQB/BqAK4BkATqp+l79PAPgQgCaAutz+GgA7ZZ++E8D3teN9AMApAM+Sn/85gG+W6R/+m+8/1hCYInyEg8ZjiMgRQtwthLgza2chxD8LIe4UIV8D8EUAz8w5/i8DeIsQ4rAQYgPhAP+yAU0f3xZC/IMQIhBCdIQQ1wkhviOE6Akh7kYoZJ5tfOePhBCrQoh7AXwFwBPkdg/A+QgFYFcI8c0B2qHzcgC/J4Q4IoQ4CuDtAF6hfe7Jzz0hxOcArCEUvgmI6DwATwHwViGEK9vzjynnu0oIsS6E6ACAEOL9QohTWp8+noh2a/v/sxDi6/LztwC4jIjO1T7P6h9mjmGBwOQihLgD4Sz+KgBHiOhjRHRW1v5E9GIi+g4RLRPRKsLZ796cU5wP4O+leWcV4WzeB3D6AM28z2jDI4jos9JkchLAH6a04SHtdRvADvn6TQi1m2sojBp6zQDt0DkLwD3a+3vkNsVxIUQvow3mcZaFEG1t230p+0XbiMgmoj+SZriTCDUXINkH0f5CiDUAy0b7svqHmWNYIDCFCCH+TgjxDISDtwDwx+ojfT9pI/80QvPG6UKIRQCfQzjA9u0vuQ/Ai4UQi9pfTQhxf1pTsppovH83gNsAXCKE2IXQ11HKBi6EeEgI8VohxFkItZf/n9JDMIvKBD+AsL8U58ltg/IggCUiamjbzk3ZT2/PfwLwkwCuALAboVkJSPZBdAwi2gFgacj2MXMECwQmFyJ6JBE9Vw72XQAdhDN4AHgYwH4iUvdRFaF56SiAHhG9GMALtMM9DGCPYbp4D4A/UE5NItpHRD+Z0Zy076exE6GdfY2IHgXgdWV+qzz/zxLROfLtCsKB1k/Z9WEAmTkJCH0Rvyt/z14AbwVQOkRXIYS4B8BBAFcRUZWILgPwHwq+thPABoDjABoINSSTlxDRM6Sz/fcBXC2ESNM8mG0ECwSmiAUAfwTgGEIzwmkIZ9wA8En5/zgRfVcIcQrA6wF8AuFg+p+g2buFELchHCgPSRPRWQgdmv8I4ItEdAqhI/bH0hqS8f00fkue+xSAvwTw8Yz90ngKgKuJaE226w1CiLtS9vtrhH6VVSL6h5TP34FwIL8BwI0Aviu3DcPLAVyGcIB/B8Lfs5Gz/4cQmqjuB3ALwj41+TsAb0NoKnqyPAezzSEheIEchtlKyFDY24QQbxvy+x9AGK31u2NtGLPlYQ2BYWYcInoKEV1EYc7GixD6B9K0EoYZCc5qZJjZ5wwAnwGwB8BhAK8TQnxvuk1i5hE2GTEMwzAA2GTEMAzDSGbKZLR3716xf//+aTeDYRhmy3DdddcdE0LsG8exZkog7N+/HwcPHpx2MxiGYbYMRHRP8V7lYJMRwzAMA4AFAsMwDCNhgcAwDMMAYIHAMAzDSFggMAzDMAAmHGVERHcjLDDmA+gJIQ5M8nwMwzDM8GxG2OlzhBDHNuE8DMMwzAiwyYhhGKYk3/zhMdx9bH3azZgYkxYIAmGd++uI6Mq0HYjoSiI6SEQHjx49OuHmMAzDDM9vffJ6/OU3Dk27GRNj0gLh6UKIJwF4MYBfI6JnmTsIId4nhDgghDiwb99Ysq8ZhmEmQrfnw+0F027GxJioQBBCPCD/HwHw9wAuneT5GIZhJknPF/CD+a0QPTGBQERNItqpXiNcW/emSZ2PYRhm0rh+gN4cC4RJRhmdDuDviUid5++EEP8ywfMxDMNMlJ4fzLWGMDGBIIQ4BODxkzo+wzDMZhIEAoEAegH7EBiGYbY1nhQE86whsEBgGIYpQc8PBQELBIZhmG2OEgjz7FRmgcAwDFMCNhkxDMMwAFhDYBiGYSSezxoCwzAMg1gzYIHAMAyzzemxhsAwDMMAgMc+BIZhGAaIM5R9zlRmGIbZ3rCGwDAMwwDgKCOGYRhGwqUrGIZhGACcqcwwDMNIOFOZYRiGAcB5CAzDMIzEk4JACYZ5hAUCwzBMCZQgmGMFgQUCwzBMGWIfAmsIDMMw2xqOMmIYhmEAcJQRwzAMI1GZykIAwZwKBRYIDMMwJdA1g3nVElggMAzDlEAPNw0ECwSGYZhti+uzhsAwDMMgqSH4PgsEhmGYbUvShzCfuQgsEBiGYUrg6RoCm4wYhmG2Lz3NTOSzU5lhGGb7opuJeuxDYBiG2b54uobAJqPhICKbiL5HRJ+d9LkYhmEmhR5lxGGnw/MGALduwnkYhmEmhhewhjASRHQOgB8H8FeTPA/DMEm+dMvDeNm7vzW3NXemQVJD4LDTYXgngDcByOw9IrqSiA4S0cGjR49OuDkMsz244fAqDt6zAneOV/fabHRH8pzKg8kJBCJ6KYAjQojr8vYTQrxPCHFACHFg3759k2oOw2wrOq4PYH5NG9PA48S0kXg6gJ8gorsBfAzAc4nowxM8H8Mwko4XCoR5DY+cBj1OTBseIcTvCCHOEULsB/DzAP5NCPGfJ3U+hmFilEDw5nQmOw16XNyOYZityIYXCgLWEMaH6wdYqIRD5rxqCJXNOIkQ4qsAvroZ52IYRtMQ2Kk8NnpBKBA2esHcCgTWEBhmDlFO5Xk1bUyDni9Qc2wA86shsEBgmDkkdiqzhjAuPD+IBMK8CloWCAwzh3Q91hDGTS8QqEcawnwKWhYIDDOHcNjp+AlNRuGQOa+ClgUCw8whyofAYafjw/MDLFTYh8AwzBaDNYTx0wsEFpz5DjtlgcAwc0iXncpjh53KDMNsOXp+EC3m4s3pwDUNOOyUYZgtR7enL/XIGsK46AUBahV2KjMMs4VQDmVgfgeuzUYIAc/XfAhzKmhZIDDMnKH8BwA7lceFMhHVVJTRnHYrCwSGmTM6ukDgsNOxoDStGiemMQyzldBNRt68TmU3GbXyHCemMQyzpUhoCHNq695slOkt0hDmVNCyQGCYOUMXCBx2Oh6UYF3gKCOGYbYSG6whjB0lWB3bgkWch8AwzBahw1FGY0cJVse2ULEs+GI++5UFAsPMGR1XS0yb05nsZqOc8xWbYFvEGgLDMFsDdiqPHxW+G2oINLeaFwsEhpkzuuxUHjtKAFQsgmXR3OYhVKbdAIZhxkvH9WFbBAJrCOPC8w0NYU4FLWsIDDNndDwfdcdGxZ7fgWuzUf2ofAgBO5UZhtkKdDwfNceGY1nRzJYZDdWPFWu+fQhsMmKYOaPr+ag5FvwgmNuBa7NR/ejYBNvmKCOGYbYI3chkZLHJaEyoKKOKbcGm+TXFsUBgmDmj4/qoV204FrFTeUx4WpQR5yEwDLNlUD4Em53KY8MzMpXntaw4CwSGmTM6XoA6O5XHSq8vU3nKDZoQLBAYZs7oulrYKTuVx0KkIVgWKvb8JqaxQGCYOaPjhT6EeTZtbDZmHsK8muImJhCIqEZE1xDR9UR0MxG9fVLnYhgmpiPDTh2beMW0MaFXO7Vpfp3Kk8xD2ADwXCHEGhE5AL5JRJ8XQnxngudkmG1P1w2dymHYKWsI48DT8xDmWEOYmEAQQggAa/KtI//msxfHyKGja/jK7Uej9xfta+LyR542xRZl89Xbj+Ap+5fQXNh++Y1HTnVx33IbTz5/adpN6aPbkz6EOc6o3Wz0PISKTdjw5lPQTvRJJiIbwHUALgbwLiHE1Sn7XAngSgA477zzJtmcLcH/+NIP8NkbHozeVysWbv/9F4GIptiqfo6c7OJVf3Mt/vhnHoefe8r2u25//Y278LFr78P1b3vBtJuSwPMDeL4Io4xsK1EKmxmeZB6ChV4wn/06UaeyEMIXQjwBwDkALiWiH0nZ531CiANCiAP79u2bZHO2BB3Xx6PO2Inr3/YC/LfnPwJuL8DaRm/azerj6NoGAGBtYz4fjCJOdns41fUgZqzImSp9Xa/aoWljXuMjN5m4dEVYy2hefQibEmUkhFgF8FUAL9qM821lXD9AvWpjd93BGbtqAIDVtjflVvWzsh62abvGuXc9H4EIr9csoTSCmmOzU3mM9IIARIBtESwuXTE4RLSPiBbl6zqAKwDcNqnzzQtuL4Bjh5dld8MBMKMCoe0CALzebA2Im0XHDQferjtbv1+1J/QhsFN5XHi+gGOFz2WFF8gZijMBfFD6ESwAnxBCfHaC55sLPD+InLStRhUAsNpxp9mkVJRAmLUZ8mbR7fnR/91wptyamI5mMuLEtPHR8wNU7NCPN88lQSYZZXQDgCdO6vjziusHWJQawqLUEFZmUENYXt/eAkFpCOr/rBCbjCw4tgVvTmeym00vEKhYoUCoWIRgTgUCZyrPGF5PoGoIhBPt2dMQlBnL3aYmI+W8nbUoHiWgajLs1GcNYSx4fmzKnec8BBYIM4brB3AqUiDUQ5PRLGsI29epHMj/syUQlClLrYfgzenAtdl4msmIo4yYTcPtBZGGUK1YaFbtmXYqb1cNoTOjGkLX1XwIHHY6Nnq+QMViDYHZZFw/QLUSJ6EtNqpYnUGTURRltE1NEkoQzJqGEDmVudrpWPECAUc5lbe7hkBEDSL6f4noL+X7S4jopZNt2vbE82MNAQj9CKudGdQQ1tmHEP6frd+vCwR2Ko+PMMpIhZ1a21sgAPgbhMXqLpPvDwN4x0RatM3R8xCAUCCszKCGsN2jjCKn8qxFGSmncpVrGY0TzxcJp/J2FwgXCSH+BIAHAEKIDoDZKq4zJ3h+gGpFFwhVnJgxH0LX86OZ6HbUEHqyXhAwgz4EFXZaUdVOxcyV19iK9IIgMhlVLJrbhL+yAsGV2cYCAIjoIoQaAzNGgkAkZiIAsFifPQ1Bb892jDLqakJwFn0ItkVwbIIj4+bn1QG6mYRO5bA/rTnWEMompl0F4F8AnEtEHwHwdACvnlSjtivK3qtrCK1GFSc6HoJAwLJmQylT5iJgewoE3Uw0cwLBDddTJqLI5u0HAo495YZtcbyED2F+o4xKCQQhxBeJ6DoAT0VoKnqDEOLYRFu2DVHmF9OpHAjgVLcX1TaaNioMtu7Y29JkpAuBWXMqd3vh4jgAIhOH5wfRNmY4eoFAzYl9CEJgpiZp46JslNG/CiGOCyH+WQjxWSHEMSL610k3bruh7NKmDwGYrXpGSkM4Y3cN7jZ0Wup+g5nzIbg+6tV44ALAjuUx0PODKA9BmY78OfTN5AoEuS7yEoC9RNQioiX5tx/AWZvRwO2Emm3rPoTWDFY8VXkRp+1cgNubrQFxM+jOsEDoeOFqaQAiEweHno5O6NtTeQixKW7eKDIZ/TKANyIc/K9DHFl0EsC7JtiubYmyxyc1BFXgbpY0hFA47du5gAdPdKfcms1npn0ImkBwWEMYG70g0DKV1bb569dcgSCE+HMAf05EvyGE+ItNatO2ZSPSEJKZygBwYoaS01baLnbVKmhUt6cPoePNsEBwfSwYGgILhNHp+SIuf600hDns17JO5b+Qy18+BkBN2/6hSTVsO6I0hIVKMuwUAFbWZ0dDWGm7aDWrYSbsNowyUo7kikUzl5jW9fxoEhE5ldlkNDJeECeMVqJw3vnr11ICgYjeBuByhALhcwBeDOCbAFggjJE0H8JuKRBmqXzF8rqLVqOKasXalhqC0gpazerMRRl1PB9nKg1hjm3dm43Xi/MQ7O3qVNZ4GYDnAXhICPFqAI8HsDCxVm1T0nwIFdvCzlplppzKK20XrYaDqm1ty9IVymTUajiz6VSuKpNRHHbKjEYvSOYhAPMpaMsKhI4QIgDQI6JdAI4AuHByzdqepGkIQJicNksVT1fWPbSaUkPwg21XGkFpCIuN6sz5ELpe0JeHwD6E0UlGGc1vv5bNVD5IRIsA/hJhtNEagGsm1qptipuiIQCqwN1saQhLjdCHIEQ4U6rY85Wgk4fSCpYaVRw5OVtRVl03jjJSzs95tHVvNj1jxTRgPjWEsk7lX5Uv30NE/wJgl1wzmRkjaZnKgFwTYUZ8CF3PR9v10WpWowfD1dL6twNd1wdR6N+ZTZNReC1U2Ol2XbNinHiBHmU0vzWiSmcqq9dCiLuFEDdwpvL4SctUBsJIo1kxGSlfRktqCEDocNtOdHsBahUb9ao9U1FGnh+EJRYqHHY6bnp+ACfKVJ5fZ32uhkBENQANyExlxIlpu8CZymPH9cPBpd+H4MyMU1mVrVhqOji6Fj4Q282x3HFDx23NsROVT6dNtDiO6VRmk9FIBIFAINCnIWw7gYD+TGXFKXCm8thRM21TQ9jdqOJk14MfiOhmnBYqY3pRVmEFtqFA8HzUKhZqThh2OwvXBYjXU46cyhZrCONACVQzD2EeBUKRyehbAJ4G4LeEEBcCeDuAmwB8DcDfTbht244Nvz9TGQg1BCGAkzPgR1ACYUlGGQGAN0Oz5M2g6/moVe3IebsxI/Wc9OUzgXhG67OGMBJKoJp5CPPorC8SCO8FsCEzlZ8F4L8D+CCAEwDeN+nGbTfUwLpgJ0sVz1I9I5UxrfsQtpuG0JX1gpRpZlb8CKbJKC5/PX8z2c0kEggcZQRbCLEsX/8cgPcJIT4N4NNE9P3JNm37oQZWp5LUEOIS2LOgIYRtWGw4sUDYZhpCxwvXHFDO21mJNFKCicNOx0tsMoqX0AS2Z5SRTURKaDwPwL9pn5XNYWBK4mWFnaryFTOgISyvu9hZq8CxrchktP00hHBVspqcic9K+QrVjlpUuoI1hHEQm4ySGkIwhwKhaFD/KICvEdExAB0A3wAAIroYodmIGSOuH4AIfQ7KltIQZiDSaKXtYqkZtqdqb08fQsf10WpUo5n4rGQrq3aolb0cDjsdC6r0h/LJqP/zqCEUlb/+A5lvcCaAL4q4RoEF4Dcm3bjthusHqNoWiEyTkfIhzIJA8CIT1vbVEHzUHCsaeGfGZJQRdsomo9HwjGAPi7avDwFCiO+kbPvBZJqzvXF7QZ+5CAB21RwQASdmwGS0su5i7w5VXllqCNtQINQde+Y0BNOHoMJO2WQ0GkoTqBiJafOoIWyfegNbAM8P+nIQAMCyCLvrs1HPSJW+BmKT0XZ0KqvENGAGo4w47HSsmBpCHGU0f/06MYFAROcS0VeI6FYiupmI3jCpc80Lbi/oy1JWtGakntGqXBwHAKoVVcto/mZKeURRRs5sRRlFPoS+8tfb6/qMG9OpvG19CCPSA/CbQojvEtFOANcR0ZeEELdM8JxbGs8XqRoCEBZSm3aUUdfzse76mlM5HHi2k1NZCBGVmFa2+o0ZiTLKMhmxU3k0lA/GqXAewtAIIR4E8KB8fYqIbgVwNoCZFwjrGz38zLu/heMyCcsmwh/9zONw+SNPi/bxA4HXfOBavPaZF+IZl+wtddw3fOx7+Nadx6P3P/Okc/DmFz8qeh9qCOklEFoNB9/44TE85Q++DAA4p1XHJ3/5skSV0QdPdPC6D38X73vFk3HarmilUxw52cXPv+87OLXRAxAOGB949VNw4b4dhW0+2fXwsnd/CyttL1r3QDm5nUhDiAfEjZ6PV73/WrzpRY/EE89rpR6z6/l45fuvwe+85NF4wrmLhW0IAoH/673fxj3L7b7PbCL84U//CJ77qNMzv//ur96JtQ0Pv/3CR2Xuk8VXbj+CD3/7HvzVKw+AiKJ1r+uOjVql36n8pVsexicO3oe//C8HEse5/aFTeMvf34gPvuZSNBfix04IgV/84EG88mn78exH7Bu4fTodz4dtUaRlWhaBKOlU/v59q/i1j3w3NRCgUbXxt6/5MZy3p1F4LiEEXvuhg3j5U8/Hc7TnIo9f+dvr8NNPOhsveOwZmfv86RduQ3Ohgl+9/OLMfd7+Tzfjszc8mPrZjz/uTFz1E49NbLv1wZN46/++CR94dbLvbzi8inf886340GsujbS9NJSGparHDlK64uPX3ovrD5/AH/7U4wr3nQU2xYdARPsBPBHA1SmfXUlEB4no4NGjRzejOYXcfXwdtz10Co8+cxeuePTpOHKqi4N3ryT2WWm7+NoPjuJbdx4rfdx/u/UI9jSruOLRp8OxCN82vuv6AaqV9BvzymddhJ89cC6uePTpuGBvE9+7dzUSWIrr7zuB79+3ipseSEYE3/7wKRw6to4nnLuISy9Ywr3Lbdzy4MlSbb5vuY0fPLyGR52xE89/zBl41dP24/mPDgffaopT+diai28fOo5r7lpOPR4A3L/awdV3LePqQ8cz99FZc3s4eM8Kzm3VccWjT0/8PXyqi+/fu5r7/a/94Ai+cPPDpc5lcu1dy/jX245gTQrTjhuHdkaZyppA+Pc7juFLtzzc51e45u5lHLxnBfevdhLb266Pf7vtCL51R/n7KIsTHS9aclXhWFbCZHTT/Sdw/2oHz7x4b6Ifn3xeC/ccb+OHR06VOtdGL8CXby3f7q7n419ufghX59wXAPCFmx/Gl2/Jv1b/fscx1Byr717YsVDB137QP4Zce/cyrr17BXcfX09sv/rQMq65axmHV/onGjpmprKKMipjMvrKbUfxpYLfM0tMPLmMiHYA+DSANwoh+kYhIcT7IMtgHDhwYCZ0MFXR89efczEuvWAJX7j5ISwb5hplvilbTsLzA5za6OEljzsTr3/eJfiVv70Oh46tJfYJo4zSNYTLLtqDyy7aAwD4/I0P4pq7lrG87uJ0TRNQbVleT/oa1O950wsfid11B/98w4OlcxqUXfq1z7wQzzJmsEqF1p3KaiDMc4CvRu0s13eqaNtPP+kc/Oennp/47B+/fz/aBU7djhcMbW5Tx15Z97Cz5qDbi80yKlNZjzJSv2m57eLsaj3aviq3m21VwqRsX+QROvyTAqFiE3qawFZt/b3/+CPYoc2W7zhyCv9y80Ol/SHqOGUDHdT9VnT81bZbmPDV8XwcOH8J//2nk7Puq/7xZnzmu4f79j++5ibaoFDP9PE1FxfnKDkqU9nMQyijIay03ej+3QpMVEMgIgehMPiIEOIzkzzXONFLPAOq/LRr7OMl9i1CDdbKIVuv2n0PR1aUkYnKA1hZN9vkpm5fjcpNVLFLziBPlHRQd9wgaq9JNaWWUTRQ5PSL6jtTw8nCzMDVSevHvu+7fsLkNQjxwBe2NbLTV21YFqFasRLnj4TymnFtjO8rygjQsiyvu9jTTC51XrEoMZONNBzjPhs0Ykpdk7KCVvVL3uAohMBK2+ubfKWdW+WA6LQaVZzs9vrCoNW5TYGg7tGiSV0vMhklfQhlNISVthtNIrYCk4wyIgB/DeBWIcT/mNR5JkEsEMKHq9WoYiVj1j3oDEnN4GqOHQ22irwoIx3l1DUfnBVtdprYHpWsdkJnqGOXfpDNUEadtFpG5gCaRtTOsgKhl8zA1ak5xQKh4/nwA4GT3V6p85nfBTSBIN8vSO2g7tiJQS4SCH0aZXj9zZwF8/ijsLzuotU0TEa2lRggO56Pqm31rXBXHzBialDNRl3zvOOf7PbgBwInOl7u7Lvr+VH/6yzJ/BizL49nDPyq7UUTk56ZqSwFQ5nSFSttD54vElraLDNJDeHpAF4B4LlE9H3595IJnm9sLK+70RKJQDizNm+myGRU8oGIhIyc3dcdu29wKKshqIe+T0PIaNNq24vqD4W/p/yCO+oBThuMbYtgW9Q34KhzZrE8oMnIjJ7RSevHvu+X0FqKzh0P6EmNKTx//PvVxCFLezMHxNgkNbpACMuKGBqCTYnBte36qddy0MqtZr8Ut63YZKSeKSHyNdiuzAMxUc+WeV+pvjUnQSsln2FP9p+Zh1CkIQghomPP0kJKeUwyyuibiFdY21Ko5Ct14VsNBzfdn25/LDuzUzeGMvfUq6GpQQgRlarYKKkhtKIbP10F7nsg2nEyGYCBktzMRVdMqraV4UPI0RAG9SF42W0os4ylrrXsR7PUORXmDL5raEw1J2kyypp1mian6PjyfZGZpIggCM0te5rVxPaK4VTOGkwHrdyqtLay7c76/Tr6/bC8HtfM0vEDAc+PlwnViTTnDGHc50MYVEOwzLDT/EF+baMXCY2u5yd8NrMKZyqnYN6MrWaoIeg26Nj+6JVWHYH4pq07dnRzK8pqCI5tYWetkqkCm9tX2l7C2bjYcHCiM7rJKGwLJX5DPIBmC5zBTUbKhzCcyaiMGav4u3KG6yY1ppqmoXQ9P20++k8AACAASURBVFMbMU1O5vFPdLyRzArKzNIyBYJNibDTjiy7YWJZhAXDH5KHmigUmXcU0Uw5V0PwtNfp18os4KdTJBDMe1K9L9IQ4iijwcpf62bmWSlvUgQLhBSW191I/QTCAXSjFxgzwfBi+4HAqRK2ad2ODyA1y1UVtyvDUrPaf+NnzLxX226kmQDAYr06sMkobVYJANWKHcXmA3GS1qohQBPtlH23ttErtdpYJ0dLqTs2OjmJYZ4fRALL9AOVITaBKdU/KSB1p7YucPr9O+k+BGUyKjKTFKHO168hUCIxreP6mdpevWqXjohR/VB2JT81+OZFhJkaQhp596MSCPoAHzqqw/f6JMgPRHRNizQEcwnNSEMoSPjT7wcWCFsYU0NQwmElYwZTRm1eWXfR0OrfxFmu8Y3i9URpgdBK8WtE9uu+mZDbpyGULYMRR6VkmYzSfQi9QESx+yaJgbOElrDRyxcIeYOY/iAOoyFEJh3lAzCEU60SawiJAU2LMnJ7QV8eQ3T8RPtGEAhqJTtDIKQ5lbOEe72EthUdRwuIKHP/r2ZoSDr69cm6VpGGkHI/qntcH+BPbfTiCYHWvyc7HtQEv+i+UJn40RKaJfMQlhMCYWv4EFggpLCi1esB0sM8l9su1LIFZQa1ZcOOnxbV4fpB32ppWZgawkbPx9pGD45NWG27CTV+Zd1LaAi7Gw5OlAzDDCM6LFgZi8hXK4YPQfs9WVrIStvFgjSNHV8r7jvTbq9TFHaaFhI6CHF4ZXKGrwt2tY8SyI5NicFAnzyYbe24sdAcJdJI3Qt9GoKdDDvtZpiMgGJtS6ebuM7l7n/zeyZJgZB+76i+XkgxGVVsK/SPZQjmNA2uWrEKNceo2qmWAW4REBQ8P3nXfVZhgWCQ5pxTMw99gFtZd3HeUiN6XcSKERKYJhC8XhDVByqi1agaNtfw9f49TQSaGq9mpy3DZOT6QambNG9GCaTMQFNCME1W1l1csDd07pYRpnkmoyIfQlebyQ4zA+93KiejjHSnshpkzt/TTNwT+nn7BEKKQ3oY4lDpNKdysQ8BUKHQg4WdAuVMcaYPJmufPc2qHKTzNYSs37CnWU1oCOqanL1Yxwnj+QWAC/c2cXx9I7ftUekKLWnUNvI70lhmH8LWJ805p16bMxhVC6iUyajtJQbltESgjYE0BCfV5nqR0abVjhvtr1hMEXBZdNzsAQQIZ1j6gKMn4aQNwH4gsNrxcNFpOxLtziPPqVxoMtLbM4aw0ygMtxI7lTtG6OhF+5qZ9vAsH8Kw7TPPYQoExwg77bh+VA3VJNR2yiamaYJsQJNRlmYaTpqqWGr0+8fM82b5QVqG5qw0hAv3NbHaibXi6Hk5bQe6XpArqMwoIyAUCEXO9FU2GW190pxz8QAaftbzA5zoeLhoXzjLLaUhGCYj06kshIDnB1go60NoVtHx/P7B6LRkm/Qs5ej31AcQCDkzSiDUEHSnsj44p5kSTnQ8CAFcLAVXmWzlPD+GCt8t+i4wuElGCJGamObYFJkP6o4d+TjUIHPB3h1Y0Uow6OdN8yGoBfJGCT1dNnxUioplJZzKaj3oNAbxIQxsMpJ9E4jsFfZW2mEwx2IjOyw6L2sd6Delqj69cG8zkZyorkl8H2ZrCWYeAtDfr2no7WCT0RYlbaZlxv0rh+zZi/VQvS0xsK4Yjup6tEB77IQVAqXyEADd0S0ji9qGhmCUsWgZPoTwd5Sz3+dVgjQ1hI70Oejn1lHt3b+3AduiUsK02/NRzfBj1B0bvUBkLtKjHsSdtcrAUUZq8NmxUEHb9bHR8/v6Q9cQVtsudtcdnLZzITTbdZWTX10Dpy/Kpuv62FVzUHOyzSRl0Bcu0qnYFEXJAOM1GVkU5qGY+TAmPT/AqW4vSvTMOsdq28Niw8FSsz9gQj8vkG0yMrULXRMAEJmNVJvV9rx7o+cHqFiUWNo21BDyZ/2rbS/SatlktEVJEwiObWHnQhz3r2ZESzsWsNSoFmc6+gFOdnvpTmVp41aDapk8BCA2Y5kDvxIIURZmpCFoJqN6+N0TZTWEHB9CX2KaF+D0XTUQpZuMVDv3NBfQajilNIQNL+irvaMoWqRGbT9rd31gDSH67mJYQHC17fU5ZdWsWgiBZZnvYcbDq9985u56aqZyo2rLgWy0KKM9O1IEQkrYaWaU0UAmo1DTWEyp82WiJlBn7g77MetaKaHWyhEIeXkIQFi+Qs8ZWlkPAxjO3B0WGoyfCxc1x8LZi+H2PA2hF4goB0FRsQh+gVN5ed2NzrvBAmFrkmWLXWzGN756cFsNJ7RZFj0Qqo5RjlNZDaqlNQTDr6EG1gv2KWet1GaMonqAZgIrEXpazoeQjGJpLlSwq5Y+UOj9G6r3+Q69qA0Zg1ijWonOm4YyYZ25WOtLLiw8rzymPpiYcfz1qo1AhI5HZQPvE9ZtD82qjd11J7WWUd2x0WpWR1oAKVtDiDU4ZQLLzENwyiemqYlCWviziRKIavBN0xCEEFhte2H/NZxCp3KmyahRheeLaO2P41Izbxn3vMo1Mp+jNDw/iArbKcr4EFbabqEQnDVYIBhkCYTwxk9WOG01qlhqZt+8imhQ1n0I1eTiKu6gGkKjfxa6q1aJzQ+GKUnPQ1DfLedDCHJNRo5NfcXtao4VPtRpGoImoNKS69Lo9vKSqWQ/Zpgh9EHd8wXWByhFrI6pNITldbfPBr+gLZKjBpk9KRrCYqOKRkqIrBJ2SyUmFnmElU77BYKjhZ3qi/ukMZAPQQrGVtMpNMWp++DMxezBcd314fpBqGHJ5WLTBtxCgWAkpylT7WJ0z8fb1T0I5Ic/9/x+DcE2NK80VtouzpACgZ3KW5Qs59xiI57BqUFN3WhFD7IuQBTqoVQzWDWoDpKpDGhlINoe9uwIC5vpdtTVtodqxUoMAjXHQrVilfYhDBp2GpoS0meOanBoNRzsaS6UMhkVxc4D2TMw9SCeLQejQez0XUNDWG174Qy7mtQQ1L4qf8WMSguLzlVRS6m7FGkIJUyPeWTV/gmdn2EfxEUCM8xvJepCKZSQLqUhyM/PkhpCmjYX+brkM5WVAR07lTNMRmqA1+oU6RqCXhF4qVnFrlol9GXl/IZeEPRVhy3SEFQp7307F+DYxBrCViVL9dZnvCvajL+MD2ElxWSkBI56OJTZpayGsLvugCgUBICc8cibvtWsJmdIjWrCIUZEWKw75XwIrp85gKj2bvT6nZatjIqqyqZbd+zSGkLHC7CQ4whV503/br/Zpyzxd2vRdzuen/Bn1LXrqAZl5fA/rgnrxYbTVxkVCH0IkYYwpEDouGENJTNLGZBOZXlvFZUhCSOmglK1uYoEv06/yah/tpx4pjLKuwNDaAhy4FcObWUyWpHPORGhlRPmCoTPpmP1+xDy8hDarg+3F6DVqCay2WcdFggGWc45fQa3su5GSyi2mtnqrULXKBSObSVmDoP6EGwrHNT1QnHq+Lr5YUUORiZlS2AXhZ1WDQ2hK2fQWTNH1U4iknbz4qJuXWMQ1jE1rbTvAvHsdJBBNzYZxRrChqExqYFped3FhhwA6tVwzQnTbJFmklHaz2LDSV3cpQxZdYyAcFEXdW92CgbTqC9L1JdSi9QsNcP7KM83E5mMpFBuu/0lTeLij06sYaVcq47no6KtG21iagjLa+HAX5EFIVc1s6/ad0+BMO75g2sI0TPfCDVDNhltUbJU71ajKuuiBFhe96JZ4FLDKSxMlmYyApJZtoNGGQFIOLT1PAddeK220zWexXq1lMnINJGYVCuWsWJagFrFlia2dB+CsufuiUwr+YIpz2yVtq5xov1uGB55+q7QnFa2qJ9+TDW7X1l3+wSkev3AahdA7KsJZ/yxRqkEhWmS0TWEQdunUMlXWRqCqnaat64EMNiaCMo53WpU0QtiJ24aq20X1YoVJUemXSu9PHxk3knpi26BT0uP8HJ74bK1aluY3+D2Rf21jCRPEy81ysjKFwjrcXRfzbFYQ9iqmJVOFcrcs9r2wkG2qW6m9JK7OqttN1yD17iR9dLJG5GGUH4JCWWuEkL0awiaymyuogWEuQhFg48v4/uLEtM802RUDZ3Kaxu9vvyAlbYXDQxZ5YpNQg1hOB+CGrjKnss8rzqHMhmaUTrq9YMnOgCg3RcOltc3wrW05eCjJgD6bLrjhWGnanAapp5RnoZQsWKTUbfAZFRkftNRmk3U7px+Vc9UPSciTJ9R5x0z9F1kD1uNqo0FWfrC1MxVuZdVTRsBwhDoIg3BjDKyCkxG+rnZZLSFydIQ9CgFvVBdmbC15fXkegSKupYINLSGIGetG70gHoy0tWXDZJ80DcEpLLdcVDdGtdcz49wdG4vNuL90lO0WQF80TnY70tfQBYrXAlYD166aA4vKrwGsH7NetaOggo4bGAIhbNf9q6FAiIXyApaNwUf1o5nZrQusYRzLKnQ31alsa07lguup+0OKUAl60ep9eetfaD4UIP1arayHqxTuqju5z1TXzU+UJCIsyXpGZsSguobm+uaFGkJKlFGlIDFNP8cg+R3ThgWCRp5zTldjVchauL2chpB2TN2mPGiUESA1hLZ240dCKo6oWO2kC6MyPoQiJyQQagiuH0AIASEEur3YqQz0DxQJYbqjnEDILdlczR/E1IzesqhURJj5XUBqCE0Hy2039CEYeQgA8GBkMopNifosdbFRjZzzakAUQqA9Dg1hPbn4kk6YqSwS5y3yIaQ5fU2U6WyxRLuV2TLW5tKcyh521x3YFqFRtcMM6DSBkBOCrFiSQRV9AkGuFNj/vCzk+gG9DB9CroagmYlrlfLhvNOGBYJGnuqtP7ChCpw0e+TNPFWIm0mtGpcbHlZDWFn3+m58JXzuWW6HhfrSNIRGWAspb+ZSNIAAcRy+5wts9AIIASw46QOcWkC9Fc2ilUDIT07LWlQdKBN2GguTvBo5abS13698Ih0vabKoGyYj/Rosr7vRwLAkZ4p6Wz1fwA9EFHEF9C+LWobl9Q3YFmFXrV/wO1aKhjCkP0YnciqXMRmpsFsnztkwUXWMAMiAg/T8nk7GmtA6WRpCS2ZV6yGuQLEfsJcRZZTnQ1hue9G67OxU3qLkOefUtqOnNkKHVJ+GkL9kZJrZpu5YfXkIZaOMgFATcP0A9y0n7dfqwTp0dA0AUs+twvDyVrsqYzJSPg/XDxL7mwUBgbiwnRKmqu+KchE2chyJanvWSlx6pnWZEGGdriw8t1AJfSLHTm2gJwdw8/z3r3ajAQAIJxVrGz08dDLUHELnYrKtnYTAURrVMBpCGHKcVuupYhMCEZZ1L7qeg/gQVLBBLPiz7yNVo4iIpJk0LcrITUTD6YmgOnnF+RSqFlKayehkt4eja0kT25LM38mamIR5CP2JaXkCQdW1si1CrcJO5S1JvoYQ3qx3H1uX78N9ohDDnAd5pe1Fg6BOwmQ0jIYg23CnHPhNDeHQ0XW5Pd1kpNqWRZHNGYhNXJ62xGg9Y6AwV/Vy5IImeSYjPxBw/exBwLYI1Zy1gPWooMWMQSaLjuuj4dhRrLqKpEkLOz22thENAPpvjK9Btc9Gr9rcqFZQc2w0qvZQuQhZfi8gnmD0AlEcZVTgj1EEWrDBzloFFmVrCIFcqlJ/XtKjjLy+gpLZTuV8gdBqVLG8FgsEVd13MXqG24n3SwWTOs8XfRO1IoGgB6fowSOzDgsEjTznXN2xUa1YuEsJBG2fvKQiVSo7VUOojuhDaBoCwXB03ykHo3SncrGpS3eqZuFIARZqCHFphDSTkZ58pDAXNDEpKmamzpeVhxAmtYXfLVNmJPnd2Nykt3khoSHE7dKj05YMYa3CTtVxgTgeX5XfGDZbOSuZEtAXhA8i8+Sw/hhFV1vS1LIoN1v5ZDdcqlI9L6GGkJ6Ypt+nWaU88taEVuxphsL74ZNdLDacyP6v+ujQsTXsWKhEZsgi02UvCFIFQp4PQWlFAFITEmcVFggaec65cJbo4JAUCPrDn2XvBOLMyFQfgjZzGDRTOTxvPOjYFmFnLQzrUzfioWPxYGRSpsBdUSITEM9A3V6gmUDCpL2FipVwXKeWFm9WE8scZrUhTyjl1eDRncBq4Cpb4E4PMU0rTAiEAlxZasxJAhAKZRVy3DDi/GMNrBJ9Z1iTUVoyJRAv++j58doOCwVJflnmN4Uu+IH8AAW9VAmAzJj8FcPP1mqmH3Ojl5+HEH43fi7051SVfT90dD1xPYv8Nz1fRIJVUehD0LS22gBFA6cNCwSNPOccEA4o9y4n1U21PSt6Ja3aqKKuCQRXzroG8iGoWeiR9YQNeaFiY8dCBfceb8v2peQhSDU6r3xFGR/CgqYhmALEnPGm9UVRyYa8RdUVjWr2WsCJWX6zio1euaVD1bl1c5NC7w9lFweSglcNBncdW9MGhqSGYOYFhImGwziVs01GkYbgB+i4PdSlCSyNIge9omNobXnX0NQK00xGHddH1wv6nqlVbZEhRV7WukKZfO84stZnhgKAwyvtvgkdkK0heH6ahmAVaAixxsMmoy2KyhfIWlC+1ahGs4Iy9k51zHCfojyE/lWZilADa8fz+7SAVtNBLxCwCKkCLtYQRpudO9EMNOgTIGZUT1pfFJuMshdVV+Qt7KI7lfMyYDO/m2IyMs1XNeP4QHxt9IHO9CGombjSHFSo6iCoJUnTkikBRM7QXiAK17ZQFXiLBi8z+iyvnlGcgRz2QcOp9F2rNFNiq1FNLDIUnbvgNwB6bpCXmHyo66ObsIB4ApWpIWSth5CTh6BHFtYGqBE1bVggaCyvb2TOtICk2cCcDZadIemo2ZIQYiinsqrUqNqgowaI3fV0AbdjIfxuXi6CsvWWcSq7vaB/xttI1vhfaceF7aJ2NvPNOGW0lLzEn4TZp0SIZPZ3001GQDww6tdgURYf1LebpSFMJ69elLAsq203jNzKcipbscDuuPkROlXbgm1RoVPZLDDXkiUh0ohrFMnBMUVDSCvtEs/ak8ctWsFPPxdgLIVb7/fx6N/J9CH4IrGeMgDYdrbJyNR4VHs3Mlb1myVYIGisrHuZzjkgNhvUHTsxS9Ezg/uPmW0yqjl2tMZsFHZqlb8kKvoFSFm/wQiLTfvuYt0p5UPIDTutaAOOOVA0kwOFXthOsdQMta6TnfRaOEXVLVX7yuQhmGWpi+hoIY66ycis7WQeHwht97uj6Jb4vgmPa/gQIg0hdIZmLQeahhm5ZRJpCL6I1qrIIgoLLXIqeymCLKPAnZoQxH1g9Qmc1Xa/5pgWlCCEkLWyivMQouNor1VElLldvc8y14UmIyPslLIFgl6GA4jLjW8FsxELBI3j6xuZzjkgvmFN88+SVufIZNm4OXSiEthuAFfedFnmqizUuc0bXJ0vrdKpYncjvwR2NBhXc8pfSw1hQ3Mq64Oo3ier7f58DNXfWUsYFi2qrj5Lm9X2/ACeLyL/Q5ms8sS5NXOTro2Z/ozIlp4y6wy3J2eKkUAw+6tEkqOJ+i17mgupn1eisNOglLmlVkogJKOVWo1qGFSQ8r3ldVf65ULHeZrASasGHJfyiO8fNcPOK7YIhPedmnPoGoLKVjfPpfbLjjJKNxll+RBWDCE4SH7HtGGBoJHnnAPiASVtdgGkzzxX215oJkm5ifUZo9frd1yVIS6VkN6mvN8Tagj5YadqIfUsqpXwQfF8kWpKWO3EM8ewf01huhB9ltqGEUxGXbVCWBTWmS24s86trpuKMlPnS5zf6dcQgPiaqO0LFQtEcanutiEQonj4IQRCtslI8yEULIcKhH2VFcKriDTBSnyd9bborMh1ppVWmOZUNgfQ8Jj9fVEmwACIS8Prx1FkbQ/9gNkaQp/JKCfKaMWIVjTXPpllWCBIipxzADLNM3np+3lCpq4to+n6wUD+g+jczXQhpWdnZpFVolqhkrqyolIAoGqHN7vXC1Jnjn4gcLIbmoPCwaF/ZgZkZyuXy0NID+szZ+DRokIlNYS2EfNumn4UaX4GoN9sZ5pkTJNRVChugPIVyymza51IQ/Dz11NWlDEZxe1Oxven3UumVlhLyRnRS0Ur0tZEKKMtmt9fMjT+KBnNmJjs2VHN1FJ7vugzGVXsYg2hZWiGWyEXYWICgYjeT0RHiOimSZ1jnBQ554D4gTUH2XwNITtpSM8MTQttK0OsCRiDkZqd5piMigrclTExOJW4dEU8c7Tk8ZMmkLQEqqKy1KP4EMzvVmwLu2pOaZOMuXSnHkuvU8vQEPakCOu6Y0eagV4aAyhXOdckLreSfp2VqUNFgRVqCAP4EFRiV979r8pqKNLWlV5pu9hZqyTu/6YscKdHhJmCKA/V9+YELzYZJU1srUYVXS9INT1mFbfLihoyq6nm1XCaNSapIXwAwIsmePyxYl7ENKKbqc+HoAa1FB/Cevp6BIA2c+iF5asHyVKOzt1IzkLjNqULL53FejW3BHZRqWFAK10hBYJjk5YZGod59vwAJ7tepiYzikDIWgs4LWy21XBKxfoLocI042ti2oQVpsknOlezXyjrNvq2VhpD//4g5SuW224i69akYpUPO43aVzLKyAzJTWv3qqEV1h0bvUAkAjDMpDRABj0YYbhlTUZ6m8zjZmoI6j40hJoQAr2gv7idTTkagtJ46kaU0RYQCJVJHVgI8XUi2j+p44+b42v5zjkg24eQV5hspe1FSzCa6Ms/er4YymSU5StoZQgKs91rGz188uB9sIiwtKOK5zzytOjzouUzgTgPQTmV00wsn7/xQdx4/4lEYTuFyuA9npGtXLSoOhCvBewHInL8AukrhLWaVfzw4VP49HWH+46zf28DTz5/CUBcibRRjR+R0Bben+lbd+ww36NuTBRSroHu7zAHaNVf37rzWGG/K248fCJXq60kwk7L+BDsPtPnHUdOYVfdwWk7w7Wl+6KM5DX92g+OoucnB8mHTnbxhHMXo/d6gb/d9bBtyxnFH83yFWUmBwoVrFD2uYgWulpzo7WfAUSDfr+GkL1i2krbxa5aJfpO2tKk9y23cc1dy9H7BcfCS3/0rMLfNWkmJhDKQkRXArgSAM4777yptSPWELJNLKfvWsCuWgWXnLYzsX2hYqOZUZhsJc9kpNW2cXv+UBrCJaftQM2xcE6rkdh+/p4mqraFS07fkfnd8/eE3/ntT90QbfvGm56Dc5caUbuKZpQLWtjpRi854JzbqqNiEd779UPxOfc2+46RFwNepnyGnvDVXIhv6bTvXrCnic9873785ievTz3OLb/3QhBR6ncfcfpOnLfU6POpnLengYv27UgIIwC45PQdqDt2YoBpaNqMnvgGhDkoZy/W8bkbH8Lnbnwo8/eaPOsR+zI/c7Sw07I+hAeMmewvffAgLr1gCX/yssfLdidt+YuNKhYbDj7z3fvxme/e33fM8/fG96ZeLyla+L7tFa5jDpS7FxSXnLYTZy/Wo6Q/xcWn7cDeHQvRuRVKYzD9CCqyyZwE6EuTmiyvJ9c/iRdxivd/+z/dgi/f+nD0fu+OBRYIACCEeB+A9wHAgQMHppbKl1fHSNGoVnDt716ROnCnZWv25IplmU5lPcrIF5E9fhCeeclefP+tL+h7SM7YXcMNV/Vv1/nJJ5yNp+xfQs8XuPqu4/jtT92AI6c2YoFQwmRk1jLSB7jTdtVwzVuuwJp0KlcrFs7YXes7xlIzuwpp18tfVB1IClZdIKTNKP/4ZT+KN17xiL5jfPzgvXjXV+7E2kYPO2tOakLca55+Af7LZfv7vvu6Z1+E1z7zwr7tz37EPnzvrc/vW3JTDzs1Z+xf+K/Pyq3tlMZpu7K1Wj3sNG9taoXpQxBC4IETXTx8Mh4ou3LyogSgbRG+/qbnYDXFZEoEnNOKBWJaRdXldReXnNY/cWk1Hdz+0Kno/UYJbVHxqqftx8ufel6f8P65A+fip590dt+MX/kUzGc4KrdiTOryoozMQJK0leiOnuri0guW8GdSyA6QfjRRpi4QZoW8jGKdLFttWjGuvMJ2QHL5R3dIHwIRFa4VkIcyZ6nw02XDZpvngwD6E9NM++5Ss5orZIG46FwaRYuqA9nLaKYN6o5t4bw9SW0KCDUqIJyt7qw5WqXX+JpYFqGakieStT3t2tQdO7ov2ikmuR0LFexYGN9jqXwIXS/MySgyGYX+mHgmuy7vTd0R33H9vlIiu2pOZg0wnbR6SSsZKwqaayIMYjKyLMKC1b9f1vasEtjKH2C2Ly8P4fi6i7MX44lPmlN5pe3hSec1U+/FaTIjcmn6rKy7aFTtUjdbGov1al/0SlGMuB6f7A4ZZTQu0so6lPEh6KUrOl5QmDSURl6Vz7JmDqA/zrtMLSaFmRlbJv9hGPRS3V23eMY+KuqeUlpamSgjvR+VtmIOzMP2i1kCvOv5aLt+6jOy1EwWuJvUNQGAXfUw8dA0XS4bIaQK2yIIufCQyYqhISyk3J8r6+lCcNpMMuz0owC+DeCRRHSYiH5xUucaB8s5tv4ypIVwFgmEpA9huDyEcbGUEmVRKuw0WjFNyIFiuOS6rNj7jYJyC0BcHK6/imZxLaa4DcnCd4PYqwdBT8wqI3BHRYWdqiJxRQJbmYxUMqGyqa8azt1h+0UPpAiPq8pWpGsIeoG7QfIQBkWVgenXENKjD21pivKNch1CCGkyis146jcrf4TbC3BqozfSeDMpJhll9AuTOvYkWMkJDy1DmtmjUEOoqEXXA3j+cCajcdGQ6xckNAS32FxDRKjaVlTcLm21uSJaMtppo9e/dnKZFbKyVvoaZEZpLvnZTYlQGgfJsNMeGtXJmgyUyehUWQ2hasMPhIx6o+gePtntoSfj8UcRZGqCoXIx4sJ2/c+eHpK82KiWSlIchbQFlMy6RApbClo/ENC7Ym2jB9cPEs+BYxMsiu9PZZ7dVhrCVmM5JYt2EBYbDk50vIQKWSQQKraFqnzApq0hEFFf1daypgHHpiiscZjZm3ow0pLkroIPnwAAGalJREFUyoRK1jI0hCiBqsQAsmiYzAYxNw1CwmRUwj8yKsp5WlYgmHV39AxylbPSHdI0qJ9fHT8v/8dMeJuU1qZIq1q8su6mhhTr+R06ac+88iWp+3ElpyT+tGGBIFlNSY4ZhEWp3qoHD9DUzRxBo1aQGjZTeZzoWk5aYlYW1YoVOZWHmTnmZeiGg2Z+G7J8CGYmcB4qiUg5fCfmQ6hahoYwYR9CpCF40flz22f0pT5A6ua0ooqjWZgCJ62wnSKukRSed2OA6zkMact2LrfdxFrZCluGBflG3oUSoGbJjLpjR3kIWVrHLMACQZK3Lm0ZWinJacfXXexcqOTO/FWi0rQ1BCB8INQNrRKzymkIymQ03MwxL9N1EKdy2kpcRbWYFBXbws5aJdJSlElj7D4ELVO3jI9mVEwNobAvVX0tt18gKHPaxgjtblSTAsdcQEfHDHTo9gJZIHDw8OwyhD4EU0Poz64HYg3B9CEsRwmuye+EGeCBPCabjGYazw9wqjuakydtjeKVtts3UzBRTjzXF9PXELQFWgZRz6sVC27JWjnp580u6lbGgRkvPJNMFBpUY9E1JLM8w7hQv2V9o4eul79gzThQTuVTG1JDGFC46hnkqwkNYcQoo0jgZDuVTc2xjPlwFPbIqCY9v2Cl7abO5K3IZJS857LMxDXHijSE5ZIh7tOABQJ0tXV4m95io9/sUUbrULVj3J4/MVW4LEsNJ7qhBxkQlVO5UyIiKP286XVkgHKLqmfVmy+jXei0tCU/1YA1bpOO6k91nomHnVpG2GmJWkaA7kPYiLJ6dVv+sO1WgkQ3GZmF7RSNqo1qxYrui1Gim8rQasqoJm1SZ2YdKyINIcgwGaVoCGZkVd5aJdOCBQK00LcRVLi4BHBSIBRF3agwRC+lxO5ms9RciFZ+S6sDlEW1YqHt+qVNTCZRVdQ0k5FbHMoaRxklV10rk5lrtmPVdGAOORPOQgkYFe++aRrCAHkIQBxltbzu4qJ9cdIeMJoz3LIIC5V41bS0wnYKIsJSI6m1TlIgqHYcT/hN3FTnr/IpmLWbltc3UHOsRA0sQAqEXtyno+Q8TRIWCNDUvFFMRlpdFv24RUJGJQINux7CONFXfhvEqerYVjSrGuYmr1Ys7FyopGoIZcJOHZvCtYD7nMrFDmkdfW3gjhdqbIOuYFeE6k9lKpm4D8FSeQjlNAQzcez4movz9zRhWxSb09zhNEH9HOr4WYXtFK1mnBsw6ags00QlhCj2IaRoCGkFMsPnPIiOP4vmIoAFAgDdsTX8RdolF19RpoA4QaVYIKxt+DKeefo+BCB8SKMZchmTUcWKkoeGHeAWU2LAgXJmgmjhGdOHMKDNebFRjerxTCqLWP0W9VsnrSEQhcIyijIa0IegtNzFemxO6/ZGs+XXHTupIeSYTnQhvdEbTRAVoQZp5Tdpu+FELW2iaOeEnaY98zVtLWkzk3mWYIGA4lWnymBbhN11J2Fy2OgFhcesOXY0u566hqBF+wySmOXYFMWoD2tiWWr0L3IeLapeog1pawEPamJYbDg4tdEbKYS2CHVMZZaYdNgpEM5mo/WIS/pj2q6Pjuuj4/lY2lGVmfguPLVO9SgCQdMQsmbgipZW1qTjDu/MLoOquKrOt5wTDaQEQmBGGWUKBM1k1PZm0n8AsEAAMD4njz6LUrOMIjNUzbGj2fU0M5WBZCLQ4CajciaJvHObtaDiQay4X+pVKzUPYdAoI0CZzCYTARQ7lTdHQwDiekaOnV81FkiWp1ZlK/Y0q2jJ5VbTCgYOil4vKSuKR6H7ELq9yYbpmuHPefkClQwfwvG1dL9hzbGjaq2j5jxNEhYICG+A5hicPLpTMi/hRqdetSKH37Q1hD0pJqMyD+BCxRpIgKSxlBIDPsjgo5shFINGw6gJwYmOO3TWdRGRhiAnDJP2IQCxY7nM79HLgMQhlAtRefeontAI7Vah1qqwXZGGsNrx4AdiYJ/QoKjFmtTvjjWENKeyTEwbxGSkmeHYhzDDrBQ4tsrS0grcHc9RN3XM0szTRC/fMEiUkd7uoavFGouhAIPlQqStBTyoDyGueOoNHKFUlsiH0N5EgWAlV+7KQw871UMoVfHGeBnL0ZzKbdcvVXJ+qeFAiLBsxqRNRoCsvCt/d17hvUpKHoJuYjNRWtE4cp4mCQsEhD6Ecahw+iI56qYqDDvVHtJpm4xUtM/xdXeg2vO6ZjPsDG6p6WDd9bGhLTNYZvnM+Lz9AmHQuHU9M3ZiPoRq0oewOSYjSpw7D9siVKXGp2fdKufuOBL2VO5NrIHkOJU1rXWj54+kmZRBz9bPq0Vmp0QZKV9klsmo4/ljyXmaJCwQEM4Ix5FGvthwcELOKvIcUjr6De5M2WQExE68QUxGuoYwig8BSIbtDqKlNKp2ig9hsDDFuOKph/aETUZRlNEmmozKCh9VgC8aEHdUsdioYqMXRD6yUWbqarYc++5yTEZawmfXCyauIejZ6ittWdguZeGftLBTJUCXUsJOa44NIYAjcuW5cVgkJsH0R6AZYGU9P/StLK1GNYpSWV53UbEIu2r5FcZnSUMAVNy3G6+bW0JI6RrCKD4EwFixraeqlZYwGVWTPoSeH8D1B3MM6071iZuMNlFDUCajsgJOmd+Or7uo2qHWqAbmB090wn1G6JuGjDIqqgasf6a0tkn6EIBwdq/7EBYb1dRcFCtFICgnfFaUEQA8eKKbuc8sMP0RaAYYlw9Bn2GqpLSiQlwJgTDEmsrjZkmaBjpeuG6uufZsGrogG3aAM8tPA9qSiSVmhabJqCsjlMpUa1U0qzYqFmGl7ZXKkB4GZZI5tVEuc3gcqNlsaQ2haqPjBVhe38CSvIdVtq4a0EYZmJXJqIwPQQnpI6c2hs6EH4SWJhBW215mieq08tfLOWZi1V8PrIYClcNOZxS1etG4fAhAGKVSpmwFkJxpVe3pp7IvNRewsu5J+3u520PXEMrM5tPP21/PaBB7tbn04yDmJgURRZFiHc/vKz8wLlSbygrcUVHnKDurVwP28bXYt7ZbDmAPRQJhtDyErhdEA2je4LhkaCaTLvew1Kyi7YYRUHmJpak+hBwzcZ01hK3BOFcv0pdhLJuert/g065lBITOruPrG+EMueQAord72BlcVPE0sXbvAHkIjh2VrA6/O1z56pamIU1q8FHJaJvhPwA0p3LJ36P8McfX3ShZS93LaoY7ah6C6wc4traRWdgu2rdqo+ZYeGBVCqJNcCoD4eCe9wwrM1wvYTJy4djpZuLYZBT2H0cZzSiq5PI4FqtYrMdmj+PrxaWvgeSANe08BCAUjF0vwHLbLf3QK83GtmhooWbWvgcGMxmp7Fe1FvCw0TCqJr7bm1xpanXczTAXAbF5Y1Afgj5DVtfnoZNj0BDkdx9Y7ZaaKS81qrGGMOFnRE9Oy8sXiDWEOOx0eS3cP81MrCY1D652UXdms7AdwAIhd03XQdHXRAgd1SVMRjOUhwDEgvGB1U7pm9apxDPQYRcvcaTzUncqDxLppKI4VHbzsNVKFxtOpNYP4n8YBNWvm1G2AtBNRuV+jx4WqgZsdW+rmfoowkzN8h9Y7ZSaKS82qprvYrJ9pjSi5XU39CEUmozibcdzTEyq3Q+c6MysuQhggRBlFo/FZKTK5665WO14pS68/mBNez0EIP4ND6x2Ss+ulVN51AgQs3xFZDIqoyE4cckFQPMhDKEhPHxy9EEvD9WmzZolDmoyqldtnOh4WNvoRX6wmmOj7tg4trYRvR+Whvzu/avlBselZhUPjkEQlUEJqPtW2mFhu4x8ATslMW15fSMSKCaqvx460Z1ZhzLAAmEshe0UKkrl7mPrEKLcMfVZ2yxoCGoAWGl75U1GlcHCGrNoNZMF7iKTUYmZrVm2edgF2RebDjxZn2ZSA3ZkMtosDWGATOVwPysyDekx9boWPcrkRf3uU91eqcGx1azC9csV5xsV9czeeWQ9PHemDyHdqZyWgwDEk5peIFhDmGXy1nQdFBWlcujYGoByWscs+hAUZQcQJchGnb0tNZw+HwJRufwMvQaP+u4wbdIHgEkN2JtuMlI+hJLnqzt2NNDpg9fuhtIWRlsnQr8mZcyqeo7QpPMQdtcdWATceTR8houijEynclZkoX4vzapDGWCBgOV1D82qjYUxZUC2Gg4OHQ1nF6XCTmfUhwCUH0CqA4Y1ZmEuct6Va/eW8UuYSz8O4n/QUQsdAfNjMho0U1m/7roJRGkIo7Zb/36ZSZO+z6T7zLbCSZ0SCFn5SapPleB0e2GNomwfQvxsj8NfOSmmPwJNmdV28apmg7DYcOLCdmWcytXZ0hB2yRkSMICGoExGIwpV04cwSLVSvWxz+F8mpg1qMmoMriENikp423Sn8gBRRgp9gFP386j9MuhseWkTBQIQDtj3y/DaTA2BkgKhqLpxfUAhOC2mPwJNmXEVtlPoA0qWg0lHH0RnoXSFmiEB5R/8yKk84gC31KxiXSYFAbIWUUkhGZuMZJSRq3wIg/WpPnubVMz7ZoedOipTeQCTkULXchfHpCHogrBMkTf9mZq0yQgA9jQXoNa9yTJpmYlpeVnKgKEVsclodlkZc21yfUAp45dQi44Ds6EhAPEsp3SUURR2OmKUkbZADTBYtVJz6cdhncrD+FAGRQmaTXMq24M5/VW7bIsShd1akQ9hRA1hwMFxaRO0Nh2VJGlbhJ0ZtcjMxLSiuky6E541hBlmUhrCzoVKab9EvWrDonjWMW2WBnzwVWLayE5l+SCqh2sggWBEGSmH9KDRMLoQn5fEtEHDTlWfLzWThd1U34wq+GsZJqks9AVqNsNkFOVe1J1M57ltJxPTlJk4yypARJF2wz6EGWZ1fbzrm6pjDTILqDv2zGgHQPwAlo8yGiwTNovFSENQAqH8ClmRD8GN8xCGSZRT2ebA5Gz8UemKTctUHszpr9plmj8Wx6UhaO0oU1Rys30IUXZ2zjNsFrdblvkZeRqPajubjGaUqLDdWE1G4bEG0Tpqjj0TEUaKyGQ0YHG7UR9Ws8DdICWo1SDWdsMqosMucFOtWNixEJoJJu5D2DST0eDrIQD993CrMdhEIQvdL1QqD0E+U1Xb2hQtOnqGc8aFyIfgxyYjonwBp/yFnIcwo4wzS1mhwhYHFQizkKWsiKJJShe3G1/YKRDnhnS88ksmxj6EUIUfdHEcncUxDXxZ1DZbIAxay6iaPnCNS0Oo2BaqtoVdBYXtFGqt44VNcCgDsdknbS1lRRRlJL3Px6UvMk9gqX7dthoCEb2IiG4nojuI6M2TPNcwjDNLWbE4hIZQd6yZ1BDKPvhKmI06gKqBeHl9cKeyaoPuQxg2ImWx4aBi0cSuiRoYNr+W0WA+BNNkNK48BNWWQZ6RVqO6aXkbZbR8yyIQJaOMin7PQsVCzbE2bSIwDBMbhYjIBvAuAC8G8BgAv0BEj5nU+YZBVTodpw9BzSoGEgjV2fIhxCajwTSEUUMCHdvCzloliukeZJZvWaHTrqtFGQ29nGejOlH7/tTCTgc2GSXLMMQawuj3at2xB1qUqtV0Nq2/4oJ++e2rWBT5EPIK2ynqVXumtQMAmMwKICGXArhDCHEIAIjoYwB+EsAt4z7Rf/iLb/atp1uGdblq1TgvknJKDnLM+oz5EFoDagjj8iEA4cP46e8exr/fcQwPnewONPjUHRsfv/Y+fOW2Izi80sFjz9o1VBsWG9WhF/opgxrYNq+43WACOzIZGREzu+sOiManIQwSbdNqVOH2guIdx0Bc8ju/fbZF+Og19+LLtzyMe5bbeN6jTsvdv1YZTAhOg0kKhLMB3Ke9Pwzgx8ydiOhKAFcCwHnnnTfUiS7a14yKXw3KcxpVXHzajqG+m8bpuxbwxisuwUt/9MzS33nl0/YnSjZMmx+7YAmvfeYFOLC/VWr/M3fX8PrnXYLnP+b0kc/9umdfhK//8CgA4BGn78RPPfHs0t/9tedcjO/euwIAuOT0HXjJ48pfA51XPPV8PP2iPUN9twxPOr+FK591IS69YGli59B58ePORMW2Sq8Ad/5SA796+UV4oXE9bYvwuz/+GDz1wtHb/frnXYx9O2ql9/+lZ16YyGKfJGcv1vH6512CF/9I/v3zG8+9BDc/cAJAeL/9wqX549cvPfOCqDz7rEJqQZGxH5joZwG8UAjxS/L9KwBcKoT4jazvHDhwQBw8eHAi7WEYhplHiOg6IcSBcRxrknaKwwDO1d6fA+CBCZ6PYRiGGYFJCoRrAVxCRBcQURXAzwP4xwmej2EYhhmBifkQhBA9Ivp1AF8AYAN4vxDi5kmdj2EYhhmNSTqVIYT4HIDPTfIcDMMwzHiYnVhHhmEYZqqwQGAYhmEAsEBgGIZhJCwQGIZhGAATTEwbBiI6CuCeIb++F8CxMTZn0my19gLc5s1gq7UX4DZvFlltPl8IsW8cJ5gpgTAKRHRwXNl6m8FWay/Abd4Mtlp7AW7zZrEZbWaTEcMwDAOABQLDMAwjmSeB8L5pN2BAtlp7AW7zZrDV2gtwmzeLibd5bnwIDMMwzGjMk4bAMAzDjAALBIZhGAbAHAgEInoREd1ORHcQ0ZuncP5ziegrRHQrEd1MRG+Q25eI6EtE9EP5vyW3ExH9T9neG4joSdqxXin3/yERvVLb/mQiulF+538SEY2h3TYRfY+IPivfX0BEV8tzf1yWLAcRLcj3d8jP92vH+B25/XYieqG2fezXhIgWiehTRHSb7OvLZrmPiei/yvvhJiL6KBHVZq2Piej9RHSEiG7Stk28T7POMUKb/1TeFzcQ0d8T0eKw/TfMNRqmzdpnv0VEgoj2yvfT7WchxJb9Q1hW+04AFwKoArgewGM2uQ1nAniSfL0TwA8APAbAnwB4s9z+ZgB/LF+/BMDnARCApwK4Wm5fAnBI/m/J1y352TUALpPf+TyAF4+h3f8NwN8B+Kx8/wkAPy9fvwfA6+TrXwXwHvn65wF8XL5+jOzvBQAXyOtgT+qaAPgggF+Sr6sAFme1jxEuH3sXgLrWt6+atT4G8CwATwJwk7Zt4n2adY4R2vwCABX5+o+1Ng/cf4Neo2HbLLefi3B5gHsA7J2Fft60gXMSf7ITvqC9/x0AvzPlNv1vAM8HcDuAM+W2MwHcLl+/F8AvaPvfLj//BQDv1ba/V247E8Bt2vbEfkO28RwA/wrguQA+K2+kY9pDFfWrvGEvk68rcj8y+1rtN4lrAmAXwgGWjO0z2ceI1xNfkn32WQAvnMU+BrAfycF14n2adY5h22x89lMAPpLWL0X9N8xzMEqbAXwKwOMB3I1YIEy1n7e6yUg9eIrDcttUkGrkEwFcDeB0IcSDACD/nyZ3y2pz3vbDKdtH4Z0A3gRArfi9B8CqEKKXco6oXfLzE3L/QX/HKFwI4CiAv6HQzPVXRNTEjPaxEOJ+AH8G4F4ADyLss+sw232s2Iw+zTrHOHgNwlnyMG0e5jkYCiL6CQD3CyGuNz6aaj9vdYGQZuedShwtEe0A8GkAbxRCnMzbNWWbGGL7UBDRSwEcEUJcV6JNeZ9tSnslFYQq97uFEE8EsI5QBc5i2n3cAvCTCM0UZwFoAnhxzjlmoY+LmPk2EtFbAPQAfERtymjDMG0e2+8hogaAtwB4a9rHA7ZtrP281QXCYYR2OMU5AB7Y7EYQkYNQGHxECPEZuflhIjpTfn4mgCNye1ab87afk7J9WJ4O4CeI6G4AH0NoNnongEUiUivo6eeI2iU/3w1geYjfMQqHARwWQlwt338KoYCY1T6+AsBdQoijQggPwGcAPA2z3ceKzejTrHMMjXSyvhTAy4W0kQzR5mMY/BoNw0UIJwvXy+fwHADfJaIzhmjzePt5WDvpLPwhnDkekp2rnEOP3eQ2EIAPAXinsf1PkXTo/Il8/eNIOo2ukduXENrJW/LvLgBL8rNr5b7KafSSMbX9csRO5U8i6Uz7Vfn615B0pn1Cvn4skg67QwiddRO5JgC+AeCR8vVVsn9nso8B/BiAmwE05PE+COA3ZrGP0e9DmHifZp1jhDa/CMAtAPYZ+w3cf4Neo2HbbHx2N2IfwlT7edMGzkn9IfTK/wBh1MBbpnD+ZyBU0W4A8H359xKE9sV/BfBD+V9dPALwLtneGwEc0I71GgB3yL9Xa9sPALhJfud/YQBnVkHbL0csEC5EGK1wh3woFuT2mnx/h/z8Qu37b5Ftuh1aVM4krgmAJwA4KPv5H+RDMbN9DODtAG6Tx/xbhIPSTPUxgI8i9HF4CGeav7gZfZp1jhHafAdC+7p6/t4zbP8Nc42GabPx+d2IBcJU+5lLVzAMwzAAtr4PgWEYhhkTLBAYhmEYACwQGIZhGAkLBIZhGAYACwSGYRhGUinehWHmDyJSIXkAcAYAH2F5DABoCyGeNpWGMcwU4bBTZttDRFcBWBNC/Nm028Iw04RNRgxjQERr8v/lRPQ1IvoEEf2AiP6IiF5ORNfI+vMXyf32EdGnieha+ff06f4ChhkOFggMk8/jAbwBwOMAvALAI4QQlwL4K4TlKADgzwH8f0KIpwD4GfkZw2w52IfAMPlcK2QJYSK6E8AX5fYbATxHvr4CwGMoXmRtFxHtFEKc2tSWMsyIsEBgmHw2tNeB9j5A/PxYCBdP6Wxmwxhm3LDJiGFG54sAfl29IaInTLEtDDM0LBAYZnReD+CAXBT9FgC/Mu0GMcwwcNgpwzAMA4A1BIZhGEbCAoFhGIYBwAKBYRiGkbBAYBiGYQCwQGAYhmEkLBAYhmEYACwQGIZhGMn/AUYIkDsDueSJAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# plotting \n",
    "plt.plot(times[::1000],states[::1000],'-')\n",
    "plt.title(f'state transition graph')\n",
    "plt.xlabel('Time')\n",
    "plt.ylabel('State')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 74,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.31596684, 0.31425686, 0.31722683, 0.01778982, 0.01747983,\n",
       "       0.01727983])"
      ]
     },
     "execution_count": 74,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "times, states, state_matrix = kmc(P,epochs=100000,init_state=0)\n",
    "probabilities(state_matrix)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Unit Test:\n",
    "\n",
    "We now check the correctness of our KMC simulation by comparing the simulated probabilities to the equilibrium probabilities predicted by the Boltzmann distribution.\n",
    "\n",
    "For a Hamiltonian $\\mathcal{H}$, the probability of a state $\\omega$ is given by:\n",
    "\n",
    "$$ p(\\omega) \\propto \\exp\\{-\\beta \\cdot \\mathcal{H}(\\omega)\\} $$\n",
    "\n",
    "Here, we have $E_A = \\beta^{-1} = .01$, $E_B = 4 \\cdot \\beta^{-1}$.\n",
    "\n",
    "Thus, in the absence of noisy microstates, the equilibrium probabilities should come out to:\n",
    "\n",
    "$$p(A) = \\frac{e^{-1}}{e^{-1}+ e^{-4}} \\approx 0.953$$, $$p(B) = \\frac{e^{-4}}{e^{-1}+ e^{-4}} \\approx 0.047$$.\n",
    "\n",
    "with a relative ratio of:\n",
    "\n",
    "$$ \\frac{p(A)}{p(B)} = \\frac{e^4}{e^1} = e^3 \\approx 20.085$$\n",
    "\n",
    "We now test our simulation with $\\beta = 10$ and 50 trials to see whether we indeed recover these equilibrium probabilities:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "num_trials = 50\n",
    "\n",
    "avg_prob = 0\n",
    "for _ in range(num_trials):\n",
    "    times, states, state_matrix = kmc(P,epochs=100000,init_state=0)\n",
    "    prob_array = probabilities(state_matrix)\n",
    "    avg_prob += prob_array / num_trials\n",
    "    avg_prob_A = prob_array[0] + prob_array[1] + prob_array[2]\n",
    "    avg_prob_B = prob_array[3] + prob_array[4] +prob_array[5]\n",
    "    \n",
    "print(\"Average equilibrium probability distribution (including microstates ):\", avg_prob)\n",
    "print(\"Average equilibrium probability of state A: \", avg_prob_A)\n",
    "print(\"Average equilibrium probability of state B: \", avg_prob_B)"
   ]
  }
 ],
 "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.6.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}