\n",
"\n",
"# Estimating Probabilities with Simulations\n",
"\n",
"In [another notebook](http://nbviewer.jupyter.org/url/norvig.com/ipython/Probability.ipynb) I showed how to solve problems by computing probabilities. The computations are simple: count the frequency of the \"favorable\" outcomes and divide by the frequency of all possible outcomes—the \"sample space.\" (In [yet another notebook](http://nbviewer.jupyter.org/url/norvig.com/ipython/Probability.ipynb) I tackle some probability paradoxes.)\n",
"\n",
"\n",
"But sometimes it is inconvenient, difficult, or even impossible to explicitly enumerate all possible outcomes. Perhaps the sample space is infinite, or perhaps it is just very large and complicated, with a bunch \n",
"of low-probability outcomes that don't seem very important. In that case, we might feel more confident in writing a program to *simulate* a random outcome. *Random sampling* from such a simulation\n",
"can give an accurate estimate of probability.\n",
"\n",
"# Simulating Monopoly\n",
"\n",
"\n",
"\n",
"Consider [problem 84](https://projecteuler.net/problem=84) from the excellent [Project Euler](https://projecteuler.net), which asks for the probability that a player in the game Monopoly ends a roll on each of the squares on the board. To answer this we need to take into account die rolls, chance and community chest cards, and going to jail (from the \"go to jail\" space, from a card, or from rolling doubles three times in a row). We do not need to take into account anything about acquiring properties or exchanging money or winning or losing the game, because these events don't change a player's location. \n",
"\n",
"A game of Monopoly can go on forever, so the sample space is infinite. Even if we limit the sample space to say, 1000 rolls, there are $21^{1000}$ such sequences of rolls. So it is infeasible to explicitly represent the sample space. There are techniques for representing the problem as\n",
"a Markov decision problem (MDP) and solving it, but the math is complex (a [paper](https://faculty.math.illinois.edu/~bishop/monopoly.pdf) on the subject runs 15 pages).\n",
"\n",
"The simplest approach is to implement a simulation and run it for, say, a million rolls. Below is the code for a simulation. Squares are represented by integers from 0 to 39, and we define a global variable for each square: `GO`, `A1` (for the first property in the first monopoly), `CC1` (the first community chest square), and so on. Wiithin the function `monopoly` the variable `loc` keeps track of where we are, and dice rolls and cards can alter the location. We use `visits[square]` to count how many times we end a roll on the square.\n",
"\n",
"The trickiest part of the simulation is the cards: chance and community chest. We'll implement a deck of cards as a double-ended queue (so we can take cards from the top and put them on the bottom). Each card can be:\n",
"- A square, meaning to advance to that square (e.g., `R1` (square 5) means \"take a ride on the Reading\").\n",
"- A set of cards (e.g., `{R1, R2, R3, R4}` means \"advance to nearest railroad\").\n",
"- The number -3, which means \"go back 3 squares\".\n",
"- `'$'`, meaning the card has no effect on location, but involves money.\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"scrolled": false
},
"outputs": [],
"source": [
"%matplotlib inline \n",
"import matplotlib.pyplot as plt\n",
"import random\n",
"from collections import Counter, deque\n",
"\n",
"# The Monopoly board, as specified by https://projecteuler.net/problem=84\n",
"board = \"\"\"\n",
" GO A1 CC1 A2 T1 R1 B1 CH1 B2 B3 \n",
" JAIL C1 U1 C2 C3 R2 D1 CC2 D2 D3 \n",
" FP E1 CH2 E2 E3 R3 F1 F2 U2 F3 \n",
" G2J G1 G2 CC3 G3 R4 CH3 H1 T2 H2\"\"\".split()\n",
"\n",
"for i, sq in enumerate(board): # Make the square names be global variables\n",
" globals()[sq] = i\n",
" \n",
"def Deck(cards):\n",
" \"\"\"A deck of cards; draw from the top and put it on bottom.\"\"\"\n",
" random.shuffle(cards)\n",
" return deque(cards)\n",
" \n",
"CC_cards = Deck([GO, JAIL] + 14 * ['$'])\n",
"CH_cards = Deck([GO, JAIL, C1, E3, H2, R1, -3, {U1, U2}]\n",
" + 2 * [{R1, R2, R3, R4}] + 6 * ['$'])\n",
"\n",
"def roll() -> int: return random.randint(1, 6)\n",
"\n",
"def monopoly(rolls):\n",
" \"\"\"Simulate a number of dice rolls of a Monopoly game, \n",
" and return the counts of how often each square is visited.\"\"\"\n",
" visits = len(board) * [0] # Counts of how many times each square is visited\n",
" doubles = 0 # Number of consecutive doubles rolled\n",
" loc = GO # Location on board\n",
" for _ in range(rolls):\n",
" d1, d2 = roll(), roll()\n",
" doubles = ((doubles + 1) if d1 == d2 else 0)\n",
" loc = (loc + d1 + d2) % len(board) # Roll, move ahead, maybe pass Go\n",
" if loc == G2J or doubles == 3:\n",
" loc = JAIL\n",
" doubles = 0\n",
" elif loc in (CC1, CC2, CC3):\n",
" loc = do_card(CC_cards, loc)\n",
" elif loc in (CH1, CH2, CH3):\n",
" loc = do_card(CH_cards, loc)\n",
" visits[loc] += 1\n",
" return visits\n",
"\n",
"def do_card(deck, loc):\n",
" \"Take the top card from deck and do what it says; return new location.\"\n",
" card = deck[0] # The top card\n",
" deck.rotate(1) # Move top card to bottom of deck\n",
" return (loc if card is '$' else # Don't move\n",
" loc - 3 if card == -3 else # Go back 3 spaces\n",
" card if isinstance(card, int) # Go to destination named on card\n",
" else min({s for s in card if s > loc} or card)) # Advance to nearest"
]
},
{
"cell_type": "markdown",
"metadata": {
"button": false,
"new_sheet": false,
"run_control": {
"read_only": false
}
},
"source": [
"Let's run the simulation for a million dice rolls and print a bar chart of probabilities for each square:"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"button": false,
"new_sheet": false,
"run_control": {
"read_only": false
}
},
"outputs": [],
"source": [
"N = 10**6\n",
"\n",
"visits = monopoly(N)"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAlsAAAG4CAYAAABy/SLAAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3Xu0ZGlZH/7vw8wQhnvIDA3h1ggEI0xQaFBEjAMkQREligES0Chh9BfAG0pGUCAmGuLSJIBoooImiAyR21JQLoIDzs8w0MMdhpswioByVUEMMPjkj10NZ86cOl21+7ynq09/Pmuddar23m+979773VXf2req7g4AAGNc42Q3AADgIBO2AAAGErYAAAYStgAABhK2AAAGErYAAAYStgAABhK2AAAGErYAAAYStgAABjrzZDdgq3POOacPHz58spsBAHBcl1122ce6+9zjTbdRYevw4cM5evToyW4GAMBxVdUfrzKdw4gAAAMJWwAAAwlbAAADCVsAAAMJWwAAAwlbAAADCVsAAAMJWwAAAwlbAAADCVsAAAMJWwAAAwlbAAADCVsAAAMJWwAAAwlbAAADCVsAAAOdebIbAPvl8IUvWXnaK558v4EtAeB0Ys8WAMBAwhYAwEDCFgDAQMIWAMBAwhYAwEDCFgDAQEPDVlXdsKqeV1XvrKrLq+ruI+sDANg0o++z9ZQkL+3uB1bVNZNce3B9AAAbZVjYqqrrJ/n6JP86Sbr7c0k+N6o+AIBNNPIw4pcl+WiSX62qN1bVr1TVdQbWBwCwcUYeRjwzyZ2TPLq7L62qpyS5MMlPbJ2oqi5IckGSHDp0KBdffPHAJnE6e8x5V648rX4IwF6p7h7zwlU3SfLa7j68eH7PJBd299IfnTty5EgfPXp0SHvAbyMCsJeq6rLuPnK86YYdRuzuP0vygaq6/WLQvZO8Y1R9AACbaPTViI9O8uzFlYjvS/Ldg+sDANgoQ8NWd78pyXF3rwEAHFTuIA8AMJCwBQAwkLAFADCQsAUAMJCwBQAwkLAFADCQsAUAMJCwBQAwkLAFADCQsAUAMJCwBQAwkLAFADCQsAUAMJCwBQAwkLAFADCQsAUAMJCwBQAwkLAFADCQsAUAMJCwBQAwkLAFADCQsAUAMJCwBQAwkLAFADCQsAUAMJCwBQAwkLAFADCQsAUAMJCwBQAwkLAFADCQsAUAMJCwBQAwkLAFADCQsAUAMJCwBQAwkLAFADCQsAUAMJCwBQAwkLAFADCQsAUAMJCwBQAwkLAFADCQsAUAMJCwBQAwkLAFADCQsAUAMJCwBQAwkLAFADCQsAUAMJCwBQAw0JkjX7yqrkjyqSRfSHJldx8ZWR8AwKYZGrYWzu/uj+1DPQAAG8dhRACAgUaHrU7y8qq6rKouGFwXAMDGGX0Y8R7d/aGqunGSV1TVO7v7NVsnWISwC5Lk0KFDufjiiwc3idPVY867cuVp9UMA9kp19/5UVPWkJJ/u7p9dNs2RI0f66NGj+9IeTj+HL3zJytNe8eT7DWwJAAdBVV22ysV/ww4jVtV1qup6xx4n+adJ3jaqPgCATTTyMOKhJC+sqmP1/EZ3v3RgfQAAG2dY2Oru9yW506jXBwA4Fbj1AwDAQMIWAMBAwhYAwEDCFgDAQMIWAMBAwhYAwEDCFgDAQMIWAMBAwhYAwEDCFgDAQMIWAMBAwhYAwEDCFgDAQMIWAMBAwhYAwEDCFgDAQMIWAMBAwhYAwEDCFgDAQMIWAMBAwhYAwEDCFgDAQMIWAMBAwhYAwEDCFgDAQMIWAMBAwhYAwEDCFgDAQMIWAMBAwhYAwEDCFgDAQMIWAMBAwhYAwEDCFgDAQMIWAMBAwhYAwEDCFgDAQMIWAMBAwhYAwEDCFgDAQMIWAMBAwhYAwEDCFgDAQMIWAMBAwhYAwEDCFgDAQMIWAMBAwhYAwEDCFgDAQMPDVlWdUVVvrKoXj64LAGDT7MeerR9Icvk+1AMAsHGGhq2qunmS+yX5lZH1AABsqtF7tv5bkscm+dvB9QAAbKQzR71wVX1zko9092VV9Q27THdBkguS5NChQ7n44otHNYnT3GPOu3LlafVDAPZKdfeYF676T0keluTKJNdKcv0kL+juhy4rc+TIkT569OiQ9sDhC1+y8rRXPPl+A1sCwEFQVZd195HjTTfsMGJ3/1h337y7Dyd5cJJX7Ra0AAAOIvfZAgAYaNg5W1t198VJLt6PugAANok9WwAAAwlbAAADCVsAAAMJWwAAAwlbAAADCVsAAAMJWwAAAwlbAAADCVsAAAMJWwAAAwlbAAADCVsAAAMJWwAAAwlbAAADCVsAAAMJWwAAAwlbAAADCVsAAAMJWwAAAwlbAAADCVsAAAMJWwAAA60UtqrqlasMAwDgqs7cbWRVXSvJtZOcU1V/N0ktRl0/yd8f3DYAgFPermEryfcm+cFMweqyfCls/VWSpw9sFwDAgbBr2OrupyR5SlU9uruftk9tAgA4MI63ZytJ0t1Pq6qvTXJ4a5nu/l+D2gUAcCCsFLaq6llJbpPkTUm+sBjcSYQtAIBdrBS2khxJ8hXd3SMbAwBw0Kx6n623JbnJyIYAABxEq+7ZOifJO6rqdUk+e2xgd3/LkFYBABwQq4atJ41sBADAQbXq1YivrqpbJbldd/9eVV07yRljmwYAcOpb9ed6HpHkeUn+x2LQzZK8aFSjAAAOilVPkH9kkntkunN8uvs9SW48qlEAAAfFqmHrs939uWNPqurMTPfZAgBgF6uGrVdX1eOSnF1V/yTJbyb57XHNAgA4GFYNWxcm+WiSt2b6cerfSfLjoxoFAHBQrHrrh7OTPLO7fzlJquqMxbDPjGoYAMBBsOqerVdmClfHnJ3k9/a+OQAAB8uqYeta3f3pY08Wj689pkkAAAfHqmHrr6vqzseeVNVdkvzNmCYBABwcq56z9QNJfrOqPrR4ftMkDxrTJACAg+O4YauqrpHkmkm+PMntk1SSd3b35we3DYAdHL7wJStNd8WT7ze4JcAqjhu2uvtvq+rnuvvuSd62D20CADgwVj1n6+VV9e1VVUNbAwBwwKx6ztYPJ7lOki9U1d9kOpTY3X39YS0DADgAVgpb3X290Q0BADiIVjqMWJOHVtVPLJ7foqruNrZpAACnvlUPI/5Ckr9Ncq8k/yHJp5M8PcldB7VrGFfxwKnPdrw+ywxOnlXD1ld3952r6o1J0t2frKpr7lagqq6V5DVJ/s6inud19xNPqLUAAKeYVcPW5xc/Pt1JUlXnZtrTtZvPJrlXd3+6qs5KcklV/W53v3Z+cwEATi2r3vrhqUlemOTGVfVTSS5J8tO7FejJsd9TPGvx13MbCgBwKlr1asRnV9VlSe6d6bYPD+juy49XbrE37LIkt03y9O6+9EQaCwBwqqnu5TubFuddfV+msPTWJM/o7ivXrqTqhpn2jD26u9+2bdwFSS5IkkOHDt3loosuWvfl1/LWD/7lStOdd7MbDG0H+2/VdZ9Y/5vudN+O58z/6b7MYITzzz//su4+crzpjhe2npvk80n+IMk3Jrmiu39wToOq6olJ/rq7f3bZNEeOHOmjR4/OefmVuSLn9LXquk+s/013um/Hc+b/dF9mMEJVrRS2jncY8Su6+7zFCz4jyevWaMC5ST7f3X9RVWcnuU+S/7xqeQCAg+B4Yevzxx5095Vr/jTiTZP8z8V5W9dI8r+7+8XrN5G95Nvt6cu6Bzg5jhe27lRVf7V4XEnOXjw/7m8jdvdbknzV3jQTOKgc3gUOul3DVnefsV8NAQA4iFa9zxYAADMIWwAAA636cz0ApzS3SwBOFnu2AAAGsmcLgJPG3kNOB/ZsAQAMZM8WsJS9DgAnzp4tAICBhC0AgIEcRgQ4iRyqhYPPni0AgIHs2YJd2OsAwIkStgC4mlW/aCS+bMDxCFsA7ImDtCf4IM0LJ59ztgAABrJnCwA4JZyqexyFLdhjp+qbwV453eef05e+zzIOIwIADGTPFgCnFHuQONXYswUAMJCwBQAwkLAFADCQsAUAMJAT5AHgFOHigFOTPVsAAAMJWwAAAzmMCKegVQ8lJA4nAJxswhYAHGDO8zr5HEYEABjIni3YAL55AhxcwhYAnCS+aJ0ehK1TmI0UADafc7YAAAYStgAABnIYkVOWw6gAnAqELXa1XzfPFJwAOKgcRgQAGEjYAgAYSNgCABhI2AIAGEjYAgAYyNWIAMAJc1X5csLWhtBJATidnE6few4jAgAMJGwBAAzkMOIKTqddnQDA3rJnCwBgIGELAGCgYWGrqm5RVb9fVZdX1dur6gdG1QUAsKlGnrN1ZZLHdPcbqup6SS6rqld09zsG1gkAsFGG7dnq7g939xsWjz+V5PIkNxtVHwDAJtqXc7aq6nCSr0py6X7UBwCwKaq7x1ZQdd0kr07yU939gh3GX5DkgiQ5dOjQXS666KKh7XnrB/9ypenOu9kNTqjMuvarXeuWWXX67fWsa1Pn5SCtF/Ni21+njPWymfOyqctrTj2bPC+rOv/88y/r7iPHm25o2Kqqs5K8OMnLuvu/HG/6I0eO9NGjR4e1J5l3z6z9uM/WfrVr3TKrTr+9nnVt6rwcpPViXmz765SxXjZzXjZ1ec2pZ5PnZVVVtVLYGnk1YiV5RpLLVwlaAAAH0chztu6R5GFJ7lVVb1r8fdPA+gAANs6wWz909yVJatTrAwCcCtxBHgBgIGELAGCgkXeQP23t1xV8AMDms2cLAGAgYQsAYCBhCwBgIGELAGAgYQsAYCBhCwBgIGELAGAgYQsAYCBhCwBgIGELAGAgYQsAYCBhCwBgIGELAGAgYQsAYCBhCwBgIGELAGAgYQsAYCBhCwBgIGELAGAgYQsAYCBhCwBgIGELAGAgYQsAYCBhCwBgIGELAGAgYQsAYCBhCwBgIGELAGAgYQsAYCBhCwBgIGELAGAgYQsAYCBhCwBgIGELAGAgYQsAYCBhCwBgIGELAGAgYQsAYCBhCwBgIGELAGAgYQsAYCBhCwBgIGELAGAgYQsAYCBhCwBgIGELAGAgYQsAYKBhYauqnllVH6mqt42qAwBg043cs/VrSe478PUBADbesLDV3a9J8olRrw8AcCo482Q3AOCgOHzhS1ae9oon329gS4BNUt097sWrDid5cXffcZdpLkhyQZIcOnToLhdddNGw9iTJWz/4lytNd97NbjC7zKrTzylzIu2aU2bOvMyxqfNykNaLeTEvI9o1p8x+z8scmzovm7q85tSzyfOyqvPPP/+y7j5yvOlOetja6siRI3306NFh7UlW/+a59VvnumXmfLvdj3bNKbNf39Q3dV4O0noxL+ZlRLvmlNnveZljU+dlU5fXnHo2eV5WVVUrhS23fgAAGGjkrR+ek+T/JLl9Vf1pVT18VF0AAJtq2Any3f2QUa8NAHCqcDUiAHAVm3Zu1KnOOVsAAAMJWwAAAzmMyJ6z+xkAvsSeLQCAgYQtAICBhC0AgIGELQCAgYQtAICBhC0AgIGELQCAgYQtAICBhC0AgIGELQCAgYQtAICBhC0AgIGELQCAgYQtAICBhC0AgIGELQCAgYQtAICBhC0AgIGELQCAgYQtAICBhC0AgIGELQCAgYQtAICBhC0AgIGELQCAgYQtAICBhC0AgIGELQCAgYQtAICBhC0AgIGELQCAgYQtAICBhC0AgIGELQCAgYQtAICBhC0AgIGELQCAgYQtAICBhC0AgIGELQCAgYQtAICBhC0AgIGELQCAgYQtAICBhC0AgIGGhq2qum9Vvauq3ltVF46sCwBgEw0LW1V1RpKnJ/nGJF+R5CFV9RWj6gMA2EQj92zdLcl7u/t93f25JBcl+daB9QEAbJyRYetmST6w5fmfLoYBAJw2qrvHvHDVdyT5Z939bxbPH5bkbt396G3TXZDkgsXT2yd515AG7e6cJB/bwDKb2q45ZbTLvJiXU7/MprZrTplNbdecMpvarjll9qtde+VW3X3ucafq7iF/Se6e5GVbnv9Ykh8bVd8JtvXoJpbZ1HYdpHnZ1HaZl81sl3nZzHaZl81s1ybPy37/jTyM+Pokt6uqW1fVNZM8OMlvDawPAGDjnDnqhbv7yqp6VJKXJTkjyTO7++2j6gMA2ETDwlaSdPfvJPmdkXXskV/a0DKb2q45ZbRrfJlNbdecMpvarjllNrVdc8psarvmlNnUds0ps6ntmlNmv9q1r4adIA8AgJ/rAQAYStjitFdVf/9ktwGAg0vYguS1J7sBABxcwtYJqKp/MqPMly8Zfv2qus0Ow//RjDq+e90yi3Jz5mftMru81trLoKpuUlU3WTw+t6q+rarusG7Va05/rO5Tbt73uJ/t2JcX487aYdg5x3m9tcpU1TWq6hqLx9esqjtX1Y2OU8dXVdUDq+of7jbdQbVH28uu634xfu31P6MNK6//qrrRLn/XmVH3nvX9mf147TI7vMZPrzP9CFX15VV176q67rbh9z1ZbRrmZN/o62T8Zfpx7NdkuuPsR5O8Osk3zXidP9mLMkn+RZIPJXlTkrcnueuWcW/Yj3bt5fwcZ/rfXTJ87WWQ5HuTvD/JFUn+vySXJnlmpl8hePiGLq+rzf9+zft+9LMk52f6aa6PJnl5ksMrzMucMg9I8udJPpzpN1cvTfKqxevcf0mZJyR5d5LnJHlfkkfMWe9bXu+Xdhh2/ST/KcmzkvzLbeN+YcnrnLFYn/8hyT22jfvxJWWuneSxSX40ybWS/OtM9zH8mSTXHbm97Nb3Z67LOctsrfW/mO/3Lf5v//vA4u9fnYS+P6cfzynz1G1/T0vyF8ee7zD9TZL8YpKnJ/l7SZ6U5K1J/neSm87cXt667fn3L/reixZ98luPt7wW426R6beW/yDJ45KctWXci+a0bT/+ht76YRNV1SMyvek8NsnRxeAjSZ5cVTfv7l/aNv2yG7FWpk64Ux1P3aXMDXcY/rgkd+nuD1fV3ZI8q6oe190vyJK9LlX1ll3qOLRk3Nz5WatMVd15l+m/csm4tZdBkkcluUOSs5P8cZLbdvefVdXfTfL7SZ6xpU1PS7LTpbfL1smxcnOW17rzP3TeT6SeGX35ZzL9TNfbq+qBSV5RVQ/r7tfuMi9zyjwxyZ0yzf+bMwXHd1XVrZI8P8lv71DmQUm+srs/U1V/L8lLk/zyktefZnL53oJK8k07DP/VJO9ZtOF7qurbMwWIzyb5miWv9T8yhafXJXlqVb26u394Me7bkvzHHcr8WqaAcHaSlyS5PMnPJrl/pg/Jh+1QZq0+M2PdJ/PW5Zxlttb67+5bL3mdaYaqzs30pfvZW4btR9+f04/nlPm2JBdnCoHH2vLgJJctadevZepX18nUN56d5H6Zwt1/X/y/mqr6tiWvV5kC3FaPyPR+9OmqOpzkeVV1uLufkuXLK5m+IDw/0+kfD0/y6qq6f3d/PMmtdil3Up12YSvJDyX5uu7+xJZhr6qqb0xySa5+v457Jnlokk9vG15J7rakju9O8pgkn91h3EN2GHZGd384Sbr7dVV1fpIXV9XNs3NASKZA9c+SfHKHdv3hkjLJvPlZt8zrM71x7bTBLHuTnrMMPt/dn0nymar6o+7+s0X5T1bV9jJHr158pXFzlte683/mjHm/co15P2bOMl63L1+zFzcv7u7nVdXlSV5QVRfuUsecMjk2z1X1J939rsWwPz52eGUH/3exzNLdH99luq0+mimYbF2XvXh+4x2mv013f/vi8Yuq6vGZ3l++ZZc67tbd/2gxLz+f5Beq6gWZlu+yD51/0N3/oqoq0x6O+3R3V9UfZPoA3sk620uy/rpP5q3LOctszvpPVZ3V3Z/fNuyc7v5oVf27bZPvR9+fNR8zyvzDTHtO75vkR7v7g1X1xO7+n0umP9TdT1vU8W+7+z8vhj+tqh6+rF1JnpspmO00v9fa9vyM7v70ou1XVNU3ZApct8ruYevc7v7vi8ePrqqHJnnNor9s7r2sTvautf3+S3L5OuOS/G6S85dM/5olw1+V5GuXjHv/DsP+MNMbztZh10/ye0k+u+R1npEpNO407jd2mcc587NWmSRvS3K7JdN/YMnwnZbB9ZK8cpdlcDSLXchJbr5l+LWSvHmP+suc5bXW/O8y77ut/7XnfeYyXrcvH01yk23Dbp7p0OWndpmXdcu8Mck1Fo/vtmX4GUnetqTMX2Q61PZbmb79b33+W0vKvCfJLddYl5cfa9eWYd+V6bDtHy95nXfuMOwJSf7/JO9ZUuZNWx4/c9u4Zet/rT6z7ro/gXU5Z5mttf5z9UN8t94ybtkhvv3o+3P68dpltkxzl0x7qn4kyRW7TPfmLY//47Zxb9ml3GVJ7rhk3Ae2PX9Vpj3NW4edmeR/JfnCLnW8Pcm1tg27T5L3JvnwbvN/Mv9OegP2fYan49t32mH4nZK8bo3XuUeSpy8Zd6MkZ6/xWnfKDh/OSb4+yfuWlFn7HJt9XMYPTHL7JeMesOYyOCtLzqVIcstsOV6/ZfjNMn3T3zrst7Plg3X73y7zcrtsO49mMfye2RZc5s5/krsmueeS9b/snJ2v36XMsnO25izjGyW59hrr/j5Ltq8bJHn8mmVuuEuZu25/w10MP5zkoUvK/OPd/paUeeRObVuMe/QOw35me99bDL9vlgenX09y3x2G/5tMe6N2KvMr2eHcrCS3SXLJkjK3zLQX9bjby5x1fwLrf84yW2v9Z9rbfIfF4wdmCtFfs3j+xiV1rPs+Pmfe5/TjZWVutazMtulq0a9/fZdpfnJJ/7ptkuftUu6eWf7l5Mi25+/NtAdtp2mv9p67ZdwP7bS9JvmqJK9Yp7/u599pdwf5qvq6TLs5fzVTCu9Mnfe7MnXUS3Yp+5VJ/mWmE43fn+T53f3zO0x3bqZdne/YNvwOST7S3R9do44X9GJ37rbp3tDdy84NWltVnZHkwd397ONOnC9eXfPxHtiB5tSxrExV/ePdynX3q5e83ouTPK6737Jt+JEkT+zu++9Q5jG56u7sznQxxiXd/f49qmPtMsusuu4X/TrH6b+z+/6qdcwtU1W37O4/WfV1T1VVVWtuMzuu/znLa2aZtbaXOfVU1Zu7+05bnt8hyQuSXJjkJ3Z6L62qb820B/Dpi+eXJjl3Mfqx3f28Nab/d939mzvUsfb2MqeeZdZ9399Le/0ZtulOu3O2uvuSxcnBj8x0BU8y7Zb86u7+8+3TV9U/yHQi4UOSfDzTMenq7vN3qeZpmU5S3e7mSR6fKUydaB03rqofXjayu//LTsOr6vqZ5v1mmfbqvCLTibM/kml399U2uqr6miRPTvKJTMf9n5XknCTXqKrv7O6Xbpt+abuWtW3dOtYts0uYukWmZb/j+ExXFF3tYoTuPro4qXMn191h2OEkj6+qJ3X3RXtQx9plZq77ynRC7qMyfSO+RlVdmeRp3f2TO1SzVt/fUs8Tkzx6xTrmtCuZrnq686L88/tL5wktVVWP7e6fWTz+jq0fZFX10939uG3Tz+n7a9VxvDJJfirTxRDby6y7/tdeXjPLrLu9zKnn81V1k/7SeWpvr6p7J3lxpr2BO3lspveGY/5Opi/m18n0Zf15a06/Uwias72sXc+Sdf/ITFezXm3dz+nHi3LLLkQ6Vu77tzyd+xm27MKFnerYGKdd2NryreAJi+evy7TRPnTxBrZ9A3pnpktM79/d712U+aHjVHPeTh/u3f2yqvq5HaafU8cZmd6k1r1H1LMynVT/fzIdpvjRJNfMdNntm5aU+flMb943yHSc/Ru7+7U13WvmOZmu6trqelsef2+mq62OZ9065pY5tvfrOzKF25sleeEu7dp+UudWZ+80sLv//ZJ6b5TpPKztHx5r1zGzzJx1/4OZDpnf9dhehqr6siS/WFU/1N3/ddv06/b9Y33969aoY067kqtuK1+2ZH63e3Cmw1xJ8mO56gfZfXP1UDOn769bx9wy667/Octr7TIztpc59VyY6aKiP9tS758u9ng/akmZa3b3B7Y8v6SnK94+Xjvfm2vd6ZMZ28vMetZd93P6cXLVi43+faYvRMvM/QzbegXl8erYHDsdWzzIf5lOOr3FludvynRs/pZJXrnD9P88056mD2S6TPzeWXJy6JYy795l3Lv2qI5Z52xly71OMnX2Tya53nHKbD0Z9/Jt43Y832HV8SdSxzplMr15fGemAPa+JD+X5E9XaNdzssP9mDJdcvzcGcv/avMyp46ZZeas+zcmOWeH4ecumZe1+v6cOk6gzBt2erzq+tqhT+1V31+7jpll1lr/M5fX2mXmLMN168l0q4JHbnl+6eJ94H1JvmNJmffu8np/dKLTL4bP2V7m1LP2tn+8dXCi5Ub2j038O+32bGXnbwWfSPKJnb4VdPcLk7xwMe4BmU7OO1RVv5jkhd398h3qeE9VfVN3/87WgTXdXuJ9e1THrLueJ/niZc/d/YWqen93f+o4Zf52y+O/2TbueOeGrHruyJw61inzkUz3MfrxTOu8q+qfr9CuH8y0bv5VvvSN6kimb4WrlP+iqrpXrn6rjrl1zCkzZ92f1d0f2z6wp0vlr3an7KzZ92fWMbfMnarqrzJtO2cvHmfxvLv7+juU2X4u0bJxO1m178+pY06Zddf/nOU1p8yOdtle5tQz5xDfpVX1iO6+yr3Yqup7M72XnOj0ybztZU49c7b9LxZZcbp1y839DFunjo1xOp4g/97uvu2ScX/U3cuO32+d7kaZDkM9qLvvtcP422W6Idwf5qofhHdP8s3d/e49qONGfdV7ha2kqr6Q5K+PPc10yOkz2eXNcEuZrdMfK3+t7l724bbySZBz6linzOJQ1YMzvbn+RqY9ia/o7pUOddR0T6o7Lp6+vbtftcu0b83V3wRulOnu7d/Z3e880TpmtmvOul+6/nYaN6fvr1vH3DJzHIS+v0OZZMX1vx/mbi9r1vH67r7rluc/392PWjx+bXdf7eapVXXjTOeGfTbJGxaD75IpqD2gt53ju+70izJztpc59cxe93O3p+OVm/sZthdtOxlOx7D17CQXL/lW8A3dvexmfevUcdtMd8u9XbZ8EGa61PWD3f1HJ1rHJtv25nnbTPP9Rb24gePJUNPvAj62cWTtAAAEZUlEQVR48Xe7TMf7X7hKAF6jjlttG9SZrpD8652m32Tb3qSvMio7fKjP6fvr1jG3zH7Y5L6/qfZjezmRL9mLPWzHfj/yuF+C1pn+RD4r1m3XOub246r61JZy185VvwTsSaDfjzpGOB3D1trfCmbUsWeX5Z+KFt/WDmU6B22rWyX5UC8uAjjZquq8TCfJP2iVPZocn75/avT9081+fMme2a6N3F7047132oWtYwZ/K3hbd99xybi3dvd5e1XXJtrEN5Bt34auMipT8H5vphsPvnJfG3bA6Pub1/fZny/ZM9u1kduLfrz3TscT5JMki3C1ZwFrmzmX5R8kh7dvpMlx7xs1VHdfb9m4mm7sd8dM95rZ8Y2Plen7G9b3Sbr7I0m+dtuX7Jfs5ZfsmTZ1e9GP99gqP8TK+l5fVY/YPrCmH/Bc9ivrB8mmvoHsqLu/0N1vznSDQU6Mvr/cxvX90013v6q7n7b4O9lBK9nc7UU/3mOn7WHEkarqUKYbZX4uO1yW34u7GB9UVfWcJK/a4fyIhyf5p939oJPTMkbT9/V9Vrep24t+vPeErYHmXMp/EGzqGwj7R9/X91ndpm0v+vHeE7YYZtPeQGC/6PscBPrx3hG2AAAGcoI8AMBAwhYAwEDCFrCRqurxVfX2qnpLVb2pqr76ZLcJYI7T9qamwOaqqrsn+eYkd+7uz1bVOZmuhBpV3xnd/YVRrw+c3uzZAjbRTZN8rLs/myTd/bHu/lBV3beq3llVl1TVUxc/K5KqelJV/cixwlX1tmN3uq6qF1XVZYu9ZBdsmebTVfWTVXVpkrtX1V2q6tWLaV9WVTddTPf9VfWOxR62i/ZvEQAHhT1bwCZ6eZInVNW7k/xekucmuTTJLye5V6bfsnzuiq/1Pd39iao6O9Mdu5/f3R9Pcp0kb+vuJ1TVWUleneRbu/ujVfWgJD+V5HuSXJjk1os9bDfcy5kETg/2bAEbp7s/nemHgi9I8tFMwer7kry/u9/T0z1rfn3Fl/v+qnpzktcmuUWS2y2GfyHJ8xePb5/pfkKvqKo3JfnxJDdfjHtLkmdX1UOTXHlCMwacluzZAjbS4hyqi5NcXFVvTfJdSZbdGPDKXPXL47WSpKq+Icl9kty9uz9TVRfnS7/79n+3nKdVmW7aePcdXvt+Sb4+ybck+YmqukN3C13AyuzZAjZOVd2+qm63ZdBXJvnzJLeuqtsshj1ky/grktx5UfbOSW69GH6DJJ9cBK0vT/I1S6p8V5JzFyfmp6rOqqo7VNU1ktyiu38/yWOT3DDJdU94BoHTij1bwCa6bpKnLc6RujLTOVoXJHlekpdU1ceSXJIv/ZTI85N85+IQ4OuTvHsx/KVJvq+q3pIpUL12p8q6+3NV9cAkT62qG2R6b/xvi9f59cWwSvJfu/sv9nxugQPNz/UAp6TFIcIf6e5vPtltAdiNw4gAAAPZswUAMJA9WwAAAwlbAAADCVsAAAMJWwAAAwlbAAADCVsAAAP9P2QDRzt6LwinAAAAAElFTkSuQmCC\n",
"text/plain": [
"