CS166 Session 2.2 Pre-Class Work

CS166 Session 2.2 Pre-Class Work.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 11.3\n",
    "Modify Code 11.5 to implement a simulator of the Game of Life CA. Simulate the dynamics from a random initial configuration. Measure the density of state 1’s in the configuration at each time step, and plot how the density changes over time. This can be done by creating an empty list in the \"initialize\" function, and then making the measurement and appending the result to the list in the \"observe\" function. The results stored in the list can be plotted manually after the simulation, or they could be plotted next to the visualization using pylab’s \"subplot\" function during the simulation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 123,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPsAAAD7CAYAAACscuKmAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAP9klEQVR4nO3dXYhc533H8e+vUhTXDkayvTaKZFcyiCQmkDq7uHZcSrASmrgh9oUDDqGIItBN2jgvkMjtRehdDSF2LkpAWA2imMSpYiohQoJRnIveqN6NTWNbdqTaRt5YsdZgJyU3rci/F3NGGq1mtWdmzstzzvP7wLA7bzv/88w+83vOc15GEYGZ9d8ftV2AmTXDnd0sE+7sZplwZzfLhDu7WSbc2c0yMVNnl/QpSa9IOi1pf1VFmVn1NO12dkkbgF8BnwSWgWeBz0fES9WVZ2ZV2TjDc+8ATkfEqwCSfgDcB6zZ2W+44YbYsWPHDC+ZpqWlpQu/z8/Pt1iJ5e7111/n7bff1rj7Zuns24A3Rq4vA3+2+kGS9gH7AG655RYWFxdneMk0SRfbto/LZ92xsLCw5n2zrLOP+/S4bJ0gIg5ExEJELMzNzc3wcsWLSpdd2hYRFy7T1NT0cqTSbtasWTr7MnDzyPXtwJuzlWNmdZmlsz8L7JK0U9Im4EHgaDVlrW00RYeXlExS0zBhm1qOpl+vq/o68pl6nT0izkv6W+CnwAbgXyLixcoqM7NKzTJBR0T8GPhxRbWYWY1m6ux9NDp8q3uo2/RQ2kP3cobtNPxf6Eu7eXdZs0z0Itmr+ATuyqf46omj1Ovtsr61rZPdLBO9SPac1J02XRnh9ElT80ROdrNMdDrZq0yhnJOsyS0Q1h4nu1kmOp3suaRQ3cmbSzuW1fS8RVOv42Q3y0Snk92sDintQVflfhVOdrNMuLObZcLD+A640tAthaFmX6UwfK+yBie7WSac7B2X0mSSVaeO99HJbpYJJ3sFxp2vzCemsNQ42c0y4WSvQB9T1XMA/eNkN8uEk92sQinM36zFyW6WCXd2s0x4GG9jpTL07JqU283JbpYJJ7t1hjcHXjTN2Yuc7GaZcLJbZ/ign9k42c0y4WS3zvHJPKZbPie7WSac7NYrXq9fm5PdLBNOduulLiV6Mt/iKulmSc9IOinpRUkPFbdfJ+lpSaeKn1tqq9LMZlZmGH8e+FpEfAi4E/iipNuA/cDxiNgFHC+um1mi1u3sEXE2In5R/P4/wElgG3AfcKh42CHg/rqKNBsl6ZJL10XEhUudyzTRBJ2kHcDtwAngpog4C4MPBODGNZ6zT9KipMWVlZXZqjWzqZXu7JLeB/wI+HJE/K7s8yLiQEQsRMTC3NzcNDWaXWI0CUfTsE8pX4dSnV3Sexh09Cci4qni5rckbS3u3wqcq6VCM6tEmdl4AQeBkxHx7ZG7jgJ7it/3AEeqL89sfaMpb2srs539buCvgV9Ker647e+BfwJ+KGkvcAb4XD0lmlkV1u3sEfEfwForQ7urLcfakPIZUa063l3WLBPeXbYnZjnwwymeBye7WSac7AWvt1rfOdnNMuHObpYJD+MLXR+yd7V+n1GmOU52s0y4s/dEXw4Esfq4s5tlwuvs1gqvqzfPyW6WCSd7x6yViF1LyK7V2wdOdrNMJJ/sq2eYU0iEptc3mzqveGrWa2ev90/GyW6WCXd2s0wkP4zPcYi2eniaehtUOZwus8ri4ft0nOxmmUg+2atS5URfXYkySWKlOHFZhTLL0ZdlbZqT3SwT2SR7F9JgkhpTWp6UarG1OdnNMpFssqc849rX9WXrNye7WSaSS/aUE30o5dpS14X3t6+c7GaZSC7Zh5/4dSeAE6Y5bRzI43mVyznZzTLhzm6WieSG8XXz8L15Vbd1mSH6egfRlHls3zjZzTKRbLLX9Wmby6d4H1UxKkvp/W96lOFkN8tE6c4uaYOk5yQdK67vlHRC0ilJT0raVF+ZZoPUSymZZzVcntFLnd/sM0myPwScHLn+CPBoROwC3gH2VlmYmVWrVGeXtB34K+Dx4rqAe4DDxUMOAfeXfdGufy/ZsP6uL4elp87RS9lkfwz4OvCH4vr1wLsRcb64vgxsG/dESfskLUpaXFlZmalYM5veurPxkj4DnIuIJUkfH9485qFjP44i4gBwoPhbIanz611dr9+6o8r9Qspsersb+Kyke4GrgGsZJP1mSRuLdN8OvDlzNWZWm3WH8RHxcERsj4gdwIPAzyLiC8AzwAPFw/YAR2qr0sxmNst29m8AX5V0msE6/MH1njA/P+8hsFkJw8nfKifsJtqDLiJ+Dvy8+P1V4I5KqjCz2iW7u6xZnVL/ssw6avLusmaZcLKvkvonvvVbnYdgO9nNMuFkL/ikFmloamSV4/vsZDfLhJO90OYnfV9GFVWfXKIv7bKeWUczZdvJyW6WCSe7Vabqc/73PdGHZknzSZ7vZDfLhDu7WSY8jE9A34arfVueWdUx0TjN33Kym2XCyd4z4yZu/CWH7UrlfXCym2XCyd4zOe6U0hVtvw9OdrNMONl7rO0ksbQ42c0y4c5ulgl3drNMuLObZSKZzu4vSUyL34/+Saazm1m9ktn0VvWx0DYdt39/OdnNMpFMsg85Udrl9u8vJ7tZJpJL9iqMm0VeK7HaPuzQrClOdrNM9DLZJ0lnJ7nlwslulolOd3bv5WVWXqc7u5mV585ulolSnV3SZkmHJb0s6aSkuyRdJ+lpSaeKn1vqLrao5cIlIqaeYPMqgKWsjv/Pssn+HeAnEfFB4CPASWA/cDwidgHHi+tmlqh1O7uka4G/AA4CRMT/RsS7wH3AoeJhh4D76yqySlWMCszqNvz/rDLhyyT7rcAK8D1Jz0l6XNI1wE0RcbYo7Cxw47gnS9onaVHS4srKSiVFm9nkynT2jcBHge9GxO3A75lgyB4RByJiISIW5ubmpizzkr83Uyp3IdE9n2BDVf6/lunsy8ByRJworh9m0PnfkrQVoPh5rpKKzKwW63b2iPgN8IakDxQ37QZeAo4Ce4rb9gBHaqkwQ3Wsr5mV3Tf+74AnJG0CXgX+hsEHxQ8l7QXOAJ+rp0Qzq0Kpzh4RzwMLY+7aXW05zenCoa2z1pTzKaZmWfa+tpv3oDPLhDu7WSYaPZ59aWnpkuFzm8Okvg3RbHZ9Hb4POdnNMtFoZ5+fn595pxgrL+d2nmbzZd/by8lulolenoNumnWvSc5Ia93h9/AiJ7tZJlpN9qpn5meZTa0qAZqa0e3CTkGWFie7WSZaTfbRNKoiEXNKt5yWdbW+bQ9vat8TJ7tZJpKZje/Lp3RfliNFfU30WQ7WmeT5TnazTLizm2WilWF834ZjTapihyG3exqafh+c7GaZaOUQVyfLZGbdNNOX9u7LclRhmrZwsptlopVDXG0yPiy4e0a/kzAVTnazTLizm2XCnd0sE8luZ/e24e7LeX+KFJfZyW6WiVaSffipd6Xtx1WezKLqv2tXlkOid3Hk6WQ3y4Q7u1kmkjlTTRf/fgpSHE62WUNTqxAptPOknOxmmUjmTDU2mRwmwaaxevLX7XORk90sE072juprYlV1ptW+ts8snOxmmSjV2SV9RdKLkl6Q9H1JV0naKemEpFOSnpS0qe5irb+Gh4P25XDeTh7iKmkb8CVgISI+DGwAHgQeAR6NiF3AO8DeOgs1s9mUHcZvBP5Y0kbgauAscA9wuLj/EHB/9eXZaqmlRRmjKbdW/X1I89St29kj4tfAt4AzDDr5b4El4N2IOF88bBnYNu75kvZJWpS0uLKyUk3VZjaxMsP4LcB9wE7g/cA1wKfHPHTsx3JEHIiIhYhYmJubm6VWM5tBmU1vnwBei4gVAElPAR8DNkvaWKT7duDN+srMW1Nf/FeXcTW3vdNLXUdEtr1cV1Jmnf0McKekqzVYkt3AS8AzwAPFY/YAR+op0cyqUGad/QSDibhfAL8snnMA+AbwVUmngeuBgzXWaVNIeTKv7Qm50U18VdbS9nJdSak96CLim8A3V938KnBH5RWZWS28u2wHTJsUPijERnl3WbNMONlbVEXiljl5hRO9el0cLTnZzTLhZG9Y1dvMu5QsfdDFRB9ysptlwsnesC4mgl3U5ffPyW6WCXd2s0w0OoxfWlq6cDaSVE0zgdbGgSqTnC++y5NKfZBK+zvZzTLRaLLPz8+zuLjY5EtObPTTd71P5DY/sdtOCSuv7t2Wyx7s5GQ3y4Q3vV3Bep/AXUnXrtTZd3W9D6N/d2FhYc3HOdnNMuFkt9rUdeonm46T3SwTTnarTQopXvcWk1S2oZfhZDfLhDu7WSbc2W1qKZ+91i7nzm6WiVYn6Lr+TSeWvq5NzNU54edkN8tEK8nepc0Vtrac378ql72pEa6T3SwTrSR7Cong0UV13JazaardnOxmmchud1mnUPWa+k65SU7FZZdzsptlIrtkdxrUp+629Xs3Gye7WSbc2c0ykWxn90EWZtVKtrObWbWSnaBranOOWS6c7GaZUJOJKWkF+D3wdmMvOpsb6E6t0K16u1QrdKfeP4mIuXF3NNrZASQtRsTaZ7JPSJdqhW7V26VaoXv1juNhvFkm3NnNMtFGZz/QwmtOq0u1Qrfq7VKt0L16L9P4OruZtcPDeLNMuLObZaKxzi7pU5JekXRa0v6mXrcsSTdLekbSSUkvSnqouP06SU9LOlX83NJ2rUOSNkh6TtKx4vpOSSeKWp+UtKntGockbZZ0WNLLRRvflWrbSvpK8T/wgqTvS7oq5bYtq5HOLmkD8M/Ap4HbgM9Luq2J157AeeBrEfEh4E7gi0WN+4HjEbELOF5cT8VDwMmR648Ajxa1vgPsbaWq8b4D/CQiPgh8hEHdybWtpG3Al4CFiPgwsAF4kLTbtpyIqP0C3AX8dOT6w8DDTbz2DDUfAT4JvAJsLW7bCrzSdm1FLdsZdJB7gGOAGOzhtXFcm7dc67XAaxQTwiO3J9e2wDbgDeA6BseOHAP+MtW2neTS1DB+2IBDy8VtSZK0A7gdOAHcFBFnAYqfN7ZX2SUeA74O/KG4fj3wbkScL66n1Ma3AivA94rVjsclXUOCbRsRvwa+BZwBzgK/BZZIt21La6qzjzswPcltfpLeB/wI+HJE/K7tesaR9BngXEQsjd485qGptPFG4KPAdyPidgbHR7Q+ZB+nmDe4D9gJvB+4hsHq52qptG1pTXX2ZeDmkevbgTcbeu3SJL2HQUd/IiKeKm5+S9LW4v6twLm26htxN/BZSa8DP2AwlH8M2CxpeNhySm28DCxHxIni+mEGnT/Ftv0E8FpErETE/wFPAR8j3bYtranO/iywq5jR3MRgwuNoQ69digYHzh8ETkbEt0fuOgrsKX7fw2BdvlUR8XBEbI+IHQza8mcR8QXgGeCB4mFJ1AoQEb8B3pD0geKm3cBLJNi2DIbvd0q6uvifGNaaZNtOpMGJj3uBXwH/DfxD25MVY+r7cwZDs/8Cni8u9zJYFz4OnCp+Xtd2ravq/jhwrPj9VuA/gdPAvwHvbbu+kTr/FFgs2vffgS2pti3wj8DLwAvAvwLvTblty168u6xZJrwHnVkm3NnNMuHObpYJd3azTLizm2XCnd0sE+7sZpn4f/2umd/fke8VAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib\n",
    "import pylab\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "n = 100 # size of space: n x n\n",
    "p = 0.1 # probability of initially dead cells\n",
    "\n",
    "def initialize():\n",
    "    global config, nextconfig\n",
    "    config = pylab.zeros([n, n])\n",
    "    for x in range(n):\n",
    "        for y in range(n):\n",
    "            config[x, y] = 1 if pylab.random() < p else 0\n",
    "    nextconfig = pylab.zeros([n, n])\n",
    "    \n",
    "denslist=[]\n",
    "\n",
    "def observe():\n",
    "    global config, nextconfig\n",
    "    pylab.cla()\n",
    "    pylab.imshow(config, vmin = 0, vmax = 1, cmap = pylab.cm.binary)\n",
    "    density =[]\n",
    "    # Here, I obtained the density of the plot by taking the sum of all the plots that were filled and dividing\n",
    "    # that by the total number of plots. Then, I took the average of those numbers and appended all of them\n",
    "    # into an array named \"denslist\". This will be used to show the change in density with respect to time,\n",
    "    # which will be displayed in a plot below.\n",
    "    density.append(sum((sum(config==1)/(sum(config==1)+sum(config==0))))/len((sum(config==1)/(sum(config==1)+sum(config==0)))))\n",
    "    denslist.append(density)\n",
    "\n",
    "def update():\n",
    "    density = []\n",
    "    global config, nextconfig\n",
    "    for x in range(n):\n",
    "        for y in range(n):\n",
    "            count = 0\n",
    "            # According to game of life, if the cell is dead, it'll become alive \n",
    "            # if there are exactly 3 alive cells around it\n",
    "            if config[x, y] == 0: \n",
    "                for dx in [-1, 0, 1]:\n",
    "                    for dy in [-1, 0, 1]:\n",
    "                        count += config[(x + dx) % n, (y + dy) % n]\n",
    "                if count == 3:\n",
    "                    nextconfig[x, y] = 1 \n",
    "                else:\n",
    "                    nextconfig[x, y] = 0\n",
    "            # On the other hand, if the cell is alive, it'll stay alive \n",
    "            # if there are 2 or 3 alive cells around it\n",
    "            else:\n",
    "                for dx in [-1, 0, 1]:\n",
    "                    for dy in [-1, 0, 1]:\n",
    "                        count += config[(x + dx) % n, (y + dy) % n]\n",
    "                if count == (2 or 3):\n",
    "                    nextconfig[x, y] = 1 \n",
    "                else:\n",
    "                    nextconfig[x, y] = 0\n",
    "    config, nextconfig = nextconfig, config\n",
    "\n",
    "%matplotlib inline\n",
    "\n",
    "initialize()\n",
    "for i in range(10):\n",
    "    update()\n",
    "    observe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 124,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x11984c400>]"
      ]
     },
     "execution_count": 124,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEWCAYAAAB8LwAVAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO3deXxU9b3/8dcnOxB2QqKsQREhiKgBF3Bp9Uq0rq1tpVqt0sW2aL3e3tba1lpvN3tbl1ZvtT+X2rpv7aVeBbVUXLBIwI2wKKvsa9gCIST5/P44JzDESRggk5PMvJ+PxzyYOcucT07CvOd8v+d8j7k7IiIijWVEXYCIiLRNCggREYlLASEiInEpIEREJC4FhIiIxKWAEBGRuBQQ7ZSZvWpmX22lbX3TzNaa2XYz69lo3kAzczPLao1amqjvDDNbkeRt3GRm9ydzG+F2lprZWcneTjpord9ZKlNAtGHhh8XO8IN5rZk9ZGb5B/geh/QBbmbZwO3A2e6e7+4bD+Z92jt3/4W7fxXaRijKvuJ9SYj9ncnBUUC0fee7ez5wPDAK+FErb78QyAMqWnm7InEpmFuPAqKdcPeVwIvA8MbzzCzDzH5kZsvMbJ2Z/dnMuoazXwv/3RweiZwcZ/1cM7vTzFaFjzvDaUcBC2LWn9pMiZeZ2cdmtsHMfri/9w7nTTOzz4XPx4bfys8NX59lZu/G25CZdTCzP5lZpZnNJQjO2PmHm9mzZrbezJaY2XUx824xs6fCfbTNzCrMrDRm/vfNbGU4b4GZnRmz3iNN7NPTzWyTmR0T8z69w6O/giZ+hq+Z2bxwO3PN7PiY2SPN7H0z22JmT5pZXrhOdzN7Pvy5KsPnfWPe81Uz+y8zezN835fMrFfM/CvCv5GNZvbj2Oas8G/oRjNbFM5/ysx6hPPyzOyRcPpmM5tpZoVN/FxDwzo2h/v2gnD6SWa2xswyY5a92MzeT2D7DUdsE8zsY2Bqo212Ivi/cXj4+9ge/g3s+Z3FvMdVZrY83H/XmNmocF9vNrO7G73v1eHvqNLMppjZgHg/cypTQLQTZtYPOBd4J87sr4SPTwGDgHyg4Y/9tPDfbmET0Vtx1v8hcBIwEjgWGA38yN0/BEpi1v90MyWOBYYAZwI3m9nQ5t47nDcNOCOmzsXA6TGvpzWxrZ8AR4SPccCVDTPMLAP4O/Ae0Ces53ozGxez/gXAE0A3YBLhvjKzIcBEYJS7dw7fe2mc7Tfep9PC97s8ZpnxwCvuvr7xymb2eeAW4AqgS1hPbNPdF4AyoBgYQfC7heD/60PAAKA/sJO9v+cGXwKuAnoDOcB3w20OA/4HuAw4DOga7p8G1wEXEez/w4FK4J5w3pXh8v2AnsA14bYb/1zZBPv+pXD71wKPmtkQd/8XUAXE/g19CXgsge03OB0YSvB72cPdq4BzgFXh7yPf3Vc1ri90IjAY+CJwJ8Hf51kEf+dfMLPTw5/lIuAm4LNAAfA68HgT75m63F2PNvog+HDaDmwGlhH8B+8QznsV+Gr4/B/At2LWGwLsBrKAgYADWc1sZxFwbszrccDS8Hmz68fM7xsz7W3g0gTe+0zg/fD5ZOCrwL/C19OAzzaxzcVAWczrrwMrwucnAh83Wv4HwEPh81sIPrgb5g0DdobPjwTWEXxgZDd6j1uAR5raJ+F2lwMZ4ety4AtN1D8F+E4zv/PLY17/Gri3iWVHApUxr18lCPaG198CJofPbwYej5nXEagBzgpfzwPOjJl/WMzf0NXAdGDEfv5eTwXWNOyDcNrjwC3h858BD4bPOxMExoAEtt+wvwc1s+0zGv4G9vM76xMzfyPwxZjXzwLXh89fBCbEzMsAdjTUmy4PHUG0fRe5ezd3H+Du33L3T3xzI/jGtSzm9TKC/1hxmwESXP/wA6xzTczzHQRHMft777eAo8LmipHAn4F+YbPIaPY25cSrd3mj92wwgKCpYXPDg+CbYOy+aFxrnplluftC4HqCD5Z1ZvaEmSW0H9x9BsEH3ulmdjRB2ExqYvF+BMHZlLj70sw6mtl9YTPRVoL90y222aapdWm0z9x9B/setQwA/hqzz+YBdQT77S8EofaEBc2Evw6PFho7HFju7vUx05ax90jlMeCzFjQxfhaY7e4Nv7vmtt8g9nd+sNbGPN8Z53XD/hoA3BVTzybA2PeoK+UpIFLDKoI/6Ab9gVqCP/5EhuuNt35Th+gtUdsq2PMhNQv4DjDH3WsIvqneACxy9w1NvOdqgg/Z2PdssBxYEoZqw6Ozu5+bSLHu/pi7jw1rduC2eIs1sfrDBM1MXwaecffqJpZbTtA8dqD+g+Do8ER378Lepi5LYN3VQGx/RQeC5qLYms5ptN/y3H2lu+9295+6+zDgFOA8guaxxlYRBHzs50p/YCWAu88lCIxz2Ld5qdntxyzT3N9ySw9LvRz4RqN6Orj79BbeTpumgEgNjwP/bmbFFpwG+wvgSXevBdYD9QR9E82t/yMzKwi/vd8MPNLM8gdaW3PvPY2g3b+hv+HVRq/jeQr4Qdhp25egrbvB28BWCzqbO5hZppkNN7NR8d9qLzMbYmafDr/hVhN8o6yLs2hT+/QvwMUEIfHnZjZ1P/BdMzvBAkcm2AHaOaxpc9iB+5ME1mnwDHC+mZ1iZjnAT9k3WO4Fft5QR/j7ujB8/ikzOyY8UtlK0PQTb780HEV9z8yyzewM4HyC/pkGjxH0N5wGPJ3I9hO0Fuhpe0/OOFT3EvyNlYT1dA37jtKKAiI1PEjw4fQasITgw+1a2PMt/efAm+Hh8klx1v8ZQZv5+8AHwOxwWkvY33tPI/jge62J1/H8lOCb6BKCDtG/NMxw9zqCD6WR4fwNBB/IiXxw5AK/CtdZQ9DRelPjhZrap+6+Ivz5nKBTMy53fzpc/zFgG/A3oEcC9d0JdAjr+xdBv01C3L2C4G/iCYKjiW0E/S27wkXuImgSe8nMtoXvf2I4r4ggYLYSNP1MI84XiPAI8AKCI4QNBH1mV7j7/JjFHifoL5ja6Aixue0n8vPND997cfg7OdAm0sbv91eCo8cnwua8OQQ/V1qxsANGRFqAmT1IcDZNa1+vckDCI83NwGB3XxJ1PdI26YITkRZiZgMJOl+Pi7aS+MzsfIIz3gz4DcER3dIoa5K2TU1MIi3AzP6LoBniv9vwN/ILCTqSVxFcC3CpqwlBmqEmJhERiUtHECIiElfK9EH06tXLBw4cGHUZIiLtyqxZsza4e9wxw1ImIAYOHEh5eXnUZYiItCtmtqypeUltYjKzMgtGxFxoZjfGmZ9rwWiVC81sRngWCGZ2mZm9G/OoN7ORyaxVRET2lbSACK+6vIfg4pJhwPhwRMlYEwgGGzsSuINwWAN3f9TdR7r7SIJhC5a6e9yhn0VEJDmSeQQxGljo7ovDKyyfIDjNLtaFBOPXQHCl5plm1nhcmfGk4zC7IiIRS2ZA9GHf0RdX8MmREPcsE44btIV9BxCDYNz2uAFhZl83s3IzK1+//hPD7ouIyCFIZkDEG2Gy8UUXzS5jZicCO9x9TrwNuPsf3b3U3UsLCuJ2wouIyEFKZkCsYN8hmfvyySGk9yxjwX1muxKMu97gUtS8JCISiWQGxExgcDgEdQ7Bh33jG6hMYu/tIi8hGOHRYc+tIz/PvkMFi4hIK0laQIR9ChMJ7kQ1D3jK3SvM7FYLb2QOPEAwhvtCgpvExJ4KexrBLQQXJ6tGgFWbd/KrF+ezdmtT93YREUlPSb1Qzt1fAF5oNO3mmOfVBEcJ8dZ9leBm90lVtauWe6ctok+3PL588sBkb05EpN1I+7GYjuydz6CCTkypWLv/hUVE0kjaB4SZMa6kiLcWb2TzjpqoyxERaTPSPiAAykqKqKt3Xpm3LupSRETaDAUEMKJvVw7rmseUijVRlyIi0mYoINjbzPTah+up2lUbdTkiIm2CAiI0rqSIXbX1TPtQQ3aIiIACYo9RA7vTo1OOmplEREIKiFBWZgZnDe3N1HnrqKmtj7ocEZHIKSBilA0vYtuuWqYv2hB1KSIikVNAxDjliF7k52apmUlEBAXEPvKyMzljSAEvVaylrr7xyOQiIulFAdFI2fAiNlbVMGtZZdSliIhESgHRyBlDepOTlcHkOWpmEpH0poBoJD83i1OP7MWUijWEt6YQEUlLCog4xg0vYuXmnVSs2hp1KSIikVFAxHHW0EIyDDUziUhaU0DE0aNTDicW92SyTncVkTSmgGhC2fAiFq7bzsJ126MuRUQkEgqIJpxdUgigi+ZEJG0pIJpwWNcOHNuvmwJCRNKWAqIZZSVFvL9iCys374y6FBGRVqeAaMa4sJnpJR1FiEgaUkA0Y1BBPkcV5ut0VxFJSwqI/SgrKWLm0k1s3L4r6lJERFqVAmI/zi4pot7hlXlroy5FRKRVKSD2o+TwLvTt3kHNTCKSdhQQ+2FmlJUU8ebCjWyr3h11OSIirUYBkYBxw4uoqavnnwvWR12KiEirUUAk4Pj+3emVn8sUNTOJSBpRQCQgM8M4u6SQfy5YR/XuuqjLERFpFQqIBI0rKWJHTR1vfLQh6lJERFqFAiJBJw/qSee8LA0BLiJpQwGRoJysDM4aWsgr89ZSW1cfdTkiIkmngDgA40oK2bxjN28v2RR1KSIiSaeAOACnHVVAXnaGhgAXkbSggDgAHXOyOP2oAqZUrKW+3qMuR0QkqZIaEGZWZmYLzGyhmd0YZ36umT0Zzp9hZgNj5o0ws7fMrMLMPjCzvGTWmqhxJUWs2VrNeys2R12KiEhSJS0gzCwTuAc4BxgGjDezYY0WmwBUuvuRwB3AbeG6WcAjwDXuXgKcAbSJcS7OPLqQrAxjSoUG7xOR1JbMI4jRwEJ3X+zuNcATwIWNlrkQeDh8/gxwppkZcDbwvru/B+DuG929TVyh1rVjNicf0ZPJc1bjrmYmEUldyQyIPsDymNcrwmlxl3H3WmAL0BM4CnAzm2Jms83se/E2YGZfN7NyMytfv771xkkaV1LE0o07+HDt9lbbpohIa0tmQFicaY2/cje1TBYwFrgs/PdiMzvzEwu6/9HdS929tKCg4FDrTdjZwwoxQ2cziUhKS2ZArAD6xbzuC6xqapmw36ErsCmcPs3dN7j7DuAF4Pgk1npAenfJ44T+3XWPCBFJackMiJnAYDMrNrMc4FJgUqNlJgFXhs8vAaZ60LA/BRhhZh3D4DgdmJvEWg/YuJIi5q7eyvJNO6IuRUQkKZIWEGGfwkSCD/t5wFPuXmFmt5rZBeFiDwA9zWwhcANwY7huJXA7Qci8C8x29/9LVq0HY1xJEaBmJhFJXZYqZ+KUlpZ6eXl5q27z3Ltep2NOJs9885RW3a6ISEsxs1nuXhpvnq6kPgTjSoqY9XEl67ZVR12KiEiLU0AcgrLhRbjDy3N10ZyIpB4FxCE4qjCf4l6ddDaTiKQkBcQhMAtuRfrWoo1s2dEmRgIREWkxCohDVFZSRG29M3WBmplEJLUoIA7RsX27UdglV81MIpJyFBCHKCPDGFdSxLQP17Ozpk2MJygi0iIUEC2grKSI6t31TPuw9QYMFBFJNgVECxhd3INuHbN1VbWIpBQFRAvIyszgrKGFvDJvLTW19VGXIyLSIhQQLaSspIht1bX8a/HGqEsREWkRCogWMnZwLzrmZDJZzUwikiIUEC0kLzuTTw3pzUsVa6mrT40BEEUkvSkgWtC44UVs2L6Ldz6ujLoUEZFDpoBoQZ8aUkBOZoYumhORlKCAaEGd87IZc2RPpsxdQ6rcZ0NE0pcCooWVDS9i+aadzF29NepSREQOiQKihZ01tJAMgylqZhKRdk4B0cJ65ucyamAPplRodFcRad8UEElQNryIBWu3sXj99qhLERE5aAqIJDi7pAhARxEi0q4pIJKgT7cOjOjbVYP3iUi7poBIknElRby7fDOrt+yMuhQRkYOigEiScWEz00tqZhKRdkoBkSRH9s7nyN75amYSkXZLAZFEZSVFzFiyiU1VNVGXIiJywBQQSTSupIi6eueVeWpmEpH2RwGRRMP7dKFPtw68pGYmEWmHFBBJZGaMKynitY82sH1XbdTliIgcEAVEko0rKaSmtp5XF6yLuhQRkQOigEiy0oE96NkpR1dVi0i7o4BIsswM4+ySQqbOW0v17rqoyxERSZgCohWcXVJEVU0d0xdtiLoUEZGEKSBawSlH9KRzbhZT5qiZSUTaDwVEK8jNyuTTQ3vz8ry11NbVR12OiEhCFBCtZFxJEZuqapi5tDLqUkREEpLUgDCzMjNbYGYLzezGOPNzzezJcP4MMxsYTh9oZjvN7N3wcW8y62wNpx9VQG5WhsZmEpF2I6GAMLNnzewzZpZwoJhZJnAPcA4wDBhvZsMaLTYBqHT3I4E7gNti5i1y95Hh45pEt9tWdcrN4rSjCnipYg3uHnU5IiL7legH/h+ALwEfmdmvzOzoBNYZDSx098XuXgM8AVzYaJkLgYfD588AZ5qZJVhTuzOupIhVW6r5YOWWqEsREdmvhALC3V9x98uA44GlwMtmNt3MrjKz7CZW6wMsj3m9IpwWdxl3rwW2AD3DecVm9o6ZTTOzU+NtwMy+bmblZla+fv36RH6USJ01tDeZGcbkOWpmEpG270CajHoCXwG+CrwD3EUQGC83tUqcaY3bVppaZjXQ392PA24AHjOzLp9Y0P2P7l7q7qUFBQUJ/RxR6tYxh5MG9WCy+iFEpB1ItA/iOeB1oCNwvrtf4O5Puvu1QH4Tq60A+sW87gusamoZM8sCugKb3H2Xu28EcPdZwCLgqMR+pLatrKSIxeurWLhuW9SliIg0K9EjiPvdfZi7/9LdV0NwBhKAu5c2sc5MYLCZFZtZDnApMKnRMpOAK8PnlwBT3d3NrCDs5MbMBgGDgcUJ/1Rt2NnhrUjVzCQibV2iAfGzONPeam6FsE9hIjAFmAc85e4VZnarmV0QLvYA0NPMFhI0JTWcCnsa8L6ZvUfQeX2Nu29KsNY2rbBLHsf176ZmJhFp87Kam2lmRQQdyR3M7Dj29hl0IWhuapa7vwC80GjazTHPq4HPx1nvWeDZ/b1/e1VWUsQvX5zPisod9O2+390oIhKJ/R1BjAN+Q9B/cDvw2/BxA3BTcktLXePCZiYNAS4ibVmzRxDu/jDwsJl9LvxWLy1gYK9OHF3UmSlz1jBhbHHU5YiIxLW/JqbL3f0RYKCZ3dB4vrvfnrTKUty4kiJ+N/Uj1m/bRUHn3KjLERH5hP01MXUK/80HOsd5yEEqG16EO7wyT81MItI27a+J6b7w35+2Tjnp4+iizvTv0ZHJc9YwfnT/qMsREfmERC+U+7WZdTGzbDP7h5ltMLPLk11cKjMzyoYXMX3RBrZW7466HBGRT0j0Ooiz3X0rcB7B1c9HAf+ZtKrSxLiSInbXOf+cvy7qUkREPiHRgGgYkO9c4PFUuWgtasf160bvzrm6qlpE2qREA+LvZjYfKAX+YWYFQHXyykoPGRnG2SWFvLpgPdW766IuR0RkH4kO930jcDJQ6u67gSo+eW8HOQhlJYexc3cdr33Y9ocrF5H00uxZTI0MJbgeInadP7dwPWnnxEE96Nohm8kVa/YM5Cci0hYkFBBm9hfgCOBdoKEtxFFAHLLszAzOHNqbf8xbx+66erIzk3qbcBGRhCV6BFEKDHPdTDkpykqKeG72SmYs3sTYwb2iLkdEBEi8k3oOoPaPJDntqAI6ZGcyuWJ11KWIiOyRaED0Auaa2RQzm9TwSGZh6SQvO5MzhhTwUsVa6ut1kCYibUOiTUy3JLMIgQtHHs6Lc9Zw72uL+NYZR0ZdjohIwqe5TgOWAtnh85nA7CTWlXbGlRRx3ojD+O1LHzJzqa5DFJHoJToW09cIbv15XzipD/C3ZBWVjsyMX372GPp278B1j79DZVVN1CWJSJpLtA/i28AYYCuAu38E9E5WUemqc14293zpeDZur+E/nn5P/REiEqlEA2KXu+/5ShteLKdPryQY3qcrP/zMUKbOX8cDbyyJuhwRSWOJBsQ0M7sJ6GBm/wY8Dfw9eWWltytOHkBZSRG3TZ7P7I8roy5HRNJUogFxI7Ae+AD4BvAC8KNkFZXuzIzbLhnBYd3yuPaxd9iyQ/eLEJHWl+hZTPUEndLfcvdL3P3/6arq5OraIZu7xx/Pum3VfPeZ99DuFpHW1mxAWOAWM9sAzAcWmNl6M7u5dcpLb8f268b3y47m5blreejNpVGXIyJpZn9HENcTnL00yt17unsP4ERgjJn9e9KrEyaMLeasoYX88sV5vL9ic9TliEga2V9AXAGMd/c9p9O4+2Lg8nCeJJmZ8ZvPj6B35zwmPvaO7l8tIq1mfwGR7e4bGk909/XsvQ2pJFm3jjn8bvxxrNq8kxuffV/9ESLSKvYXEM1dzqtLfVvRCQO685/jhvDCB2t45F/Loi5HRNLA/gbrO9bMtsaZbkBeEuqRZnzt1EG8tXgj//X8PI7r353hfbpGXZKIpLBmjyDcPdPdu8R5dHZ3NTG1sowM4/YvjKRHpxwmPjab7btqoy5JRFKY7m/ZzvToFPRHLK/cyU3PfaD+CBFJGgVEOzS6uAc3/NtRTHpvFU/MXB51OSKSohQQ7dQ3Tz+CUwf34pZJFcxbHa+bSETk0Cgg2qmG/oguHbKZ+NhsqtQfISItTAHRjhV0zuWuS0eyZEMVP/7fOVGXIyIpJqkBYWZlZrbAzBaa2Y1x5uea2ZPh/BlmNrDR/P5mtt3MvpvMOtuzU47oxXVnDua52St5ulz9ESLScpIWEGaWCdwDnAMMA8ab2bBGi00AKt39SOAO4LZG8+8AXkxWjani2k8P5uRBPbn5fyv4aO22qMsRkRSRzCOI0cBCd18c3o3uCeDCRstcCDwcPn8GONPMDMDMLgIWAxVJrDElZGYYd106kk65mXz7sdnsrKmLuiQRSQHJDIg+QGybx4pwWtxl3L0W2AL0NLNOwPeBnyaxvpTSu0sed3xxJB+t284tk5SpInLokhkQFmda46u6mlrmp8Ad7r692Q2Yfd3Mys2sfP369QdZZuo4dXAB3z7jSJ4sX87f3lkZdTki0s4lMyBWAP1iXvcFVjW1jJllAV2BTQT3nPi1mS0luCfFTWY2sfEG3P2P7l7q7qUFBQUt/xO0Q9efNZjRA3tw018/YNH6ZvNVRKRZyQyImcBgMys2sxzgUmBSo2UmAVeGzy8BpnrgVHcf6O4DgTuBX7j73UmsNWVkZWbwu/HHkZuVwbcfnU31bvVHiMjBSVpAhH0KE4EpwDzgKXevMLNbzeyCcLEHCPocFgI3AJ84FVYOXFHXPG7/4kjmr9nGrc/PjbocEWmnLFUGeystLfXy8vKoy2hTfvniPO6btpjfjz+O8489POpyRKQNMrNZ7l4ab56upE5h3z17CCcM6M4PnvuApRuqoi5HRNoZBUQKyw77IzIzjImPz2ZXrfojRCRxCogU16dbB377+WOZs3Irv/i/eVGXIyLtiAIiDZw1rJAJY4t5+K1lvPjB6qjLEZF2QgGRJr5fdjTH9uvG9559n+WbdkRdjoi0AwqINJGTlcHd448DYOJjs6mprY+4IhFp6xQQaaRfj4789yUjeG/FFm6bPD/qckSkjVNApJmy4Ydx5ckDeOCNJbw8d23U5YhIG6aASEM3fWYow/t04btPv8eKSvVHiEh8Cog0lJuVyd3jj6eu3rnu8XfYXaf+CBH5JAVEmhrYqxO/+twxzP54M795aUHU5YhIG6SASGPnjTicy07sz33TFvPP+euiLkdE2hgFRJr78XnDOLqoMzc89S6rt+yMuhwRaUMUEGkuLzuTey47nl219Vz3+DvUqj9CREIKCOGIgnx+cfExzFxayZ2vfBR1OSLSRiggBICLjuvDF0v7cc+rC3n9I93fW0QUEBLjlgtKGNw7n289Opvn3298+3ARSTcKCNmjQ04mD35lFEcU5DPxsXf4z6ffo2pXbdRliUhEFBCyj77dO/L0NScz8VNH8szsFZz3+zd4f8XmqMsSkQgoIOQTsjMz+O64ITz+tZOo3l3HZ/9nOn94dRH19alx/3IRSYwCQpp00qCeTP7OaZxdUshtk+dz+QMzWLOlOuqyRKSVKCCkWV07ZnPPl47n158bwTsfb6bsrteYPGdN1GWJSCtQQMh+mRlfGNWP568bS9/uHbjmkVn84LkP2FGjDmyRVKaAkIQdUZDPc98cwzdOH8Tjb3/M+b9/g4pVW6IuS0SSRAEhByQnK4MfnDOURyacyLbqWi6+Zzr3v75YHdgiKUgBIQdl7OBeTL7+NE4fUsDP/m8eVz70Nuu2qQNbJJUoIOSg9eiUwx+/fAI/u2g4M5du4pw7X2fqfN3GVCRVKCDkkJgZl580gL9PHEvvLnlc/adyfvK/c6jeXRd1aSJyiBQQ0iIGF3bmb98+hQlji3n4rWVcePebLFizLeqyROQQKCCkxeRmZfLj84bxp6tGsbGqhvPvfoOHpy/FXR3YIu2RAkJa3BlDejP5+lMZc0RPfjKpggkPl7Nh+66oyxKRA6SAkKTolZ/Lg18ZxS3nD+ONhRsou/N1pn2o+0yItCcKCEkaM+MrY4qZNHEMPTplc+WDb/Oz5+eyq1Yd2CLtgQJCku7ooi5MmjiWK04ewP1vLOHie6azcJ06sEXaOgWEtIq87ExuvXA4919Rypqt1Zz3+zd4dMYydWCLtGEKCGlVZw0rZPJ3TmXUwB788K9z+MZfZlFZVRN1WSISR1IDwszKzGyBmS00sxvjzM81syfD+TPMbGA4fbSZvRs+3jOzi5NZp7Su3l3yePiq0fzoM0P554J1lN31GtMXboi6LBFpJGkBYWaZwD3AOcAwYLyZDWu02ASg0t2PBO4AbgunzwFK3X0kUAbcZ2ZZyapVWl9GhvHVUwfx12+NoVNuFpc9MINfvTifmtr6qEsTkVAyjyBGAwvdfbG71wBPABc2WuZC4OHw+TPAmWZm7r7D3RtuNpAHqKE6RQ3v05Xnrx3LpaP6c++0RXzuD9NZvH571GWJCMkNiD7A8pjXK8JpcZcJA2EL0BPAzE40swrgA+CamMCQFNMxJ4tffvYY7r38BJZX7uC837/BU7yAQcUAAAuVSURBVDOXqwNbJGLJDAiLM63x//gml3H3Ge5eAowCfmBmeZ/YgNnXzazczMrXr9dFWO1d2fAiXvzOqRzbtxvfe/Z9rnxoJm8u3KCgEIlIMgNiBdAv5nVfYFVTy4R9DF2BTbELuPs8oAoY3ngD7v5Hdy9199KCgoIWLF2icljXDjzy1RP54blDqVi5hcvun8E5d73OUzOXa4RYkVaWzICYCQw2s2IzywEuBSY1WmYScGX4/BJgqrt7uE4WgJkNAIYAS5NYq7QhmRnG104bxJs3fppfXzICgO89+z5jfjWV21/+UDcmEmkllszDdzM7F7gTyAQedPefm9mtQLm7Twqbjf4CHEdw5HCpuy82sy8DNwK7gXrgVnf/W3PbKi0t9fLy8qT9LBIdd+etRRt58M0l/GP+OrIyjPOPPZyrxxQzvE/XqMsTadfMbJa7l8adlyrtuwqI9LBkQxUPT1/KU+XL2VFTx+jiHkwYW8xZQwvJzIjXpSUizVFASMrZsnM3T81czp+mL2Xl5p3069GBr5xSzBdK+9I5Lzvq8kTaDQWEpKzaunpenruWB95YQvmySvJzs/h8aV+uOqWY/j07Rl2eSJungJC08N7yzTz05hKef381de7829BCrh5bzInFPTBT85NIPAoISStrt1bzl7eW8eiMZVTu2E3J4V24ekwx5x17GLlZmVGXJ9KmKCAkLe2sqeNv767kwTeW8NG67fTKz+XLJw3gspP60ys/N+ryRNoEBYSkNXfnjYUbeOCNJby6YD05WRlcNPJwrhpTzNDDukRdnkikmgsIjZAqKc/MOHVwAacOLmDhuu38afoSnp21kqfKV3DKET2ZMLaYTw3pTYZOkxXZh44gJC1t3lHDEzOX8/D0pazeUs3Anh25akwxl5zQl065+t4k6UNNTCJN2F1Xz5SKNTzwxhLe+XgznfOyGD+6P1ecPIC+3XWarKQ+BYRIAmZ/XMlDby7lhQ9W4+6UDS/i6jHFnDCgu06TlZSlgBA5AKs27+TPby3j8bc/ZsvO3eRkZdCjYw7dO+XQs1Pwb4+O2fTolEuPTsG/3Ttl06NTDj065dC9Yw7Zmbrdu7QPCgiRg7Cjppbn31vNog3bqayqYVOjx9bqpu9h1SUvKwiLhlDpmLM3QPYJmhx65OfQOTdLRykSCZ3FJHIQOuZk8YVR/Zqcv7uunsodNVRW7WZj1S4qq3azaUcNm7bXULmjho1VNVRW1bBqczVzVm5lU1UNNXXx77mdnWl7QqR7GBqNj1oKO+dS1DWPwi555GXrgj9JPgWEyEHKzsygd+c8enfOAzrvd3l3p6qmjsqqveGx54gkDJZNO4Lp81ZtZdOOGjbv2B33vbp3zKawSx5FXfMo6hKExmFd8ygMXxd1yaNbx2wdlcghUUCItBIzIz83i/zcLPr1SOwMqdq6ejbv3M2mqhrWbq1mzZbwsbU6eL01ODrZsH3XJ9bNzcrYc8RRFBMme6Z1zaN351z1l0iTFBAibVhWZga98nPplZ/LUYVNH6XU1NazblsYGlt2sWZrNWu27GTN1l2s3VLNu8s3s6aimprafZu4zKBXfu6eo5Cirrkc1rVDTKjkUtglT0OopykFhEgKyMnKoG/3js1eu+HuVO7YzZote48+Yo9IVlTuoHzZprjNWvm5WRR2yaV35zzysjPIycogJyuT7EwjNyuDnMyGaRlkNzzPzAjmNZrWsNx+n2dmkKWjm0gpIETShJntOZNq2OFNj0FVvbuOtVurWd0QJDFNWuu37WJjVS01tfXU1Nazq7aemrrg+e7w39r6ljszMsPYEzCxQZSVmYF6V/Y6Y0gBP/zMsBZ/XwWEiOwjLzuTAT07MaBnp4Nav67e2V0XhEdDaNTEBElT02viTQuX3dVoWm19/LPB0lVhl7ykvK8CQkRaVGaGkZmRqVNxU4Aa+EREJC4FhIiIxKWAEBGRuBQQIiISlwJCRETiUkCIiEhcCggREYlLASEiInGlzA2DzGw9sOwQ3qIXsKGFymnvtC/2pf2xl/bFvlJhfwxw94J4M1ImIA6VmZU3dVeldKN9sS/tj720L/aV6vtDTUwiIhKXAkJEROJSQOz1x6gLaEO0L/al/bGX9sW+Unp/qA9CRETi0hGEiIjEpYAQEZG40j4gzKzMzBaY2UIzuzHqeqJkZv3M7J9mNs/MKszsO1HXFDUzyzSzd8zs+ahriZqZdTOzZ8xsfvg3cnLUNUXJzP49/H8yx8weN7Pk3NYtQmkdEGaWCdwDnAMMA8abWcvf2LX9qAX+w92HAicB307z/QHwHWBe1EW0EXcBk939aOBY0ni/mFkf4Dqg1N2HA5nApdFW1fLSOiCA0cBCd1/s7jXAE8CFEdcUGXdf7e6zw+fbCD4A+kRbVXTMrC/wGeD+qGuJmpl1AU4DHgBw9xp33xxtVZHLAjqYWRbQEVgVcT0tLt0Dog+wPOb1CtL4AzGWmQ0EjgNmRFtJpO4EvgfUR11IGzAIWA88FDa53W9mnaIuKiruvhL4DfAxsBrY4u4vRVtVy0v3gLA409L+vF8zyweeBa53961R1xMFMzsPWOfus6KupY3IAo4H/uDuxwFVQNr22ZlZd4LWhmLgcKCTmV0ebVUtL90DYgXQL+Z1X1LwMPFAmFk2QTg86u7PRV1PhMYAF5jZUoKmx0+b2SPRlhSpFcAKd284onyGIDDS1VnAEndf7+67geeAUyKuqcWle0DMBAabWbGZ5RB0Mk2KuKbImJkRtDHPc/fbo64nSu7+A3fv6+4DCf4uprp7yn1DTJS7rwGWm9mQcNKZwNwIS4rax8BJZtYx/H9zJinYaZ8VdQFRcvdaM5sITCE4C+FBd6+IuKwojQG+DHxgZu+G025y9xcirEnajmuBR8MvU4uBqyKuJzLuPsPMngFmE5z99w4pOOyGhtoQEZG40r2JSUREmqCAEBGRuBQQIiISlwJCRETiUkCIiEhcCghJO2ZWaGaPmdliM5tlZm+Z2cUR1XKGmZ0S8/oaM7siilpEGkvr6yAk/YQXNf0NeNjdvxROGwBckMRtZrl7bROzzwC2A9MB3P3eZNUhcqB0HYSkFTM7E7jZ3U+PMy8T+BXBh3YucI+732dmZwC3ABuA4cAs4HJ3dzM7AbgdyA/nf8XdV5vZqwQf+mMIrs7/EPgRkANsBC4DOgD/AuoIBsK7luCK3O3u/hszGwncSzBS6CLganevDN97BvApoBswwd1fN7MS4KFwGxnA59z9o5bZc5KO1MQk6aaE4OrXeCYQjMo5ChgFfM3MisN5xwHXE9w3ZBAwJhy36vfAJe5+AvAg8POY9+vm7qe7+2+BN4CTwoHungC+5+5LCQLgDncf6e6vN6rnz8D33X0E8AHwk5h5We4+OqypYfo1wF3uPhIoJRg/SeSgqYlJ0pqZ3QOMBWqAZcAIM7sknN0VGBzOe9vdV4TrvAsMBDYTHFG8HLRckUkw9HODJ2Oe9wWeNLPDCL7hL9lPXV0JAmZaOOlh4OmYRRoGUpwV1gLwFvDD8D4Wz+noQQ6VjiAk3VQQMwqpu3+boFmngGD492vDb/Mj3b04Zoz/XTHvUUfw5cqAipjlj3H3s2OWq4p5/nvgbnc/BvgGcKi3p2yop6EW3P0xgr6UncAUM/v0IW5D0pwCQtLNVCDPzL4ZM61j+O8U4Jth0xFmdtR+boqzAChouDezmWWH/QDxdAVWhs+vjJm+DejceGF33wJUmtmp4aQvA9MaLxfLzAYBi939dwT9HiOaW15kfxQQklY8OCvjIuB0M1tiZm8TNN98n+DWonOB2WY2B7iPZpphw9vUXgLcZmbvAe/S9D0BbgGeNrPXCTqzG/wduNjM3o0JgwZXAv9tZu8DI4Fb9/PjfRGYEzaBHU3QhyFy0HQWk4iIxKUjCBERiUsBISIicSkgREQkLgWEiIjEpYAQEZG4FBAiIhKXAkJEROL6/y1zzERrtvNdAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Here is the plot of how density changes over time, as aforementioned.\n",
    "\n",
    "plt.title('Plot of how density changes over time')\n",
    "plt.xlabel(\"Generations\")\n",
    "plt.ylabel(\"Density\")\n",
    "plt.plot(denslist)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 11.5\n",
    "Implement an interactive parameter setter for $p$ in Code 11.5. Then conduct systematic simulations with varying $p$, and identify its critical value below which isolated clusters are formed but above which the whole space is filled with panic."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 126,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPsAAAD7CAYAAACscuKmAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAARwElEQVR4nO3dX8hl1X3G8e/TmUysBvFPR5k4to4wmEggNe+L1SSUoAk1NkQvDBhCGYplbpLGJoFE24vQi14IIRpoCQzaIEUS04l0xIsEmWjo1dSZGhp1NDPVohMnOqWaQm7awV8vzn6b19fzZ5+z19577b2eDxze9/xfe529zrP22uvsrYjAzMbvt/ougJl1w43drBBu7GaFcGM3K4Qbu1kh3NjNCtGosUu6SdILkk5KuitVocwsPa26n13SNuDnwCeAU8BTwGcj4rl0xTOzVLY3eO61wMmIeBFA0veAW4CZjV2SZ/D0aG1tbertx44dW/o5y7zOMq+xzOvNK3dfVlnW1MsREZp2e5Nkvw24KSL+rLr+J8AfRMQXtjxuP7C/urrap25JzPqspanrxtznLPM6Ddaxua83r9x9WWVZUy/HrMbeJNmnveA7ljQiDgAHwMmemzor2cZj5q3Ei15n8/3LNIZFj82p0efQyBdpMkB3Crh80/XdwKvNimNmbWnS2J8C9kraI2kHcDvwaJpiWRskve3SZxn6FhHvuHSlrzpYuRsfEWclfQH4EbAN+PuIeDZZycwsqZUH6FZ6M2+zD1aTbfZlX29Zqd6/SdqmGMhMpY0BOivItIG6HLrjKa0y4Dek40F4uqxZIZzstpSxpPm83YE57dJLycluVggnuw1eW+MI8xJ+SNvqG5zsZoUYTLKn3lVi/aozDbcri8qSuox9jQk42c0KkX2yt/XN3/W3q3sm06WYyOJ6rMfJblaIbJN9ld9eN3ndPjiZmsmp3nIqyyxOdrNCuLGbFSK7bnxb3fetr7PxPu5K2yqGuL442c0KkU2y9z1w1lbC971c1twQU3waJ7tZIbJJ9llTFtv6Vp31fql+VFEn0ceSGGM1ts/HyW5WiGySfUNO36ZtbG/ntHw20VVvsm9OdrNCZJfsXWv7p5ZjTYkxGvtn5WQ3K4Qbu1khiu/Gb8jpyCnWnmmTp8befd/gZDcrhJN9i1VPL5yjUnYpzTP0zzAlJ7tZIZzsc2xNwlVSIodj3c26b4xJv+gzGuMy1+VkNyuEk30JTUbs2z6lb9tlmqevtPT2+HKc7GaFUJffjpJG9VXcVd3NS84c062tpPeU5noiYuoCOdnNCpFNsg95pHjVOpy1jDmmdRNtHes/9dlVh7TOzbNysku6XNITko5LelbSndXtF0l6XNKJ6u+FqQttZunU6cafBb4SEe8HrgM+L+lq4C7gcETsBQ5X180sU0t34yUdAv62unwsIk5L2gU8GRFXLXhu7W58zbIs/Zy2pZzUMbbufCp16rDk7vysbvxS+9klXQFcAxwBLo2I09WLn5Z0yYzn7Af2L/M+ZpZe7WSX9B7gJ8DfRMQjkt6MiAs23f9GRMzdbq+z663tySddSXFmmz537eXYq+iq7nJer+potOtN0ruAHwAPRcQj1c2vVd13qr+vpyiombWjzmi8gAeA4xHxzU13PQrsq/7fBxxKUSBJc79ZN+7P/dt3CGUcilXq0vX/Tgu78ZI+Cvwz8DPgrermv2Sy3f594HeBl4HPRMR/LXit2n2rtk/w2JUm8wfcjZ9o+2QdKd8vB7O68dlMqhmrefXb9qSaFD/R7VOKRufG/hueLmtWCP/EtUddn3FmCAfVTJ2qQ1jmrjjZzQrhZE/MCWK5crKbFcKN3awQ7sY3UFKXfWhHy9nKA3VOdrNiONmXkDoVhnCkmjq7woaUmnXK2vWx/rviZDcrhJO9hpSJ1WdSjuXnw10ZW8I72c0KkU2yj/0cXauUP1XCO9HbM6SjIjvZzQqRTbIvMi2dcv4WTWmZc8b38ZPOIYzCb1WnTpuccy/HddPJblYIN3azQvTajU818NR2l2lIk0aWkareUtbP5tfoqiucY5e7DU52s0Jkcwy6IR3nO4fjmuVcX11NK85JTgNzPgadWeGy2fWWwzfikCyzndx13aY+F1uT1Oxq0ssQ1l8nu1khskn2IckpVYeQKH2Y9dnktG3dNSe7WSGc7A2Mdf97CfrYn983J7tZIZzsCZSSDCnkMEdhq1K2453sZoVwYzcrhLvx1rpVTlvd1CqDp2PvzjvZzQrhZLfW5LBLcpVj84814Z3sZoWo3dglbZP0tKTHqut7JB2RdELSw5J2tFdMG5KIyCLV55H0/5dZhrAcy1gm2e8Ejm+6fg9wb0TsBd4A7khZMDNLq1Zjl7Qb+GPg/uq6gBuAg9VDHgRubVKQjW/RJhfrx9A/h7oJX2fZcq6Dusl+H/BV4K3q+sXAmxFxtrp+Crhs2hMl7Zd0VNLRRiU1s0YWNnZJnwJej4hjm2+e8tCpX2cRcSAi1iNifdNtraTB0BNmCFJ9dovStA91ylS3d5njelhn19tHgE9Luhk4BzifSdJfIGl7le67gVfbK6aZNbUw2SPi7ojYHRFXALcDP46IzwFPALdVD9sHHGqtlGbWWJP97F8DvizpJJNt+AcWPWFtba3Tbk2OXakcdD0Iunk3V47d962GUMZpIoK1tbWZ9y81gy4ingSerP5/Ebi2QdnMrEPFTJcd6xTIOvrq2ZRY132o+/l6uqxZIXo9I0wfiVNi2nRdz2Or41T111a9bO21+owwZoUbzDb7tG/FVb5xS9x2X+VnnnVfY9nXKVnKdW+VOneymxWil2RPlSg+bnszrr96Nq+DKeqqrfqOCNbX12fe72Q3K0SnjX2ZGXTLzGJaZcaTZ9b9xlBmtg1BzrMFnexmhXBjNytEdrvemnR7PODU3Kw6dJ3O12d3ve57O9nNCtHLdNmuzhBSZ9lyGkDJSZP1Ysx1Oqtemkw6qlNfy0zI8XRZs8L1ss3ubev8+TNq3zI9oBS9JSe7WSGyG4234RrztvoYONnNCuFkt7nqbLs70ZfTV3052c0K4cZuVgh3462WkrvqdU/omDsnu1khskv2sR8jbnMCjHUZbbq+P28nu1khskv2DW0nfNffstO26cbeixm6IWyHL8PJblaI7I8uu0r6De0b2Qk/brl8rk52s0Jke3TZrVKdQ7zrI376KLaWCye7WSGyHY1PIZdtJRuWrntiKd/PZ4QxMzd2s1LUauySLpB0UNLzko5Lul7SRZIel3Si+nthigJNO33Oosus5w7F0Mo7Vm0NpjYZUE6pbrJ/C/hhRLwP+CBwHLgLOBwRe4HD1XUzy9TCxi7pfOAPgQcAIuJ/IuJN4BbgwephDwK3NilIk3RzMloTfe0erdNbTalOsl8JnAG+I+lpSfdLOg+4NCJOA1R/L5n2ZEn7JR2VdPTMmTPJCm5my6nT2LcDHwK+HRHXAL9miS57RByIiPWIWN+5c+eKxRwn90j6k2qbeZUxpjqfexvrRp3Gfgo4FRFHqusHmTT+1yTtqgq2C3g9acnMLKmFjT0ifgm8Iumq6qYbgeeAR4F91W37gEOtlNCsYCkTvu4Muj8HHpK0A3gR+FMmXxTfl3QH8DLwmSQlMrNW1GrsEfFTYNo8vBubvHmJ26slLnOOpn0Oy2y/d/05pjj3nmfQmRXCjd2sEKP+1VsO3G0fjlmfVU5HBG7SnXeymxXCyW62QN9pPs3mMtVNeSe7WSHc2M0Gru7EGzd2s0J4m92sY1u3sbsaE3CymxXCyZ5YX9/aZovWNSe7WSF6TfacZiY1MW8/p8/jZlv1tS442c0K4cZuVohsGvvYT4A49uWz/GXT2M2sXdntehv7gNa8dB/rMlsenOxmhcgu2TcMKeFX+bmhWdec7GaFyDbZNwwp4Vc15mXrwrTelOv0nZzsZoXIPtmHZpkDAnadPrPKNLQU9LjIapzsZoUYTLKPZdu9y/LXTcDcf5DUdZKP9WfKTnazQrixmxViMN34oUlxIr5VNH2/rjaX2qqXNsqd+2ZOXU52s0I42TvWVkqkTsqSd28NOb3ncbKbFcLJ3qOU28d9jRGsqs4ZUy0tJ7tZIWo1dklfkvSspGckfVfSOZL2SDoi6YSkhyXtWPbNN85RtflSopSHrJpWp6tc+pJDGYYqIlhbW5t5/8LGLuky4IvAekR8ANgG3A7cA9wbEXuBN4A7kpTYzFpRtxu/HfhtSduBc4HTwA3Awer+B4Fb677pvG/ukr/VNxJ+66UPOXwOOZShb7PWiVXWk4WNPSJ+AXwDeJlJI/8VcAx4MyLOVg87BVw27fmS9ks6KunomTNnai6imaVWpxt/IXALsAd4L3Ae8MkpD5361RIRByJiPSLWd+7c2aSsZtZAnV1vHwdeiogzAJIeAT4MXCBpe5Xuu4FX2ytm2cYyXXNVQ9ut2ESby1hnm/1l4DpJ52pS6zcCzwFPALdVj9kHHGqniGaWQp1t9iNMBuL+FfhZ9ZwDwNeAL0s6CVwMPJCiQIsGG/oeuKqrrTJ2vew5DZLlVJaUulqna82gi4ivA1/fcvOLwLXJS2Rmrehluuy0aaK5J3XphrLdPJYjGrXB02XNCtHrD2GGcqCFZbSdfH0vax/v39cBQPqu69Sc7GaFGPRPXHP85k25bZvj8o1R3+MQXY1dOdnNCqGO99ku/WZjO4b3MvU9bVnHuj25jBRn25n3Gn3WbZP2uKlXOXUBnOxmhXBjNytE9t34MfMgXhqp1+Ec6nSVZXI33syAge96G7qhTEHN3byBzFXkMCjcxrrhZDcrhJM9A0749FZJ41n1P+32rtI+5brhZDcrhJM9I074fm1N60UHUZn2nLakWDec7GaFcLJnqM63eA77gseuzucwpOnLTnazQrixmxXCjd0sgSEc8diN3awQHqAzW2CZ6bhtDdh5Uo2Z1eZkN1vBogk4qc7P5x/CmNnSnOxmCcybgJPLKL2T3awQTvaMdfXDGE/LTaetz2yZH+nM4mQ3K4STvWB10qGr7U33IKabVS+r1JeT3awQbuxmhXA3fgBynHqZ2pB+F96FNurByW5WCCe7vY2Pg9dM03prs2fjZDcrRNfJ/p/Ar6u/Q/A7DKesULO8ddKjg23nYuu2Zb83645OT+wIIOloRKx3+qYrGlJZYVjlHVJZYXjlncbdeLNCuLGbFaKPxn6gh/dc1ZDKCsMq75DKCsMr7zt0vs1uZv1wN96sEG7sZoXorLFLuknSC5JOSrqrq/etS9Llkp6QdFzSs5LurG6/SNLjkk5Ufy/su6wbJG2T9LSkx6rreyQdqcr6sKQdfZdxg6QLJB2U9HxVx9fnWreSvlStA89I+q6kc3Ku27o6aeyStgF/B3wSuBr4rKSru3jvJZwFvhIR7weuAz5flfEu4HBE7AUOV9dzcSdwfNP1e4B7q7K+AdzRS6mm+xbww4h4H/BBJuXOrm4lXQZ8EViPiA8A24Dbybtu69k4bU2bF+B64Eebrt8N3N3Fezco8yHgE8ALwK7qtl3AC32XrSrLbiYN5AbgMUBMZnhtn1bnPZf1fOAlqgHhTbdnV7fAZcArwEVMZpg+BvxRrnW7zKWrbvxGBW44Vd2WJUlXANcAR4BLI+I0QPX3kv5K9jb3AV8F3qquXwy8GRFnq+s51fGVwBngO9Vmx/2SziPDuo2IXwDfAF4GTgO/Ao6Rb93W1lVjnzZhOMt9fpLeA/wA+IuI+O++yzONpE8Br0fEsc03T3loLnW8HfgQ8O2IuIbJ7yN677JPU40b3ALsAd4LnMdk83OrXOq2tq4a+yng8k3XdwOvdvTetUl6F5OG/lBEPFLd/JqkXdX9u4DX+yrfJh8BPi3pP4DvMenK3wdcIGnjx0051fEp4FREHKmuH2TS+HOs248DL0XEmYj4X+AR4MPkW7e1ddXYnwL2ViOaO5gMeDza0XvXosnPlR4AjkfENzfd9Siwr/p/H5Nt+V5FxN0RsTsirmBSlz+OiM8BTwC3VQ/LoqwAEfFL4BVJV1U33Qg8R4Z1y6T7fp2kc6t1YqOsWdbtUjoc+LgZ+Dnw78Bf9T1YMaV8H2XSNfs34KfV5WYm28KHgRPV34v6LuuWcn8MeKz6/0rgX4CTwD8C7+67fJvK+fvA0ap+/wm4MNe6Bf4aeB54BvgH4N05123di6fLmhXCM+jMCuHGblYIN3azQrixmxXCjd2sEG7sZoVwYzcrxP8BRTDl9toNZNAAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib\n",
    "import pylab\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import matplotlib.animation as ani\n",
    "\n",
    "n = 100 # size of space: n x n\n",
    "p = 0.35\n",
    "# probability of initially panicky individuals\n",
    "\n",
    "def initialize():\n",
    "    global config, nextconfig\n",
    "    config = pylab.zeros([n, n])\n",
    "    for x in range(n):\n",
    "        for y in range(n):\n",
    "            config[x, y] = 1 if pylab.random() < p else 0\n",
    "    nextconfig = pylab.zeros([n, n])\n",
    "\n",
    "\n",
    "def observe():\n",
    "    global config, nextconfig\n",
    "    pylab.cla()\n",
    "    pylab.imshow(config, vmin = 0, vmax = 1, cmap = pylab.cm.binary)\n",
    "\n",
    "def update():\n",
    "    global config, nextconfig\n",
    "    for x in range(n):\n",
    "        for y in range(n):\n",
    "            count = 0\n",
    "            for dx in [-1, 0, 1]:\n",
    "                for dy in [-1, 0, 1]:\n",
    "                    count += config[(x + dx) % n, (y + dy) % n]\n",
    "            nextconfig[x, y] = 1 if count >= 4 else 0\n",
    "    config, nextconfig = nextconfig, config\n",
    "\n",
    "%matplotlib inline\n",
    "\n",
    "initialize()\n",
    "for j in range(10):\n",
    "    update()\n",
    "    observe()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Around which value of the initial density parameter does the phase transition (the whole space getting filled with panic) happen?\n",
    "\n",
    "Around the values of $p = 0.35$ to $p = 0.45$. This is known as the critical density of the system that divides a calm state from a fully-panic state."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}