happiness threshold-Agent_based.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import math\n",
    "import random\n",
    "import matplotlib.pyplot as plt\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.close()\n",
    "\n",
    "N=40\n",
    "M=40\n",
    "e=20 #empty percent\n",
    "n_race=3 #number of races\n",
    "empties=[]\n",
    "\n",
    "agents=np.random.random((N,M))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def race_up(n,agents):\n",
    "    \n",
    "    for i in range(N):\n",
    "        for j in range(M):\n",
    "            agents[i,j]=math.floor(n*agents[i,j])+1      \n",
    "    return agents\n",
    "\n",
    "def empty_up(e,agents):\n",
    "    \n",
    "    for i in range(int(round(e*N*M/100.0))):\n",
    "        a,b=random.choice(range(N)),random.choice(range(M))\n",
    "        empties.append([a,b])\n",
    "        agents[a,b]=0\n",
    "        \n",
    "    return agents,empties"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAsUAAAFTCAYAAADRHHHbAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xu8HWV97/Hvz+xwCZBwDSdAEBUQsKfEsgup2taCGsmxDR60x0sRtJ6kHi+10dZLq5haD9Wj0kNbBRQLFa0XRPHY1Agobe2RW0pALmqUi8COIJFEEA828Dt/zEQW271+v7XWs2etvTOf9+u1X1l7zZqZZ5555pkns2e+y9xdAAAAQJs9YdQFAAAAAEaNQTEAAABaj0ExAAAAWo9BMQAAAFqPQTEAAABaj0ExAAAAWo9BMYCRMjM3s0NHXY6pmNmDZvbkUZcDANA8BsXALGFmt5vZc0ZdjmEysyvM7NWT3nu2md01jPW7++7ufusw1pUxsxVmtsHMfmxm95nZ5WZ2yKjL1ckqf2Fmd5vZ1nr/Pa1j+s5m9rF6G35gZquDZZ1qZuvrz95lZu8zs7GO6Rea2aZ6+ncmt5OOz51e/8frOR3vvc/M7qznvcPM/rRj2q/X/xnq/HEzO7mefpqZPTJp+rMLqw7ADMCgGMDjdA48MDPUV9L/XtKbJC2Q9CRJH5L06DSuw8ys9JzwYkmvkvTrkvaW9A1JH++Y/i5Jh0l6oqTfkvQnZvb8LsuaJ+mNkvaVdJykEyS9uWP6GZIOcff5kn5H0l+Y2TGTtukpkl4kadOkZZ8n6Yh63mdIepmZ/VdJcvd/rf8ztLu77y7pBZIelPTljvm/0fkZd78irhYAswGDYmAWqq9W/ZuZnWlmW8zsVjN7Rv3+nWZ2r5md2vH5883sbDO71MweMLN/NrMndkx3M3utmW2UtLF+7xlmdk19xe8aM3tG/f5LzOzaSeX5IzP7Yv16ZzN7v5l938zuqde7a8dn/7i+wjdhZq+ahrp4pZndUm/XrWa2qmPas+urjG+vr67ebmYv77NeDu347N+a2T/Wn72qHnRt/+wR9XJ+ZGbfNrPf7Zi23Mxurue728zeXL+/r5l9qd6HPzKzf+0yMF0i6TZ3v9wrD7j759z9+/Vy5tTb+L16HevNbHE9bcr9WE+7wszeY2b/JukhSU82swVmdl69j+626srvnB53x5Mkfd3db3X3RyRdKOmojumvkPRud7/f3W+R9BFJp021IHf/cD1A/Zm73y3pE5Ke2TH9Jnd/ePuv9c9TJi3mbyS9RdLPJi372+7+k463HpXU7RaeUyVdNOnzAHZADIqB2es4STdI2kfSJyV9StKvqjq5/56kvzGz3Ts+/3JJ71Z15W2DqkFGp5PqZR5lZntL+kdJZ9XL/6CkfzSzfSR9UdJTzeywjnlfVpdBkt4r6XBVA7lDJR0o6Z2SVF8VfLOk56q6Yjgdt4Pcq+pq3nxJr5R0ppn9Ssf0/6Rqmw9UNcA518ye2jE9q5dOL5W0RtJekr4r6T2SZGa7SbpUVR0srD/3IXvs1oHzJK1y9z0k/ZKkr9bvv0nSXZL2k7S/pLerGtxN9u+Sjqj/E/Rbk/arJK2u17m8rodXSXoo2Y/bnSJppaQ9JN0h6QJJ21Ttu6dLep6kV9fbeXA9gD+4S/18StKhZna4mc1VVd9frufdS9IBkq7v+Pz1kp72C0uZ2m9IuqnzDTP7kJk9JOlbqq4Gr+2Y9mJJP3P3tZqCmb3VzB5UVf+76bH22/mZeaquNF8wadLT6/9kfcfM3mH8dQXYMbg7P/zwMwt+JN0u6Tn169MkbeyY9p9VDab273hvs6Ql9evzJX2qY9rukh6RtLj+3SUd3zH9FElXT1r/NySdVr++UNI769eHSXpA1Z+7TdJPJD2lY75fU3WVU5I+JukvO6YdXq/70C7bfIWqK5hbOn4elHRXUE9fkPSH9etnqxrg7dYx/TOS3tFHvRza8dmPdnx2uaRv1a//m6R/nVSOcySdXr/+vqRVkuZP+syfS7qk2/ZP+uzSuuw/lPT/6vLsXk/7tqQVU8yT7ccrJP15x7T9JT0sadeO914q6Ws9ttGdJP3vut62SbpN0pPqaYvr93fp+PxzJd3ew3JfqWrwuu8U0+ZIepakP5M0t2M/buxY9+2qj51J85qqgf8aSXt0qb/bJFnHe09WdUX8CaqOu5slvW26jnN++OFndD9cKQZmr3s6Xv9Uktx98nudVxTv3P7C3R+U9CNVV+5+YXr9/h2T1neHqqutUnVV7aX165dJ+oK7P6Tqiuc8SevrK4pbVF0p3K9juZ3rmbyOqbzB3ffc/qPqqvDPmdmJZnZlffvBFlWD1X07PnK/P/5P33eoy3Z3qZdOP+h4/ZAeq98nSjpu+zbX5Xi5qqvUknRyXa476ls0fq1+/3+puuL8lfrWj7d2qwR3v9Ldf9fd91N1z+5vSNr+gNhiSd+bYrZsP0qP3x9PlDRX0qaO7ThH1dXvXpyu6q8ViyXtomqw+dX6iuuD9Wfmd3x+vqr/UHVlZidJ+ktJJ7r7fZOnu/sj7v51SQdJek399hpJH3f326Jle+U6VcfKmik+cqqkv3d375jnVne/zd0fdfdvqvqPzYui9QCYHRgUA+2xePuL+s/ve0ua6Jje+Wf7CVUDpE4HS7q7fv0VSfua2RJVg+Ptf3q+T9UA42kdA9kFXj2wJFV/4l48aZkDM7OdJX1O0vtVXSXfU9Wf0K3jY3vVtzd0rrNzu7N66cWdkv65c/Du1QNYr5Ekd7/G3VeoGlx+QdUVX3l1b/Cb3P3Jkn5b0mozOyFbmbtfI+liVbdibF//5PtppXw/So/f73equlK8b8d2zHf3Xm9xOFrSp939Lnff5u7nq7rV5Ch3v1/V/j960udv+sXFVOrbbT4i6bfrAWhkTI/VwQmS3mBVwsUPVO3jz5jZW3qYd/u6F6v6S8PfJ+t1Pb69AZilGBQD7bHczJ5lZjupuof2Kne/s8tn10o63MxeZmZjZvbfVD0w9SVJcvdtki5SdaVzb1X308rdH1U1iDnTzBZKkpkdaGbL6uV+RtJpZnZUffXw9MJt2knSzqpuKdhmZiequgd2sjVmtpOZ/bqqK82f7ZjWT7108yVV9XWKmc2tf37VzI6s1/tyM1vg7v8h6ceqbtGQmb3AzA41M+t4/5HJC6/L99876vQIVYkLV9Yf+aikd5vZYVb55fq+4XA/Tubum1T9h+cDZjbfzJ5gZk8xs9/ssR6ukfRiM9u/nvcUVVeev1tP/3tJf2Zme9Xb8N9V3QbyC8zseFX3d5/s7ldPmrbQqgc+d7fqIcNlqv5ztv1e7RNU/YdhSf0zoer2lb+ty7WqLoOZ2bGSXivp8klFOEXS/3X3x12Br/8ysX/9+ghJ71B1CwyAWY5BMdAen1Q1CP2RpGNU/Xl/Su6+WdXg8U2q7k3+E0kvmPTn60+qelDus/Ugebu3qBoEXWlmP5Z0maSn1sv9J0l/pWrw8l09NogZiLs/IOkNqgbb96u6leOLkz72g3rahKpB1h+4+7cmbUdP9ZKU43mSXlKv5weqHjjcuf7IKZJur+vjD1Q9CClV92NfpurWgm9I+pBPHe+1RdUg+Jv1w2FflvR5Se+rp39QVR18RdXg+jxV9wX3sh8ne4Wq/2zcrKreLpK0SPr5g3YPBg/avVfVw3Mb6jL/kapB7ZZ6+umqbvO4Q9I/S/pf7r79QbzJy36Hqvi5tfZYHvA/1dNc1a0Sd9VlfL+kN7r7JVLVft39B9t/VP1H4/769hhJemFdjgdU3R//1/XP5HqY/ICdVA24bzCzn6j6T8fFkv5nl/oAMItYx61SAHZQZna+qofT/mzUZRkmq75U4UJ3P6jL9PPVwnoBAPwirhQDAACg9RgUAwAAoPW4fQIAAACtx5ViAAAAtB6DYgAAALQeg2IAAAC0HoNiAAAAtB6DYgAAALQeg2IAAAC0HoNizApm9nYz++h0f7aHZbmZHTodywIA9K7jq7/njLosaAdyijESZnaapDdJeoqkH0v6vKS3ufuWUZZrMjNzSYe5+3dHXRYAGAUzu13Sq939so73Tqvfe9aoygVMN64UY+jM7E2S3ivpjyUtkLRU0hMlXWpmO03x+bHhlhAAMJtx3sAgGBRjqMxsvqQ1kl7v7l929/9w99sl/a6qgfHvmdm7zOwiM7vQzH4s6bT6vQs7lvMKM7vDzDab2TvM7HYze0497eefNbND6lsgTjWz75vZfWb2px3LOdbMvmFmW8xsk5n9zVQDcwDA1MzsrWb2PTN7wMxuNrMXdkw7zcz+zcz+2sy2mtm3zOyEjulXmNkZZnZ1Pf0SM9u7nra9/x6rf19gZufVffXdZvYX22+t6FjPmWb2I0nvGm4tYEfAoBjD9gxJu0i6uPNNd39Q0j9Jem791gpJF0naU9InOj9rZkdJ+pCkl0tapOpq84HJep8l6amSTpD0TjM7sn7/EUl/JGlfSb9WT/8fA2wXALTV9yT9uqq+eI2kC81sUcf04yTdqqqfPV3SxdsHvrVXSHqVpAMkbZN0Vpf1XFBPP1TS0yU9T9Krp1jPQknvKdsktBGDYgzbvpLuc/dtU0zbVE+XpG+4+xfc/VF3/+mkz71I0v9x96+7+88kvVNSdnP8Gnf/qbtfL+l6SUdLkruvd/cr3X1bfcX6HEm/OdimAcAO6wv1X9S2mNkWVRcmJEnu/ll3n6j7609L2ijp2I5575X0V/VfBj8t6duS/kvH9I+7+43u/hNJ75D0u5MfrjOz/SWdKOmN7v4Td79X0pmSXtLxsQl3/+u6P5983gBS3HODYbtP0r5mNjbFwHhRPV2S7gyWcUDndHd/yMw2J+v9QcfrhyTtLklmdrikD0oalzRP1TGxPtsIAGiZk6Z60K5+/QpJqyUdUk/eXY9d4JCku/3xT/Xfoaof3+7OSdPmTppfqm6vmytpk5ltf+8Jk+aNzhtAiivFGLZvSHpY0n/tfNPMdlN1FeDy+q3oyu8mSQd1zLurpH0GLM+HJX1LVcLEfElvl2TxLAAASTKzJ0r6iKTXSdrH3feUdKMe348eaB0jWUkHS5ro+H3xpGn/occukGx3p6pzx77uvmf9M9/dn9bxGeK0UIRBMYbK3bequufsr83s+WY218wOkfRZSXdJ+ngPi7lI0m+b2TPqh+LWaPCB7B6qIuEeNLMjJL1mwOUAQBvtpmow+kNJMrNXSvqlSZ9ZKOkNdX//YklHSlrbMf33zOwoM5sn6c8lXeTuj3QuwN03SfqKpA+Y2Xwze4KZPcXMuN0N04ZBMYbO3d+n6ors+1UNSK9SdRXgBHd/uIf5b5L0ekmfUnXV+AFV96yl807hzZJeVi/jI5I+PcAyAKCV3P1mSR9Q9VfAeyT9Z0n/NuljV0k6TNXV3/dIepG7d97y9nFJ56u6zW0XSW/osrpXSNpJ0s2S7ld1gWRRl88CfePLOzDrmdnukraougXitlGXBwBQyb7kw8yukHShu0/Lt5ACJbhSjFnJzH7bzObV9yK/X9I3Jd0+2lIBAIDZqudBsZnNMbPrzOxL9e9PMrOrzGyjmX2aLzzAkK1Q9aDGhKo/y73E+bMHAAAYUM+3T5jZalWxVfPd/QVm9hlJF7v7p8zsbEnXu/uHGywrAAAA0IierhSb2UGqgrY/Wv9uko5XdZO7VH3LzElNFBAAAABoWq+3T/yVpD+R9Gj9+z6StnR8+cJdyr9mFwAAAJiR0m+0M7MXSLrX3deb2bO3vz3FR6e8D8PMVkpaKUlP2GXuMbssnvo7FnzL4F+ud9SiHw487y0/3WvgeSVpzsZBUsAqjxy2czg9qpNsm6Ptyup67Ic/Cacf/ssPhdMHdfOm/cLptudU3wzdmybbV0m5s/azbb/dwulR2UrKdeSu94fzZsuO2lC2TYMuV2qubUrS+hsevs/d4w3fAXT22XM055h5mj/iEvUnawMlfX52XHznhnldp2XlKjmmmpSVu2SbB11uL6J+pmTMIDXXhkq3OVKyH0tlY51IyRhLkh7Q/X312+k9xWZ2hqRTJG1TlR84X9LnJS2T9J/cfZuZ/Zqkd7n7smhZux2+yJ921iunnLbtksHPNdeePvitzEs3vGjgeSVpwfKNA8+7de1h4fSoTrJtjrYrq+t9zvm/4fR1E9eH0wc1vib+3oyxFYN3ZE22r5JyZ+1n86pnhNOjspWU68olF3Wd1suyozaUbdOgy5Waa5uSNGfRxvXuPt7YCmag+ba3H2cnjLoYfcnaQEmfnx0Xyw44uuu0rFwlx1STsnKXbPOgy+1F1M+UjBmk5tpQ6TZHSvZjqWysEykZY0nSZX5RX/12evuEu7/N3Q9y90MkvUTSV9395ZK+Jml7yzhV0iUDlBcAAAAYuZKc4rdIWm1m31V1j/F501MkAAAAYLj6utHS3a+QdEX9+lZJx05/kQAAAIDh4hvtAAAA0HoMigEAANB6PX+j3XSYt3CxH3Hy6imnZU+DljwZWfIEbPZEcFTurMwlTyNnSQzRE66liRslsqe3I02WO0qnyOq6JM0jS2LInjiP5s/KHT3Vm5WrpP1lx0VJOkWmJAmk36eYdwRR+kSTbbc0HSAyqif8m0xTyET7osnElpJzYHbuLTHK9Imon2lyX2RK6rvkWM/mjZIrekmmmPb0CQAAAGBHx6AYAAAArcegGAAAAK3HoBgAAACtx6AYAAAArcegGAAAAK3X1zfalbI9t3WNRMoiTrat6h57lUV6NCkq95UTg8ePZUqizbIYkygCJVt3FsGzTN2nNxlHUxKbtnVFXB9pnOA53bc5izbbrMGj0Ur285gGj1yTyuKnom0qafeZpWuTmKUTG1v1DqnJaL0So4q9Kl3vsgMKln16tNyyqLho3aOMGJupon53fE18zDQZVViiJJ4xnXd5cG5OxiqS+u63uVIMAACA1mNQDAAAgNZjUAwAAIDWY1AMAACA1mNQDAAAgNZjUAwAAIDWG2ok25G73t81UimLZIuiSMb1mnDepRu6x3ZksVVja8uiqSLja+Jyh1FzQcROtuwmtymL4In2cxQ5JMWxfFIc5ZXF9kWxMNcuKYvBiWJjsva3OdnmaF+Nr4rbV7RdWdtcqvh4jeIIswioaJuyfiITRfNlcUdzitaMyUripaI2VBIpWaok3izrO6PppbFqTcn6kZI2kMVZNrmfS/qoqM/P+vusPrM6GXTe/DzVXPxitOworm1QXCkGAABA6zEoBgAAQOsxKAYAAEDrMSgGAABA6zEoBgAAQOsxKAYAAEDrmbsPbWXzFi72I05ePeW0kiiRLHoljCdL1htFOEnNxQpJebRQJKqT0vWGMWANRsZk+6Jk2ZEsjiYT1WeT7SuLBippI1kET5PHc6Qk5jDbpuvOXr3e3ccHKtgsNX70Ln71uoOnnFYaA9ZUxFgWbVbatiNZ9GOJbLt2NCXRZlJZFGZJXTfZvjJNnQMz2XmsKb2cH+cs2thXv82VYgAAALReOig2s13M7Gozu97MbjKzNfX755vZbWa2of5Z0nxxAQAAgOnXyzfaPSzpeHd/0MzmSvq6mf1TPe2P3b25r40BAAAAhiAdFHt10/GD9a9z65/h3YgMAAAANKyne4rNbI6ZbZB0r6RL3f2qetJ7zOwGMzvTzHZurJQAAABAg/pKnzCzPSV9XtLrJW2W9ANJO0k6V9L33P3Pp5hnpaSVkrTTwvnHHH3BawcqaPS0aElaQiZ7krTJJ4KjJ+lLUgkyJU/PZuUa1TZloieds/aTPSUdaXLZpeuOlCSYlKR5lB7rUbmzY7nfp5hnq84+++ADx4657donNbKekn1Rkk6RPf1fklBSMm/W/5XUV4kmk0AiTfZ9M1WTKQ4lyRRZn91kYkvpueQyv6i59Al33yLpCknPd/dNXnlY0t9JOrbLPOe6+7i7j48tmNfP6gAAQ9bZZ++3z5xRFwcAhqaX9In96ivEMrNdJT1H0rfMbFH9nkk6SdKNTRYUAAAAaEov6ROLJF1gZnNUDaI/4+5fMrOvmtl+kkzSBkl/0GA5AQAAgMb0kj5xg6SnT/H+8Y2UCAAAABgyvtEOAAAArcegGAAAAK3Xyz3F0+bIXe/vGpkURXVJ0rqJ7lFLSzfEMU3ZsiP7KI4LiWJjsniVLCIliujJ4mqidWfrvTKo60r3mJ0svieqz3HF+2lUcW/ZNmX11WRkW7TN2X4uOS6unYjrMzoms/ieqNwl9SFJWtV9UtaPSGck03c837lhXtf2XxoxVhLXFc1b0q4zWXxUVCclkVhSvM0lMZrl54PBNRm7VhL9mMmiISNRHzZ+SXNtNxPVV0kUYTZ/SV02gSvFAAAAaD0GxQAAAGg9BsUAAABoPQbFAAAAaD0GxQAAAGg9BsUAAABovaFGskWySI+lKwaPPouigUqibCRJl3SflEXdZJExS9d23+Zs3s2rutdJFmuVCaPPkqiuSB6JNbgspilqI+OrspicwaOlsvaXGVs7eMxT1D6z9lUUn7cknDVeb7Ifs2Mu6iuyfqSNDv/lh7RuXbf2Hbf7kradxTRFfVh2Lkmj4E6PJ0eiPiyNCJuIJ4dtP4gazKTtvuB4LYlcK41UK4ntK4nPKzm/lsSOSmXljsY6JZFr2fybVTBvNj6TpLP72x9cKQYAAEDrMSgGAABA6zEoBgAAQOsxKAYAAEDrMSgGAABA6zEoBgAAQOsNNZLtOzfM6xrTUxLBM35JHFNSEiE2vmbwuJB1p8fRP1nU17VLupctimuTpDF1j2bJYnJKoqmyGKYoJidbb7afs9iYcNlBxFhJzI0Ub9fWtfG8JbFEJfsxivSTyqKDSuqzNLIojN5L5m2jW366V9f+ojQyK1IaG1ki6sOyOLew3EnkWiY6brJjvSTWqskYsKgNZWOCNNI0OM+V9umjap9N1bWUnIvOGXi1knqMTmtg3kFwpRgAAACtx6AYAAAArcegGAAAAK3HoBgAAACtx6AYAAAArcegGAAAAK031Ei2Rw7bWVvPmjpmJY2PWtJ9UhbTFCmOmwmiSrJ4siwOrkQUGZOVqySOa+mKOO5t6YbuMTtZpFpJNEsWpSRl07vLIu4iWbtP90UU26fBo/dKo89K4qPC9SaxfFm5i2Kazo4nt00WmVUi6iekOF4qP9Zj0fxZuy+JPsv6vysngj59edanDzfWahiabH8l45GSWNKSc4lUFnEXKW0/JXGp0bpLltsNV4oBAADQeumg2Mx2MbOrzex6M7vJzNbU7z/JzK4ys41m9mkz26n54gIAAADTr5crxQ9LOt7dj1b1R4Pnm9lSSe+VdKa7Hybpfkm/31wxAQAAgOakg2KvPFj/Orf+cUnHS9p+k9MFkk5qpIQAAABAw3q6p9jM5pjZBkn3SrpU0vckbXH3bfVH7pJ0YDNFBAAAAJrVU/qEuz8iaYmZ7Snp85KOnOpjU81rZislrZSknRbO77qOkgSJTPTEcMnT6lL8ZGT2NPv4mviJzmj+6ClTSVq6tnu5FyieN3uiM0qYyMoVyZ4az55I36zu9ZnNGz1tnO3HKOlDUvikcvY0eyZqn1m5snSKEtG6s/WG8ybHY0l9Zn3BnIGXPLt09tm7aF7XYzo7XrOn8KOn4bN+JH6Svix9IlLytHvWvpadE9dXVJ/ZvojONU2ee1MT0cTsfNBc6k8mHFM0mCqVKUmYKKmTkuSoNPXnkkFKNLi+0ifcfYukKyQtlbSnmW0fVB+kLs3b3c9193F3Hx9bMK+krACAhnX22XO186iLAwBD00v6xH71FWKZ2a6SniPpFklfk35+yedUDX08DwAAAEyPXm6fWCTpAjObo2oQ/Rl3/5KZ3SzpU2b2F5Kuk3Reg+UEAAAAGpMOit39BklPn+L9WyUd20ShAAAAgGHiG+0AAADQegyKAQAA0HrmPmWSWiPm295+nJ0w5bSSKJFMFPFUGokVRZFEMSTZvJkml53F+0SymJwoaqlkvVIeARWJ1l0S/ZPJoqdKYq9K2kg2byaK2cmif5qat5f5I9edvXq9u48PvIBZKOqzMyV9elkkWyyNUCwwqr6itO9sSlN9ci9Kz+2DajLibpQxdJGSSLYm1ytJl/lFffXbXCkGAABA6zEoBgAAQOsxKAYAAEDrMSgGAABA6zEoBgAAQOsxKAYAAEDrDTWSbbfDF/nTznrlQPOWxKplMU2RLMYkigQZZaRMtM2lMWCRLIKnJJqlJG6ryaiuTJMRY5GS6KnS/dhU+8siiZo8Xucs2kgkW4eSNiDl7WBQpfGe0XGTtb8m496aijfLlpsdF9F5Kosni9bd5PmzJLJUKjuPNRnZVjJmaDJeNpq/6Ti3fqM0uVIMAACA1mNQDAAAgNZjUAwAAIDWY1AMAACA1mNQDAAAgNZjUAwAAIDWGxvmyuZsfLhrDE9J/EoW/dNknNu60wePj8rie0riuCJZxMnSDXGkURh9tjYu85gGjyfL6mvZ8iDeJ9hPkrRU8b6KFJV7SbzsrH2WxPuEUUoT2XIL4pImBp8Vw7Vtv920+eSySKQdSZMxh5mo3x5l9FkYqXV60aJDJfFjmdIYsKakbWhV90lNxqZlwrHMOfG8pevuF1eKAQAA0HoMigEAANB6DIoBAADQegyKAQAA0HoMigEAANB6DIoBAADQekONZIvifcbXxBEoJdFTWTRaiSgWJosSWbo2Lle3+DpJ2ro2jk2LooNK4scyUZmlslihTEm0UEldZ9F5URvJ5h1Vu8+if7K6jtadRdhFsmNq3UQcmTWu7vsiiyKUzkim73iOWvRDXd2lDZb2q9lx1ZSs/Y1fEkQVFvZRJUrWHdV11u6zGLrSSLdRKI1ci/qhrF2X9I37aLjxZNsV19fy7uUuWXZPsbUn9rdMrhQDAACg9dJBsZktNrOvmdktZnaTmf1h/f67zOxuM9tQ/yxvvrgAAADA9Ovl9oltkt7k7v9uZntIWm9ml9bTznT39zdXPAAAAKB56aDY3TdJ2lS/fsDMbpF0YNMFAwAAAIalr3uKzewQSU+XdFX91uvM7AYz+5iZ7TXNZQMAAACGwty9tw+a7S7pnyW9x90vNrP9Jd0nySW9W9Iid3/VFPOtlLRSknbRvGMbY05zAAAcjElEQVSe1eXW4+wJxJ6eMuwieqKz9GniKFkgkz1J31RSQ/bUeEmCRLafsieZI1m5S/ZzyX7MtrmpFBEpTokoSYgoFe2LknZfKmoHWRu47uzV6919fLrLNNP02mdn7StLMClJnwjTdQr6iVKjSqdo8ljOlPTpJf1Xpsk+PZImSBQkV5SuO1JSrpJtzvr70vq65sQz+uq3e7pSbGZzJX1O0ifc/WJJcvd73P0Rd39U0kckHTvVvO5+rruPu/v4XO3ca7kAACNAnw2grXpJnzBJ50m6xd0/2PH+oo6PvVDSjdNfPAAAAKB5vaRPPFPSKZK+aWYb6vfeLumlZrZE1e0Tt0ta1UgJAQAAgIb1kj7xdUk2xaS1018cAAAAYPj4RjsAAAC0HoNiAAAAtF7PkWzTYd7CxX7EyaunnFYSa1US4ZTFQ5XEWjUZ/VMSfVYaWVQSaxXVd2kET2SU5YrqO4szKil3yXGRRUuVtKEmowizcpXER81ZtLEVkWydxo/exa9ed/DQ19tknFumqX67yfjPkgixLIKzqei8TNYGmjxfjDLirkRTkWyZJmM0I72U+TK/aPoj2QAAAIAdGYNiAAAAtB6DYgAAALQeg2IAAAC0HoNiAAAAtB6DYgAAALTeUCPZonifLHoqUhJPlsW+ZKJYmNJYoSgqJ5s32uaSupbiaKGSSKwmo+JKlJarJC4pizSKypZFLUVtd5RtpKn4RSneF9l+vO7s1a2LZItiNEuPt6h/zOK2onlnY1xb00r6oFH2b5Hs/FoS7ZgZVWTbqNpfaXRtSZxgSbynRCQbAAAA0DcGxQAAAGg9BsUAAABoPQbFAAAAaD0GxQAAAGg9BsUAAABovaFGss23vf04O2HKaVmMThTLURKvMspItpLYoUxUnyVRXVJzMWClMTdRGxllJFFJPF6T5W5qvdm6s+O1qXZf6poTz2hdJFsUo5kpOZ5LoroyJeXKjqleIqIGFUWMlcRaNXmsZ6J1Z+eSJmNJS8YUTbavUcnqI6vP6LjIYjZLjyki2QAAAIA+MSgGAABA6zEoBgAAQOsxKAYAAEDrMSgGAABA6zEoBgAAQOuNjboAvYqiW5qMcMoiUsbXBHEiq+J1Lzsgnp5Fo0WWbui+Xdk2R/NKcVROFq8idd+mLIIn133+rSsG38+lkUWRNG7mnHhy1EbGLxk8Jierrzx6b/BIrZLoqUy0n5uM05qtbvnpXsVRiTNNSdxbdkzl/V8zSqK8ZmoMWGm8YjQuuHZi8Mi1TJPtK9NUPO24yspVovh8cGJ/6+NKMQAAAFovHRSb2WIz+5qZ3WJmN5nZH9bv721ml5rZxvrfvZovLgAAADD9erlSvE3Sm9z9SElLJb3WzI6S9FZJl7v7YZIur38HAAAAZp10UOzum9z93+vXD0i6RdKBklZIuqD+2AWSTmqqkAAAAECT+rqn2MwOkfR0SVdJ2t/dN0nVwFnSwukuHAAAADAMPadPmNnukj4n6Y3u/mMz63W+lZJWStJOC+dr6wVTP1GaP83e/Sn7ZYrTJ0qeYs2eSC9JiMie6IxSIEqeGC55QlWSNNF9UpaoEW1z9kTw+JrB0xTyZIsgUUNlT5RHTxRn25w9+R/V99ja+MncaLv2WR63+/yYClIxkv0Ytc8saSZ7+p+EidzkPntQ2TEX7avsKfySFJKSdICs78zadqTJbYqOm1ElZkhlaTBZHxRND1OjCpWMGYrPzQXC+rykbNlhSlKyL6I6KTneuunpSrGZzVU1IP6Eu19cv32PmS2qpy+SdO9U87r7ue4+7u7jYwvmTUeZAQANoc8G0Fa9pE+YpPMk3eLuH+yY9EVJp9avT1Xx/yUAAACA0ejl9olnSjpF0jfNbEP93tsl/aWkz5jZ70v6vqQXN1NEAAAAoFnpoNjdvy6p2w3EJ0xvcQAAAIDh4xvtAAAA0HoMigEAANB6PUeyNa0k2iyLlBlT96ibLMpm6do4EiuKE0njn1bFk6O4mpJooCzmKxPV2SjjfaJImZIIu9IYr6hcWcTYlRNx+xxfFUTSJI++RhFQW1cMHmNYKmq7WcxcyeO+aczciYMve7Y6ctf7B44wC9vmCGX9X1ORbdmxvu70+BwYlTtbdoksKi6qryYis6ZDadxbtM1LVwweo5kpGSdlSs6RmagdjDKGbipcKQYAAEDrMSgGAABA6zEoBgAAQOsxKAYAAEDrMSgGAABA6zEoBgAAQOuZuw9tZfMWLvYjTl497csdZWRMFCeSLbs0FiayYPnGrtOyWJcs3qckdq0ksihbb8m+iObNIpyyKJusfTalpFxZLFVWnyXLjtpBSRuQ4n0ZHTOSdJlftN7dx8MP7WCa6rMzozpmpLLzRZPxUtG6s3NJdNyUxHv2Mn+kNB40Eh3P2bk16wsipefXkmU3td4mo1abjmSbs2hjX/02V4oBAADQegyKAQAA0HoMigEAANB6DIoBAADQegyKAQAA0HoMigEAANB6Q41k2+3wRf60s1457cstiU8pFUW7ZOVqKtpMKou6yaK8SiJUSmKFsjiaaJuzGLASWV1H7SBrAyVRS1msVbafI021gdJlZ0qOi2tOPKN1kWzzbW8/zk6YclppTNOoYtea7AtKlMaDRkYVFZcpKVdJxFhJ3KnUXBsqjXGNlByvpfGxUX01HVFHJBsAAADQJwbFAAAAaD0GxQAAAGg9BsUAAABoPQbFAAAAaD0GxQAAAGi9sWGuzLeMdY2ByuJ5olirLJYjivwojhVS93I3GVmUxZhsW9U9biuPPovjZqLYmJKInTwmJ97PUcTYsuVxfTUbbda97WbLXnd6FkPXvc5GFbkmlUWfRW27NEopqpNsm+cUrXl22rbfbtp88tTHRmlMU1OyuKwmI8RK2n2JJiPXmlRyrC8I+tVe5g+XncW8Tgy86Eaj9yJZGykZJ127ZPD2VzJ+660uV/dVHq4UAwAAoPXSQbGZfczM7jWzGzvee5eZ3W1mG+qf5c0WEwAAAGhOL1eKz5f0/CneP9Pdl9Q/a6e3WAAAAMDwpINid/8XST8aQlkAAACAkSi5p/h1ZnZDfXvFXtNWIgAAAGDIBk2f+LCkd0vy+t8PSHrVVB80s5WSVkrSTgvnd33aPnvaM3tCcdB5lx0Qz5s9dRk9zZ49nZ3ZrGjdcX1F687qsuQJ6uxp0OgJ2Gze8Uvi6VGSQ1yXyb46J5xVm4OkDymr73hflNRnSeJGJmsj8fE8+FPj6VPhiWhf5U8y9/cU82zV2WfP3X001zuyBIlRaTJdIkuiaapOSrcpK3dk64rRJJRkSpKQStIl0rpMz0Xd+/SSsU46ljk9nhy1sZI+vZfUlTln97fMga4Uu/s97v6Iuz8q6SOSjg0+e667j7v7+NiCeYOsDgAwJI/rs3fdbdTFAYChGWhQbGaLOn59oaQbu30WAAAAmOnS2yfM7B8kPVvSvmZ2l6oL5c82syWqbp+4XdKqBssIAAAANCodFLv7S6d4+7wGygIAAACMBN9oBwAAgNZjUAwAAIDWM3cf2srGj97Fr1538NTTmowxCUSRar0II7EOiCOxSmLmsvpqsk6ieJYs5isqVxY5VBIdlG1TSbkyUblLlx21gyxGJ4odKqkvqWy7SuJ7SmLmsnifOYs2rnf38YFXMAtFfXbWvzW5L6J2XzJvJmv3JfFSeQzYaJRsU9YGoj5qptaHVD5uGFRJv1vS7ptUEl3bSx9z3dmr++q3uVIMAACA1mNQDAAAgNZjUAwAAIDWY1AMAACA1mNQDAAAgNZjUAwAAIDWG2ok23zb24+zE6ac1mT8ShQpUxIbJMURKVmUTRbJFkWolEQDlUTBZbLYlybjV6IopizOrSRCrCQerzRirGSbS9pItuySyKIsUqtESYTdZX5R6yLZoj67pP8qVdJGRhX/mSmNQRxUk1GYTSopd9bvZuORaF9l/UjUp48q6lIqGyeVnNczpeOVfqM0uVIMAACA1mNQDAAAgNZjUAwAAIDWY1AMAACA1mNQDAAAgNZjUAwAAIDWG2ok2/jRu/jV6w6ectqyA44O540iUkpirUriU6Rm46MiWX1FsriZrD4jTdZXk5EyURROSQxOJtuPJeUuiXjKtrnkeM2WXRJFWLLN2bzXnb26dZFsUZ9dGss3qrabKYlsazKaquQc2GTkaaTJ2LRRxVVmmuyzS8crgypt103FwPYyDuo3SpMrxQAAAGg9BsUAAABoPQbFAAAAaD0GxQAAAGg9BsUAAABoPQbFAAAAaL2hRrLNW7jYjzh59UDzRpEgJXEfWfxOFgFVoiRyJos2K4kVKtnmbJuifZXF5GSiOJusvkrXHYliibL91OR+LjmmsiicaP6SbS7ZXolItn6V9NklSo7XJiMUS9puaYRddLxmkWvRsrO4rVHFuWVK4t4y2b6INNlnlyiJs8zmbXIsE51remmb15x4BpFsAAAAQD/SQbGZfczM7jWzGzve29vMLjWzjfW/ezVbTAAAAKA5vVwpPl/S8ye991ZJl7v7YZIur38HAAAAZqV0UOzu/yLpR5PeXiHpgvr1BZJOmuZyAQAAAEMz6D3F+7v7Jkmq/104fUUCAAAAhmus6RWY2UpJKyVpp4Xzuz7FmD29GCl5yrTkicxs/vQp6Il48via7k9sZk//76PBn1TOyh2tO3rKtDJ4UkjJ09njGk0aR9PLLml/UZ0s3RC3kQWKj9eSJ6yj9nXtRPykcomlSvqRsxtb9YzS2WcffOBY16fDS9MUIqWpLCVK2m6T56LN6t63jmnwfiTtsy+JJzfZP5YoaX9ZIkd+nuuuJF2nJNmipD6aHMtk4iSkfDwxp8/1DXql+B4zWyRJ9b/3dvugu5/r7uPuPj62YN6AqwMADENnn73fPv2eUgBg9hp0UPxFSafWr09V+n9JAAAAYObqJZLtHyR9Q9JTzewuM/t9SX8p6blmtlHSc+vfAQAAgFkpvafY3V/aZdIJ01wWAAAAYCT4RjsAAAC0HoNiAAAAtF7jkWyd5mx8uGv0WknESRZHk8aJBLKYppJlZ7FqWtV9UhxTEsvitrJyhZFuySOXUWRRHssXT28qJqeoPpL5s2izZYrXHc2fzRu1r5KIxEwWdxTV57ID4mVn+6IklqiNbvnpXl2P2dLjNZIdy1H0VJMRYVkkVkkk22wVHVMlkadZP5HFgJXE9pVEeKbLLoj8K5k3U1RfBVFy6TgoUDIO6oYrxQAAAGg9BsUAAABoPQbFAAAAaD0GxQAAAGg9BsUAAABoPQbFAAAAaL2hRrI9ctjO2npWl8ikJMorksWiRXEhJTEkUlkEz7ZVcTxUVLaSaJZ0m5OoGymIQVkSz1lSX3n8Svfp2Xqj6VdOZLF7cbnGVw2+r9JYokAWT7bP8u7LzubNXLskaLtJ3FEUsZi1gfE1cZRXFBFFXNsv8i1jXetlc9J/ZW23KIazwdi1aNkl/VcWYZcer0F9bl0RzztT231Urs2K20cWzxj1FU1Gm5XEkzVZriaPmexYX3ZO99i1rE+PjrmsDVT6i4bkSjEAAABaj0ExAAAAWo9BMQAAAFqPQTEAAABaj0ExAAAAWo9BMQAAAFrP3H1oK9vt8EX+tLNeOeW0LFZt2QHdIz2yKJsoCqckDkSK42yyCJQsCieaP6uvqNzZektiwDIlMUyZKOomaj9S3A6yeUviy0ragJS3g0gU/5O1gTwarZlll0RiSWXxU9edvXq9u48XFWCWmbdwsR9x8uopp83UyLUs+iwrV7RdpeeL2ajkmGkyBqxJM3WbS88Xgyo5z2RKzq/ZsS5Jl/lFffXbXCkGAABA6zEoBgAAQOsxKAYAAEDrMSgGAABA6zEoBgAAQOsxKAYAAEDrjY26AL2KYnT2Wd5chFgmikDJ4kK2rh18vSXRP1lsy7rT49ihKEIlizuKYtNK4u8yWbmWbghi1ZL9lJUr2ubxS7pHl0nNRu+F5VJcrrC+FMdaZRF20bJLtkmStCSeHJlz9uDzzla257au/UXWT0gzM74s6/82q3tfseyAeNkL1L3PL4luzPQSTdVNVq4mYyFnaoRdyTY3uU1ZuaJ2UNL+Go3CXBXPO6bu29zTNp2Yf6QTV4oBAADQekVXis3sdkkPSHpE0ra2BdsDAABgxzAdt0/8lrvfNw3LAQAAAEaC2ycAAADQeqWDYpf0FTNbb2Yrp6NAAAAAwLCZuw8+s9kB7j5hZgslXSrp9e7+L5M+s1LSSknaRfOOeZYtn3JZWTpAJHrSPVOSliDFT2WWPJUrSeNruicANFmuaL29rHtQUapFL6InUUv2RVaurA1FTwxnaQolbTt7Mjd6Ujnbpqxc0fxZ+8naX1OyJ7uvOfGM9W14bqKzz95p4fxjjr7gtQMtJ2vbWX3PRCUJOKX9ZtQPNZlsUXoeG1STKQ5ZWse6iSxZZXAzNXGjSSXHTaSXPqTffrvoSrG7T9T/3ivp85KOneIz57r7uLuPz9XOJasDADSss88eWzBv1MUBgKEZeFBsZruZ2R7bX0t6nqQbp6tgAAAAwLCUpE/sL+nzZrZ9OZ909y9PS6kAAACAIRp4UOzut0oquwkUAAAAmAGIZAMAAEDrMSgGAABA603HN9r17JHDdtbWs6aOjrl2yeAxTWXxKWXRK2HUyJJ43iyapSTCJ4rRydabxZw0FRW3bVVz8WTZNkf1lcUd7bM8Kdc53SdtXRvPGs0rxdFnY4r3Y7Rd2fG4dEXSdoP5S9p9aVxb1IY2a/BYyDbKYq02J8fzqGQRY1H7HGWMXHyeG12EWNQOSs7NWfvK+uVwP08MUqLHRPF42TaXnJszJdFnJW27qci1TC9xgXP6XCZXigEAANB6DIoBAADQegyKAQAA0HoMigEAANB6DIoBAADQegyKAQAA0HpDjWQ7ctf7u0ZoRBEnkqRV3SeVxJhkkR5ZBFRJbFqmNH6qm9Iyj18yeLmimJ3SeLJIFhkTbdMoY5hKlEQaZcfjAsXLHl81eGxfiSy2L4pLGl9DJNtkUZ+dx1rF0VRRG8vitiJZu1+6dvDzRckxVRILWSo9vwaunEjKVRBvFp7jgnO+lMdGlsjOvWNro3PC4DF0pW2g5NzcZKxaGIUZxIpKzZ4vpsKVYgAAALQeg2IAAAC0HoNiAAAAtB6DYgAAALQeg2IAAAC0HoNiAAAAtJ65+9BWNn70Ln71uoMHmzeISGkyMiuLSImibkqjRqIInyw+JYumikSxVZksdiiKNMrWm8UKRfOXRBKV1IdUFhmYKYnwKTmmSqKpMlHbzsqV1Ue0L7JtuswvWu/u4+GHdjBRn11yPGZK+pGs7WX7OZItu6TtZkqWXVJf2bkmOo+VxIqW9kFNtr+onyk51zQZT1ayL7L1NrnsSC913W+/zZViAAAAtB6DYgAAALQeg2IAAAC0HoNiAAAAtB6DYgAAALQeg2IAAAC03tgwV3bzpv26RnekEWKruk/KImOiaJcs1mXp2iSaZaJ7NMvSDYPHUmWyuJrN6h7tkkWglETKbFsV74t1QX1lsS7XTgxe7izqJqrPbD9mMWBR+8zafVbuZcu7b3MWtRSte+uKeN7SGKdISURPFqUUl2vwqK4dVdRnR31yZfBIrKz9RP1IdrxmbTc6nkv69JLzlBQfF1m7L4lIbDKOq0QWudZkjGvJOTKSjoNOb2S16bqXnRNvb3ZujvZFSbRjL7F7cxalH3kcrhQDAACg9YoGxWb2fDP7tpl918zeOl2FAgAAAIZp4EGxmc2R9LeSTpR0lKSXmtlR01UwAAAAYFhKrhQfK+m77n6ru/9M0qckrZieYgEAAADDUzIoPlDSnR2/31W/BwAAAMwq5u6DzWj2YknL3P3V9e+nSDrW3V8/6XMrJa2sf/0lSTcOXtzG7CvpvlEXYgqUqz+Uqz+Uqz9Pdfc9Rl2IptFnF6Fc/aFc/Zmp5ZJmbtn66rdLItnukrS44/eDJE1M/pC7nyvpXEkys2vdfbxgnY2gXP2hXP2hXP2ZyeUadRmGgT57cJSrP5SrPzO1XNLMLVu//XbJ7RPXSDrMzJ5kZjtJeomkLxYsDwAAABiJga8Uu/s2M3udpHWS5kj6mLvfNG0lAwAAAIak6Bvt3H2tpLV9zHJuyfoaRLn6Q7n6Q7n6Q7lmjpm6zZSrP5SrP5SrfzO1bH2Va+AH7QAAAIAdBV/zDAAAgNYbyqB4Jn8dtJndbmbfNLMNo3y63Mw+Zmb3mtmNHe/tbWaXmtnG+t+9Zki53mVmd9d1tsHMlo+gXIvN7GtmdouZ3WRmf1i/P9I6C8o10jozs13M7Gozu74u15r6/SeZ2VV1fX26fmh2JpTrfDO7raO+lgyzXB3lm2Nm15nZl+rfR1pfwzRT+2367IHKRZ/df7nos/sr147RZ7t7oz+qHsL7nqQnS9pJ0vWSjmp6vX2U73ZJ+86AcvyGpF+RdGPHe++T9Nb69VslvXeGlOtdkt484vpaJOlX6td7SPqOqq8bH2mdBeUaaZ1JMkm716/nSrpK0lJJn5H0kvr9syW9ZoaU63xJLxplG6vLtFrSJyV9qf59pPU1xO2esf02ffZA5aLP7r9c9Nn9lWuH6LOHcaWYr4Pugbv/i6QfTXp7haQL6tcXSDppqIVS13KNnLtvcvd/r18/IOkWVd+oONI6C8o1Ul55sP51bv3jko6XdFH9/ijqq1u5Rs7MDpL0XyR9tP7dNOL6GiL67QR9dn/os/tDn92/6eizhzEonulfB+2SvmJm6636JqeZZH933yRVB66khSMuT6fXmdkN9Z/qhv4nwk5mdoikp6v6H+uMqbNJ5ZJGXGf1n5U2SLpX0qWqrgRucfdt9UdGcmxOLpe7b6+v99T1daaZ7Tzsckn6K0l/IunR+vd9NAPqa0hmcr9Nnz0Y+uz+yiXRZ/dUrh2pzx7GoNimeG9G/K+i9kx3/xVJJ0p6rZn9xqgLNAt8WNJTJC2RtEnSB0ZVEDPbXdLnJL3R3X88qnJMNkW5Rl5n7v6Iuy9R9e2Tx0o6cqqPDbdUv1guM/slSW+TdISkX5W0t6S3DLNMZvYCSfe6+/rOt6f46Ezqy6bTTN5W+uz+jbz/2Y4+u3f02b2brj57GIPinr4OelTcfaL+915Jn1fV8GaKe8xskSTV/9474vJIktz9nvqgeFTSRzSiOjOzuao6sU+4+8X12yOvs6nKNVPqrC7LFklXqLoPbE8z255XPtJjs6Ncz6//pOnu/rCkv9Pw6+uZkn7HzG5XdevA8aquQsyY+mrYjO236bP7N1P6H/rswdBn92Ra+uxhDIpn7NdBm9luZrbH9teSnifpxniuofqipFPr16dKumSEZfm57R1Y7YUaQZ3V9wqdJ+kWd/9gx6SR1lm3co26zsxsPzPbs369q6TnqLp37muSXlR/bBT1NVW5vtVxkjRV94ANtb7c/W3ufpC7H6Kqz/qqu79cI66vIZqR/TZ99mBG3f/UZaDP7q9c9Nl9mLY+O3oKb7p+JC1X9UTn9yT96TDW2WO5nqzqqerrJd00yrJJ+gdVf6L5D1VXaX5f1f0wl0vaWP+79wwp18clfVPSDao6tEUjKNezVP0Z5AZJG+qf5aOus6BcI60zSb8s6bp6/TdKemf9/pMlXS3pu5I+K2nnGVKur9b1daOkC1U/7TyKH0nP1mNPMo+0voa83TOu36bPHrhc9Nn9l4s+u79y7RB9Nt9oBwAAgNbjG+0AAADQegyKAQAA0HoMigEAANB6DIoBAADQegyKAQAA0HoMigEAANB6DIoBAADQegyKAQAA0Hr/HwW4FaOyfaROAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 864x360 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "\n",
    "#___________________________________________\n",
    "neixy=np.array([[-1,1],[0,1],[1,1],[-1,0],[1,0],[-1,-1],[0,-1],[1,-1]])\n",
    "\n",
    "def neibs(i,j):\n",
    "    neibs=[]\n",
    "    neib=neixy+np.array([i,j])\n",
    "    for i in neib :\n",
    "        if -1<i[0]<N and -1<i[1]<M : neibs.append(list(i))\n",
    "    return list(neibs)\n",
    "    \n",
    "def is_happy(i,j):\n",
    "    neighbours=[int(agents[i,j])]\n",
    "    for k in neibs(i,j):\n",
    "        neighbours.append(int(agents[k[0],k[1]]))\n",
    "    #print(neighbours)\n",
    "    race_slice=(neighbours.count(int(agents[i,j])))/float(len(neighbours)-neighbours.count(0))\n",
    "    \n",
    "    if neighbours.count(0)>0:\n",
    "        dominancy=1.0/(len(set(neighbours))-1)#+1-neighbours.count(0)) # 1 over number of unique races in the neibourhood minus the empty houses\n",
    "    else: \n",
    "        dominancy=1.0/(len(set(neighbours)))\n",
    "    if race_slice>dominancy : return True\n",
    "    else: return False\n",
    "\n",
    "\n",
    "\n",
    "race_up(n_race,agents)\n",
    "empty_up(e,agents)\n",
    "\n",
    "fig, axs = plt.subplots(1, 2, figsize=(12, 5), sharey=True)\n",
    "axs[0].pcolor(agents)\n",
    "axs[0].set_title('Original')\n",
    "\n",
    "\n",
    "def happiness_score():\n",
    "    trues=0\n",
    "    for i in range(N):\n",
    "        for j in range(M):\n",
    "            #print (i,j,agents[i,j], is_happy(i,j))\n",
    "            \n",
    "            if int(agents[i,j])!=0 and is_happy(i,j) : trues+=1#;print(trues)\n",
    "    return 100*float(trues)/(M*N-round(e*N*M/100.0))\n",
    "        \n",
    "        \n",
    "def move(i,j):\n",
    "    new_address=random.choice(empties)\n",
    "    agents[new_address[0],new_address[1]]=agents[i,j]\n",
    "    agents[i,j]=0\n",
    "    empties.append([i,j])\n",
    "    empties.remove(new_address)\n",
    "    \n",
    "\n",
    "\n",
    "while happiness_score() <80.0:\n",
    "    #print agents, '  |happiness score',\n",
    "    #print(happiness_score())\n",
    "    for i in range(N):\n",
    "        for j in range(M):\n",
    "            if agents[i,j]!=0:\n",
    "                if is_happy(i,j): pass\n",
    "                else : move(i,j)\n",
    "\n",
    "axs[1].pcolor(agents)\n",
    "axs[1].set_title('Happier')\n",
    "fig.suptitle('Improved Happiness Score: '+str(happiness_score()))\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}