{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# *Introduction to Monte Carlo Algorithms*\n", "`Doruk Efe Gökmen -- 20/07/2018 -- Ankara`\n", "## Calculation of $\\pi$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Direct Sampling:\n", "Calculation of $\\pi$ with Monte Carlo algorithm using direct sampling and the calculation of rms deviation." ] }, { "cell_type": "code", "execution_count": 23, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Populating the interactive namespace from numpy and matplotlib\n", "3.138\n", "3.153\n", "3.139\n", "3.152\n", "3.168\n", "3.115\n", "3.16\n", "3.176\n", "3.185\n", "3.161\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZwAAAEaCAYAAAAlqOH8AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xm8VfP+x/HX55xOdaTBkKGQKRG55ZJSKFOmSFzzEF2k+F0XIdN1TUUuriuSzPOUiqgoiQiRK1PGyIlrLMlJ0+f3x3cddsfe56zTae+1zznv5+OxH+fsNe3P3t+112d/1/qu79fcHRERkWwrSDoAERGpG5RwREQkJ5RwREQkJ5RwREQkJ5RwREQkJ5RwREQkJ2pkwjGz4WZ2cdJxJMXMLjWz+6L/NzGzn82sMIE4rjCz78zs61VYt06XYV2Wuv9WczvV2oei783m1Y2jpjOzY8xsYsxlq1V29VZ1xWwxsznA+sAyYDnwHnAPMMLdVwC4e78svXYf4K/u3jUb288Gd/8CWDPXr2tmGwNnA63c/Zuqrp+tMowj2sf+6u7PJRWDVF9V9iEzmwLc5+4jU9bP+fcml8xsU+AzoMjdl2Vazt3vB+7PRUz5WsPp6e6NgVbAEOA84PY4K5pZ3iXRWqoV8P2qJJvqUhknR5997ZLz8nT3vHoAc4C9yk3rCKwAtoue3wVcEf3fDfiSkJS+Bu6Nph8IvAXMB14Gtk/Z3sbAKOBb4HvgJmAbYDGhVvUzMD9DfH2AT4GFhF8Px0TTtwAmR9v7jvCLoVm59zUQeBtYREig6wPPRNt6DlgrWnZTwIFTgHnAV8DZKdu6lPBrLXXZetHzKcDlwLRouxOBdVPWPR74PIrz4nSfd8qyTQm1y2+jdS4i/EjZCyiNyuRn4K4065aVywXR5zGn7LMqX4Yx9ok+0fu5HvgBuCKK46Iorm+iOJumrHMQ8G5U/lOAbaLp90Zxl0axn1tB7OdG2/4K6AXsD3wYxXBByvIFwPnAJ9Hn+giwdsr8Rwn75gJgKrBtuc9hGDAuKq9XgS0yfA4Ngfui15gPvA6sH83bDHgh2sazhH36vtT3k+l7Rvh+vRJt86to3fopyzowAPgI+CyatnX0Oj8As4HDKyi/jLFF8zsRvqPzgf8C3aLpRwIzym3r78DYNMeBtYCnCPvqj9H/G0XzriR8rxdHZX5TyvvasqJ9PWX/ewm4Ntr2Z8B+Gd7r+cBj5ab9G7ixouNHjO/AFCr4XmdY54voPf4cPTqT/rvUB3ipXLxzgZ+AN4BdMxx7Mu6PGWOK82Zz+SDDATD68E5Ls6N1I5x+uxpoABQDOxAOFDsDhcAJ0XYbRM//G33gjaIPrWvqjlVBbI2iQmgTPd+Q6OABbAnsHb1Gc8KB5YZy72s6Icm0jOJ7E+gQrTMZ+Ee07KbRjvJg9JrtCF+EsgNEaqGXLZuacD4Btoo+iynAkGhe22jH6wrUJ3yBlqb7vKPl7wHGAI2j1/kQ6JvpIFZu3bJyuS56f7sTEm2b8mUYPZ9fVg5pttUn2tYZhNPAxcBJwMfA5oRTiqP4/cfGVtFr7Q0UERLHx0QHUSpIsuVivyRa/+To838g+iy2JRy8No+WPzMq242i93or8GDK9k6K1msA3AC8lTLvLsIXv2P03u4HHsoQ16nAk8AahP34z0CTaN4rKZ/1boSDUtyE82fCQb9eVM7vA2emLOuERLF29Nk3IhyQTozW2YHwo2LbDHFXFFtLwgFrf0Li3jt63jx6nwuB1inbeh04Ms1xYB3g0GidxoQkPzplvSmE06iUe19lCaeifb0P4XtycvS5n0b4IWhp3msr4JeUcikkJPFOVHD8iHFcnEKG73UF62xKyrGhgu9SH1ZOOMdGn2c9wmnzr4GGaY49GffHjDHFebO5fJA54UwHLkyzo3UDlpR9ING0W4DLy60/m3DQ60w4eNRL8xorffBp5jciHBgPBYoreR+9gJnl3lfqL/zHgVtSnp9R9gVJ2VG2Tpl/DXB7mkJfaaeKdsSLUtbrD4yP/r+ElQ+Ea0SfXbrPuxD4FWibMu1UYErK5x4n4TRKmfYIcHH5MoyxT/QBvig3bRLQP+V5G8JBoR6h5vZIyrwCoITffzmn3cfKxV4KFEbPG0ef8c4py7wB9Ir+fx/YM2XehmWxpNl2s2hbTVM+h5Ep8/cHPsgQ10mUq61H0zdJ81k/QMyEk+Z1zgSeSHnuwB4pz48AXiy3zq1EP5iqGNt5RD8UUuZPAE6I/r8PuCT6vzUhAa1R2T4EtAd+THk+hQwJh8r39T7Ax+W+Nw5skOG1XwKOj/7fG/gk+j/28SPNNqeQ4XtdwTqbkj7hlP8u9aHi496PwJ+i/y9NKbu0+2NFj3y9hpNOS8IvwXS+dffFKc9bAWeb2fyyB+E0Wovo7+dewUW0TNx9EeHL1g/4yszGmdnWAGa2npk9ZGYlZvYT4YuybrlN/C/l/9I0z8tfxJyb8v/nUfxxpLYa+yVluy1St+nuvxB+TaazLqEW9Hm5GFrGjAHCF35RufXjvofy5pZ73iJNbPUINciV5nlobDKXqsX+vbsvj/4vjf5mKq9WwBMp+9r7hFM465tZoZkNMbNPov1iTrRO6r6RqbzKu5dwMH7IzOaZ2TVmVkR4v+k+61jMbCsze8rMvo5ivIo/7rupn38rYOdy369jgA3SbL6y2FoBfym3ra6EpA0hOR0V/X804UfZL2newxpmdquZfR69h6lAs5itN+Ps67+VUcrrZyqn8jE/EK2X8fgRU9z9pDLlv0srMbOzzex9M1sQlUdT/rg/QOb9MaMakXDMbCdC4b+UYREv93wucKW7N0t5rOHuD0bzNslwsaz8dv64gPsEd9+b8IX4ALgtmjU4Wn97d29CqJZaZdurxMYp/29CqMZXx1eE0z4AmFkxoeqczneEX+mtysVQUoXXW8vMGpVbf1XfQ/mymZcmtmWEpLDSPDMzwmdZFnul5VxFcwnn9FP3t4buXkI44BxMuO7VlPCrE1Zh33D3pe7+T3dvC+xCuE55PKFc033WZRYRfpWHFw4H4eYp828h7Muto333gjTxpX5mc4EXyr3fNd39tDRhVxbbXEINJ3Vbjdx9SDR/IrCumbUnHMQfSPMaEE79tCHUQpsQTt2R8j4qKvPVsa+nehToZmYbAYekxlzB8SMbMr3njJ+Fme1KqHUeTrim3Ixw7fEP+2sF+2NGeZ1wzKyJmR0IPESoxs2KueptQD8z29mCRmZ2gJk1Bl4jfAmGRNMbmlmXaL3/ARuZWf0M8axvZgdFX55fCddDyn4FN46ezzezloQGAtV1cfTLbVvC+fKHq7m9x4CeZrZL9B7/SYYDX/Tr/hHgSjNrbGatgLMINbeq+KeZ1Y925AMJX8bV4UHg72a2mZmtSfhV/nBUc30EOMDM9ox+cZ1NKK+Xo3X/R7j2s7oMJ3xOrQDMrLmZHRzNaxy99veEg/5Vq/oiZtbdzNpFCeMnwkFyubt/Dszg98+6K9AzZdUPgYbRd6CIcEG8Qcr8xtH2fo5+cadLHKmeArYys+PMrCh67GRm25RfMEZs9xH2yR5RbbChmZUdrInK8zFgKOEa0rMZYmpMqHXON7O1gX+Um5+xzFfjvl62vW8Jp8DuJDSyeB8qPX5kw7eEBjJV2dcbE364fQvUM7NLgCbpFsy0P1a08XxNOE+a2ULCr58LCRccT4y7srvPIFzgu4lw/vFjwnnKsp2rJ+Hc7ReE1khHRKtOJrRs+trMvkuz6QLCwWse4fTe7oRzqRAO3jsQfg2MI1zErq4XotgnAde6e6ybszJx93cJ14oeIiTdhYTGC79mWOUMwq/jTwm1yweAO6rwkl8TPv95hIvh/dz9g3QLWrgJb9cqbPsOQpV+KqG1z+IoXtx9NqGG+R/Cr9eehKb2S6J1BwMXRadwzqnCa2byb2AsMDHab6cTGqxAuBj9OeHX8nvRvFW1AeHg+xPhtN0L/H5QPDp6zR8IB9t7ylZy9wWE/XRkFMciwn5f5pxo/YWEH2sV/rBx94XAPoRWZPMI5VzWaCedimKbS6gBXkA4yM0l/FhLPTY9QKghPlrBqfAbCBfAvyN8xuPLzf83cJiZ/WhmN6ZZv7r7enllMafWyDIeP8xsVzP7uRqv9wfRqb8rgWnRvt4pxmoTCC1nPyTst4vJfAquov0xLYsu/kgesZg3bK2G11mTcBGztbt/tpq33Y1QK92osmVl9TOzSwktsI5NOhaRMvlaw5EsMbOe0Wm6RoRm0bP4/UK2iEjWKOHUPQcTqvTzCM1Mj3RVc0VqJAv9oP2c5vFu0rGlo1NqIiKSE6rhiIhITijhiIhITtTKnl/XXXdd33TTTZMOQ0SkxnjjjTe+c/fmlS+56mpVwjGznkDPLbfckhkzZiQdjohIjWFmsbtDWlW16pSauz/p7qc0bdo06VBERKScWpVwREQkfynhiIhITijhiIhITijhiIhITijhiIhITijhpPpkMsz/IukoRERqJSWcMsuXwpgz4ObO8NptsGJF0hGJiNQqSjhlCovgxKdh447w9Dlw1/7w3UdJRyUiUmso4aRaqxUcOwp63QLfvA+3dIEXrwu1HxERqZZalXCiwcVGLFiwoDobgfZHw4DXYKseMOmfcNse8NV/V1+gIiJ1UK1KOKu1a5vG68MR98Lh98DCr2FEd5h0GSxdXP1ti4jUQbUq4WRF24Ph9NfgT0fBi/+C4V3hi+lJRyUiUuMo4cRRvBb0Ghau7yz7Fe7YF54eCL8uTDoyEZEaQwmnKrbcE/q/AjufGppO39wZPn4u6ahERGoEJZyqarAm7Hc1nDQeiorhvkPhidPglx+SjkxEJK8p4ayqTTrBqS/CrufA2w/DsJ3hvTFJRyUikreUcKqjqCHseTGcMgWabAiPHA8PHxtatYmIyEqUcFaHDbeHv06GvS6FDyfCsI4w8z5wTzoyEZG8oYSzuhTWg65/h9NehvW2hTED4N5D4MesDxMuIlIjKOGsbutuCX3Gwf7Xwpevh5Zs04fDiuVJRyYikiglnGwoKICOJzNh99G8snwrGH8eb1++C89NnZp0ZCIiialVCWe19KW2moyeWcKZ47/nqF/O4cwl/dl4xZfsNukQ3nvoYnUGKiJ1Uq1KOKu1L7VqGjphNqVLlwPG6BVd2evXoUxcsSNtP7gx9Ms2b2bSIYqI5FStSjj5ZN780pWef09TTl/6f5yy5O+w6Fu4bU949h+wtDTDFkREahclnCxp0aw47fR3m+wGA14NQyBMuyGMuTNnWo6jExHJPSWcLBnYow3FRYUrTSsuKmRgjzZQ3AwOvgmOHwMrloXRRcedDYt/SihaEZHsU8LJkl4dWjK4dztaNivGgJbNihncux29OrT8faHNu4XOQDv1h9dvD02oP3o2oYhFRLLLvBbeDb/jjjv6jBkzkg6jaua+DmNPh28/gO2PhH0HwxprJx2ViNQRZvaGu++YzddQDSdfbLwTnDoVdj8P3nkMbtoJ3hml7nFEpNZQwskn9RpA9wvglBeg2cbw2Inw0DHw01dJRyYiUm1KOPlog+2g73Ow9+XwyaQw9MGb96i2IyI1mhJOviqsB13+L3QGukE7GHsG3HMQ/PBZ0pGJiKwSJZx8t84WcMKTcOD1UDITbtkFXrlZnYGKSI2jhFMTFBTAjieFG0Y33RUmDILb94Fv3k86MhGR2JRwapKmLeHoh+HQ2+HHz2D4rjDlali2ZKXFRs8socuQyWx2/ji6DJnM6JklCQUsIvI7JZyaxgzaHQYDXoO2B8OUq2BENyh5AwjJZtCoWZTML8WBkvmlDBo1S0lHRBKnhFNTNVoXDrsdjnoISn+EkXvBxIu4cfzbUS/VvytdupyhE2YnFKiISKCEU9O12Q8GTIcdjoeX/8PtpX+jU8F7f1isfO/VIiK5VmnCMbPmZnaBmY0wszvKHrkIrqryaQC2nGrYFHr+G054kqJC46H6V3BlvdtpzC+/LZKp92oRkVyJU8MZAzQFngPGpTzyTj4NwJaIzXZj5v7juGPFgRxZOJmJDc5lj4I3f++lWkQkQfViLLOGu5+X9Uhktei505aMrncNpz6zK+csvok76l/L3JYHsPFWNyYdmojUcZX2Fm1mVwAvu/vTuQmp+mpkb9HZsGwJvHQ9TB0KDZvAftfAdoeGlm4iIinypbfovwFPmdliM1sYPTRSWE1Qrz50Oy/0Qr3WpvB4X3jwSFigJtIiknuVJhx3b+zuBe7eMPq/sbs3yUVwspqs3xb6Pgs9roJPX4CbO8GMO2HFiqQjE5E6JFazaDM7yMyujR4HZjsoyYKCQug8APq/DC3aw1Nnhs5Av/8k6chEpI6I0yx6COG02nvR42/RNKmJ1t4cjh8LPW+Er/4bOgOddiMsX5Z0ZCJSy8VpNPA20N7dV0TPC4GZ7r59DuJbJWo0ENNP82Dc2TD7aWixAxx8E6y/bdJRiUgC8qXRAECzlP/r6E0utVCTFnDkA3DYHTD/C7h1N3j+Klj2a9KRiUgtFCfhDAZmmtldZnY38AZwVXbDkpwxC02lB7wW/r5wNdy6O3ypGqKIrF5xWqk9CHQCRkWPzu7+ULYDkxxrtA70HgFHPwK//hQ6Ax1/ASxZlHRkIlJLZEw4ZrZ19HcHYEPgS2Au0CKaJrXRVj2g//Qw4Nv0YaFRwacvJB2ViNQCFXVtcxZwCvCvNPMc2CMrEUnyGjaBA68Lp9jGnhGaT+9wPOx9ORQ3q3DV0TNLGDphNvPml9KiWTEDe7ShV4eWOQpcRPJZnFZqDd19cWXT8olaqa1GS0thymB4+T/QaL2QiLY+IO2iZYO/pY7HU1xUyODe7ZR0RPJcvrRSeznmNKmNioph78vgr5PCoG8PHQ2Pngg/f/uHRYdOmK3B30Qko4yn1MxsA6AlUGxmHYCyHh+bAGvkIDbJJy13gFOmwEs3wNRr4NPnYd+rYfvDf+sMNNMgbxr8TUSg4ms4PYA+wEbAdSnTFwIXZDEmyVeFRbD7QNimJ4w9HZ44Bd55DA68HppuRItmxZSkSS4a/E1EoIJTau5+t7t3B/q4e/eUx0HuPiqHMUq+WW9rOGlCqOHMeQmG7Qyvj2TgPq0pLipcaVEN/iYiZSptNABgZgcA2wINy6a5+2VZjKta1Gggh36cA0/+DT6dApvswrOtL+LSab+qlZpIDZOLRgOVjvhpZsMJ12y6AyOBw4DXshmU1CBrbQrHjYa37ocJF7D3vN7s3W0QdD4dCuMMKCsidUWcVmq7uPvxwI/u/k+gM7BxdsOSGsUMOhwbusfZci947h8wcg/4elbSkYlIHomTcMquAv9iZi2ApcBm2QtJaqzGG8AR98Ff7g49UY/oBpOvUGegIgLESzhPmVkzYCjwJjAHyMu+1Mysp5mNWLBgQdKh1F1msG2vUNtpdzhMHQrDd4UvXk06MhFJWKxGA78tbNYAaOjueX1EV6OBPPLxc/DkmbDgS9j5VNjjYmiwZtJRiUg5iTYaMLM93H2ymfVOMw81jZZYttwL+r8Cky6DV4eHwd56/hu2UFd8InVNRafUdo/+9kzzODDLcUlt0qAx7D8UThwPhQ3g3kNg9AAo/THpyEQkh+J03lno7ssrXCjP6JRaHlu6OAzyNu3foW+2A/4Vei4QkUTlS+edn5nZCDPb08ys8sVFKlDUEPb6B5w8GdZcDx4+Fh45Hhb+L+nIRCTL4iScNsBzwABC8rnJzLpmNyyp9Vq0h5Ofhz0vgdnjYVhHeOsBqEIjFhGpWeIMMV3q7o+4e2+gA6G3aA0BKdVXWAS7ng39XoLmW8Po0+C+Q2H+F1Xe1OiZJXQZMpnNzh9HlyGTGT2zJAsBi0h1xKnhYGa7m9nNhPtwGgKHZzUqqVuabwUnPgP7DYUvpsOwTvDqCFixItbqZQO/lcwvxYGS+aUMGjVLSUckz1SacMzsM+BM4EVgO3c/3N0fz3pkUrcUFMDOp8CA6bBJJ3hmINy5H3z3UaWrauA3kZohTg3nT+5+iLs/6O6Lsh6R1G3NNoFjH4det8C3H8AtXeDFf8HypRlX0cBvIjVDnISzgZlNMrN3AMxsezO7KMtxSV1mBu2PDt3jtNk33DR6W3f46r9pF880wJsGfhPJL3ESzm3AIEKnnbj728CR2QxKBIDG68Ph98Dh98LP38CI7vDcpeFenhQDe7TRwG8iNUCchLOGu5cf/2ZZNoIRSavtQTDgVfjTUfDS9TC8C3z+ym+ze3VoyeDe7WjZrBgDWjYrZnDvdhr4TSTPxBkh6zsz2wJwADM7DPgqq1GJlFe8FvQaBu0ODSOM3rkv7HRyuIm0QWN6dWipBCOS5+LUcAYAtwJbm1kJocXaaVmNSiSTLfaA016BnfvB6yPh5s6hR2oRyXtxbvz81N33ApoDW7t7V3efk/XIRDJpsCbsdzWcNAGKisPNok/0g19+SDoyEalARcMTnJVhOgDufl2WYhKJZ5Od4dQXwyBv024INZ39r4W2B4eWbiKSVyqq4TSOHjsSTqG1jB79gLbZD00khqKGsOfFoV+2Ji3g0RNCh6ALv046MhEpJ87wBBOBQ919YfS8MfCou++bg/hWiYYnqKOWL4NXboIpg6FeA+hxFbQ/RrUdkRjyZXiCTYAlKc+XAJtmJRqR6iisB13PhH7TYL1tYcwAuLcX/Dgn6chEhHgJ517gNTO71Mz+AbwK3J3dsESqYd0toc+4MLjblzNCS7bpw2FFjRpHUKTWidNK7UrgROBHYD5worsPznZgItVSUAA7/RX6T4dWXWD8eXDHvvDNB0lHJlJnVXoNpybSNRxZiTu8/QiMPx+W/Ay7nRtOvRUWJR2ZSN7Il2s4IjWbGfzpiNAZ6NYHwvNXwIhuMG9m0pGJ1ClKOFJ3rNkc/nInHPkALPoObtsDnr0ElmoYA5FciDviZysz2yv6vzhqGi1SM219QOgMtMOxMO3fYcydOdOSjkqk1osz4ufJwGOE/tQANgJGZzMokawrbgYH/QeOHwMrlsFd+8NTZ8Hin5KOTKTWitt5ZxfgJwB3/whYL5tBieTM5t2g/yvQaQDMuANu7gQfTqzSJkbPLKHLkMlsdv44ugyZzOiZJVkJVaSmi5NwfnX33278NLN6REMViNQK9RvBvldB32ehQWN44C/w+Mmw6PtKVx09s4RBo2ZRMr8UB0rmlzJo1CwlHZE04iScF8zsAqDYzPYGHgWezG5YIgnYeCc4dSrsfj68OwqGdYR3Hg/NqjMYOmE2pUtXvqG0dOlyhk6Yne1oRWqcOAnnfOBbYBZwKvA0cFE2gxJJTL0G0H1QSDzNNobHToKHjoaf0o85OG9++hZumaaL1GUVJhwzKwTucffb3P0v7n5Y9L9OqUnttv620Pc52OcK+GQyDNsZ3rj7D7WdFs2K066eabpIXVZhwnH35UBzM6ufo3hE8kdhPdjlDDjtZdigHTz5f3DPQfDDp78tMrBHG4qLCldarbiokIE92uQ6WpG8l3EAthRzgGlmNhZYVDZRA7BJnbHOFnDCk/Dm3TDxYrh5F9jjIuh0Gr06tATCtZx580tp0ayYgT3a/DZdRH4XJ+HMix4FhAHZROqeggLY8URovQ+MOwsmXhgaFhx0E706tFWCEYkh7zvvNLPNgQuBpu5+WJx11HmnZJV7aL32zLnhRtHdzoGuZ0E9nXmWmisvOu80s+ZmNtTMnjazyWWPOBs3szvM7Bsze6fc9H3NbLaZfWxm51e0DXf/1N37xnk9kZwwg3aHhc5At+0VRhgdsTuUvJF0ZCJ5LU6z6PuBD4DNgH8Srum8HnP7dwErDUUdtXwbBuwHtAWOMrO2ZtbOzJ4q91CPBpK/Gq0Lh46Eox6C0vkwci+YcCEs+SXpyETyUpyEs4673w4sdfcX3P0koFOcjbv7VOCHcpM7Ah9HNZclwEPAwe4+y90PLPf4pipvRiQRbfaDAdNhhxPglZvgll3gs6lJRyWSd+IknKXR36/M7AAz60DowHNVtQTmpjz/MpqWlpmtY2bDgQ5mNqiC5U4xsxlmNuPbb7+tRngiq6BhU+h5Q2jNBnB3T3jyb7B4QbJxieSROAnnCjNrCpwNnAOMBP5ejde0NNMytlxw9+/dvZ+7b1HR0NbuPsLdd3T3HZs3b16N8ESqYbPdwn07u5wBb94Tbhid/UzSUYnkhUoTjrs/5e4L3P0dd+/u7n9297HVeM0vgY1Tnm9EaHYtUjvUXyP0UPDX56B4bXjwSHisbxj0TaQOq/Q+HDO7kzQ1kOhazqp4HWhtZpsBJcCRwNGruC2R/NXyz3DKFHjpepg6NHSRs981oYWbpavoi9RucU6pPQWMix6TgCbAz3E2bmYPAq8AbczsSzPr6+7LgNOBCcD7wCPu/u6qBC+S9+rVh27nQb8XYe3NYdRfQ41ngYYvkLqnyjd+mlkB8Jy775GdkKpPN35KXlqxHF4dDpMuh4J6sM9lsEOf0IuBSMLy4sbPNFoDm6zuQFYHM+tpZiMWLFDLIMlDBYXQeUAYYbRlB3jq76E12/efJB2ZSE7E6WlgoZn9VPaXMPjaedkPrerc/Ul3P6Vp06ZJhyKS2dqbwfFj4aD/wNezwn07026E5cuSjkwkqyptNODu6rBTZHUzgx2Ohy33hnFnw7MX/9YZKBtsl3R0IlkRp5XaDhXNd/c3V184InVMkw3hyPvh3Sfg6YGhT7Zdzw6Peg2Sjk5ktYozPMHNwA7A24SbNrcHXiX0QOBA3jYeEKkRzGC73rB5Nxh/PrxwNbw3JtR2Nt4p6ehEVps4jQbmAH+O7uL/M9CB0Bda93xuqSZS46yxNvQeAUc/Cr8uhNv3hvGDYMmiytcVqQHiJJyt3X1W2RN3fwdon72QVp1aqUmtsNU+0H867NQXpt8MN3eGT6ckHZVItcVJOO+b2Ugz62Zmu5vZbYQbNvOOWqlJrdGwCRzwL+jzdLhn556DYczpYRgEkRoqTsI5EXgX+BtwJvBeNE1Esm3TLnDaNOhyJrz1QOgM9INxSUclskridN652N2vd/dDgL7AJHdfnP3QRASAomLY+59w8iRo1BweOhoe7QM/a7goqVni3Pg5xcyamNnawFvAnWZ2XfZDE5H3ySL9AAAUNUlEQVSVtOgApzwPe1wUajnDOsJ/H4Iqdk8lkpQ4p9SauvtPQG/gzqil2l7ZDUtE0iosgt0GQr+XYJ3W8MSpcP9fYP7cytcVSVichFPPzDYEDif0HC0iSWveBk4aD/teDZ9Pg5s7wWu3wYoVSUcmklGchHMZYSiBj939dTPbHPgou2GJSKUKCqFTv9CEeqOd4Olz4K4D4LuPk45MJK0qD0+Qz8ysJ9Bzyy23PPmjj5QTpQ5xD63YJgyCpYuh+yDofAYUxulMRCR/hyfIW7oPR+osM+hwDAx4DVrvDc9dCiP3gK/eTjoykd/UqoQjUuc13iB0Bnr4PfDTVzCiG0y6LNR6RBKmhCNSG7U9GAa8CtsfAS/+C27dFb54NemopI6LMzxBA+BQYNPU5d39suyFJSLVtsbacMgt0O5QePJMuKMHdDwF9rwEGqyZdHRSB8Wp4YwBDgaWAYtSHiJSE2y5VxjWuuPJ8NqI0Bnox5OSjkrqoDhNWDZy932zHomIZE+DxrD/UNi2N4w9A+7rDe2PgR5XQvFaSUcndUScGs7LZtYu65GISJWMnllClyGT2ez8cXQZMpnRM0sqX6lV59BLQdezQrc4w3aG98ZmP1gR4iWcrsAbZjbbzN42s1lmpraWIgkaPbOEQaNmUTK/FAdK5pcyaNSseEmnqCHs9Y/QL9ua68Ejx8HDx8HC/2U9bqnbKr3x08xapZvu7p9nJaJq0I2fUld0GTKZkvmlf5jeslkx086vwkC8y5fCy/+BKUNCr9Q9roL2R4f7eqROyYsbP9398yi5lAKe8sg7uvFT6op5aZJNRdMzKiyCXc8KY+6stw2M6R+u7/yYd78npRaIMzzBQWb2EfAZ8AIwB3gmy3GJSAVaNCuu0vRKrds6jC66/7Uw97XQku3VW9UZqKxWca7hXA50Aj50982APYFpWY1KRCo0sEcbiosKV5pWXFTIwB5tVn2jBQWh6XT/V2CTTvDMuXDnfvDth9WMViSIk3CWuvv3QIGZFbj780D7LMclIhXo1aElg3u3o2WzYoxw7WZw73b06tCy+htvtgkc+zj0Gg7ffgDDu8DUa8P1nsgqtZCTOi9Oo4HngF7AEGAd4BtgJ3ffJfvhrZodd9zRZ8yYkXQYIjXfz9+EYQ/eGwMbtIODbmL0/5ozaNQsSpcu/22x4qLC1ZfwJBF50WiA0MvAL8CZwHjgE6BnNoMSkTyx5nqhI9Aj7gvJ57Y9WPjURaxYunLjhNKlyxk6YXZCQUpNEaeV2iJgY6Cbu98NjASWZDswEckj2/QMnYG2P4rjlo/i6fqD2NE+WGmRKreQkzonTiu1k4HHgFujSS2B0dkMSkTyUPFacPAw/lZ0KfVZxmMNLuOf9e6kESHRrHILOakz4pxSGwB0AX4CcPePgPWyGZSI5K/u+x9OL7+WO5bty3GFzzGxwbnsU/R29VrISZ0Qp/POX919iUV3HptZPfL0xs+UngaSDkWk1iprGDB0QhOeWtCJfzUcyQgbAp/NgTaDw7AIImnEaaV2DTAfOB44A+gPvOfuF2Y/vFWjVmoiObTsV5g6FF66Ppx2238otO2l7nFqmHxppXY+8C0wCzgVeBq4KJtBiUgNUq8B7HERnDIFmrSER/vAw8fCwq8TDkzyTaU1nJpINRyRhCxfBtOHwfNXQWGDMN5Oh2NV26kB8qKGY2YHmtlMM/vBzH4ys4Vm9lM2gxKRGqqwHnT5G/SbBhtsB2NPh3t7wY9zko5M8kCcU2o3ACcA67h7E3dv7O5NshyXiNRk624JJzwFB1wHX74ROgOdfgusWF75ulJrxUk4c4F3vDaeexOR7CkogJ36woDp0KoLjD8f7ugB33xQ+bpSK8VpFn0u8LSZvQD8WjbR3a/LWlQiUns03QiOeRRmPQrPnAe37gq7DYQuZ0K9+klHJzkUp4ZzJaEvtYZA45SHiEg8ZrD94TDgtdBNzvNXwm3doeTNpCOTHIpTw1nb3ffJeiQiUvut2RwOuwO2OwzGnQUj94TOp0P3C8IQ11KrxanhPGdmSjgisvpsvT/0nw4djoOXb4RbdoE5LyUdlWRZ3L7UxptZqZpFi8hqU9wMDroRjh8LvgLuOgCe+jss1uGltoozPEFjdy9w92I1ixaR1W7z3eG0l8OptTfugps7wYcTko5KsiBODafGMLOeZjZiwYIFSYciIlVRv1HolaDvs9CgCTxwODx+Miz6PunIZDVS1zYikl+WLYEX/xUeDZvAftfAdoeCGaNnljB0wmzmzS+lRbNiBvZoo2GtV5O86NpGRCSn6tWH7oPg1BegWSt4vC88dDTjX36TQaNmUTK/FAdK5pcyaNQsRs8sSTpiiSlOX2r3xpkmIrJarb8t/PU52OcK+OR5uk48gINXPEvqcFylS5czdMLs5GKUKolTw9k29YmZFQJ/zk44IiIpCgphlzPgtGnMWr4pQ4pG8kDRlWxi//ttkXnzSxMMUKoiY8Ixs0FmthDYPqU59ELgG2BMziIUEVlnCwaucTmDlvZlu4LPmFD/PPoWjqOAFbRophtGa4qMCcfdB7t7Y2BoSnPoxu6+jrsPymGMIiKcs+82jC7Yh31+vYZpK7bl4qL7eaLBpVzWWZeia4pKu7Zx90FmdhCwWzRpirs/ld2wRERWVtYabeiE2Zw8/xyOa/gGF9qdNHjhMPBzoOtZ6gw0z1XaLNrMBgMdgfujSUcBM/K5lqNm0SJ1xKLvYfx5oSfq9drCQTfBRrrEvCrypVn0AcDe7n6Hu98B7BtNExFJVqN14NCRcNTDUDofbt8LJlwIS35JOjJJI+7Jz2Yp/zfNRiAiIquszb5hoLcdToBXboJbOsNnU5OOSsqJk3AGAzPN7C4zuxt4A7gqu2GJiFRRw6bQ84YwtDUGd/eEsf8Hi9XVVb6I1bWNmW0I7AQY8Kq7f53twKpD13BE6rglv8CUwaG2s+b6cOD10Ga/pKPKa/lyDQdCstkN2DX6X0Qkf9VfA/a5HP46CYrXhgePhMdOgkXfJR1ZnRana5shwN+A96LH/0Ut10RE8lvLHeCUKdD9QnhvLNy0E7z9CNTCTotrgjjNot8G2rv7iuh5ITDT3bfPQXyrRKfUROQPvnkfxpwOJTOgdQ848DpoulHSUeWNfDqlplZqIlKzrbcN9J0IPQbDnBdhWCd4/XZYsSLpyOqMWtVKTQOwiUiFCgqhc/8wwmjLHWDcWaE12/efJB1ZnVDhKTUzM2AjYBlqpSYitYk7zLwv3Ci6/FfofgF0GgCFlfb4VSslfkrNQzYa7e5fuftYdx+T78lGRCQWM9jhOBjwKmyxJzx7CYzcE76elXRktVacU2rTzUxNoUWkdmqyIRx5P/zlLvipBEZ0g8lXwLJfk46s1omTcLoDr5jZJ2b2tpnNilquiYjUDmaw7SEw4DXY7jCYOhSG7wpzX0s6slolzslK3Z4rInXDGmtD71uh3WHw5Jlw+z6wcz/Y82Ko3yjp6Gq8Sms47v55ukcughMRSUTrvaH/K7BTX3j1Fri5E3zyfNJR1XgaKk9EJJ2GTeCAf8GJz0BBEdzbC8YMCMMgyCpRwhERqUirXeC0adD17/DWgzBsZ3hfgx6vCiUcEZHKFBXDXpfCyZOgUXN4+Bh45AT4+ZukI6tRlHBERGIYPbOELvf8SOvPB3JrvWNY/sG40BnoWw+qM9CYlHBERCoxemYJg0bNomR+KUupx+CfD+CgpUP4vnhTGN0P7j8M5s9NOsy8p4QjIlKJoRNmU7p0+UrT3l26Ib1+uQj2uwY+fyW0ZHvtNnUGWgElHBGRSsybX5p2+pcLlsDOp4Ym1BvtBE+fA3ftD999lOMIawYlHBGRSrRoVlzx9LVawXFPwME3wzfvwS1d4MXrYPnSHEaZ/5RwREQqMbBHG4qLCleaVlxUyMAebX6fYAYdjoEBr8NW+8Ckf8Jte8BX/81xtPlLCUdEpBK9OrRkcO92tGxWjAEtmxUzuHc7enVo+ceFG68PR9wHh98DC7+GEd1h0mWwdHHO4843lQ4xXRNpPBwRyQu//BDG2/nvA7BOazj4JtikU9JRpZX4eDgiIlINa6wNh9wCxz4ehju4Y194+lz49eekI0uEEo6ISLZtuVdoydbxFHhtBNzcGT5+Lumock4JR0QkFxqsCftfAyeNh3oN4L5D4YnTwmm3OkIJR0QklzbpBP1egl3PhrcfDp2Bvjcm6ahyQglHRCTXihrCnpfAKVOg8QbwyPHw8LGhVVstVqsSjpn1NLMRCxYsSDoUEZHKbbg9nDw59ET94UQY1hFm3l9rOwOtVQnH3Z9091OaNm2adCgiIvEUFoWxdk6bBuu1hTH94d5D4MfaN7ByrUo4IiI11rqtoc/TsP+18OXroSXbq7fCiuWVr1tDKOGIiOSLggLoeHJoQt2qMzxzLty5H3w7O+nIVgslHBGRfNNsEzjmMTjkVvjuQxjeFaYOrfGdgSrhiIjkIzP405Ew4DVosz9MviL0yzbvraQjW2VKOCIi+WzN9eDwu0OHoIu+CT1QP/sPWJp+jJ58poQjIlITbNMTBrwK7Y+GaTeE02yfv5x0VFWihCMiUlMUrxV6nD5uNCxfEhoUjDsbFv+UdGSxKOGIiNQ0W3SH/tOhU394/fbQhPqjZ5OOqlJKOCIiNVH9RrDvYOg7Mfx//2Ew6lQonZ90ZBkp4YiI1GQbd4R+L8Ju54YbRgvqJR1RRko4IiI1Xb0GsMeF4YbRBmsmHU1GSjgiIrVFvQZJR1AhJRwREckJJRwREckJJRwREckJJRwREckJJRwREckJJRwREckJJRwREcmJ/L0ldRWZWU/gOzMrPyB4U2BBjGnrAt9lKbyKpIslV9uJu05ly1U0P9O8OOWSVJmkiyVX28n3MgF9V6qzXFXLJW5ZVadMWq3ievG5e616ACPiTs8wbUY+xZ2L7cRdp7LlKppfnXJJqkySLJd8L5Mky6UuflfillWS35U4j9p4Su3JKkzPtGwSVlcsq7KduOtUtlxF81Uu2VleZZLb7SRVLlUpq7xlUVaUiJnNcPcdk45DfqcyyU8ql/yT72VSG2s41TUi6QDkD1Qm+Unlkn/yukxUwxERkZxQDUdERHJCCUdERHJCCUdERHJCCacSZra5md1uZo8lHYsEZtbLzG4zszFmtk/S8QiY2TZmNtzMHjOz05KOR35nZo3M7A0zOzDpWOpkwjGzO8zsGzN7p9z0fc1stpl9bGbnA7j7p+7eN5lI644qlslodz8Z6AMckUC4dUIVy+R9d+8HHA7kbbPc2qAq5RI5D3gkt1GmVycTDnAXsG/qBDMrBIYB+wFtgaPMrG3uQ6uz7qLqZXJRNF+y4y6qUCZmdhDwEjApt2HWOXcRs1zMbC/gPeB/uQ4ynTqZcNx9KvBDuckdgY+jGs0S4CHg4JwHV0dVpUwsuBp4xt3fzHWsdUVVvyfuPtbddwGOyW2kdUsVy6U70Ak4GjjZzBI95te6zjuroSUwN+X5l8DOZrYOcCXQwcwGufvgRKKrm9KWCXAGsBfQ1My2dPfhSQRXR2X6nnQDegMNgKcTiKuuS1su7n46gJn1Ab5z9xUJxPYbJZzfWZpp7u7fA/1yHYwAmcvkRuDGXAcjQOYymQJMyW0okiJtufz2j/tduQslszp5Si2DL4GNU55vBMxLKBYJVCb5R2WSn2pEuSjh/O51oLWZbWZm9YEjgbEJx1TXqUzyj8okP9WIcqmTCcfMHgReAdqY2Zdm1tfdlwGnAxOA94FH3P3dJOOsS1Qm+Udlkp9qcrmo804REcmJOlnDERGR3FPCERGRnFDCERGRnFDCERGRnFDCERGRnFDCERGRnFDCEYnJzKaYWda73jez/zOz983s/nLT25vZ/hWst6OZVdjlj5l1M7OnVlesIlWhvtREcsDM6kU358XRH9jP3T8rN709YayZP3SOGW1/BjCjepGKZI9qOFKrmNmmUe3gNjN718wmmllxNO+3GoqZrWtmc6L/+5jZaDN70sw+M7PTzewsM5tpZtPNbO2UlzjWzF42s3fMrGO0fqNoUKzXo3UOTtnuo2b2JDAxTaxnRdt5x8zOjKYNBzYHxprZ31OWrQ9cBhxhZm+Z2RFmdqmZjTCzicA9qbUXM+sYxTkz+tsmzevvHm3rrWi5xtUvAZHMlHCkNmoNDHP3bYH5wKEx1tmOMGZIR8JwFL+4ewdCFyLHpyzXKBrzpT9wRzTtQmCyu+9EGH9kqJk1iuZ1Bk5w9z1SX8zM/gycSBhuoRNhrJIO0aiZ84Du7n592fLRGCeXAA+7e3t3fzia9WfgYHc/utz7+QDYLXoPlwBXpXnP5wAD3L09sCtQWumnJFINOqUmtdFn7v5W9P8bwKYx1nne3RcCC81sAfBkNH0WsH3Kcg9CGATLzJqYWTNgH+AgMzsnWqYhsEn0/7PuXn6wLICuwBPuvgjAzEYRDvoz47zBFGPdPV2iaArcbWatCd3UF6VZZhpwXXStaJS7f1nF1xapEtVwpDb6NeX/5fz+w2oZv+/zDStYZ0XK8xWs/MOsfOeDThiL5NCo5tHe3Tdx9/ej+YsyxJhu/JJVkWn7lxOS6HZAT/74fnH3IcBfgWJgupltvZpiEklLCUfqkjmEU1AAh63iNo4AMLOuwAJ3X0DoofcMM7NoXocY25kK9DKzNaLTb4cAL1ayzkIg7nWWpkBJ9H+fdAuY2RbuPsvdryY0NlDCkaxSwpG65FrgNDN7GVh3FbfxY7T+cKBvNO1ywimrt83sneh5hdz9TeAu4DXgVWCku1d2Ou15oG1Zo4FKlr0GGGxm04DCDMucGTVY+C/h+s0zlcUtUh0ankBERHJCNRwREckJJRwREckJJRwREckJJRwREckJJRwREckJJRwREckJJRwREckJJRwREcmJ/wdRC3LfTe0atAAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "%pylab inline\n", "import random, math\n", " \n", "def direct_pi(N):\n", " n_hits = 0\n", " for i in range(N):\n", " x, y = random.uniform(-1.0, 1.0), random.uniform(-1.0, 1.0)\n", " if x ** 2 + y ** 2 < 1.0:\n", " n_hits += 1\n", " return n_hits\n", " \n", "n_runs = 10\n", "n_trials = 4000\n", "for run in range(n_runs):\n", " print 4.0 * direct_pi(n_trials) / float(n_trials)\n", "\n", "#Calculation of rms deviation for 500 runs:\n", "n_runs = 100\n", "n_trials_list = []\n", "sigmasqs = []\n", "for poweroftwo in range(4, 13):\n", " n_trials = 2 ** poweroftwo\n", " sigmasq = 0.0\n", " for run in range(n_runs):\n", " pi_est = 4.0 * direct_pi(n_trials) / float(n_trials)\n", " sigmasq += (pi_est - math.pi) ** 2\n", " sigmasqs.append(math.sqrt(sigmasq/(n_runs)))\n", " n_trials_list.append(n_trials)\n", "\n", "pylab.plot(n_trials_list, sigmasqs, 'o')\n", "pylab.plot([10.0, 10000.0], [1.642 / math.sqrt(10.0), 1.642 / math.sqrt(10000.0)])\n", "pylab.xscale('log')\n", "pylab.yscale('log')\n", "pylab.xlabel('number of trials')\n", "pylab.ylabel('root mean square deviation')\n", "pylab.title('Direct sampling of pi: root mean square deviation vs. n_trials')\n", "pylab.savefig('direct_sampling_rms_deviation.png')\n", "pylab.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Bunching method to calculate the error without knowing the value of $\\pi$." ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "3.14122 1.64244240401 1.64244236172\n" ] } ], "source": [ "import random, math\n", "n_trials = 400000\n", "n_hits = 0\n", "var = 0.0\n", "var_est = 0.0\n", "Obs_exp = 0.0\n", "Obs2_exp = 0.0\n", "for iter in range(n_trials):\n", " x, y = random.uniform(-1.0, 1.0), random.uniform(-1.0, 1.0)\n", " Obs = 0.0\n", " if x**2 + y**2 < 1.0:\n", " n_hits += 1\n", " Obs = 4.0\n", " Obs_exp+=Obs/n_trials #expectation value of the observable\n", " Obs2_exp+=(Obs)**2/n_trials #expectation value of the square of the observable\n", " var+=(Obs-math.pi)**2/n_trials #calculation of the variance if the value of pi had been known beforehand \n", "var_est = Obs2_exp - Obs_exp**2 #best estimate to the variance\n", "print 4.0 * n_hits / float(n_trials), math.sqrt(var), math.sqrt(var_est)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Markov-chain Sampling:\n", "Calculation of $\\pi$ with Monte Carlo algorithm using Markov chain sampling, multiple runs." ] }, { "cell_type": "code", "execution_count": 28, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "3.109\n", "3.178\n", "2.955\n", "3.107\n", "3.225\n", "3.161\n", "2.85\n", "3.397\n", "2.841\n", "2.97\n" ] } ], "source": [ "import random\n", "\n", "def markov_pi(N, delta): \n", " x, y = 1.0, 1.0\n", " n_hits = 0\n", " n_accept = 0\n", " for i in range(N):\n", " del_x, del_y = random.uniform(-delta, delta), random.uniform(-delta, delta)\n", " if abs(x + del_x) < 1.0 and abs(y + del_y) < 1.0:\n", " n_accept += 1\n", " x, y = x + del_x, y + del_y\n", " if x**2 + y**2 < 1.0: n_hits += 1\n", " return n_hits\n", "\n", "n_runs = 10\n", "n_trials = 4000\n", "delta = 0.1\n", "for run in range(n_runs):\n", " print 4.0 * markov_pi(n_trials, delta) / float(n_trials)" ] }, { "cell_type": "code", "execution_count": 52, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "0.970458984375\n", "0.96240234375\n", "0.873046875\n", "0.747314453125\n", "0.548583984375\n", "0.24609375\n", "0.068359375\n" ] } ], "source": [ "import random\n", "\n", "def markov_pi_acceptance(N, delta): \n", " x, y = 1.0, 1.0\n", " n_hits = 0\n", " n_accept = 0\n", " for i in range(N):\n", " del_x, del_y = random.uniform(-delta, delta), random.uniform(-delta, delta)\n", " if abs(x + del_x) < 1.0 and abs(y + del_y) < 1.0:\n", " n_accept += 1\n", " x, y = x + del_x, y + del_y\n", " if x**2 + y**2 < 1.0: n_hits += 1\n", " return n_accept\n", "\n", "n_runs = 1\n", "n_trials = 2**12\n", "for delta in [0.062,0.125,0.25,0.5,1.0,2.0,4.0]:\n", " print markov_pi(n_trials, delta)/ float(n_trials)" ] }, { "cell_type": "code", "execution_count": 54, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Populating the interactive namespace from numpy and matplotlib\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAbIAAAEaCAYAAAB0PNKfAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzsnXl8VOX1/99nJsmELCRhkyXsSwgQlkASFo11wQWEFrXiVpXFpdqKtV+7ffut1rbqr1WrtdqqKHUD3CsuRVBZFDVhCSJIEVmEgOwJZF9mnt8f9yZMkpnJZLm5yczzfr3ygrnbc+567nOec89HlFJoNBqNRtNRcdhtgEaj0Wg0LUE7Mo1Go9F0aLQj02g0Gk2HRjsyjUaj0XRotCPTaDQaTYdGOzKNRqPRdGjanSMTkQEiokQkwkYb9orI+c1c9z8icn1r29RWiMi/ROSP5v/PEpEdNtggIrJIRApEJLcZ63foc6BpPt7Xbwu30+xrSET6iUixiDhbakdHR0R+IyILg1y22eeu2c5CRPYCvYHeSqljXtM3A2OAgUqpvc3dfkdFKXWx3Ta0Fkqpj4EUG5o+E5gKJCulSpq6sp3nQEQUMFQp9Y1dNmhaTlOuIfNZOF8p9YG57j4gziLT2gUi8j3gRaVUcqDllFL3tYU9Le2R7QGuqvkhImlAp+ZuzM5emKZd0R/Y2xwn1lL0NWgf+tiHFm15PlvqyF4ArvP6fT3wvPcCIjJdRPJE5JSI7BeRe7zm1YQR54nIPuCj+g2IyGVmqG+U+XumiGwTkUIRWS0iqeb0X4nIa/XWfVRE/ubPeBG5UUS2i0iRiHwlIules8eKyBYROSkiL4tItLlOkoi8IyJHzdDXOyKS7LXN1SIy3/z/DSLyiYg8aC67R0T8vumJyC9F5IBpzw4ROc+cnikin5n7/J2I/F1EorzWUyJyq4jsNNf9g4gMNtc5JSKv1CwvIt8TkXyzy3/MPLbX+LHneyKS7/V7r4j8j6/jYs7/hWnfQRGZb9o1xM+2e4vIMhE5ISLfiMiN5vR5wEJgkhme+b2PdW8QkXUi8phpx39rjlX9c9AYInKPiLwmIi+KyCngBhFxicgj5n4cNP/v8lrnRtPmE+Y+9DanrzUX+cK0fXYA2/9qns/dIjLZnL5fRI6IV0jLtOVBEdknIodF5J8i0smcF8y1+AezvSIRWSEi3fwch27m+oXmfn0sIg5z3jgR2WRu42URWSqnw883iMgn9bZVe96lGfe/iEwUkU9NW74Q4+3f3/mrYxsQXW/+JSKy2dzWpyIy2pwe8Hkhde/jwSLykYgcN++Zl0Qk0Zz3AtAPeNs857+QesMj/q51c949Ytyfz5v7sE1EJvjZ13+KyIP1pr0lInea//f5/GgMaeS+9rF8LPAfoLe5z8XmPvq6l+4RkRe91n1VRA6Z7awVkZF+2vB7PfpEKdWsP2AvcD6wA0gFnMB+jLdpBQwwl/sekIbhNEcDh4EfmPMGmMs+D8Ri9OZqpkUAc4BvgCHm8sOAEoywUyTwC3N+lNluKdDZXNYJfAdM9GP/D4EDQAYgwBCgv9e+5WKETrsA24FbzHldgcuAGCAeeBX4t9d2V2OEGQBuAKqAG017fgwcBMSHPSnm8evtdWwGm/8fD0w0j8kA0547vNZVwDKgMzASqAA+BAYBCcBXwPVe56MaeBhwAWebxzTFnP8v4I9ey+bXO+f+jstFwCGz/RiMlxxVc+587O8a4AmMB89Y4Chwntdx+yTAtXeDuQ8/M6+D2cBJoIuPc9APKAT6+dnWPeY5+gHGNdoJuBf4HOgBdAc+Bf5gLn8ucAxIN4/fY8DaeufC5z7Xs32OeU38EdgHPG5u7wKgCIgzl3/EPLddMK63t4H7m3At7sK4bzqZvx/wY9f9wD/N4xkJnIVxX0QB33od68vN4/VHf+fK+xjQ9Pu/D3AcmGauM9X83d2HzY3Zlg4cAbLMY309xjXsopHnBXWvoSGmHS7zelgLPFL/Wej1u2a/IoK41u8Bys39dZrn4XM/5ygb4xkh5u8koAzjfvT7/AjyWe7zvg6wzvfwejYEuJfuwQhB1iwzF+NadWFc25u95v3L69z5vB792hPMjgbY+fOB35qNXgSsxHjY1joyH+s9Avy13gkf5OMi+B+MB3Cy17z/A17x+u3AcEbfM39/Alxn/n8qsCuA/e8DCwLs27Vev/8M/NPPsmOBgnoPD29H9o3XvBhz33r62M4QjJvufCCykWN/B/BmvQfHFK/fG4Ffev1+CPPG47Qji/Wa/wrwfz4upjoXa6DjAjyL+YD12h+fD3WgL+AG4r2m3Q/8y+u4NebI6rwQYNyIP6p/DoK4ju/ByxGZ03YB07x+X4gR6gR4Bviz17w4jJt3gNe5aMyR7fT6nWauc4bXtOPmdSUYLxmDveZNAvY04Vr8rdfvW4Hlfta9F3irvu0YD8/6x/pTgnRkPtpp7P7/JfCCj3v1eh/basy2f2C+gHjN3wGcbf7f7/Mi0DWE8aDOq3df+HRkNH6t3wN84DVvBFDmp13BeOnJNn/fCHzkdb8F9fzwsd29BPm881rme/h2ZPXvpXvwcmT15iWaxynB/P0vr3Pn83r099caWYsvAFdjXNDP158pIlkissoMf5wEbgHqhzf2+9juXcDjSql8r2m9Md7AAFBKecx1+5iTFnN6zO5q83dN9l1NF3ibOb8vxgPLH4e8/l+KOXgrIjEi8qSIfGt2n9cCieI/Q6l2O0qpUvO/DQaClZEccAfGiT9ihm9qQlbDzG72IbPN+2h4DA97/b/Mx2/vNgtU3fGnbzGObTD4PC7m+t7n0dc5raE3cEIpVVTPhj5+lvfFAWVe8V7rB7sP9alva53rrN6261+DxRiOpym21z83KKV8na/uGC8/G80QSyGw3Jwe7LXo73zV5y8Y0Y0VYoQ7f+W1v76OdVA04/7vD/ywZn/NfT4T6OVj843Z1h/4eb1t9eX0ufT5vPCxDz3M+/GAeZxf9LEP/gjmWq9/jqLFx/iSuZ9L69n8kjnP7/MjSIK9ThrD730vIk4ReUBEdpnHca85y9ex9Hc9+qTFjkwp9S1G0sc04A0fiyzGCI30VUolYHQXpf5mfKx3AfBbEbnMa9pBjIsTMNK0MS7MA+akV4HviTFOMMtsG6XUx0qpOPOvJia7Hxgc9I6e5ucY3fgspVRnjLdCfOxTk1FKLVZKncnp8Oz/M2f9A/gvRjZcZ+A3LWwvyYxz19AP49i2hO8A7wymvgGWPQh0EZH4ejYc8LO8L/qY5997/ebuQ/3rr851Vm/b9a/BWIwQX1NsD5ZjGE5tpFIq0fxLUErVPGRa7VpUShUppX6ulBoEzADuNMdYvsP3sa6hBMPZGg2L9Ky36abe//sxemSJXn+xSqkHfJjdmG37gT/V21aMUmqJOd/n88IH95s2jjaP87X19sHX86uG1rjWvVkCXC4i/TFCpq/XGuH/+WEF/vY50LG4Gvg+Rq8xAaPnCj6u1wDXo09a6zuyecC5yneWWTzGG0m5iGRi7EwwbMMIVz4uIjPNaa8A00XkPBGJxLiRKzDCCSiljmKEBBZhhF+2B9j+QuB/RGS8GAwxL47GiMd4uBSKSBfg7iD3JyAikiIi54qRVFButuH2avMUUCwiwzHG2lrK70UkSkTOAi7BuKlbwivAHBFJFZEY4Hf+FlRK7cc4Z/eLSLQYA/DzMN8ug6QHcLuIRIrIDzHGad9rvvl1WILxEtVdjOSI32G8hYPxsJsjImPNc3UfkKNOf2pyGGNsssWYEYengb+KSA8AEekjIheai7TatShGUsQQ0ymcwrj23MBnGKHo20UkQkQuBTK9Vv0CGGkej2iMHoE3Tb3/XwRmiMiF5ht8tBhJR77SvBuz7WngFrNXKCISK0bySTw06XkRDxRjHOc+GNEib/ye81a61r23l4cxxrYQeF8pVQiNPj+s4DDQVUQSmrBOPMbz+jjGy4/f1PwA16NPWsWRKaV2KaU2+Jl9K3CviBRhPBBeacJ2v8B4yD4tIhcrpXZgvA09hvG2OgOYoZSq9FptMYbH9/d2VbPtV4E/mcsVAf/GGOhsjEcwBjGPYSQELA92fxrBBTxgbvcQxoP6N+a8/8F4ABRh3Jwvt7CtQ0ABxtviSxgDu/9tyQaVUv8B/gaswggJfGbOqvCzylUYb2QHgTeBu5VSK5vQZA4wFON4/Qm4XCl1vP5Ccvrj1H715wXgj8AGYAvwJbDJnIZS6kOMsdrXMXoEg4Ervda9B3jODGVd0YQ2/fFLjOP5uRmO+YDT3/a15rU41Nx2Mca5e0Iptdq8ty7FGDoowEisqY28KKW+xhjP+ADYiTHu5E2T7n/zwf99jGv/KEav6i58PKuCsG0DxjjS383535jLehPM8+L3GIkjJ4F3aRh5uh/jxadQRP7Hx/otvdbrs8SHzX6fHyJyjZweUmkVzOfFEmC3ud/BhDGfxwirHsDIf/g8wLI+r0d/C9dkv2jCBAnyQ8ZWaCcV2Aq4lFLVrbztGzAG4s9sze1qgkNE/oUx0P9bu23RaKAdlqjSdFxEZJYZrkzCiM+/3dpOTKPRaOqjHZmmNbkZIxy0CyOe3RpjeRqNxgbEKJpQ7OPvP3bbVh8dWtRoNBpNh0b3yDQajUbTodGOTKPRaDQdmpCqNi0iM4AZ8fHxNw4bNsxuczQajabDsHHjxmNKqe5229EcQnKMbMKECWrDBn+ftWk0Go2mPiKyUSnls/J+e0eHFjUajUbTodGOTKPRaDQdGu3INBqNRtOhCalkj/ZESUU1T67dxYuffUtBaRVJMZFcO6k/N2cPJtalD7tG0xSqqqrIz8+nvLzcblM6PNHR0SQnJxMZGWm3Ka2GfqJaQElFNbOeWMe3x0upqPYAcKK0iifX7Gb51kO8eesUy5yZdqCaUCQ/P5/4+HgGDBhAXdUWTVNQSnH8+HHy8/MZOHCg3ea0Gjq0aAFPrt3F0eMn+LF6hY2um9ntupqNrpv5sXqFo8dP8OTaQHqezafGgT65ZjcnSqtQnHags55YR0mFLnuo6ZiUl5fTtWtX7cRaiIjQtWvXkOvZ6ld0C3j90/+y1PFb+sthoqUKgK4UcUvE21yscpn72f3cOTWlka00nSfX7qrTC6yhotrDt8dLeXLtLkvarUH3BjVWop1Y6xCKx1H3yCzgisp/13FiNURLFf3lMFdUvmlJuy9+9m0DJ1ZDRbWHFz/fZ0m7oHuDGo3GPrQjs4DrIlc2cGI1REsVP4r4wJJ2C0p9t3l6fmXA+S0hmN6gRqPRWEFIOTIRmSEiT508edJWOxIpDjg/iSJL2k2KCZyFlBQTZUm7YG9vUKPxpqSimodX7iD93hUM/NW7pN+7godX7miVqMDJkyeZNWsW48ePJy0tjYULFzZrO8uXLyclJYUhQ4bwwAMPNHmZwsJCLr/8coYPH05qaiqfffYZ+/fv55xzziE1NZWRI0fy6KOPNsu2jkhIDVwopd4G3p4wYcKNthrSqQuUHfc7W8V0wYoo9bWT+vPCmq+4Xi3jRxErSaKIAuJ5oXoqz8lMrp04xIJWDezsDWo0NVidMfz6668THx/Pxo0bASgrK2vyNtxuN7fddhsrV64kOTmZjIwMZs6cyYgRI4JeZsGCBVx00UW89tprVFZWUlpaSllZGQ899BDp6ekUFRUxfvx4pk6dWme7oUpI9cjaC5I5HxXh8jlPRbhwZMy3pN2bJ57Bm1G/45aIt+kqRTgEuoqRZPJm1O+4eeIZlrQL9vYGNZoarA5xp6ens2bNGiZMmMDdd9+Ny+X7Pg9Ebm4uQ4YMYdCgQURFRXHllVfy1ltvBb3MqVOnWLt2LfPmzQMgKiqKxMREevXqRXp6OgDx8fGkpqZy4MCBFu1vR0E7MiuYfDuSNAgioutOj4g2pk++3ZJmYzc8QX+HnyQTx2FiNzxhSbtg9AZdEb4vJ1eEg2sn9rOsbStDSZqOhZUh7pMnT/KLX/yCLVu28Pnnn7Nq1aoGDuiss85i7NixDf4++OD0uPiBAwfo27dv7e/k5OQGDifQMrt376Z79+7MmTOHcePGMX/+fEpKSuqsv3fvXvLy8sjKymr2/nYkQiq02G5wxeH50fsUvfw2JbuT8Kg4HFJMbJ8C4mfPwOGKs6bd9QuR6gqfs6S6AtY/A+f8xpKmb84ezPKthxq8DbsiHPTvGsPN2YMtadfOj8817Q8rQ9xPPvkkF154IQkJCQBMmjSJQ4cO1Vnm448/bnQ7vhRH6qfEB1qmurqaTZs28dhjj5GVlcWCBQt44IEH+MMf/gBAcXExl112GY888gidO3cObuc6OLpHZgGeCjdHnvmGom8H4lGdAQce1Zmibwdy5Jlv8FS4rWm49ETg+QHG7VpKrCuCN2+dws1nD6JLbBQi0CU2ipvPHmSpM9HZkhpvrAxx5+XlMXLkyDq/09LS6iwTTI8sOTmZ/fv31/7Oz8+nd+/edbYTaJnk5GSSk5Nre1uXX345mzZtAoxSXpdddhnXXHMNl156abP3taOhX1UtoGhtPtXHy6C63ltVtYfq42UUrc0nYWr/1m84pguUBnBWnbq2fptexLoiuHNqiqUfXdcnmFBSW9qjsZdrJ/XnyTW7fV4TLQ1xJyUlkZeXx0UXXcS7777LqVOnmDx5cp1lgumRZWRksHPnTvbs2UOfPn1YunQpixcvDnqZnj170rdvX3bs2EFKSgoffvghI0aMQCnFvHnzSE1N5c4772z2fnZEdI/MAoo/O9jQidVQrYz5VpAxv+G4XA0R0ZAxz5p2bURnS2q8uTl7MP27xjQYr22NEPddd93Fm2++yZgxY3j66ad54403cDia/giNiIjg73//OxdeeCGpqalcccUVtT29adOmcfDgwYDLADz22GNcc801jB49ms2bN/Ob3/yGdevW8cILL/DRRx/V9gTfe++9Zu9vR0IrRFvA/l+tRQIk2CsUfR/Ibv2GK4ph4flQsAeqvWqpRURD0kCY/wFYNT5nE+n3ruBEAGfWJTaKTf83tQ0t0ljB9u3bSU1NDWrZ2lJpn++joLSSpJgorp3YT5dK88LX8ezICtH6rFpAkbOUzu5Yv/NPOUv8zmsRrjjDWX36NyOxo+y4EU7MmGdkSoaYEwNrQ0majokdIW6NvWhHZgFvJ67m8hNTcamGA8sVUsk7iWsYycXWNO6KMzITLcpObG/YlS2p0WjaD3qMzAI+6L2B7yKPUSF1x2cqpJLvIo/xQW/7wp6hhl3ZkhqNpv2g73IL+MGIWfyq+lFmHMlmemE2nd2xnHKW8G7iWt7usZarRlxjt4nWUFFshjUXGp8CxHQxElAsDmvqUJJGE95oR2YBc0bNYeW+lbzq/IAXe7xbO93ldJEcn8ycUXMsa9tT4aZobT4lnx3EU1qNIyaC2Em9ic9OxuFyWtauz0ST0uOw7lH4allIJpqA1mDTaNoDOrRoATGRMSyetpg5I+eQ5EpCEJJcScwZOYfF0xYTExljSbueCjdHnthM0Zp8PKVGeSZPaTVFa/I58sRm6z7EBqMnVj9bEozfBXuM+SGG1mDTaNoHIfXKKCIzgBlDhlhX5T1YYiJjuG3cbdw27rY2a9P4ELsc6mfwVXuoPl5u3YfYYIQT6zux2vbLLS2PZRd2K3JrNBqDkOqRKaXeVkrdVFMLLdwo+exgQydWQ7WHks+/s65xG8tj2YXWYNNo2gch5cjCnZpwov/5gatgtIiYLoHnW1weyw50VZF2SkUxrLoP/jwI7kk0/l11nzG9hbSlsObcuXPp0aMHo0aNqjM9kIDmgAEDSEtLY+zYsUyY0CG/bW4WIRVaDHccMREBnZmjkYKqLSJjvpHY4Su8GKLlsZJiIgNWFdEabDZgcdJRWwlrAtxwww385Cc/4brrrqszPSIiIqCA5qpVq+jWrVsz97BjontkIUTspN7gRxOMCAexE3tZ1/jk240yWD402EgaaJkGm53YqcGm8YPFSUdtJawJkJ2dTZcuDSMd4Syg6Q/tyEKI+OxkHF2icDvqjtu4HR4cXaKIz062rvGa8lhTFkBMNxAx/p2yIGRT760sUKtpJsEkHTWTthTWDJb6ApoiwgUXXMD48eN56qmnmrXNjogOLYYQ5Y4K7hjwZzLVMC46MaX2Q+zlXdaRO+Br/uV4jhisSf0Hwq48Vk1VEV2gth1hYdJRWwprBoMvAc1169bRu3dvjhw5wtSpUxk+fDjZ2RYUKG9n6DsthFi0dRG7S/ewvdt/ea7bsjrzXKUuFm1d1KafA4QDuqpIO8NCTb68vDyuv/76Or9nzJhRZ5mzzjqLoqKiBus++OCDnH/++UBwwpqN4U9As2Y7PXr0YNasWeTm5oaFI9OhRYvwlJRw5LHH+HrSZLanjuDrSZM58thjeEosqnwPvLzjZSrcFT7nVbgreHnHy5a1rdG0CyzU5KsR1gQCCmtu3ry5wV+NE4O6opmVlZUsXbqUmTNnBm2HPwHNkpKSWidaUlLCihUrGmQ8hirakVmAp6SEPbOv5MTCZ3AXFIBSuAsKOLHwGfbMvtIyZ1ZYUdii+RpNh8fCpKO2FNYEuOqqq5g0aRI7duwgOTmZZ54xxvf8CWgePnyYM888kzFjxpCZmcn06dO56KKLmr2/HQktrGkBRx57jBMLn0FVNOwdictFl/nz6PHTn7Z6u9lLsymoKPA7P8mVxNor17Z6uxqN1TRFWPN08erw0ORrDlpYU9MohYuX+HRiAKqigsIlSyxxZLNTZrNo2yKf4UWX08XslNmt3qY3thUs1mi8CbOkI412ZJbgLgwcwnMXWBPiq6m6n1+UX8eZtVXV/SNPbK5T67GmYHHZ1mP0uHWsdc7MJvkYO9FV9zWa0+gxMgtwJiYGnp8UeH5zsavqPgRXsNgSaio5rHvUzFZTpys5LDy/VcoStTd01X2Npi7akVlA4tVXIX6++BeXi8SrrrKs7Zqq+2uvXMuW67ew9sq13DbuNkudGNhYsDgM5WOCqbqv0YQTjToyEekuIr8RkadE5Nmav7YwrqPSbe5cIvv2beDMxOUism9fus2da5Nl1mFbwWILKzm0V3TVfY2mLsEE098CPgY+ACxUZgwdHLGx9HnuX6y9/1527N1JpUOI8ihSBgwl+9e/wxEba7eJrY5tBYvDUD5GV93XaOoSTGgxRin1S6XUK0qp12v+LLesA1NZXsaSP/2WbQe/pdLpABEqnQ62HfyWJX/6LZXlTa+Y3d6xrWBxGMrHJDXyUqCr7mvCjWAc2TsiMs1yS0KInDde48TB7/C46/ZQPO5qThz8jpw3XrPJMuuIz04momt0Q2cW4SCia7R1BYstrOTQXtFV9zWaugTjyBZgOLNyESky/05ZbVhHZtPyd0D5CbOpavKWv9u2BrUBDpeTHreOJf7sZByxkSDgiI0k/uxka1Pvw1A+RlfdD0xpVSmP5z1O9tJsRj83muyl2Tye9zilVaUt3nZbCWsGEs+E8BXQ9EejY2RKqfi2MCSUqK4IXIKqKgRTwsFwZglT+5MwtX/bNVojHxNGlRx01X3/lFaVcvV7V9f5lrKgooBF2xaxct/KFn+G0lbCmo2JZ0J4Cmj6I6grXkRmAjUllFcrpd6xzqQQQDqBCnCBS6e2s6UNKa0qZdHWRby842UKKwpJdCUyO2U2c0bNsTb9PwwrOeiq+75ZtHVRg4IAYBTNzi/Kb7ECRHp6Or///e+ZMGEC06dP5+67727yNryFNYFaYU1vJ9WrVy969TLGlb3FM+urSGsMgkm/fwAjvPiV+bfAnKbxgyt+HOAvlOY054cWNW/Ci7YtoqCiAIWqfRO++r2rWyWso9E0hpUKEHYJa9YXz4TwFdD0RzA9smnAWKWUB0BEngPygF9ZaVgNIhILPAFUYvQGX2qLdlvC+Itn8dlrO1DuQup+seBEnImMv3iWXaZZhtVvwhpNMFipAGGHsKYv8UwIXwFNfwRb2cO7plJCSxs1P6o+IiJb602/SER2iMg3IlLjKC8FXlNK3QgEL9pjI+OnDeWMofOJjMk8HUaUTkTGZHLG0PmMnzbUXgMtQGuhadoDia7A5d8amx+IvLy8WrmVmt9paWl1lgmmRxassKY/8UzwLaAZzgTTI7sfyBORVYBgjJX9uoXt/gv4O/B8zQQRcQKPA1OBfGC9iCwDkoEvzcU6xAfZUdER/PDXk9i8sg9frsmmvKSK6NhI0s7uw9ip/YiKDr3BeK2FpmkPWKkAUSOsedFFFwUU1mwMb2HNPn36sHTpUhYvXlxnGX/imWCIZno8HuLj42sFNH/3u981e79CgWCyFpeIyGogA8OR/VIpdSjwWo1uc62IDKg3ORP4Rim1G0BElgLfx3BqycBmAvQgReQm4CaAfv3s/44mKjqCzBmDyJwxyG5T2oREV2JALbSWvAlrNMFipQLEXXfdxezZs1m6dCkDBw5sFWFNt9vN3Llz6whrLly4kN27d/PCCy/UptgD3HfffUybNo3Dhw8za5YxPFFdXc3VV18dNgKa/vArrCkiw5VS/xWRdF/zlVKbWtSw4cjeUUqNMn9fDlyklJpv/v4RkAX8EqP3Vg58EswYmd3CmnbiKSnh2LPPUrh4Ce7CQpyJiSRefRXd5s61tDTW43mPB3wTnjNyjmVjZFoHrW2xQ0KmKcKatmXPdiDCSVjzTowezkM+5ing3Fa2xddop1JKlQDWCWmFEJ6SEvbMvpKq/ftrhT3dBQWcWPgMRe+vYODLSy1zZnZpodmqgxaG1EjIeFffr5GQWb71EG/eOsX279hqFCB0clH44LdfrJS6yfzvxUqpc7z/MDIZW5t8oK/X72TgoAXthCzHnn2Wsvx8diTGsHLkAN4bPYiVIwewIzGGsvx8jj1rnWiBXVpotumghSlaQkbTHgnm1elToH540de0lrIeGCoiA4EDwJXA1U3ZgIjMAGYMGTKklU3rGBxbsoR1/bpR6orEY8buqyKc7O6RyKGEWM5auoQeP/2pZe3b8SYcjA6apZVGwkydOhgJGf2RtqatCZQ80VNExgOdRGSciKSbf98DWvR6LSJLgM+AFBHJF5F5Sqlq4CfA+8B24BWl1LambFcp9bZS6qaa7zwLxyUdAAAgAElEQVTCjZ1R1HFiNXgcDkpdkey0SEnFTmzTQYOwVKfWEjKa9kigHtmFwA0YIb6HvaYXAS2qBaSU8imRrJR6D3ivJdsOZ/Z1S2zgxGrwOBzs6xZ6mYO26aBBcOrUIVY2KykmkhMBnJmWkNHYQaAxsufM8bAb6o2RzVRKvdGGNmqCpNIZOBW4sfkdEdt00CAs1am1hIymPRLMd2Svi8h0YCQQ7TX9XisN0zSdTnHxlBUX+Z8f39nvvI5KfHYyZVuPNUz4sFoHDcJSnfrm7MEs33qoQcKHlpDR2EkwRYP/CcwGfoqRIv9DoA11OoJHRGaIyFMnT5602xRbGHPhJTgjfYfSnJGRjLlgehtbZD226aBBWKpT10jI3Hz2ILrERiECXWKjuPnsQe0i9V4TngRz1U1WSo0WkS1Kqd+LyENAuwwtKqXeBt6eMGHCjXbbYgcZMy/l688/ofBQXXVqhzOChB49yZh5aYC1Oy626KCBkZ247lHf4cUQVaeG9i8hY2VRgJMnT3LDDTewb98+KisrWbBgAfPnz2/ydpYvX86CBQtwu93Mnz+fX/2qYQ32AQMGEB8fj9PpJCIignAt8hAMwTiyGmGtUhHpDRwHBlpnkqb5RBIVfzXOgrV4yjYbmmjSCadrLFHx2UAIpi3ayeTb4atlDRM+Qlidur1jdVGAthLWrEGLZwZHMKP/74hIIvAXYBOwF1hqpVGa5rF55T6KTrhxuiYRnfhjopPuJDrxxzhdkyg64Wbzyn12mxha1KhTT1kAMd1AxPh3ygJjegh+R9beOfbss3WcWA2qooKq/ftbXBQgPT2dNWvWMGHCBO6++25cLleTt+EtrBkVFVUrrKlpPsEke/zB/O/rIvIOEK2UCs9BqHbOl2sO4K7y/bGqu8rDl2sOhE0R4zYjDNWp2zOFi5c0cGI1qIoKCpc0vyiAt7BmbGws5557LmPHjq0t4AuGjEtRUcOEqwcffJDzzz8f8C2smZOT02CdGvFMEeHmm2/mpptuarCMxsCvIxORc5VSH4lIg4EVEaE9puCHe2WP8uLAH6uWl1j4cbBG0w5wFwaWC3IXdBxhTS2eGTyBemRnAx8BM3zMU7TDhI9wT/aIjosM6MyiY0NzjExXO9fU4ExMxF3gX07ImdQyYc3rr7++zu8ZM+o+HoPpkQUrrOlLPFM7Mt/4dWRKqbvN/85XSnUIQctwJ+3sPmxasc9neNEZ6SDt7D42WGUtpVWlXP3e1XWq7hdUFLBo2yJW7ltpacFiTfsj8eqrOLHwGZ/hRXG5SLzKZ1GhoGhLYU0tntk0gkn22CMiT4nIeeKr/6tpN4yd2o+Ebp1wRtY9rc5IBwndOjF2auhVXVi0dVED6RiACncF+UX5LNq6yCbLNHbQbe5cIvv2ReolYYjLRWTfvnSbO7fZ277rrrt48803GTNmDE8//XSrCGumpqZyxRVX1BHWPHjwIIcPH+bMM89kzJgxZGZmMn369LAXzwyEX2HN2gVEOmGEF6/EqHj/DrBUKfWJ9eY1j3AW1qwsr2bzyn18ueYA5SVVRMdGknZ2H8ZO7UdUdOh9rJq9NDugMnWSK4m1V661rH0t6tk2NEVYs/Y7siVLcBcU4kxKJPEq68VlOxLhJKwJgFKqDHgFeEVEkoBHgTWAvkvbIVHREWTOGGRLdmL58eOsvf9eduzdSaVDiPIoUgYMJfvXvyO6qzVVLgorAg/eNza/JWhRz/aJIzaWHj/9qaWSRZr2RVD9YhE5W0SewPiOLBq4wlKrmkm4l6iyk/Ljx3n+luvZ9u1OozixCJVOB9u+3cnzt1xP+XFr6g4mugIP3jc2vyXYKupZUQyr7oM/D4J7Eo1/V90XktIxNbg9ikOnyvnq4Em25Bfy1cGTHDpVjtsTOKqkCX2CqbW4B7gD+BgYpZS6Qin1uuWWNYNw1yOzk7X330sJHp9aaCV4WHu/NTWmZ6fMxuX0/VGqy+lidspsS9qF4EQ9LSEMddA8SrHraDHHiiqoNh1XtUdxrKiCXUeLtTMLc4LpkY1RSs1SSi1RSpVYbpGmQ7Jj786AWmg79u60pN05o+aQHJ/cwJm5nC6S45OZM2qOJe2CjaKeweighRjF5dVUVnvw1BvT9yhFZbWHo8W+P4LWhAfBOLKeIvKhiGwFEJHRIvJbi+3SdDAqHYETWhub31xiImNYPG0xc0bOIcmVhCAkuZKYM3KO5an3jpjAQ8yWiXqGoQ5aSUV1AydWg0cpThRrZepwJpg0tqeBu4AnAZRSW0RkMfBHKw2zg1WrVhEbG8uYMWOaVUMtnInyKCqd/p1VlIWhn5jIGG4bdxu3jbvNsjZ8ETupN0Vr8n2HF60U9QxDHbTGLp9qj58QryYsCKZHFqOUyq03LXBMpQPi8XjYvXs37733Hg8//DDLly/nuEUJCqFIyoChOPw8TBweDykDhraxRdYTn51MRNfohgrVVot6hqEOWmMd+ohmfM+lCR2COfvHRGQwRlkqRORywKJRbPtwOBzMnTuXefPmMWzYMHJzc3nsscd46aWX+Oabb/DoN76AZP/6d8TiaODMHB4PsTjI/nXoVSWwTdQzY74hFeOLENVBi3VF4PBTj8EhQpe4qDa2SNOeCOaD6EHAU8BkoADYA1yrlNpruXVNxKto8I07d7YsuaCoqIgNGzawYcMGSkpK6Nq1K5mZmYwZM4boaD8PkTCn6OARlt37Vw4X7kSpckSiOSNxKDN/9zPie/ew27zQoSZr0Z8OWghKyGz76isiu/ZtkPDhECEqwsHg7nE4zW5bnaIAxVVEx7VeUYDWEtacO3cu77zzDj169GDr1q0+lwlGfLO5hNoH0Y06stoFRWIBh1KqYUXMdkZrVvaorq7mq6++IicnhwMHDhAVFcXYsWPJzMzUgndeVJZX8/r/28jJY2V1aj3WlMe67JfjQ7KyiG1UFBvZieufMcbEOnU1emKTbw85JwbGg3dYynCOFldworiSao+HCIeDLnFRdI9z1XFiVl6Hzz77LKtXr+b5558HDGHNTp06NXk7a9euJS4ujuuuu86nI3O73QwbNqyO+OaSJUt8im82h1BzZIFkXO70Mx0ApdTDFtnUroiIiGD06NGMHj2a/Px8cnNz2bBhQ604XmZmJkOGDGlWzbVQYvPKfQ0eHmDooJ08VsbmlftCUgvNtsr7YaiD5nQIPTtH07Oz/4iI1ddheno6v//975kwYQLTp0/n7rvvbnwlH2RnZ7N3716/873FN4Fa8c3WcmShRqBXk3jz3xQgA1hm/p4BWFe8rh2TnJxMcnIyF1xwARs3bmT9+vUsXryYLl26kJmZydixY8M27BiOop52Vt7XNR59Y+V12FrCmsEQrPimxiCQjMvvAURkBZBeE1IUkXuAV9vEunZKXFwcZ599NlOmTGH79u3k5uayfPlyPvzww9qwY/fu3e02s00JR1HPYCrvW/FJgK7x6B8rr8PWEtYMhmDFNzUGwQSL+wHeXxtWAgMssaaDERERQVpaGmlpaRw8eJCcnBw2bdrE+vXrGTRoEFlZWQwdOjQswo7hKOr58o6XGzixGircFby842VLHFkwNR4TpvZv9XZrqR2fW2h80xbTxcikbAfjc1Zeh60lrBkMwYpvagyCcWQvALki8iZGCv4s4DlLreqA9O7dm1mzZtUJOy5ZsoSkpCQyMjIYN25cswaFOwp2inpWlpexftkbfLHiXcqKi+gUF8+YC6aTMfNSoqKtO+Z2Vd4PpsajZY7MV8ZkTZ3Hr5bZnjFp5XXYWsKawRCM+KbmNI12FZRSfwLmYKTeFwJzlFL3W21YRyU2Npbs7GzuuOMOfvjDHxIfH8+KFSt4+OGHefvttzly5IjdJlqCXaKeleVlLP7fn7N+2euUFZ0CpSgrOsX6Za+z+H9/TmV5mSXtgn2V922r8Qjtvs6jVdeh26O44ZbbWfzKa6SMGMXDjz3BP59bjKJ54b6rrrqKSZMmsWPHDpKTk3nmGaOsWI2wZiDxTU1Dgk6/7wi05ndkrcl3331Hbm4uX375JdXV1QwcOJDMzExSUlJCKuxoh6jnuldeYv2y13BXNXx4OyMjyZh5OVOuuMaSth/Pe5xF2xb5DC+6nC7mjJxjSWjx4L2fBXRmjthIev/fxFZvFzDkYkoDVLyJ6Qa/2NXqzTZFWLO1r0O3x6i8H8w3bB2FUEu/DylHVkN7VYguKSmpHUM7deoUCQkJZGZmMm7cOGJiLEzVDmGemHcVZcX+P23sFN+ZWxdaE5LxlbUIpyvvW5W1eHLltwFrPMafnWxdaPGeRMwiP74RgbtbP6TaFEfW2hw6Vc6xogqfRYsdInSLdwX8JKA9EmqOLHS6Ax2A2NhYzjrrLBYsWMAVV1xBYmIiK1eu5OGHH2bZsmUNMqA0jRPIiQFGuNEi7Kq8b1uNRwjLOo8nin07MdCV99sLQfWzRaQ/MFQp9YGIdAIiOkKFj/aK0+lkxIgRjBgxgkOHDpGbm8uWLVvYtGkT/fv3Jysri5SUFJzO8EyhbgpRbo+hSB1gvpXYUXm/psZj0dp8Sj7/Dk9pFY6YSGIn9rL+O7KM+UZihy8ZmRCt81jdSOl9XXnffhp1ZCJyI3AT0AUYDCQD/wTOs9a08KBnz57MnDmT888/n7y8PHJzc3nllVfo3LkzGRkZpKenExsba7eZ7ZZ+xwrZ3T3Rp6inw+Oh3zFrMgftxuFykjC1v7Vp9r6YfLuRneivzuPk29vWnjYgwiEBnZmuvG8/wfTIbgMygRwApdROEdEVYFuZmJgYpkyZUpvJlJuby4cffsiaNWtIS0sjMzOTXr0s0rfqwAythEMVVZS6Ius4M4fHQ0xFFUND7ztse3HFGSn2YVTnsUucK+AYma68bz/BOLIKpVRlzVflIhJBwNFeTUtwOBykpqaSmprK4cOHa8OOeXl59OvXj6ysLIYPH67DjibdrrqKKc88y66ETuzrlkCl00GU20O/YycZfLKMbvPm2m1i6BFmdR67x7k4VVblN2uxe5wW4bWbYBzZGhH5DdBJRKYCtwJvW2uWBuCMM85gxowZtWHH9evX8+qrrxIfH09GRgbjx48P+7Bjt7lzKVixikjPUFwJk3FExhFZVUxk+adEdt5Jt7nakWlahtMhDO4e12jlfY19BOPIfgXMA74EbgbeAxZaaZSmLp06dWLy5MlMnDiRnTt3kpOTw0cffcSaNWsYNWoUWVlZYVu+ptrpYmP6XZw8UoJHGaHFqqh4vu03lcIeP2Cg04UO/LQu4ViwOJjK+xr7COjIRMQJPKeUuhZ4um1M0vjD4XCQkpJCSkoKR48eJTc3l82bN/PFF1+QnJxMVlYWI0aMCKuw4+aV+zh1oqLWidXgUQ5OnagIWfkYu+gIBYutLFnWGsKa+/fv57rrruPQoUM4HA5uuukmFixY0GA5K4U1Q42Ajkwp5RaR7iISpZTSH0u0I7p378706dM577zz2Lx5M7m5ubz++uu8//77TJgwgQkTJhAXF3oD7/UJR/kYsE8HzfaCxY1QU7Ks8PAh3FXGI6umZNnOnHVc/aeHWuTMXn/9deLj49m4caOx7bKml0CLiIjgoYceIj09naKiIsaPH8/UqVPraI253W5uu+22OsKaM2fO1HpkfggmtLgXWCciy4CSmontUVjTq0SV3aa0GdHR0UycOJHMzEy++eYbcnNzWb16NR9//DEjR44kMzOT5GQLP5C1GbvlY+woWGynDpqtBYuDYP2yN+o4sRrcVZUUHj7E+mVvtKhkWWsIa/bq1as2Azk+Pp7U1FQOHDhQx0lpYc2mEYwjO2j+OTgtttkuUUq9Dbw9YcKEG+22pa1xOBwMGzaMYcOGcezYsdqw45YtW+jTp09t2DEiwpqah3Zhp3yM1W///rBLBw1sLljscUPxESg9Bp5qcEQYtR3jeoDDCGd+seLdBk6sBndVJV+sfK/ZjswKYc29e/eSl5dHVlZWnelaWLNpNPpUqxHY1HQcunXrxrRp0zj33HP54osvyM3N5Y033qgTdoyPb9fvJEFjp3yM1W///rBLBw3AERMRuGBxjEUvDsoDx76G6gpqv/7xVEPxYSgvhG7DwOG0tGRZawtrFhcXc9lll/HII4/QuXPnOvO0sGbTCKayR3fgF8BIoDZlRyl1roV2aVqB6OhosrKyyMjIYPfu3eTk5LBmzRo+/vhjRowYQVZWFsnJyR36Bhk7tR+7Nh3l5LGyOs7MavkYsPbtPxB26aABxE7qHbBgcexEiz7aryiC6mgafsKqDOdWfAQ696JTXHxAZ9UpvrPfeY3RmsKaVVVVXHbZZVxzzTVceumlDZbXwppNI5g400vAy8AlwC3A9cBRK43StC4Oh4MhQ4YwZMgQjh8/zvr168nLy2Pr1q307t2bzMxMRo0a1SHDjlHREVz2y/FtLh8D9hUsTnQlUlBREHC+VcRnJ1O29VjDhA+rCxZXFAP+PjxWRrixcy/GXDCd9cte9/mC4YyMYszUac02obWENZVSzJs3j9TUVO68806fy2hhzaYRzF3eVSn1jIgsUEqtwfhAeo3VhmmsoWvXrlx00UWcc845tWHHf//736xYsaI27Fg/zNHeiYqOIHPGoDbPTuwUGxdYQibOmvDt7JTZAXXQZqfMtqRdsLFgsXIHnu8xwp0ZMy9lZ866BiFfZ2QUiWf0JGNmw95PsNx1113Mnj2bpUuXMnDgQN54441m6QmuW7eOF154gbS0NMaOHQvAfffdx7Rp05g2bRoLFy6kd+/etcKabrebuXPnamHNADSqRyYinyulJorI+8DfMBI/XlNKDW4LA5tDe9Uja48opWrDjl9//XVtiaysrCz69u3bocOOVrPifxaw7dudfgsWj+w/lAsefLTV27VLB81Otn/+Aan9AkjEOCKgZxrglUm68j3Kik7RKb4zY6ZOszSTtKMRanpkwfTI/igiCcDPgceAzsDPLLVK02aICIMHD2bw4MGcOHGiNuy4bds2evbsSVZWFqNGjSIy0rrsv45K8qfr2dMt1m/B4uTP1lvSbo0Omh3fkdmGKw4QfJd5FSN70SQquhNTrrjGMmVwTftDK0RrGlBZWcmWLVvIycnh6NGjxMTEkJ6eTkZGRm3Glga2p46gWmB398QGBYsHHS0kQkHq9q/sNjMk2P7VNlK7OepmLQIgEOGqzVrUBEfY9chEZBE+XoOUUroaa4gSFRXFhAkTGD9+PHv27CE3N5d169axbt06UlNTyczMpH///mEfdnQmJkJBAcMOG38N5ndJssGqEEUcqK5DcZ84iafCCTgBNw6XG2eXBEQ7sbAmmNDiO17/jwZmYYyTaUIcEWHQoEEMGjSIgoICNmzYwMaNG/nqq68444wzyMrKIi0tLWzDjolXX8WJhc+gKhomXYjLReJVV9lgVWiilKL6WAWq2sXp92onnsoI1LEKIro7EV2FPmxpcmhRRBzAB+35OzIdWrSOyspKvvzyS3Jzczl8+DCdOnWqDTsmJlqX9t0e8ZSUsPPKH7HLM5T8MyZTZUrIJB/+lMGOnQxd+gKOEJTZsaPO47ZNXzL0jIHg63klgiM+kojOWhcsWMIutOiDoYB1X5lq2jVRUVGMHz+e9PR0vv32W3Jycvj000/59NNPSUlJISsriwEDBoRF2NFOCRk7ajyCfXUeVYXbtxMDUApPcRVoRxa2BDNGVoTRl69JGToE/NJiuzTtHBFhwIABDBgwgMLCwtqw43//+1969OhBZmYmo0ePJioqdNXA7JKQsavGI9hX57HRyJEn9JLWNMHT6Nd8Sql4pVRnr3+HKaVebwvjNB2DxMREzj//fO68806+//3v43A4eOedd3j44Yd5//33KSjwX4WiIxOMhIwVBFPj0SqCqfNoBY328PX4WFgTTI8sPdB8pdSm1jOnZYSjjEt7IjIyknHjxjF27Fj27dtHbm4un3/+OZ999hkpKSlkZmYyaNCgkAk72iUhY1eNR7CvzqO4nCDif4ws7nTCkZUK1q0hrFmD2+1mwoQJ9OnTh3feeafBfC2sGTzBjJE9AaQDWzDCi6OBHKAKI9TYbpI+wlnGpT0hIvTv35/+/ftz8uTJ2rDjjh076N69e23Y0eXq2GMadknI2FXjEeyr8yjRTiRCUNXUdWYiSITgjDNC2FYrWLeGsGYNjz76KKmpqZw61fB8aWHNphFMobC9wHil1ASl1HhgHPCNUuqc9py5qGkfJCQkcN555/Gzn/2MH/zgB0RERPDuu+/y8MMPs3z5co4fP263ic0m7ew+OCN930JWSsg0VsOxJRXeG2N2ymxcTt8vIFbWeRQRIrrH4IiPPB1GdJjZit1jalPvg1Gwbgnp6emsWbOGCRMmcPfddzf7ZSw/P593333Xb2/OW1gzKiqqVlhT45tgemTDlVJf1vxQSm0VkbEW2qQJQSIjIxk7dixjxowhPz+fnJyc2tDj0KFDycrKYtCgQc0qwmoXdknIjD5nKuvfeh2PjwitQ8Ho7zUUcGwt5oyaw8p9K/3WeZwzao5lbYtDjBT7ANmJVipYt6aw5h133MGf//xnn8uCFtZsKsE4su0ishB4ESOUeC2w3VKrNCGLiNC3b1/69u3LqVOn2LhxIxs2bODFF1+ka9euZGZmMnbs2A4RdrRLQmbg0UK2VVRRGuVsWOOx0s3Ao9bpkbX3Oo+WKVh73Dz56F+4cPJoEkp2Q1kEk9LTOHSwbm2IYGRc3nnnHXr06MH48eNZvXq1z2W0sGbTCOZOmwP8GFhg/l4L/MMyizRhQ+fOnTnnnHM466yz2LZtG7m5ufznP//hww8/ZNy4cWRmZtK1a4CK5+0AOyRkSl5+lcknC/3WeCw5+ircYV1d75jIGG4bd5tlKtQtwRIFa48bjn1N3qYNXP/DS8xp1eTlbWLGeZOM+WaJrGB6ZOvWrWPZsmW89957lJeXc+rUKa699lpefPHF2uW1sGbTaNSRKaXKgb8CfxWRLkCyOU2jaRUiIiIYM2ZMnbDj+vXrycnJYciQIWRlZTF48OAOFXa0EndhIRFK+a3x6C6wrkfW3rFEwbr4CFRXkJTQmbytO7jonCm8+8HHnCoqZvK4EbXq1BBcj+z+++/n/vvvB2D16tU8+OCDdZwYaGHNphJM+v1qYKa57GbgqIisUUr5ljbVaFpAcnIyycnJXHDBBbVhx5deeokuXbrUhh2jo6PtNtNWnImJuAN8m+dMCq9SYd5YomBdegxQ3PXj65j941+x9K33GdivD28sfAiHQ2rVqVsDLazZPIIR1sxTSo0TkflAX6XU3SKyRSk1um1MbDq61mLoUF1dzfbt28nJySE/P5+oqCjGjBlDZmYm3bt3t9s8Wzjy2GMBixV3mT+PHj/9qQ2WWYev2oD+qP2OrLUUrA/mNb5M73FN366NhGOtxQgR6QVcAfyvxfZoNHWIiIggLS2NtLQ0Dhw4QG5uLps2bWL9+vUMHjyYzMxMhg4dGlZhx25z51KwYhW7/RQr7jY3vBWWHC4nCVP7Nzs7seEGI8ATIInEYU1SjyZ4gjkD9wLvA58opdaLyCBgp7VmaTQN6dOnD7NmzWLq1Km1zmzJkiUkJSWRkZHBuHHj6NQp9KXs7SxWHJbEdIPiwwSjTq2xh2CSPV4FXvX6vRu4zEqjNJpAxMXFkZ2dzZQpU9i+fTu5ubmsWLGCVatW1YYde/ToYbeZlmFXseKwJa4Hquwk7uoYPKoztaKecgpnRCkSF7rXWkdB94k1HRan08moUaMYNWoU3333HTk5OeTl5bFhwwYGDhxIVlYWw4YNC7mwYzDFiq10ZHZJyLg9bo6XH+dE+QncHjdOh5Mu0V3oGt0Vp4UK0QoH1aovSnkfcycelYRSXYnAgf7Cy160I9OEBL169eIHP/hBnbDj0qVLSUxMJCMjg/T09JAJO9pVrBjsk5DxKA97Tu2h0l1Z+7Gw2+PmWNkxTlWeYmDngZY5M3dxJcrtOylOuRXu4kot6mkzofWqqgl7YmNjOeuss1iwYAFXXHEFCQkJrFy5koceeohly5Zx+PBhu01sMdFxgT/qtapYMdgnIVNSVVLHidWglKLSXcnxcutqdnqKqxoX9dTYSjDfkbkwxsQGeC+vlLrXOrM0mpbhdDoZMWIEI0aM4NChQ+Tm5rJlyxY2bdrEgAEDyMzMJCUlBafTupCUVaSd3YdNK/b5DC9aWawY7JOQKakqIVbF+pynlOJE+Ql6xFg0VuVHtPMPD99HbEwcRcWnOOeSqXVqKTaHwsJCFi9ezK233tqi7YQjwYQW3wJOAhsB34p6Gk07pmfPnsycOZPzzz+/Nuz4yiuvkJCQUBt2jImxt0ZgU7CrWDFAmZ8it6fnWyMh41F+CgGbuD1uS9oFjGr7ARSo7/7F/xHVO66hTW53k16UCgsLeeKJJ7QjawbBOLJkpdRFllui0VhMTEwMZ555JpMmTeLrr78mJyeHDz74gNWrV5OWlkZWVhY9e/a028xGsatYMYArIoKKav+hNFeENW07JPAoiJXJHo64SDxFRnjxgb/9hRdfX0Jy72S6d+nKuNHjmH/Xj5k56/tcfvnlDBgwgLlz57JixQp+8pOfkJGRwW233cbRo0eJiYnh6aefZvjw4Rw+fJhbbrmF3bt3A/CPf/yDv/3tb+zatYuxY8cydepU/vKXv1i2T6FGMFfdpyKS5i3lotF0ZJxOJ6mpqaSmpnL48GFyc3P54osvyMvLo1+/fmRlZTF8+PB2HXa0o1gxQL9jJ9nVObpO1f0aHB6jcLEVxEbGIiIopdj88WYKvSv8C0Q6IvnM8Vmzt9+zZ08uvvhin/OccVGosmo2bsrjlWWvk7v8E6qrq8m6OJv0sek4IupeJ9HR0XzyyScAnHfeefzzn/9k6NCh5OTkcOutt/LRRx9x++23c/bZZ/Pmm2/idrspLi7mgQceYOvWrWzevLnZ+xGuBOPIzhu732MAABrlSURBVARuEJE9GKFFAVR7LlGl0QTLGWecwYwZMzj//PPJy8sjNzeXV199lc6dOzNhwgTGjx9PbKzvsZlwZGD+YQ4OSabU5UTJ6XCbKKFTpZuB+61JpomNjCXKGUWlu974nIAgRDqsS3ARhyHque6LHL5/8QxiOsWAQ5hxySU44iKpn3s/e7YhLlpcXMynn37KD3/4w9p5FWZZsY8++ojnn38eMF6sEhISKAhQP1MTmGAcme/XFI0mhOjUqROTJ09m4sSJfP311+Tm5vLRRx+xZs0a0tLSyMzM1DIagDOxO664K6lw7KO6aiuoMpBOOKNG4YrshzPpOUvadYiDgZ0Hcrz8OOPPHt+m35GB4cyc0RFEdHYRlWwodDuinD41wmpefDweD4mJibqH1QY0mn6vlPpWKfUtUIZRo6XmT6MJORwOB8OHD+e6667j1ltvZdy4cWzbto2nnnqKZ555hq1bt+J2W5hY0M45dM4tlMf0IiL2LKITf0x00p1EJ/6YiNizKI/pxaHv3WJZ206Hkx4xPRjeZTgju41keJfh9IjpYbkTqyE7O5s333yTsrIyioqKePvttwMu37lzZwYOHMirrxqFkZRSfPHFF4ARcvzHPwxZR7fbzalTp4iPj/erGK0JTKOOTERmishOYA+wBtgL/MdiuzQa2+nRoweXXHIJd955JxdeeCHFxcW89tprPPLII6xZs4bi4mK7TWxz9pT1xuP0XcnR44xiT1no9lrT09OZPXs2Y8eO5bLLLuOss85qdJ2XXnqJZ555hjFjxjBy5EjeeustAB599FFWrVpFWloa48ePZ9u2bXTt2pUpU6YwatQo7rrrLqt3J6QIRsblC+Bc4ANTzuUc4Cql1E1tYqBRpPh/gQSl1OXBrKNlXDRW4PF4+Oabb8jJyWHXrl04nU5GjhxJVlYWffpY9+1We+LxWz4KvIDAbf84t9XbbYqMixXYVR7LKsJRxqVKKXVcRBwi4lBKrRKR/xfMxkXkWeAS4IhSapTX9IuARzGqby5USj3gbxtmkeJ5IvJaMG1qNFbhcDgYNmwYw4YN49ixY+Tm5rJ582a2bNlCcnIymZmZjBgxggiLUtDbA9FxkQFLZFlZVcTj8VBaWEDpqZN43G4cTicxnROISUyytJ6m2+O2rTyWJjiCueMKRSQO+Bh4SUSOAAHEeerwL+DvwPM1E0TECTwOTAXygfUisgzDqd1fb/25SqkjQbal0bQZ3bp1Y9q0aZx77rls3ryZ3Nxc3njjDVasWFGb7RgfH2+3ma2OXVVFlFKcOLAfd1VVrTPxuN2UFBZQXlJMlz59LXNmx8uPN1oey7KqIpqgCMaRfR8j0eMO4BogAUOjrFGUUmtFZEC9yZnAN2ZPCxFZCnxfKXU/Ru+tWYjITcBNAP36WVfZQKPxJjo6mokTJ5KZmcmuXbvIzc1l9erVrF27tjbsmJycbLeZrYZdVUUqS0txR0f5dCbuqipKCwuI69LVkrZPlJ9o0K53+5aWx9IERTB6ZCUi0h8YqpR6TkRiMHpPzaUPsN/rdz6Q5W9hEekK/AkYJyK/Nh2eLzufAp4CY4ysBfZpNE3G4XAwdOhQhg4dyvHjx2vDjl9++SW9e/cmKyuLkSNHdviwo11VRSrKSvF4OvtMd1dKUXrqpGWOrLHyV5aWx7KAxvIiOiLBFA2+EaOn0wUYjOGI/gmc18w2fUn3+D2ySqnjgHU5vRpNK9O1a1cuvvhizj33XL744gtyc3N58803WbFiBePHj2fChAl07tzZbjObjR1VRYoOf0dJUhKx0S6fzsxj4ScRToczoLPqSONjSimOHz9OdHS03aa0KsG8Pt2GEQ7MAVBK7RSRlvSj84G+Xr+TgYMt2J5G0y5xuVxkZmaSkZHB7t27ycnJYe3atXzyySekpqaSlZVF3759fT6YNXXZvWYFAHHdz/B5vMThoKAy2KH7plFUWURxZTHKx/u2IMRFxbH98HZL2raC6OjokAp3Q3COrEIpVVlz8YhIBC37IHo9MFREBgIHgCuBq1uwvVpEZAYwY8iQIa2xOY2mVRARBg8ezODBgzlx4gTr169n06ZNbNu2jV69epGZmcmoUaOIjLQu46+jM2Limax/5QU8Pny+Q0HG9y8j9czGv+tqDqVVpVz93tXkF+VT4T4tAOJyukiOT2bxtMXERHYc9YRQJJg0nzUi8hugk4hMBV4FAn/SbiIiS4DPgBQRyReReUqpauAnwPvAduAVpdS25plfF6XU20qpmxISElpjcxpNq9OlSxcuvPBC7rzzTqZPn051dTVvvfUWf/3rX/nwww85edKaorsdnYFHC4mpqMLhqZst6fB4iKmoYqB3EeFWJiYyhsXTFjNn5BySXEkIQpIriTkj52gn1k4I5oNoBzAPuABjfOt9jG+/2u2Iof4gWtNRUEqxZ88ecnJy2LFjByJSG3bs16+fDjuafD1pMhUnC9ndPZF93RKodDqIchvV9gcdLcSVmMiwTz+128wOTUh/EK2U8gBPm38ajaYVEREGDRrEoEGDKCgoqA07fvXVV5xxxhlkZWWRlpYW9mFHd2EhOKKIiJ6EKyEbR2QckVXFRBSvBfkQd4F1PTJN+yeYHtklwB+A/hiOr0bGpd2mXekemaYjU1lZyZdffklOTg5HjhyhU6dOpKenk5GRQWJiot3m2cJXU75H7sD5lHXqVqfWo8NdSaeyY2TuXciIT1bbZ2AI0JF7ZME4sm+AS4Ev23M4Eeoke/z/9u49OKrzvOP499ENxFWAwHgsBBISCAnrAuhs2qau49hM0jbETTpxnKStM65ncu3ESdrUTSbTaZombuJ0JjOeepwmdVu7cdLETSF2artO3KQOObsSQiAQGNvcZAwOxiBxFdp9+8cuaxmEtELsnj27v88MY+ns2bOP/CJ+eh+dfd+79uzZE3Q5IlPinGPfvn1Eo1F27doFQFNTE57nsWzZsqJqO/70C4+y+0jVmAsWl8SHWXnNcW762/cHUFnhCHOQZXLX4kGgL99DDJI3ewCb1q1bd1fQtYhMlZlRV1dHXV0dx48fT7cd+/v7WbRoUbrtWFEx9mr0hSS56v7Yt9cX+qr7MrFMguwvgCfM7H9J7hANgHPuG1mrSkTepKqqiltuuYUbb7wx3XbctGkTTz/9dLrtOG/evKDLzJqzp8Z/j9jZ09l5D9kFw2fPENv4GL1PPc6Zk0NUzppN2/rfo3PDe6iYXpnV15aJZRJkXwZOAtOBwv/RTySPlZeXs2bNGjo6Ojhw4AC+77N582Y2b97MihUriEQi1NXVFVzbMchV94fPnuHfP/8Zjh85TPz8MABnhgaJbfwhe/zn+MCX71OYBSyTIJvvnFuf9UpEJGNmxtKlS1m6dCknTpygq6uL7u5udu/ezcKFC/E8j7a2toJpOwa16j5AbONjHD/yCvHzbw7S+Plhjh95hdjGx/it930wa68vE8vkDdH/Y2ahCDIze5eZPag3lUoxmTt3Lm9/+9u5++67ufXWWykrK+Pxxx/nvvvu48knn+TYsWNBlzhl7bfUMre6ktLyN/+Tle1V9wF6n/zxJSF2Qfz8eXqfejxrry2ZyeSuxSFgJsnfj51Ht9+L5DXnHAcPHiQajbJz504SiQQrVqzA8zyWL18e2rbj8NmRnK+6D3DfbRPvLvWZ7/04a6+fKwV916JzrvB2BxQpYGZGbW0ttbW1DA4OptuODz/8MNXV1em247Rp04IudVKCWHUfoCKeYLj08s2rivil7U7JrXBvjiQi45ozZw433XQTN9xwAzt27MD3fZ544gmeeeYZ2tvb8TyPBQuys49Xoag9mlwaKzHGDtQliQS1R7WqSNAUZCJFoKysjLa2Ntra2hgYGMD3fWKxGL7v09DQQCQSYfny5ZSM8Y91sWschsPnznN6WvmbwuzCgsWNl7+ZUnKkoIJM27iITKympoaamhrWr19Pd3c3XV1dPPLII8yfPx/P82hvby+4jRenovr224n888NsrVnJscrTOM5iTGfeuRm0D+ym+sMfCrrEopfJzR7/5pz7o4mO5RPd7CGSuZGREXbu3Ek0GmVgYICKiop027G6ujro8gJ37tggj/75f3O6ZM4l6zzOSAzy/q+9g2nz8/bet4wV9M0eQMvoT8ysFFibnXJEJNfKyspobW2ltbWVl19+Gd/36e7uJhqNsnz5ciKRCA0NDUXbdux97ihnKheSGHnzD/2J0grOTFtI73NH8d4V/iALs8vOyMzsHuCvgErgNMnb7gGGgQedc/fkpMIroBmZyNScPHky3XYcGhpi3rx5eJ5HR0dH0bUdv/3ZX4y/qsiscu78enZ2p86lMM/IMmktfiWfQ2ssCjKRqyMej9Pf34/v+xw8eJDy8nLa2trwPI9FixYFXV5O3P+Rn45/gsHH//Gm3BSTRWEOskzeR3aPmW0AbkgdetY5F/53/4nIhEpLS1m9ejWrV6/m0KFDRKNRenp66Orqor6+Hs/zWLFiRUG3HYNc51Eyk9GMDPCAR1KHbge68nmWphmZSPacOnUq3XYcHBykqqoq3XasrCy8xXOjm14ad53HNetrc/4m7WwI84wskyDbBrQ75xKpz0uBHudcaw7qmxRtrCmSO/F4nF27duH7PgcOHKC8vJzW1lY8z+Oaa64JuryrZvjsCD+8t5sTR8+8KcwurPP43s+tzeoSWblSDEF2o3PuWOrz+STbi3kXZBdoRiaSW4cPH8b3fbZv387IyAjLli0jEomwcuXKgmg7BrXOYy4VepDdDnwV+BnJOxdvAO5xzj2a/fKujIJMJBinT59my5YtxGIxTpw4wdy5c+ns7GTNmjXMmDEj6PJkHAUdZABmdi3QSTLIfOfc4WwXNhUKMpFgxeNxnn/+eXzfZ9++fZSVlXH99dcTiURYvHhx0OXJGMIcZJnOiTt5467FBLApO+WISCEoLS1l1apVrFq1iiNHjhCNRunt7aWnp4elS5fieR5NTU2UlpYGXaoUgExai18lGWS6a1FErtjp06fp6ekhFotx/Phx5syZk247zpw5M+jyil6YZ2QFddfiBQoykfyVSCTSbce9e/dSWlqabjtee+21QZdXtMIcZJm2FquAC/ulz81SLSJSBEpKSmhqaqKpqYlXX3013XbcunUrS5YsIRKJsGrVKrUdJWO6a1FEAnfmzBm2bt1KNBrl9ddfZ/bs2axbt461a9cya9asoMsrCmGekY0bZGZmQA0wQgjuWtQbokXCLZFIsGfPHqLRKC+++GJ6iSzP87juuuuCLq+gFWyQQfqLC9W2LZqRiYTfr3/963TbcXh4mJqamnTbsaysMN6EnE8KPcjuBx5yzsVyU9LUKchECsfZs2fTbcdjx44xa9asdNtx9uzZQZdXMAo9yHYCK4D9wCmS7UWnuxZFJJcSiQQvvvgivu/zwgsvUFJSQktLC5FIhJqamqDLC70wB1km8/N3Zr0KEZEJlJSU0NjYSGNjI0ePHiUWi9HT08P27du57rrr8DyPlpYWtR2LUEZLVIWNZmQixeHcuXP09vbi+z6vvfYaM2fOTLcd58yZE3R5oRLmGZmCTERCL5FI8NJLL+H7Pnv27KGkpITm5uZ02zF5A7aMJ8xBpjm4iIReSUkJDQ0NNDQ0cOzYsfRO1n19fVx77bVEIhFaWlooL9duzoVIMzIRKUjnzp1j27Zt+L7P0aNHmTFjBmvXrqWzs1NtxzGEeUamIBORguacY+/evfi+z+7duzEzmpub8TyP2tpatR1TwhxkBdVaHLWyR9CliEieMDPq6+upr6/n9ddfJxaLsWXLFnbs2MHixYuJRCKsXr1abccQ04xMRIrO8PAw27ZtIxqN8uqrr1JZWcnatWtZt24dVVVVQZcXiDDPyBRkIlK0nHPs27cv3XYEaGpqIhKJsHTp0qJqO4Y5yAqqtSgiMhlmRl1dHXV1dRw/fjzdduzv7+eaa67B8zyuv/56Kioqgi5VxqEZmYjIKMPDw/T19eH7PkeOHGH69OmsWbOGzs5O5s2bF3R5WRPmGZmCTERkDM459u/fTzQapb+/H4CVK1fieR51dXUF13YMc5CptSgiMgYzY9myZSxbtowTJ04Qi8Xo7u5m165dLFy4kEgkQmtrq9qOeUAzMhGRDJ0/fz7ddjx8+DDTp0+no6ODzs5O5s+fH3R5UxLmGZmCTERkkpxzHDx4EN/32blzJ845VqxYQSQSob6+PpRtxzAHmVqLIiKTZGbU1tZSW1vL4OAgXV1ddHV18fzzz1NdXY3nebS1tTFt2rSgSy0KmpGJiFwF58+fZ8eOHUSjUQ4dOsS0adPSbccFCxYEXd6EwjwjU5CJiFxFzjkGBgbSbcdEIkFjY2O67VhSUhJ0iWMKc5CptSgichWZGUuWLGHJkiUMDQ2l244PP/wwCxYsSLcdp0+fHnSpBUMzMhGRLBsZGWHnzp34vs/LL79MRUUF7e3teJ5HdXV10OUB4Z6RKchERHJoYGCAaDRKX18fiUSChoYGPM+joaEh0LajgizPKMhEJN+dPHmS7u5uYrEYJ0+eZP78+XieR3t7eyBtRwVZnhi1H9lde/bsCbocEZEJjYyM0N/fTzQa5eDBg5SXl6fbjgsXLsxZHQqyPKMZmYiE0aFDh/B9n76+PuLxOPX19UQiERobG7PedlSQ5RkFmYiE2cmTJ9myZQuxWIyhoSHmzZtHZ2cnHR0dVFZWZuU1FWR5RkEmIoUgHo+za9cufN/nwIEDlJeX09raSiQSYdGiRVf1tcIcZHofmYhIniotLaWlpYWWlhZeeeUVotEovb29dHd3U1dXh+d5rFy5Mm/fZJ0rmpGJiITIqVOn0m3HwcFBqqqq0m3HGTNmXPF1wzwjU5CJiIRQPB5n9+7d+L7P/v37KSsro7W1lZtvvvmKAi3MQabWoohICJWWltLc3ExzczOHDx8mGo2yd+/eotzoU0EmIhJyixcvZsOGDcTjcUpLS4MuJ+eK+zeEIiIFpBhDDBRkIiIScgoyEREJNQWZiIiEmoJMRERCTUEmIiKhpiATEZFQU5CJiEioFdwbolObax41s/0XPTQXOJHBsWrgaJbKG89YteTqOpk+Z6Lzxnv8co9lMi5BjclYteTqOvk+JqDvlamcN9lxyXSspjImS6/wecFzzhXUH+DBTI9f5lhXPtWdi+tk+pyJzhvv8amMS1BjEuS45PuYBDkuxfi9kulYBfm9EuSfQmwtbprE8cudG4SrVcuVXCfT50x03niPa1yyc77GJLfXCWpcJjNWRacgV7+fCjPrciFdAbpQaUzyk8Yl/xTrmBTijGyqHgy6ALmExiQ/aVzyT1GOiWZkIiISapqRiYhIqCnIREQk1BRkIiISagqyCZhZvZl928x+EHQtkmRmt5rZt8zsv8xsfdD1CJjZKjN7wMx+YGYfDboeeYOZzTSzbjP7/aBryZaiDDIz+46ZvWpmfRcdf4eZ7TazF8zsLwGccy855+4MptLiMckx+ZFz7i7gDuC2AMotCpMck37n3EeA9wFFd/t3Lk1mXFI+B3w/t1XmVlEGGfAQ8I7RB8ysFLgfeCfQDNxuZs25L61oPcTkx+QLqcclOx5iEmNiZhuA/wOeyW2ZRechMhwXM7sZ2AkcyXWRuVSUQeac+zlw7KLDHvBCagY2DDwKvDvnxRWpyYyJJd0L/MQ5tyXXtRaLyX6fOOc2Oud+E/hgbistLpMcl7cBbwE+ANxlZgX5b37BLRo8BdcBB0d9PgBEzGwB8GWgw8zucc59JZDqitOYYwJ8ErgZmGtmDc65B4Iorkhd7vvkRuA9wDTgiQDqKnZjjotz7hMAZnYHcNQ5lwigtqxTkL3BxjjmnHOvAR/JdTECXH5Mvgl8M9fFCHD5MXkWeDa3pcgoY45L+gPnHspdKblXkNPMKzQALBn1eQ1wKKBaJEljkn80JvmpqMdFQfaGGNBoZnVmVgG8H9gYcE3FTmOSfzQm+amox6Uog8zMvgtsBlaa2YCZ3emcGwE+ATwJ9APfd87tCLLOYqIxyT8ak/ykcbmUFg0WEZFQK8oZmYiIFA4FmYiIhJqCTEREQk1BJiIioaYgExGRUFOQiYhIqCnIRDJkZs+aWda3KDGzPzOzfjN75KLj7Wb2u+M8b52Zjbt0l5ndaGY/vlq1iuQDrbUokgNmVpZ602omPga80zm396Lj7ST3+rpkUd7U9buArqlVKhI+mpFJQTGzZanZzLfMbIeZPWVmlanH0jMqM6s2s32pj+8wsx+Z2SYz22tmnzCzT5tZj5n9yszmj3qJD5nZL82sz8y81PNnpjY7jKWe8+5R1/0PM9sEPDVGrZ9OXafPzD6VOvYAUA9sNLO7R51bAfwNcJuZbTWz28zsr83sQTN7CvjX0bMtM/NSdfak/rtyjNf/ndS1tqbOmz31ERDJPQWZFKJG4H7nXAtwHHhvBs9ZTXLPJo/ktj2nnXMdJJcC+uNR581M7bn1MeA7qWOfB37qnOskuf/T18xsZuqx3wD+xDl30+gXM7O1wIdJbkvzFpJ7RXWkdlk+BLzNOfcPF85P7TH1ReB7zrl259z3Ug+tBd7tnPvARV/PLuCG1NfwReDvxviaPwt83DnXDvw2cGbC/0sieUitRSlEe51zW1MfdwPLMnjOz5xzQ8CQmZ0ANqWObwdaR533XUhubmhmc8ysClgPbDCzz6bOmQ7Upj5+2jl38SaIAG8F/tM5dwrAzB4jGSY9mXyBo2x0zo0VQHOBfzGzRpLbeZSPcc5zwDdSv4t7zDk3MMnXFskLmpFJITo36uM4b/zANsIbf+enj/OcxKjPE7z5B76LFyd1JPeCem9qptTunKt1zvWnHj91mRrH2j/qSlzu+l8iGc6rgXdx6deLc+6rwJ8ClcCvzKzpKtUkklMKMikm+0i24gD+8AqvcRuAmb0VOOGcO0FyxfFPmpmlHuvI4Do/B241sxmpNuQfAL+Y4DlDQKa/x5oLvJz6+I6xTjCz5c657c65e0neJKIgk1BSkEkx+TrwUTP7JVB9hdd4PfX8B4A7U8e+RLJ1t83M+lKfj8s5twV4CIgCPvBPzrmJ2oo/A5ov3Owxwbl/D3zFzJ4DSi9zzqdSN5r0kvz92E8mqlskH2kbFxERCTXNyEREJNQUZCIiEmoKMhERCTUFmYiIhJqCTEREQk1BJiIioaYgExGRUFOQiYhIqP0/boYzP7FDPlQAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "%pylab inline\n", "import random, math, pylab\n", "\n", "def markov_pi(N, delta):\n", " x, y = 1.0, 1.0\n", " n_hits = 0\n", " for i in range(N):\n", " del_x, del_y = random.uniform(-delta, delta), random.uniform(-delta, delta)\n", " if abs(x + del_x) < 1.0 and abs(y + del_y) < 1.0:\n", " x, y = x + del_x, y + del_y\n", " if x**2 + y**2 < 1.0: n_hits += 1\n", " return n_hits\n", "\n", "n_runs = 500\n", "for delta in [0.062, 0.125, 0.25, 0.5, 1.0, 2.0, 4.0]:\n", " n_trials_list = []\n", " sigmas = []\n", " for poweroftwo in range(4, 13):\n", " n_trials = 2 ** poweroftwo\n", " sigma = 0.0\n", " for run in range(n_runs):\n", " pi_est = 4.0 * markov_pi(n_trials, delta) / float(n_trials)\n", " sigma += (pi_est - math.pi) ** 2\n", " sigmas.append(math.sqrt(sigma/(n_runs)))\n", " n_trials_list.append(n_trials)\n", " pylab.plot(n_trials_list, sigmas, 'o', ms = 8, label = '$\\delta = $' + str(delta))\n", "\n", "pylab.xscale('log')\n", "pylab.yscale('log')\n", "pylab.xlabel('number of trials')\n", "pylab.ylabel('root mean square deviation')\n", "pylab.plot([10,10000],[1.642 / math.sqrt(10.0), 1.642 / math.sqrt(10000.0)], label = 'direct')\n", "pylab.title('Markov-chain sampling of pi: root mean square deviation vs. n_trials')\n", "pylab.legend(loc='upper right')\n", "pylab.savefig('markov_sampling_rms_deviation.png')\n", "pylab.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Bunching method for Markov-chain sampling. By this method, Markov-chain data are analyzed as direct-sampling ones and the error is underestimated at the beginning since the data is correlated, but that at later stages of the bunching they become more and more uncorrelated, just as direct-sampling data and the saturation (plateau) comes from the fact that the data are almost unrelated. For uncorrelated data, one has a plateau." ] }, { "cell_type": "code", "execution_count": 23, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Populating the interactive namespace from numpy and matplotlib\n", "3.15369415283 mean value, estimate of pi\n", "0.0121014992422 mean value, estimate of pi\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZIAAAEWCAYAAABMoxE0AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3XucHFWZ//HPl9wcFIiE6JoETFhiNMBKYIwo6rKCJLhKIgtrUHdRcVEXVlnXaNCfLIu4wOKuNy4aBUVEEkTIjhINrgERxZCBADFAZIRgZgIYhIRbgCQ+vz/qdOh0umdq0lPT0zPf9+vVr6k+darqqb7M01Wn6hxFBGZmZjtrl0YHYGZmzc2JxMzM6uJEYmZmdXEiMTOzujiRmJlZXZxIzMysLk4kQ5ykNZKOrDHvzZJW93dMvSXp65I+1+g4BjJJT0nat4Hbb4rPUnckTZQUkoY3OpaBRr6PZGCQtAZ4ObAV2Az8GvhIRKzth+1+KCL+r8jt2MAh6TtAZ0T8vwK3EcDkiOgoahv9TdJE4AFgRERsaWw0A4uPSAaWd0bES4BXAI8AX2twPJZTtV+pvf3l2iy/dJslzoFssL2GTiQDUEQ8C1wNTC2VSbpR0ofKnr9f0s1lz0PSRyTdJ+lxSRdKUtn8f5J0j6QnJd0t6eCyTR4k6S5JGyUtlPSitMzhkjrL1rFG0ier1U3zPyXpIUnrJH0oxbRfnn1O+/d5Sb9KMV4vaa+y+T+Q9HDa7k2S9i+b9x1JZ6fpeyS9o2zecEmPlvZX0qGSfi1pg6Q7JR3eTUzjJP1Q0npJD0j6WNm8MyVdLel7kp4A3l+jbJSkL6fXZF2aHlX++kr6tKSHgW9XbH9UivOAsrKxkjZJepmkvST9ONV5TNIvJVX9TpfeC0knA+8FPpVOd/1oJ/d1uqRb0rYfknSBpJGp/k1p0TvTNt5d5bP0mvSeb5C0StIxFe/nhZKuS5+FZZL+ssZ+lU43nSjpD+m9/mzFus4ue17tMz03faaflnSJpJdL+kna9v9JemnFZj+Y3suHJP1b2bp2kTRP0u8l/UnSVZL2rIjzJEl/AJZW25+mFRF+DIAHsAY4Mk3vClwGfLds/o1kp6BKz98P3Fz2PIAfA6OBfYD1wMw073igC3gdIGA/4JVl270VGAfsCdxDdkoN4HCyUyDkqDsTeBjYP8V/eYppvzT/PcBd3ez/jcDvgVcBLen5uWXzPwjsBowCvgzcUTbvO8DZafoM4IqyeX8L3JumxwN/At5O9iPqben52Crx7ALcltY3EtgXuB+YkeafSXYKcnaq21Kj7CzgN8DLgLFkpyw/X/b6bgHOS/vVUiWOS4EvlD0/Bfhpmj4H+DowIj3eTDpdXWU95e/Ftterjn09BDgUGA5MTJ+F06ptr/KzlGLtAD6TtvdW4ElgSll8jwHT0/qvABbU2K+JaVvfTHG9FngOeE2Nfd0WR9ln+jdkp5XHA38EbgempfdkKfDvFdu6EngxcCDZ96z0vT0trWtCWvYbwJUVy343LbvDe93MDx+RDCyLJG0AniD7J3d+L5c/NyI2RMQfgBuAg1L5h4D/iojlkemIiAfLlvtqRKyLiMeAH5UtV02tun8PfDsiVkXEM8B/lC8UEd+PiL/qIf5vR8TvImITcFV5HBFxaUQ8GRHPkf1je62kPaqs4/vAMZJ2Tc/fk8oA3gcsjojFEfHniPgZ0E6WWCq9jizBnBURz0fE/WT/rOaU1bklIhaldW2qUfZe4KyI+GNErE+vyz+UrePPZP+onitbR+X+nFD2vHx/NpOdBn1lRGyOiF9G+q/VS73e14i4LSJ+ExFbImIN2T/Nv865vUOBl5B9Xp+PiKVkP4LK9/OaiLg1sraIK+j+MwnwHymuO4E7yRJKXl+LiEciogv4JbAsIlakz9q1ZEmlcltPR8RKsqPIUtwfBj4bEZ1ln9PjtP1prDPTstXe66blRDKwzI6I0WS/Zk4FfiHpL3qx/MNl08+QfVkB9ib7td/b5XpTdxxQfmHAzlwkUHXdkoZJOjedMniC7FckwF4VyxNZ4+49wDtTMjmGF/7xvhI4Pp1O2ZCS9pvI/hlXeiUwrqLuZ8h+uXa3j5Vl44DypP1gKitZH9mpzFqWAi2SXi/plWT/UK9N884n+2V/vaT7Jc3rZj3d6fW+SnpVOq32cHpP/pMq70cN44C1EfHnsrIHyY4ISnrzmdyZ+uUeKZveVOV55brKX4vy9/OVwLVlr+E9ZBfP9PSZaXpOJANQRGyNiGvIPoRvSsVPk50yKulNglkLVD3H3IceIjukL9m7D9f9HmAWcCSwB9lpAshO01VzJdmvxFnA3fHClUNrgcsjYnTZ48URcW6VdawFHqiou1tElB+9VPv1X1m2juwfTMk+qay7dbwwM/tne1Xan/cAP46IJ9O8JyPi3yJiX+CdwCckHdHd+mpsc2f29WLgXrIrs3YnSzy13o9K64C9K9pz9iE7/drX6vne1FL+2S5/P9cCR1e8ji9KRzolg/IyWSeSAUiZWcBLyX7VANwBHCtpV2UN2Cf1YpXfAj4p6ZC07v3Sr9u+dBXwgdSIuivZ+fa+shvZee8/kf1T+M8e6i8AjgI+ygtHIwDfIztSmZGOcl6UGl8nVFnHrcATyhrCW1L9AyS9rpexXwn8P2WN5HuRvS7f6+U6vg+8m+w02bb9kfSO9F6K7HTo1vToySNk7SAlO7Ovu6VtPiXp1WSvdXfbKLeM7B/8pySNUHbBwzvJ3re+dgfwdkl7pqP70/pgnZ9L38P9gQ8AC1P514EvlL5b6T2f1QfbG/CcSAaWH0l6iuwL+gXgxIhYleZ9CXie7At6Gdl541wi4gdpfd8na9RcRNZY3mci4ifAV8naZjqAW9Ks5wAkvVfSqhqL9+S7ZKcQuoC7yRo0u4vlobT9N/LCl5zI7smZRfbreT3ZL8i5VPkeRMRWsn9uB5HdO/AoWUKu1i7TnbPJ2mHuAlaSNeSe3e0SO8ZS+sc7DvhJ2azJwP8BT5Ht70URcWOOVV4CTE2nYBbt5L5+kuwI6Umy9pSFFfPPBC5L2/j7iv15nuyU49FpWxcB/xgR9+aIvbcuJ2szWQNcXyXOnfELss/4z4EvRsT1qfwrQBvZqcYnyT6nr++D7Q14viHRCiHpNcBvgVHhm7fMBjUfkVifkfQuSSPTdffnAT9yEjEb/JxIrC99mOyU0e/JztVXnjc3s0HIp7bMzKwuPiIxM7O6DKqOw2rZa6+9YuLEiY0Ow8ysqdx2222PRsTYnuoNiUQyceJE2tvbGx2GmVlTkfRgz7V8asvMzOpUaCKRNFPSakkd1foBUtZN9sI0f5mygWOQNEbSDcq6oL6gYpkTJK1U1u3zT1XW1biZmfW/whKJpGHAhWR3r04FTpA0taLaScDjEbEf2Z3b56XyZ4HPkd09W77O4WR3j/5N6kn2LrLODc3MrEGKPCKZDnRExP2pS4QFZN1TlJtF1t0HZAM5HSFJqZvlm8kSSjmlx4tT/0K7s30HeGZm1s+KTCTj2b7L5E627yZ6uzrpDuiNwJhaK4yIzWQ3ua0kSyBTyfoN2oGkkyW1S2pfv379zu6DmZn1oMhEUq1L6cq7H/PUeaGyNIIskUwj68DuLuD0anUjYn5EtEZE69ixPV69Zma2zaIVXRx27lImzbuOw85dyqIVRfRwP3gUmUg62b7f/gnseBpqW53U/rEH2RCbtRwEEBG/TyPBXUXWw6uZWZ9YtKKL069ZSdeGTQTQtWETp1+z0smkG0UmkuXAZEmTJI0kG7azraJOG3Bimj4OWNrDUKFdZN1flw4x3sYL43WYmdXt/CWr2bR5+2FdNm3eyvlLVjcoooGvsBsSI2KLpFOBJcAw4NKIWCXpLKA9ItrI2jcul9RBdiSybYxoSWvIGtNHSpoNHBURd0v6D+AmSZvJxqh4f1H7YGZDz7oN1YdTr1VuBd/ZHhGLgcUVZWeUTT8LHF9j2Yk1yr9ONhKZmVmfGze6ha4qSWPc6JYGRNMcfGe7mVmZuTOm0DJi2HZlLSOGMXfGlAZFNPANib62zMzymj0tu0vh/CWrWbdhE+NGtzB3xpRt5bYjJxIzswqzp4134ugFn9oyM7O6OJGYmVldnEjMzKwuTiRmZlYXJxIzM6uLE4mZmdXFicTMzOriRGJmZnVxIjEzs7o4kZiZWV2cSMzMrC5OJGZmVpdCE4mkmZJWS+qQNK/K/FGSFqb5yyRNTOVjJN0g6SlJF1QsM1LSfEm/k3SvpL8rch/MzKx7hfX+K2kYcCHZcLidwHJJbRFxd1m1k4DHI2I/SXOA84B3A88CnwMOSI9ynwX+GBGvkrQLsGdR+2BmZj0r8ohkOtAREfdHxPPAAmBWRZ1ZwGVp+mrgCEmKiKcj4mayhFLpg8A5ABHx54h4tJjwzcwsjyITyXhgbdnzzlRWtU5EbAE2AmNqrVDS6DT5eUm3S/qBpJfXqHuypHZJ7evXr9/ZfTAzsx4UmUhUpSx2ok654cAE4FcRcTBwC/DFahUjYn5EtEZE69ixY/PEa2aDxKIVXRx27lImzbuOw85dyqIVXY0OaVArMpF0AnuXPZ8ArKtVR9JwYA/gsW7W+SfgGeDa9PwHwMF9EayZDQ6LVnRx+jUr6dqwiQC6Nmzi9GtWOpkUqMhEshyYLGmSpJHAHKCtok4bcGKaPg5YGhE1j0jSvB8Bh6eiI4C7a9U3s6Hn/CWr2bR563ZlmzZv5fwlqxsU0eBX2FVbEbFF0qnAEmAYcGlErJJ0FtAeEW3AJcDlkjrIjkTmlJaXtAbYHRgpaTZwVLri69NpmS8D64EPFLUPZtZ81m3Y1Ktyq19hiQQgIhYDiyvKziibfhY4vsayE2uUPwi8pe+iNLPBZNzoFrqqJI1xo1saEM3Q4DvbzWxQmTtjCi0jhm1X1jJiGHNnTGlQRINfoUckZmb9bfa07C6D85esZt2GTYwb3cLcGVO2lVvfcyIxs0Fn9rTxThz9yKe2zMysLk4kZmZWFycSMzOrixOJmZnVxYnEzMzq4kRiZmZ1cSIxM7O6+D4SM+tzi1Z0+YbAIcSJxMz6VKkb91IPvKVu3AEnk0HKp7bMrE+5G/ehx4nEzPqUu3EfepxIzKxP1equ3d24D15OJGa2g3rGPHc37kNPoYlE0kxJqyV1SJpXZf4oSQvT/GWSJqbyMZJukPSUpAtqrLtN0m+LjN9sKKp3zPPZ08ZzzrEHMn50CwLGj27hnGMPdEP7IFbYVVuShgEXAm8DOoHlktrScLklJwGPR8R+kuYA5wHvBp4FPgcckB6V6z4WeKqo2M2aXT2X33bXWJ53He7GfWgp8ohkOtAREfdHxPPAAmBWRZ1ZwGVp+mrgCEmKiKcj4mayhLIdSS8BPgGcXVzoZs2r3iMKN5ZbbxWZSMYDa8ued6ayqnUiYguwERjTw3o/D/w38Ex3lSSdLKldUvv69et7E7dZU6v38ls3lltvFZlIVKUsdqLOC5Wlg4D9IuLanjYeEfMjojUiWseOHdtTdbNBo94jioHQWF5PY7/1vyLvbO8E9i57PgFYV6NOp6ThwB7AY92s8w3AIZLWkMX+Mkk3RsThfRW0WbMbN7qFripJI+8RRaPHPPed8c2nyESyHJgsaRLQBcwB3lNRpw04EbgFOA5YGhE1j0gi4mLgYoB0hdePnUTMtjd3xpTt/hFD748oGtlY3heN/da/CkskEbFF0qnAEmAYcGlErJJ0FtAeEW3AJcDlkjrIjkTmlJZPRx27AyMlzQaOqrjiy2zQqueqq0YfUdTLjf3Np9BOGyNiMbC4ouyMsulngeNrLDuxh3WvocqlwWbNri9O7TTz5bf1npqz/uc7280KUE9j8VDv9HAgNPZb77gbebM+Vu8RxVA/tdPsp+aGIicSsz5Wb2OxT+0096m5ocintsz62GC4j8OsN5xIzPpYvXeGu9NDazY+tWXWx5r9Pg6z3nIiMetjbiy2ocaJxKwAPqKwocRtJGZmVhcnEjMzq4sTiZmZ1cWJxMzM6uJEYmZmden2qi1Jw4AlEXFkP8VjNiDU04272VDTbSKJiK2SnpG0R0Rs7K+gzBrJI/SZ9U6e+0ieBVZK+hnwdKkwIj5WWFRmDeQR+sx6J08byXXA54CbgNvKHj2SNFPSakkdkuZVmT9K0sI0f1kaPhdJYyTdIOkpSReU1d9V0nWS7pW0StK5eeIw642h3o27WW/1eEQSEZdJGgm8KhWtjojNPS2X2lcuBN4GdALLJbVVDJd7EvB4ROwnaQ5wHvBusqOgz5GNgFg5CuIXI+KGFNPPJR0dET/pKR6zvNyNu1nv9HhEIulw4D6ypHAR8DtJb8mx7ulAR0TcHxHPAwuAWRV1ZgGXpemrgSMkKSKejoibyRLKNhHxTETckKafB24HJuSIxSw3d+Nu1jt52kj+GzgqIlYDSHoVcCVwSA/LjQfWlj3vBF5fq05EbJG0ERgDPNpTUJJGA+8EvlJj/snAyQD77LNPT6sz28adLpr1Tp5EMqKURAAi4neSRuRYTlXKYifq7LhiaThZMvtqRNxfrU5EzAfmA7S2tva4TrNy7nTRLL88iaRd0iXA5en5e8nX2N4J7F32fAKwrkadzpQc9gAey7Hu+cB9EfHlHHXNzKxAea7a+iiwCvgY8HHgbuAjOZZbDkyWNCk1jM8B2irqtAEnpunjgKUR0e3Rg6SzyRLOaTliMDOzguW5s/2SiHgf8D+9WXFq8zgVWAIMAy6NiFWSzgLaI6INuAS4XFIH2ZHInLJtrwF2B0ZKmg0cBTwBfBa4F7hdEsAFEfGt3sRmZmZ9J8+d7WMljUxXSfVKRCwGFleUnVE2/SxwfI1lJ9ZYbbV2FTMza5A8bSRrgF9JamP7O9t7dYRiZmaDU55Esi49dgF2KzYcMzNrNnnaSF4SEXP7KR4zM2sy3V61FRFbgYP7KRYzM2tCeU5t3ZHaR37A9m0k1xQWlZmZNY08iWRP4E/AW8vKAnAiMTOzXL3/fqA/AjEzs+bUYyJJnTReDLw8Ig6Q9FfAMRFxduHRme0kD5VrjTTUPn95ukj5JnA6sBkgIu6i7A50s4GmNFRu14ZNBC8MlbtoRVejQ7MhYCh+/vIkkl0j4taKsi1FBGPWF7obKtesaEPx85cnkTwq6S9J3btLOg54qNCozOrgoXKtkYbi5y9PIjkF+AbwakldZL3u5un916whag2J66FyrT8Mxc9fj4kkDZV7JDAWeHVEvCkiHiw+NLOd46FyrZGG4ucvz30kAETE0z3XMms8D5VrjTQUP3/qYRypQaG1tTXa29sbHYaZWVORdFtEtPZUr8dTW5JG5SmrsexMSasldUiaV209kham+cskTUzlYyTdIOkpSRdULHOIpJVpma8qjW5lZmaNkaex/ZacZdtJPQdfCBwNTAVOkDS1otpJwOMRsR/wJeC8VP4s8Dngk1VWfTFwMjA5PWbm2AczMytIzUQi6S8kHQK0SJom6eD0OBzYNce6pwMdqbH+eWABMKuizizgsjR9NXCEJEXE0xFxM1lCKY/pFcDuEXFLGtv9u8DsHLGYmVlBumtsnwG8H5jA9uO1Pwl8Jse6xwNry553Aq+vVSeN8b4RGAM82s06OyvWWbUFS9LJZEcu7LPPPjnCNTOznVEzkUTEZcBlkv4uIn64E+uu1nZR2bKfp85O1Y+I+cB8yBrbu1mnmZnVIc/lvz+W9B5gYnn9iDirh+U6gb3Lnk8gG7K3Wp1OScOBPYDHeljnhB7WaWZm/ShPY/v/krVlbCEb2Kr06MlyYLKkSZJGknX02FZRpw04MU0fByyNbq5HjoiHgCclHZqu1vrHFJ+ZmTVIniOSCRHR6yujUpvHqcASYBhwaUSsknQW0B4RbcAlwOWSOsiORLb1KixpDbA7MFLSbOCoiLgb+CjwHaAF+El6mJlZg+RJJL+WdGBErOztyiNiMbC4ouyMsulngeNrLDuxRnk7cEBvYzEzs2LkSSRvAt4v6QHgObIG74iIvyo0MjMzawp5EsnRhUdhZmZNK0/vvw+SXVn11jT9TJ7lzMxsaMjT19a/A58mG24XYATwvSKDMjOz5pHnyOJdwDGkS34jYh2wW5FBmZlZ88jTRvJ8RISk0lC7Ly44JjMWregaUuM5mDWzPInkKknfAEZL+ifgg8A3iw3LhrJFK7o4/ZqVbNq8FYCuDZs4/Zrs6nMnE7OBJ09j+xfJeub9ITAFOCMivlZ0YDZ0nb9k9bYkUrJp81bOX7K6QRGZWXe6PSJJY4osSWO2/6x/QrKhbt2GTb0qN7PG6vaIJCK2As9I2qOf4jFj3OiWXpWbWWPluWrrWWClpEvS0LZflfTVogOzoWvujCm0jBi2XVnLiGHMnTGlQRGZWXfyNLZflx5m/aLUoO6rtsyaQ4+JJA1wZdavZk8b78Rh1iR6TCSSJgPnAFOBF5XKI2LfAuMyM7MmkaeN5NvAxWQDW/0N8F3g8iKDMjOz5pEnkbRExM8BRcSDEXEm8NZiwzIzs2aR66otSbsA90k6VdK7gJflWbmkmZJWS+qQNK/K/FGSFqb5yyRNLJt3eipfLWlGWfm/Slol6beSrpT0osr1mplZ/8mTSE4DdgU+BhwC/AMvjLNeU7qZ8UKy8UymAidImlpR7STg8YjYD/gScF5adirZsLv7AzOBiyQNkzQ+xdEaEQeQDeE7BzMza5g8V20tB0hHJR+LiCdzrns60BER96flFwCzgLvL6swCzkzTVwMXSFIqXxARzwEPpDHdpwN/SDG3SNpMluDW5YzHzMwKkGc8klZJK4G7yG5MvFPSITnWPR5YW/a8M5VVrRMRW4CNwJhay0ZEF/BFsoTyELAxIq6vEffJktolta9fvz5HuGZmA8OiFV0cdu5SJs27jsPOXcqiFV2NDqlbeU5tXQr8c0RMjIiJwClkV3L1RFXKImedquWSXkp2tDIJGAe8WNL7qm08IuZHRGtEtI4dOzZHuGZmjVfq/bprwyaCF3q/HsjJJE8ieTIifll6EhE3A3lOb3WSDdFbMoEdT0NtqyNpOLAH8Fg3yx4JPBAR6yNiM3AN8MYcsZiZNYVm7P06TyK5VdI3JB0u6a8lXQTcKOlgSQd3s9xyYLKkSZJGkjWKt1XUaeOFhvvjgKUREal8TrqqaxIwGbiV7JTWoZJ2TW0pRwD35N1ZM7OBrhl7v87T19ZB6e+/V5S/kew0VNV7SiJii6RTgSVkV1ddGhGrJJ0FtEdEG3AJcHlqTH+MdAVWqncVWcP8FuCU1BPxMklXA7en8hXA/Nx7a2Y2wI0b3UJXlaQxkHu/VnYAMLi1trZGe3t7o8MwM+tR5QihkPV+fc6xB/Z7/3OSbouI1p7q5elrawzZ0cibyI5AbgbOiog/1R2lmZltpy96v160oqtfe8/Oc2prAXAT8Hfp+XuBhWQN32Zm1sfq6f268oimdNVXab1FyNPYvmdEfD4iHkiPs4HRhURjZmZ1acRVX3mOSG6QNAe4Kj0/Dg90ZT3o70NrM8s04qqvPEckHwa+DzyfHguAT0h6UtIThUVmTasZb6gyGyxqXd1V5FVfPSaSiNgtInaJiOHpsUsq2y0idi8sMmtazXhDldlgMXfGFFpGDNuurGXEMObOmFLYNvOc2iJ1TTKZ7UdIvKmooKy5NeMNVWaDRV9c9dVbeS7//RDwcbJuSu4ADgVuwYNbWQ3NeEOV2WBSz1VfOyNPG8nHgdcBD0bE3wDTAHenazU14tDazBonz6mtZyPiWUlIGhUR90ryfwSrqRGH1mbWOHkSSaek0cAi4GeSHseDSVkP+vvQ2swaJ88Iie9Kk2dKuoGsq/efFhqVmZk1jVxXbZVExC+KCsTMzJpTnsZ2MzOzmpxIzMysLoUmEkkzJa2W1CFpXpX5oyQtTPOXSZpYNu/0VL5a0oyy8tGSrpZ0r6R7JL2hyH0wM7PuFZZIJA0DLgSOBqYCJ0iaWlHtJODxiNgP+BJwXlp2KtloifsDM4GL0voAvgL8NCJeDbwWD7VrZtZQRR6RTAc6IuL+iCh19jiros4s4LI0fTVwRBqLfRawICKei4gHgA5guqTdgbeQDdFLRDwfERsK3AczM+tBkYlkPLC27HlnKqtaJyK2ABuBMd0suy/ZXfXflrRC0rckvbiY8M3MLI8iE4mqlFUOEF+rTq3y4cDBwMURMQ14Gtih7QVA0smS2iW1r1/vHl3MzIpSZCLpBPYuez6BHe+I31ZH0nCymx0f62bZTqAzIpal8qvJEssOImJ+RLRGROvYsWPr3BUzM6ulyESyHJgsaZKkkWSN520VddqAE9P0ccDSiIhUPidd1TWJrAv7WyPiYWBtWV9fRwB3F7gPZmbWg17d2d4bEbFF0qnAEmAYcGlErJJ0FtAeEW1kjeaXS+ogOxKZk5ZdJekqsiSxBTglIkojJf0LcEVKTvcDHyhqH8zMrGfKDgAGt9bW1mhvb290GGZmTUXSbRHR2lO9wo5IrLktWtHlbuDNLBcnEtvBohVdnH7Nym3jrndt2MTp16wEcDIxsx24ry3bwflLVm9LIiWbNm/l/CWrGxSRmQ1kTiS2g3VVxlvvrtzMhjYnEtvBuNEtvSo3s6HNicR2MHfGFFpGDNuurGXEMObOmFJjCTMbytzYbjsoNaj7qi0zy8OJxKqaPW28E4eZ5eJTW2ZmVhcnEjMzq4sTiZmZ1cWJxMzM6uJEYmZmdXEiMTOzujiRmJlZXZxIzMysLoUmEkkzJa2W1CFpXpX5oyQtTPOXSZpYNu/0VL5a0oyK5YZJWiHpx0XG38wWrejisHOXMmnedRx27lIWrehqdEhmNkgVlkgkDQMuBI4GpgInSJpaUe0k4PGI2A/4EnBeWnYq2bC7+wMzgYvS+ko+DtxTVOzNrjSeSNeGTQQvjCfiZGJmRSjyiGQ60BER90fE88ACYFZFnVnAZWn6auAISUrlCyLiuYh4AOhI60PSBOBvgW8VGHtT83giZtafikwk44G1Zc87U1nVOhGxBdgIjOlh2S8DnwL+3N3GJZ0sqV1S+/r163d2H5qSxxMxs/5UZCJRlbLIWadquaR3AH+MiNt62nhEzI+I1ohoHTt2bM/RDiIkAw5FAAAKqUlEQVQeT8TM+lORiaQT2Lvs+QRgXa06koYDewCPdbPsYcAxktaQnSp7q6TvFRF8M/N4ImbWn4pMJMuByZImSRpJ1njeVlGnDTgxTR8HLI2ISOVz0lVdk4DJwK0RcXpETIiIiWl9SyPifQXuQ1OaPW085xx7IONHtyBg/OgWzjn2QHcLb2aFKGw8kojYIulUYAkwDLg0IlZJOgtoj4g24BLgckkdZEcic9KyqyRdBdwNbAFOiYitVTdkVXk8ETPrL8oOAAa31tbWaG9vb3QYZmZNRdJtEdHaUz3f2W5mZnVxIjEzs7o4kZiZWV2cSMzMrC6FXbVl9Vm0oovzl6xm3YZNjBvdwtwZU3wVlpkNSE4kA1Cp08VSf1mlThcBJxMzG3B8amsAcqeLZtZMnEgGIHe6aGbNxIlkAHKni2bWTJxIBiB3umhmzcSN7QNQqUHdV22ZWTNwIhmg3OmimTULn9oyM7O6OJGYmVldfGqrIL4z3cyGCieSAvjOdDMbSgo9tSVppqTVkjokzasyf5SkhWn+MkkTy+adnspXS5qRyvaWdIOkeyStkvTxIuPfWb4z3cyGksISiaRhwIXA0cBU4ARJUyuqnQQ8HhH7AV8CzkvLTiUbdnd/YCZwUVrfFuDfIuI1wKHAKVXW2XC+M93MhpIij0imAx0RcX9EPA8sAGZV1JkFXJamrwaOkKRUviAinouIB4AOYHpEPBQRtwNExJPAPcCAO1fkO9PNbCgpMpGMB9aWPe9kx3/62+pExBZgIzAmz7LpNNg0YFm1jUs6WVK7pPb169fv9E7sDN+ZbmZDSZGJRFXKImedbpeV9BLgh8BpEfFEtY1HxPyIaI2I1rFjx+YMuW/Mnjaec449kPGjWxAwfnQL5xx7oBvazWxQKvKqrU5g77LnE4B1Nep0ShoO7AE81t2ykkaQJZErIuKaYkKv//Jd35luZkNFkUcky4HJkiZJGknWeN5WUacNODFNHwcsjYhI5XPSVV2TgMnAran95BLgnoj4n6ICL12+27VhE8ELl+8uWtFV1CbNzJpWYYkktXmcCiwhaxS/KiJWSTpL0jGp2iXAGEkdwCeAeWnZVcBVwN3AT4FTImIrcBjwD8BbJd2RHm/v69h9+a6ZWX6F3pAYEYuBxRVlZ5RNPwscX2PZLwBfqCi7mertJ33Kl++ameXnvraq8OW7Zmb5OZFU4ct3zczyc19bVXhgKTOz/JxIavDlu2Zm+fjUlpmZ1cWJxMzM6uJEYmZmdXEiMTOzujiRmJlZXZR1bTW4SVoPPLiTi+8FPNqH4fQ1x1cfx1cfx1efgR7fKyOix+7Th0QiqYek9ohobXQctTi++ji++ji++gz0+PLyqS0zM6uLE4mZmdXFiaRn8xsdQA8cX30cX30cX30Geny5uI3EzMzq4iMSMzOrixOJmZnVxYkkkTRT0mpJHZLmVZk/StLCNH+ZpIn9GNvekm6QdI+kVZI+XqXO4ZI2lg1BfEa1dRUY4xpJK9O226vMl6SvptfvLkkH92NsU8pelzskPSHptIo6/fr6SbpU0h8l/basbE9JP5N0X/r70hrLnpjq3CfpxH6M73xJ96b371pJo2ss2+1nocD4zpTU1dMw3D191wuMb2FZbGsk3VFj2cJfvz4XEUP+AQwDfg/sC4wE7gSmVtT5Z+DraXoOsLAf43sFcHCa3g34XZX4Dgd+3MDXcA2wVzfz3w78hGyo5EOBZQ18rx8mu9GqYa8f8BbgYOC3ZWX/BcxL0/OA86ostydwf/r70jT90n6K7yhgeJo+r1p8eT4LBcZ3JvDJHO9/t9/1ouKrmP/fwBmNev36+uEjksx0oCMi7o+I54EFwKyKOrOAy9L01cARkgofPx4gIh6KiNvT9JPAPUCzDZYyC/huZH4DjJb0igbEcQTw+4jY2Z4O+kRE3AQ8VlFc/hm7DJhdZdEZwM8i4rGIeBz4GTCzP+KLiOsjYkt6+htgQl9vN68ar18eeb7rdesuvvR/4++BK/t6u43iRJIZD6wte97Jjv+ot9VJX6aNwJh+ia5MOqU2DVhWZfYbJN0p6SeS9u/XwCCA6yXdJunkKvPzvMb9YQ61v8CNfP0AXh4RD0H24wF4WZU6A+V1/CDZEWY1PX0WinRqOvV2aY1TgwPh9Xsz8EhE3FdjfiNfv53iRJKpdmRReV10njqFkvQS4IfAaRHxRMXs28lO17wW+BqwqD9jAw6LiIOBo4FTJL2lYv5AeP1GAscAP6gyu9GvX14D4XX8LLAFuKJGlZ4+C0W5GPhL4CDgIbLTR5Ua/voBJ9D90UijXr+d5kSS6QT2Lns+AVhXq46k4cAe7Nyh9U6RNIIsiVwREddUzo+IJyLiqTS9GBghaa/+ii8i1qW/fwSuJTuFUC7Pa1y0o4HbI+KRyhmNfv2SR0qn+9LfP1ap09DXMTXuvwN4b6QT+pVyfBYKERGPRMTWiPgz8M0a22306zccOBZYWKtOo16/ejiRZJYDkyVNSr9a5wBtFXXagNIVMscBS2t9kfpaOqd6CXBPRPxPjTp/UWqzkTSd7L39Uz/F92JJu5WmyRplf1tRrQ34x3T11qHAxtJpnH5U85dgI1+/MuWfsROB/61SZwlwlKSXplM3R6WywkmaCXwaOCYinqlRJ89noaj4ytvc3lVju3m+60U6Erg3IjqrzWzk61eXRrf2D5QH2VVFvyO7ouOzqewssi8NwIvITol0ALcC+/ZjbG8iO/y+C7gjPd4OfAT4SKpzKrCK7CqU3wBv7Mf49k3bvTPFUHr9yuMTcGF6fVcCrf38/u5Klhj2KCtr2OtHltAeAjaT/Uo+iazN7efAfenvnqluK/CtsmU/mD6HHcAH+jG+DrL2hdJnsHQV4zhgcXefhX6K7/L02bqLLDm8ojK+9HyH73p/xJfKv1P6zJXV7ffXr68f7iLFzMzq4lNbZmZWFycSMzOrixOJmZnVxYnEzMzq4kRiZmZ1cSIx6wVJv05/J0p6Tx+v+zPVtmU20PnyX7OdIOlwsp5m39GLZYZFxNZu5j8VES/pi/jM+pOPSMx6QdJTafJc4M1pzIh/lTQsjdexPHUa+OFU/3BlY8l8n+xmOSQtSh3yrSp1yifpXKAlre+K8m2l3gDOl/TbNE7Fu8vWfaOkq5WNE3JFf/VIbVZueKMDMGtS8yg7IkkJYWNEvE7SKOBXkq5PdacDB0TEA+n5ByPiMUktwHJJP4yIeZJOjYiDqmzrWLKOCF8L7JWWuSnNmwbsT9Zf1K+Aw4Cb+353zWrzEYlZ3ziKrC+xO8i6+B8DTE7zbi1LIgAfk1TqimXvsnq1vAm4MrIOCR8BfgG8rmzdnZF1VHgHMLFP9sasF3xEYtY3BPxLRGzXgWJqS3m64vmRwBsi4hlJN5L149bTumt5rmx6K/5OWwP4iMRs5zxJNuxxyRLgo6m7fyS9KvXeWmkP4PGURF5NNuxwyebS8hVuAt6d2mHGkg3jemuf7IVZH/CvF7OdcxewJZ2i+g7wFbLTSrenBu/1VB8q96fARyTdBawmO71VMh+4S9LtEfHesvJrgTeQ9QgbwKci4uGUiMwazpf/mplZXXxqy8zM6uJEYmZmdXEiMTOzujiRmJlZXZxIzMysLk4kZmZWFycSMzOry/8Hb1XnLhwFTOEAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "%pylab inline\n", "import random, pylab, math\n", "\n", "def markov_pi_all_data(N, delta):\n", " x, y = 1.0, 1.0\n", " data = []\n", " for i in range(N):\n", " del_x, del_y = random.uniform(-delta, delta), random.uniform(-delta, delta)\n", " if abs(x + del_x) < 1.0 and abs(y + del_y) < 1.0:\n", " x, y = x + del_x, y + del_y\n", " if x ** 2 + y ** 2 < 1.0:\n", " data.append(4.0)\n", " else:\n", " data.append(0.0)\n", " return data\n", "\n", "poweroftwo = 20\n", "n_trials = 2 ** poweroftwo\n", "delta = 0.1\n", "data = markov_pi_all_data(n_trials, delta)\n", "errors = []\n", "bunches = []\n", "for i in range(poweroftwo):\n", " new_data = []\n", " mean = 0.0\n", " mean_sq = 0.0\n", " N = len(data)\n", " while data != []:\n", " x = data.pop()\n", " y = data.pop()\n", " mean += x + y\n", " mean_sq += x ** 2 + y ** 2\n", " new_data.append((x + y) / 2.0 )\n", " errors.append(math.sqrt(mean_sq / N - (mean / N) ** 2) / math.sqrt(N))\n", " bunches.append(i)\n", " data = new_data[:]\n", "print mean / float(N), 'mean value, estimate of pi'\n", "print mean / float(N) - pi, 'mean value, estimate of pi'\n", "pylab.plot(bunches, errors, 'o')\n", "pylab.xlabel('iteration')\n", "pylab.ylabel('apparent error')\n", "pylab.title('Bunching: naive error vs iteration number')\n", "pylab.savefig('apparent_error_bunching.png')\n", "pylab.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## The Pebble Game\n", "Simulation of the 3$\\times$3 pebble game with animation." ] }, { "cell_type": "code", "execution_count": 22, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Populating the interactive namespace from numpy and matplotlib\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAdYAAAHyCAYAAABWJ+96AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAWJQAAFiUBSVIk8AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAFlRJREFUeJzt3XuwnXV97/HPkwsthEtAhUBRCbQUKraQoDN2sHDsYCm2gVKsxXY62OLUooJ0Wk9rRcHQgnWGc5AB0eMU7agDONA2dIqXeihoW4smiCOGjnJTyk2hAnJNwnP+eHYgcCDs7P1Nfmut/XrNrNmXrPXsb2at2e/9XNbzdH3fBwCoMa/1AAAwSYQVAAoJKwAUElYAKCSsAFBIWAGgkLACQCFhBYBCwgoAhYQVAAoJKwAUElYAKCSsAFBIWAGgkLDCZnRdd2LXdWd0XXdw61lmo+u65V3Xre+6rp+67bOZ+3Zd1/1O13Vf6rruvq7rHuu67tau6y7qum7ptpsaxlPneqzw/Lqu+5ckhyd5S9/3n2g7zcx0XTc/yX8kWb7Jt5f2fX/bc9x3YZLPJjlm6lvrkzyUZNepr3+c5Ji+7//vVhsYxpw1Vph878gQ1f+Yxn0/mCGq65OclmSXvu93S/LSDMHdMckVXdftuZVmhbEnrDDBuq7bO8nKJHdMfdzcfXdP8vapL8/t+/5/933/SJL0fX9HkhOSrE2yS5L3brWhYcwJKzyHqX2rfYbNwEly8Sb7J/uu625rON6WOD/JTkneleThF7jv65JsN/X5/3r2P/Z9vyHJh6e+PGFqszHwLMIKz+3RJPckWTf19YNTX2+8/aDRXNPWdd2KJMcm+Vzf95dP4yEvn/r4QN/3dz/PfW6a+rhrkmWzHBEmkrDCc+j7/tK+75ck+bepb53a9/2STW6vajnfC+m6blGGtdXHk7xzmg/beCTj5n4vLNjk81fMYDSYeMIKW9nGzcozvH1ihj92ZZKXJTmn7/vvTvMxt0993Knrupc+z31+bpPP95rhbDDRFrzwXYBZ2rhZeSYe2NIHTL3n9pQkNyc5ZwseenWSJzLsZ/2fGY4m3nS522XYV7vRTls6G8wFwgpbWd/3lya5dFv8rK7r5iX5WJL5Sd7Z9/1j031s3/f3dl13UYYon9x13QNJPpLhj4KDknwoydIM+50XJnmyeHyYCDYFw2R5e5JXJbmi7/urZvD4dye5MkmX5D1Jvp9hLXZNkl9OckGSW6bu+6NZTwsTyBorTIiu63ZJclaSx5K8t+u6HZ91l+03+XyHqX9f1/f94xu/2ff9413XHZPk+CS/k+EApfkZjgb+Pxmi++DU3b+zVf4jMOaEFbayruvelOS8GT780r7vT53mfXdNsvPU599+gfveOPXxk0lO3PQf+uE8p5+duj1D13WvztOB/uo054I5RVhh8zbuR+xmsYztk+wxw8fuMoufuzW8Zerjv/R9f2fTSWBECSts3sbNnotnuoCpk/d/omKYF/g5t2UzfwB0XXdEhiN/k+c5Cf/mdF33miQnTX159pZPCHODg5dg8zZuMj1uah/mROu67n90XXda13X7Tl0VJ13X7dp13TuTfD7DH+Mf6/v+C00HhRHmsnGwGV3XHZDkhgzv7Vyf5N4Mbze5o+/7w1rOtqWms8badd2JSS6e+nJ9hsvE7ZKn14Q/nuRtU+cNBp6DNVbYjL7vb0pyZJLPZThZw5IM59Tdu+VcW9FXMhxodX2G/+8OGa6Mc0mS1/V9/1ZRhc2zxgoAhayxAkAhYQWAQsIKAIWEFQAKCSsAFBJWACgkrABQSFgBoJCwAkAhYQWAQiWXjeu67tYMF1i+rWJ5ANDAPkke7Pt+6WwWUnU91p2333773Q488MDdipbHHLNmzdOfL1vWbg7Gl9cQs7V27do8+uijs15OyUn4u65bvWzZsmWrV6+e9bKYm7pNLs/tuhDMhNcQs7V8+fKsWbNmTd/3y2ezHPtYAaCQsAJAIWEFgELCCgCFhBUACgkrABQSVgAoJKwAUEhYAaCQsAJAIWEFgELCCgCFhBUACgkrABQSVgAoJKwAUEhYAaCQsAJAIWEFgELCCgCFhBUACgkrABQSVgAoJKwAUEhYAaCQsAJAIWEFgELCCgCFhBUACgkrABQSVgAoJKwAUEhYAaCQsAJAIWEFgELCCgCFhBUACgkrABQSVgAoJKwAUEhYAaCQsAJAIWEFgELCCgCFhBUACgkrABQSVgAoJKwAUEhYAaCQsAJAIWEFgELCCgCFhBUACgkrABQSVgAoJKwAUEhYAaCQsAJAIWEFgELCCgCFhBUACgkrABQSVgAoJKwAUEhYAaCQsAJAIWEFgELCCgCFhBUACgkrABQSVgAoJKwAUEhYAaCQsAJAIWEFgELCCgCFhBUACgkrABQSVgAoJKwAUEhYAaCQsAJAIWEFgELCCgCFhBUACgkrABQSVgAoJKwAUEhYAaCQsAJAIWEFgELCCgCFhBUACgkrABQSVgAoJKwAUEhYAaCQsAJAIWEFgELCCgCFhBUACgkrABQSVgAoJKwAUEhYAaCQsAJAIWEFgELCCgCFhBUACgkrABQSVgAoJKwAUEhYAaCQsAJAIWEFgELCCgCFhBUACgkrABQSVgAoJKwAUEhYAaCQsAJAIWEFgELCCgCFhBUACgkrABQSVgAoJKwAUEhYAaCQsAJAIWEFgELCCgCFhBUACgkrABQSVgAoJKwAUEhYAaCQsAJAIWEFgELCCgCFhBUACgkrABRaULWgNWuSrqtaGnOZ1xGz5TVES9ZYAaCQsAJAobJNwcuWJatXVy2NuWbTTXd9324OxpfXELO1fPmwW3O2rLECQCFhBYBCwgoAhYQVAAoJKwAUElYAKCSsAFBIWAGgkLACQCFhBYBCwgoAhYQVAAoJKwAUElYAKCSsAFBIWAGgkLACQCFhBYBCwgoAhYQVAAoJKwAUElYAKCSsAFBIWAGgkLACQCFhBYBCwgoAhYQVAAoJKwAUElYAKCSsAFBIWAGgkLACQCFhBYBCwgoAhYQVAAoJKwAUElYAKCSsAFBIWAGgkLACQCFhBYBCwgoAhYQVAAoJKwAUElYAKCSsAFBIWAGgkLACQCFhBYBCwgoAhYQVAAoJKwAUElYAKCSsAFBIWAGgkLACQCFhBYBCwgoAhYQVAAoJKwAUElYAKCSsAFBIWAGgkLACQCFhBYBCwgoAhYQVAAoJKwAUElYAKCSsAFBIWAGgkLACQCFhBYBCwgoAhYQVAAoJKwAUElYAKCSsAFBIWAGgkLACQCFhBYBCwgoAhYQVAAoJKwAUElYAKCSsAFBIWAGgkLACQCFhBYBCwgoAhYQVAAoJKwAUElYAKCSsAFBIWAGgkLACQCFhBYBCwgoAhYQVAAoJKwAUElYAKCSsAFBIWAGgkLACQCFhBYBCwgoAhYQVAAoJKwAUElYAKCSsAFBIWAGgkLACQCFhBYBCwgoAhYQVAAoJKwAUElYAKCSsAFBIWAGgkLACQCFhBYBCwgoAhYQVAAoJKwAUElYAKCSsAFBIWAGgkLACQCFhBYBCwgoAhYQVAAoJKwAUElYAKCSsAFBIWAGgkLACQCFhBYBCwgoAhYQVAAoJKwAUWlC1oDVrkq6rWhpzmdcRs+U1REvWWAGgkLACQKGyTcHLliWrV1ctjblm0013fd9uDsaX1xCztXz5sFtztqyxAkAhYQWAQsIKAIWEFQAKCSsAFBJWACgkrABQSFgBoJCwAkAhYQWAQsIKAIXKzhXMlCeeSO6+O7nzzuTee5PHHkvWr0+efDKZPz9ZsCBZtChZsiTZa69k992Tef6+AZgUwjpT99wzXHVg9erhrM0335zcdVdy331bdgbw+fOTPfYYIvuzPzucBfrQQ5NDDkl23HHrzQ/AViGs03XHHcmqVckXv5h87WvJf/1XzXI3bBjWbu+8M/n615NPf3r4/rx5yf77J69+dXL00cmv/mqy8841PxOArUZYN2f16uTKK4egXn/9tv3ZTz6Z3HTTcPvbv00WLkwOPzxZsWK4vfzl23YeAKZFWJ/tv/87ufji5KKLku98p/U0T1u3Lvnnfx5up5ySvPa1ycknJ7/5m0N0ARgJwrrR6tXJBRckl1ySPPpo62le2Je/PNz22CM56aTkD/8weelLW08FMOc5HPXaa5PDDhsOGLr44vGI6qbuuSf5y79Mli5Nfvd3k1tuaT0RwJw2d8N6ww3DQUGHH57867+2nmb2NmwYDnw64IDkHe8YggvANjf3wnrbbcOa3SGHJFdd1XqaeuvWDZu099svOf305KGHWk8EMKfMnbD2ffLhDyeveMWwZrcl7zUdRw8/nJx1VnLQQcNbhADYJuZGWG++OTniiOTUU5NHHmk9zbb1ve8lr3/9cHCTtVeArW6yw9r3yfnnJz//88NBSnPZxz6WvPKVw9t1ANhqJjesDz2UHHvs8J7PubaW+nxuv31Ye33/+yd/UzhAI5MZ1ptvTl7zmuGMSTxT3ycf+EBy/PHDflgASk1eWL/0peH8ujfe2HqS0XbFFckv/uJwlDQAZSYrrB/5SHLUUcn997eeZDx885vJq16VfOUrrScBmBiTE9ZzzhnOnbt+fetJxssPf5j8yq84qAmgyGSE9cwzkz//89ZTjK9HHkl+/dcn84QZANvY+If17LOTM85oPcX4e+yx5Ljjhn3UAMzYeIf1vPOS97yn9RST47HHhmu92ucKMGPjG9Z//Mfkj/+49RST55FHkt/4DUcLA8zQeIb1299O3vzm5MknW08ymX74w+SYY7zPFWAGxi+s998/bK503tut65vfTH7v95yhCWALjVdY169Pfuu3hjMrsfVdcYUDwwC20HiFdeVKR61uaytXeo8rwBYYn7Bef33yV3/Veoq5p++Tk06y6R1gmsYjrOvWJSee6KxKrdx+e/Knf9p6CoCxMB5hXblyOJiGdj76UZuEAaZh9MP6jW8MZ1eivT/4A5uEAV7A6If1Xe+yCXhUfO97yYc+1HoKgJE22mG96qrkmmtaT8Gmzj03uffe1lMAjKzRDWvfu2LNKHr44WGfNwDPaXTD+pnPJDfc0HoKnstHP5rcemvrKQBG0miGdd265H3vaz0Fz2fduuT001tPATCSRjOsl1+e3HJL6ynYnEsuSe64o/UUACNnNMN64YWtJ+CFbNgwbBIG4BlGL6zf+lby5S+3noLp+PjHh83CADxl9MJqbXV83H33cAUcAJ4yWmF96KHkU59qPQVbwh9CAM8wWmH9+793yrxxc+21w0n6AUgyamFdtar1BMzElVe2ngBgZIxOWJ94Ivn851tPwUz4gwjgKaMT1quvthl4XF1zTfLgg62nABgJoxNWaz3j64knks99rvUUACNhdML6T//UegJmw/MHkGRUwvqDHyS33dZ6CmbjuutaTwAwEkYjrKtXt56A2frP/0x+/OPWUwA0J6zUePLJ5BvfaD0FQHOjEdavf731BFTwPAKMSFitsU4GzyPACIT10UeT73+/9RRUuOmm1hMANNc+rHfe2XoCqtx1V+sJAJprH1a/jCfHPfcMBzEBzGHCSp3164f3JAPMYe3DalPwZPF8AnNc+7DefXfrCahkCwQwx7UP68MPt56ASo880noCgKbah3X9+tYTUMnzCcxxwkqtDRtaTwDQVPuwzp/fegIqzWv/kgJoqf1vwQULWk9AJc8nMMe1D+tP/mTrCaj0Ez/RegKAptqHdffdW09ApT32aD0BQFPtw7rnnq0noNJee7WeAKCp9mH1i3hydF2yZEnrKQCaah9Wa6yT48UvThYubD0FQFPCSh3PJcAIhHXx4mS33VpPQYV99209AUBz7cOaJMuWtZ6ACsuXt54AoLnRCKtfyJPB8wggrBTyPAIIK0X23tvJPgAyKmHdd18HMI27Qw9tPQHASBiNsCbJkUe2noDZeP3rW08AMBJGJ6wrVrSegNnw/AEkGaWwHn20S46Nq2XLkp/6qdZTAIyE0Qnr4sXJa1/begpmwtoqwFNGJ6yJX9DjyvMG8JTRCuvxxyfz57eegi2x//7JIYe0ngJgZIxWWPfeO/m1X2s9BVvibW9rPQHASBmtsCbJySe3noDp2mGH5C1vaT0FwEgZvbAeeWTyMz/Tegqm47d/ezjoDICnjF5Yu87mxXHx9re3ngBg5IxeWJPk938/2WWX1lOwOYcd5nJ/AM9hNMO6eHHyJ3/Sego25+yzW08AMJJGM6xJctppyZIlrafgubzhDcMaKwD/n9EN66JFyemnt56CZ5s3z9oqwGaMbliT5K1vTfbbr/UUbOrNb05e+crWUwCMrNEO68KFyV//desp2GjRomTlytZTAIy00Q5rkhx33HCqQ9o755xkn31aTwEw0kY/rEly4YXJS17Seoq57YgjvG8VYBrGI6wveUlywQWtp5i7Fi1K/uZvhpN3ALBZ4xHWJHnjG20SbuWcc5KlS1tPATAWxiesSXLRRX7Bb2srVtgEDLAFxiusL3pRsmpVsuOOrSeZGw46KPnUp2wCBtgC4xXWxC/7beVFL0r+4R+SnXZqPQnAWBm/sCbJMcckH/hA6ykm14IFyWWXJfvu23oSgLEznmFNkve+dzgLEPXOPz953etaTwEwlsY3rEnyyU8mxx7beorJ8sEPuh4uwCyMd1gXLEguvTQ5+ujWk0yGM89M3v3u1lMAjLXxDmuSbLdd8nd/Z811ts4+O3nf+1pPATD2xj+syRDXz37WPteZmDcvOe+85M/+rPUkABNhMsKaDJuFP/3p5KyzvBVnunbeeXhf8CmntJ4EYGJMTlg3+ou/GDYNO4nE5v30Tydf/Wryhje0ngRgokxeWJPhfa7//u/eh/l8jjwyue665MADW08CMHEmM6zJcIam664bTt7PYOHC5P3vT666Ktl119bTAEykyQ1rMpyW77LLhttcv57rL/zC8IfGGWck8+e3ngZgYk12WDd64xuTG2+cm2uvCxcOb6P52teSgw9uPQ3AxJsbYU2GNdbLLksuvzzZb7/W02wbRxwxrKWeeeYQWAC2urkT1o2OOy5Zuza58MJkyZLW02wdBx887Ee9+mprqQDb2NwLazKsvf3RHyXf/e7wvtedd249UY399ks+85lkzZrkqKNaTwMwJ83NsG60aNHwvtfbb0/OPTfZf//WE83ML/1Scsklw5r4CSc4QQZAQ3M7rBstXpycdlpy003JF74wnHd41I+c3Wmn5OSTk299K7nmmuRNb7IfFWAELGg9wEjpuuHkCUcemdxxx3Cg06pVybXXJuvXt55uiOlRRyUrVgzxd3YpgJEjrM9n772TU08dbj/60XAw0KpVyRe/mNx337abY+nS4bJ4K1YMR/lut922+9kAbDFhnY7Fi4d9lyecMHx9663J6tXPvN1//+x/zstelixfnhx66PBx+fLkxS+e/XIB2GaEdSaWLh1uxx//9PceeCC5667hduedw8d77kkef3zYjLxhw3AFngULkh12SPbcM9lrr+Hjxtv227f7PwFQQlir7LLLcDvggNaTANCQo4IBoJCwAkAhYQWAQsIKAIWEFQAKCSsAFBJWACgkrABQSFgBoJCwAkAhYQWAQsIKAIWqTsK/z9q1a7N8+fKixTGXeRkxW15DzMTatWuTZJ/ZLqfr+37Ww3Rdd2uSnZPcNuuFAUAb+yR5sO/7pbNZSElYAYCBfawAUEhYAaCQsAJAIWEFgELCCgCFhBUACgkrABQSVgAoJKwAUEhYAaCQsAJAIWEFgELCCgCFhBUACgkrABQSVgAo9P8ApW3N7LzTAgYAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "%pylab qt\n", "import random, pylab\n", "\n", "sigma = 0.4 # sigma and s_map are needed for the graphical output\n", "s_map = [(1.0, 1.0), (2.0, 1.0), (3.0, 1.0), \n", " (1.0, 2.0), (2.0, 2.0), (3.0, 2.0), \n", " (1.0, 3.0), (2.0, 3.0), (3.0, 3.0)] \n", "neighbor = [[1, 3, 0, 0], [2, 4, 0, 1], [2, 5, 1, 2],\n", " [4, 6, 3, 0], [5, 7, 3, 1], [5, 8, 4, 2],\n", " [7, 6, 6, 3], [8, 7, 6, 4], [8, 8, 7, 5]]\n", "site = 8\n", "N_runs = 50\n", "for run in range(N_runs):\n", " if run < 10: number_string = '0'+str(run)\n", " else: number_string = str(run)\n", " # Begin of graphical output\n", " pylab.clf() #clears the previous figure\n", " cir = pylab.Circle(s_map[site], radius=sigma, fc='r')\n", " pylab.gca().add_patch(cir)\n", " pylab.plot([0.5, 3.5], [1.5, 1.5], 'b')\n", " pylab.plot([0.5, 3.5], [2.5, 2.5], 'b')\n", " pylab.plot([1.5, 1.5], [0.5, 3.5], 'b')\n", " pylab.plot([2.5, 2.5], [0.5, 3.5], 'b')\n", " pylab.title('t = '+ number_string)\n", " pylab.axis('scaled')\n", " pylab.axis([0.5, 3.5, 0.5, 3.5])\n", " pylab.xticks([])\n", " pylab.yticks([])\n", " pylab.pause(0.05) #Pause to generate a real time animation\n", " pylab.show()\n", " # End of graphical output\n", " site = neighbor[site][ random.randint(0, 3)]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Basic pebble throw game with homogeneous steady state probabilities. Multiple runs. The probability of landing on each grid element after a specified number of runs is visualised by a two dimensional histogram." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "7\n", "7\n", "7\n", "7\n", "4\n", "2\n", "7\n", "6\n", "4\n", "6\n" ] } ], "source": [ "import random\n", "\n", "#define the grid in terms of transitions:[0=right,1=up,2=left,3=down]\n", "neighbor = [[1, 3, 0, 0], [2, 4, 0, 1], [2, 5, 1, 2],\n", " [4, 6, 3, 0], [5, 7, 3, 1], [5, 8, 4, 2],\n", " [7, 6, 6, 3], [8, 7, 6, 4], [8, 8, 7, 5]] \n", "t_max = 4 #maximum number of pebble throws at each run\n", "N_runs = 10 #number of runs\n", "for run in range(N_runs): #run many times\n", " site = 8 #initial site\n", " t = 0 #initialise time\n", " while t < t_max: #a pebble game of t_max throws\n", " t += 1\n", " site = neighbor[site][random.randint(0, 3)] #Pick a random transition from the initial site\n", " print site #prints the final site after t_max throws in each run" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[[0. 0. 2.25]\n", " [0. 0. 0. ]\n", " [0. 0. 0. ]]\n", "[[0. 0. 1.]\n", " [0. 0. 0.]\n", " [0. 0. 0.]]\n", "[[0. 0.5585625 1.1299725]\n", " [0. 0. 0.561465 ]\n", " [0. 0. 0. ]]\n", "[[0. 0.24825 0.50221]\n", " [0. 0. 0.24954]\n", " [0. 0. 0. ]]\n", "[[0.137655 0.4276125 0.841815 ]\n", " [0. 0.2833425 0.419985 ]\n", " [0. 0. 0.13959 ]]\n", "[[0.06118 0.19005 0.37414]\n", " [0. 0.12593 0.18666]\n", " [0. 0. 0.06204]]\n", "[[0.1775925 0.4267125 0.6258825]\n", " [0.10782 0.212625 0.4200525]\n", " [0. 0.1052325 0.1740825]]\n", "[[0.07893 0.18965 0.27817]\n", " [0.04792 0.0945 0.18669]\n", " [0. 0.04677 0.07737]]\n", "[[0.2203875 0.358515 0.5262525]\n", " [0.1221075 0.26874 0.360585 ]\n", " [0.052425 0.1224675 0.21852 ]]\n", "[[0.09795 0.15934 0.23389]\n", " [0.05427 0.11944 0.16026]\n", " [0.0233 0.05443 0.09712]]\n", "[[0.2317275 0.3440025 0.44424 ]\n", " [0.1649925 0.2417625 0.33723 ]\n", " [0.0882 0.1676925 0.2301525]]\n", "[[0.10299 0.15289 0.19744]\n", " [0.07333 0.10745 0.14988]\n", " [0.0392 0.07453 0.10229]]\n", "[[0.2418075 0.311715 0.39024 ]\n", " [0.18441 0.2544525 0.31446 ]\n", " [0.12519 0.18315 0.244575 ]]\n", "[[0.10747 0.13854 0.17344]\n", " [0.08196 0.11309 0.13976]\n", " [0.05564 0.0814 0.1087 ]]\n", "[[0.2434725 0.3036825 0.356175 ]\n", " [0.19872 0.2446875 0.304065 ]\n", " [0.152775 0.20106 0.2453625]]\n", "[[0.10821 0.13497 0.1583 ]\n", " [0.08832 0.10875 0.13514]\n", " [0.0679 0.08936 0.10905]]\n", "[[0.2458575 0.287865 0.3297825]\n", " [0.20898 0.248625 0.2894625]\n", " [0.18018 0.2109375 0.24831 ]]\n", "[[0.10927 0.12794 0.14657]\n", " [0.09288 0.1105 0.12865]\n", " [0.08008 0.09375 0.11036]]\n", "[[0.2522925 0.2791125 0.2991825]\n", " [0.2209725 0.24426 0.28458 ]\n", " [0.1943325 0.2231775 0.25209 ]]\n", "[[0.11213 0.12405 0.13297]\n", " [0.09821 0.10856 0.12648]\n", " [0.08637 0.09919 0.11204]]\n" ] } ], "source": [ "%matplotlib qt\n", "import random\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "\n", "xvec = {1:3, 2:2, 3:1, 4:3, 5:2, 6:1, 7:3, 8:2, 9:1} \n", "yvec = {1:1, 2:1, 3:1, 4:2, 5:2, 6:2, 7:3, 8:3, 9:3} \n", "\n", "neighbor = {1 : [2, 4, 1, 1], 2 : [3, 5, 1, 2], 3 : [3, 6, 2, 3],\n", " 4 : [5, 7, 4, 1], 5 : [6, 8, 4, 2], 6 : [6, 9, 5, 3],\n", " 7 : [8, 7, 7, 4], 8 : [9, 8, 7, 5], 9 : [9, 9, 8, 6]}\n", "\n", "N_runs = 10\n", "for run in range(N_runs):\n", " list_vec = []\n", " if run < 10: run_str= '0'+str(run)\n", " else: run_str = str(run)\n", " for n_runs in range(100000): \n", " pos = 9\n", " for iter in range(run):\n", " pos = neighbor[pos][ random.randint(0, 3)]\n", " list_vec.append(pos)\n", "\n", " x = [xvec[k] for k in list_vec]\n", " y = [yvec[k] for k in list_vec]\n", "\n", " plt.xticks([])\n", " plt.yticks([])\n", " H, xedges, yedges = np.histogram2d(x, y, bins=(3, 3), \n", " range=[[1,3],[1,3]], normed=True)\n", " print H\n", " H /= np.sum(H)\n", " print H\n", " extent = [yedges[0], yedges[-1], xedges[-1], xedges[0]]\n", " plt.clf() #clears the previous figure\n", " histo = plt.imshow(H, extent=extent, interpolation='nearest', vmin=0, vmax=1.00)\n", " histo.set_cmap('hot')\n", " plt.colorbar()\n", " plt.title('t = '+str(run),fontsize=22)\n", " plt.savefig('3x3_pebble_run_'+run_str+'.png')\n", " plt.pause(0.05)\n", " plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Transfer Matrix Method\n", "Exact simulation of the 3$\\times$3 pebble game can be done using the transfer matrix for 50 throws. The evolution of the probability vector indicates that at large times the effect of the multiplication by the transfer matrix on the probability vector is negligible compared to the changes at small times. That means that the equilibrium probability vector is an eigenvector of the transfer matrix with eigenvalue 1." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "0 ['0.00000', '0.00000', '0.00000', '0.00000', '0.00000', '0.00000', '0.00000', '0.00000', '1.00000']\n", "1 ['0.00000', '0.00000', '0.00000', '0.00000', '0.00000', '0.25000', '0.00000', '0.25000', '0.50000']\n", "2 ['0.00000', '0.00000', '0.06250', '0.00000', '0.12500', '0.18750', '0.06250', '0.18750', '0.37500']\n", "3 ['0.00000', '0.04688', '0.07812', '0.04688', '0.09375', '0.18750', '0.07812', '0.18750', '0.28125']\n", "4 ['0.02344', '0.05469', '0.09766', '0.05469', '0.11719', '0.16016', '0.09766', '0.16016', '0.23438']\n", "5 ['0.03906', '0.07324', '0.10254', '0.07324', '0.10742', '0.15234', '0.10254', '0.15234', '0.19727']\n", "6 ['0.05615', '0.08057', '0.10767', '0.08057', '0.11279', '0.13989', '0.10767', '0.13989', '0.17480']\n", "7 ['0.06836', '0.08929', '0.10895', '0.08929', '0.11023', '0.13379', '0.10895', '0.13379', '0.15735']\n", "8 ['0.07883', '0.09421', '0.11024', '0.09421', '0.11154', '0.12758', '0.11024', '0.12758', '0.14557']\n", "9 ['0.08652', '0.09871', '0.11057', '0.09871', '0.11089', '0.12373', '0.11057', '0.12373', '0.13657']\n", "10 ['0.09261', '0.10167', '0.11089', '0.10167', '0.11122', '0.12044', '0.11089', '0.12044', '0.13015']\n", "11 ['0.09714', '0.10410', '0.11098', '0.10410', '0.11106', '0.11818', '0.11098', '0.11818', '0.12530']\n", "12 ['0.10062', '0.10582', '0.11106', '0.10582', '0.11114', '0.11638', '0.11106', '0.11638', '0.12174']\n", "13 ['0.10322', '0.10716', '0.11108', '0.10716', '0.11110', '0.11508', '0.11108', '0.11508', '0.11906']\n", "14 ['0.10519', '0.10814', '0.11110', '0.10814', '0.11112', '0.11408', '0.11110', '0.11408', '0.11707']\n", "15 ['0.10666', '0.10889', '0.11110', '0.10889', '0.11111', '0.11334', '0.11110', '0.11334', '0.11557']\n", "16 ['0.10777', '0.10944', '0.11111', '0.10944', '0.11111', '0.11278', '0.11111', '0.11278', '0.11446']\n", "17 ['0.10861', '0.10986', '0.11111', '0.10986', '0.11111', '0.11236', '0.11111', '0.11236', '0.11362']\n", "18 ['0.10923', '0.11017', '0.11111', '0.11017', '0.11111', '0.11205', '0.11111', '0.11205', '0.11299']\n", "19 ['0.10970', '0.11041', '0.11111', '0.11041', '0.11111', '0.11182', '0.11111', '0.11182', '0.11252']\n", "20 ['0.11005', '0.11058', '0.11111', '0.11058', '0.11111', '0.11164', '0.11111', '0.11164', '0.11217']\n", "21 ['0.11032', '0.11071', '0.11111', '0.11071', '0.11111', '0.11151', '0.11111', '0.11151', '0.11190']\n", "22 ['0.11052', '0.11081', '0.11111', '0.11081', '0.11111', '0.11141', '0.11111', '0.11141', '0.11171']\n", "23 ['0.11067', '0.11089', '0.11111', '0.11089', '0.11111', '0.11133', '0.11111', '0.11133', '0.11156']\n", "24 ['0.11078', '0.11094', '0.11111', '0.11094', '0.11111', '0.11128', '0.11111', '0.11128', '0.11145']\n", "25 ['0.11086', '0.11099', '0.11111', '0.11099', '0.11111', '0.11124', '0.11111', '0.11124', '0.11136']\n", "26 ['0.11092', '0.11102', '0.11111', '0.11102', '0.11111', '0.11121', '0.11111', '0.11121', '0.11130']\n", "27 ['0.11097', '0.11104', '0.11111', '0.11104', '0.11111', '0.11118', '0.11111', '0.11118', '0.11125']\n", "28 ['0.11101', '0.11106', '0.11111', '0.11106', '0.11111', '0.11116', '0.11111', '0.11116', '0.11122']\n", "29 ['0.11103', '0.11107', '0.11111', '0.11107', '0.11111', '0.11115', '0.11111', '0.11115', '0.11119']\n", "30 ['0.11105', '0.11108', '0.11111', '0.11108', '0.11111', '0.11114', '0.11111', '0.11114', '0.11117']\n", "31 ['0.11107', '0.11109', '0.11111', '0.11109', '0.11111', '0.11113', '0.11111', '0.11113', '0.11116']\n", "32 ['0.11108', '0.11109', '0.11111', '0.11109', '0.11111', '0.11113', '0.11111', '0.11113', '0.11114']\n", "33 ['0.11109', '0.11110', '0.11111', '0.11110', '0.11111', '0.11112', '0.11111', '0.11112', '0.11114']\n", "34 ['0.11109', '0.11110', '0.11111', '0.11110', '0.11111', '0.11112', '0.11111', '0.11112', '0.11113']\n", "35 ['0.11110', '0.11110', '0.11111', '0.11110', '0.11111', '0.11112', '0.11111', '0.11112', '0.11113']\n", "36 ['0.11110', '0.11111', '0.11111', '0.11111', '0.11111', '0.11112', '0.11111', '0.11112', '0.11112']\n", "37 ['0.11110', '0.11111', '0.11111', '0.11111', '0.11111', '0.11112', '0.11111', '0.11112', '0.11112']\n", "38 ['0.11111', '0.11111', '0.11111', '0.11111', '0.11111', '0.11111', '0.11111', '0.11111', '0.11112']\n", "39 ['0.11111', '0.11111', '0.11111', '0.11111', '0.11111', '0.11111', '0.11111', '0.11111', '0.11112']\n", "40 ['0.11111', '0.11111', '0.11111', '0.11111', '0.11111', '0.11111', '0.11111', '0.11111', '0.11111']\n", "41 ['0.11111', '0.11111', '0.11111', '0.11111', '0.11111', '0.11111', '0.11111', '0.11111', '0.11111']\n", "42 ['0.11111', '0.11111', '0.11111', '0.11111', '0.11111', '0.11111', '0.11111', '0.11111', '0.11111']\n", "43 ['0.11111', '0.11111', '0.11111', '0.11111', '0.11111', '0.11111', '0.11111', '0.11111', '0.11111']\n", "44 ['0.11111', '0.11111', '0.11111', '0.11111', '0.11111', '0.11111', '0.11111', '0.11111', '0.11111']\n", "45 ['0.11111', '0.11111', '0.11111', '0.11111', '0.11111', '0.11111', '0.11111', '0.11111', '0.11111']\n", "46 ['0.11111', '0.11111', '0.11111', '0.11111', '0.11111', '0.11111', '0.11111', '0.11111', '0.11111']\n", "47 ['0.11111', '0.11111', '0.11111', '0.11111', '0.11111', '0.11111', '0.11111', '0.11111', '0.11111']\n", "48 ['0.11111', '0.11111', '0.11111', '0.11111', '0.11111', '0.11111', '0.11111', '0.11111', '0.11111']\n", "49 ['0.11111', '0.11111', '0.11111', '0.11111', '0.11111', '0.11111', '0.11111', '0.11111', '0.11111']\n" ] } ], "source": [ "import numpy\n", "\n", "neighbor = [[1, 3, 0, 0], [2, 4, 0, 1], [2, 5, 1, 2],\n", " [4, 6, 3, 0], [5, 7, 3, 1], [5, 8, 4, 2],\n", " [7, 6, 6, 3], [8, 7, 6, 4], [8, 8, 7, 5]]\n", "\n", "transfer = numpy.zeros((9, 9)) #initialise the transfer matrix\n", "#the transfer matrix is constructed by adding the probabilities for all of the the terms \n", "#that contribute to a specific transition to determine each element in the matrix\n", "for k in range(9):\n", " for neigh in range(4): transfer[neighbor[k][neigh], k] += 0.25 \n", "\n", "#initial probability vector: pebble is at position 8 with certainty\n", "position = numpy.zeros(9)\n", "position[8] = 1.0 \n", "for t in range(50):\n", " print t,' ',[\"%0.5f\" % i for i in position]\n", " position = numpy.dot(transfer, position) #probability vector at time t" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Find all of the eigenvalues of the transfer matrix. What do the eigenvalues mean, apart from 1, which is associated with equilibrium?" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[ 1. 0.75 0.75 0.5 0.25 0.25 -0. -0. -0.5 ]\n" ] } ], "source": [ "np.set_printoptions(suppress=True) #suppress scientific format\n", "\n", "eigenvalues, eigenvectors = numpy.linalg.eig(transfer) #eigenvalues of the transfer matrix\n", "idx = np.argsort(-eigenvalues) #sort the eigenvalues in descending order\n", "eigenvalues = eigenvalues[idx]\n", "print eigenvalues #print the sorted eigenvalues\n", "# you may print the eigenvectors by uncommenting the following lines...\n", "#for iter in range(9):\n", "# print eigenvalues[iter]\n", "# for i in range(9):\n", "# print eigenvectors[i][iter]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To obtain the deviation from the equilibrim value at each time, subtract the equilibrium probabilities from each probability vector and take the absolute value. Specifically, the evolution of the deviation in the site 0 is plotted in semi-log scale. It is found that the resulting line is approximately parallel to $0.75^t$ in semi-log scale." ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "0 ['0.11111', '0.11111', '0.11111', '0.11111', '0.11111', '0.11111', '0.11111', '0.11111', '0.88889']\n", "1 ['0.11111', '0.11111', '0.11111', '0.11111', '0.11111', '0.13889', '0.11111', '0.13889', '0.38889']\n", "2 ['0.11111', '0.11111', '0.04861', '0.11111', '0.01389', '0.07639', '0.04861', '0.07639', '0.26389']\n", "3 ['0.11111', '0.06424', '0.03299', '0.06424', '0.01736', '0.07639', '0.03299', '0.07639', '0.17014']\n", "4 ['0.08767', '0.05642', '0.01345', '0.05642', '0.00608', '0.04905', '0.01345', '0.04905', '0.12326']\n", "5 ['0.07205', '0.03787', '0.00857', '0.03787', '0.00369', '0.04123', '0.00857', '0.04123', '0.08615']\n", "6 ['0.05496', '0.03054', '0.00345', '0.03054', '0.00168', '0.02878', '0.00345', '0.02878', '0.06369']\n", "7 ['0.04275', '0.02182', '0.00216', '0.02182', '0.00088', '0.02268', '0.00216', '0.02268', '0.04624']\n", "8 ['0.03228', '0.01690', '0.00087', '0.01690', '0.00043', '0.01647', '0.00087', '0.01647', '0.03446']\n", "9 ['0.02459', '0.01241', '0.00054', '0.01241', '0.00022', '0.01262', '0.00054', '0.01262', '0.02546']\n", "10 ['0.01850', '0.00944', '0.00022', '0.00944', '0.00011', '0.00933', '0.00022', '0.00933', '0.01904']\n", "11 ['0.01397', '0.00701', '0.00014', '0.00701', '0.00005', '0.00707', '0.00014', '0.00707', '0.01419']\n", "12 ['0.01049', '0.00529', '0.00005', '0.00529', '0.00003', '0.00527', '0.00005', '0.00527', '0.01063']\n", "13 ['0.00789', '0.00395', '0.00003', '0.00395', '0.00001', '0.00397', '0.00003', '0.00397', '0.00795']\n", "14 ['0.00592', '0.00297', '0.00001', '0.00297', '0.00001', '0.00297', '0.00001', '0.00297', '0.00596']\n", "15 ['0.00445', '0.00223', '0.00001', '0.00223', '0.00000', '0.00223', '0.00001', '0.00223', '0.00446']\n", "16 ['0.00334', '0.00167', '0.00000', '0.00167', '0.00000', '0.00167', '0.00000', '0.00167', '0.00335']\n", "17 ['0.00250', '0.00125', '0.00000', '0.00125', '0.00000', '0.00125', '0.00000', '0.00125', '0.00251']\n", "18 ['0.00188', '0.00094', '0.00000', '0.00094', '0.00000', '0.00094', '0.00000', '0.00094', '0.00188']\n", "19 ['0.00141', '0.00070', '0.00000', '0.00070', '0.00000', '0.00070', '0.00000', '0.00070', '0.00141']\n", "20 ['0.00106', '0.00053', '0.00000', '0.00053', '0.00000', '0.00053', '0.00000', '0.00053', '0.00106']\n", "21 ['0.00079', '0.00040', '0.00000', '0.00040', '0.00000', '0.00040', '0.00000', '0.00040', '0.00079']\n", "22 ['0.00059', '0.00030', '0.00000', '0.00030', '0.00000', '0.00030', '0.00000', '0.00030', '0.00059']\n", "23 ['0.00045', '0.00022', '0.00000', '0.00022', '0.00000', '0.00022', '0.00000', '0.00022', '0.00045']\n", "24 ['0.00033', '0.00017', '0.00000', '0.00017', '0.00000', '0.00017', '0.00000', '0.00017', '0.00033']\n", "25 ['0.00025', '0.00013', '0.00000', '0.00013', '0.00000', '0.00013', '0.00000', '0.00013', '0.00025']\n", "26 ['0.00019', '0.00009', '0.00000', '0.00009', '0.00000', '0.00009', '0.00000', '0.00009', '0.00019']\n", "27 ['0.00014', '0.00007', '0.00000', '0.00007', '0.00000', '0.00007', '0.00000', '0.00007', '0.00014']\n", "28 ['0.00011', '0.00005', '0.00000', '0.00005', '0.00000', '0.00005', '0.00000', '0.00005', '0.00011']\n", "29 ['0.00008', '0.00004', '0.00000', '0.00004', '0.00000', '0.00004', '0.00000', '0.00004', '0.00008']\n", "30 ['0.00006', '0.00003', '0.00000', '0.00003', '0.00000', '0.00003', '0.00000', '0.00003', '0.00006']\n", "31 ['0.00004', '0.00002', '0.00000', '0.00002', '0.00000', '0.00002', '0.00000', '0.00002', '0.00004']\n", "32 ['0.00003', '0.00002', '0.00000', '0.00002', '0.00000', '0.00002', '0.00000', '0.00002', '0.00003']\n", "33 ['0.00003', '0.00001', '0.00000', '0.00001', '0.00000', '0.00001', '0.00000', '0.00001', '0.00003']\n", "34 ['0.00002', '0.00001', '0.00000', '0.00001', '0.00000', '0.00001', '0.00000', '0.00001', '0.00002']\n", "35 ['0.00001', '0.00001', '0.00000', '0.00001', '0.00000', '0.00001', '0.00000', '0.00001', '0.00001']\n", "36 ['0.00001', '0.00001', '0.00000', '0.00001', '0.00000', '0.00001', '0.00000', '0.00001', '0.00001']\n", "37 ['0.00001', '0.00000', '0.00000', '0.00000', '0.00000', '0.00000', '0.00000', '0.00000', '0.00001']\n", "38 ['0.00001', '0.00000', '0.00000', '0.00000', '0.00000', '0.00000', '0.00000', '0.00000', '0.00001']\n", "39 ['0.00000', '0.00000', '0.00000', '0.00000', '0.00000', '0.00000', '0.00000', '0.00000', '0.00000']\n", "40 ['0.00000', '0.00000', '0.00000', '0.00000', '0.00000', '0.00000', '0.00000', '0.00000', '0.00000']\n", "41 ['0.00000', '0.00000', '0.00000', '0.00000', '0.00000', '0.00000', '0.00000', '0.00000', '0.00000']\n", "42 ['0.00000', '0.00000', '0.00000', '0.00000', '0.00000', '0.00000', '0.00000', '0.00000', '0.00000']\n", "43 ['0.00000', '0.00000', '0.00000', '0.00000', '0.00000', '0.00000', '0.00000', '0.00000', '0.00000']\n", "44 ['0.00000', '0.00000', '0.00000', '0.00000', '0.00000', '0.00000', '0.00000', '0.00000', '0.00000']\n", "45 ['0.00000', '0.00000', '0.00000', '0.00000', '0.00000', '0.00000', '0.00000', '0.00000', '0.00000']\n", "46 ['0.00000', '0.00000', '0.00000', '0.00000', '0.00000', '0.00000', '0.00000', '0.00000', '0.00000']\n", "47 ['0.00000', '0.00000', '0.00000', '0.00000', '0.00000', '0.00000', '0.00000', '0.00000', '0.00000']\n", "48 ['0.00000', '0.00000', '0.00000', '0.00000', '0.00000', '0.00000', '0.00000', '0.00000', '0.00000']\n", "49 ['0.00000', '0.00000', '0.00000', '0.00000', '0.00000', '0.00000', '0.00000', '0.00000', '0.00000']\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY4AAAEKCAYAAAAFJbKyAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzs3Xd8T9f/wPHXyRIkNjGiglgpasQqMuy9qkWNGqVKVIuu77e79e1CtSi1q1Wz2tpqJLFH0NqbVsxStYnE+/fHTf00JT6J3IxP3s/H4/Ooz83n3vc5Unnn3nPO+xgRQSmllHKUS1o3QCmlVMaiiUMppVSSaOJQSimVJJo4lFJKJYkmDqWUUkmiiUMppVSSaOJQSimVJJo4lFJKJYkmDqWUUkniltYNsEO+fPnEz88vWedevXqV7Nmzp2yDMgDtd+aSWfsNmbfvjvR769at50Qk/4Ou5VSJwxjTEmjp7+9PVFRUsq4RERFBSEhIirYrI9B+Zy6Ztd+QefvuSL+NMb85ci2nelQlIgtEpE/OnDnTuilKKeW0nCpxKKWUsp8mDqWUUkniVGMcSinncevWLaKjo7lx40aKXjdnzpzs3bs3Ra+ZEdzdb09PT3x9fXF3d0/WtdJ94jDGZAe+BGKACBGZnsZNUkqlgujoaLy9vfHz88MYk2LXvXz5Mt7e3il2vYzi736LCOfPnyc6OprixYsn61pp8qjKGDPZGHPWGLMrwfEmxpj9xphDxpjX4g+3A+aKSG+gVao3VimVJm7cuEHevHlTNGkoMMaQN2/eh7qTS6sxjqlAk7sPGGNcgTFAUyAA6GSMCQB8gePxH4tLxTYqpdKYJg17POzfa5okDhFZDfyZ4HB14JCIHBGRGGAm0BqIxkoeYHd7Dy6n8IklcFvzk1JK3U96GuMowv/fWYCVMGoAXwCjjTHNgQX3O9kY0wfoA+Dj40NERESSG1Bm35eUPr2KyyOWcci/DxdzBST5GhnVlStXkvV3ltFpv9OvnDlzcvny5RS/blxcnC3XTe8S9vvGjRvJ/n8gPSWOe907iYhcBXo86GQRGQ+MBwgMDJRkrQwNDmb3nKE8Gv0dlX95HSp2gIbvgXfBpF8rg9HVtJlLRuj33r17bRnEdrbB8Vu3bvH2229z7do1YmJiGDlyJIMGDcLLy4vIyEjmz59P/vz5/9VvT09PKleunKyY6WkdRzRQ9K73vsDJpFzAGNPSGDP+4sWLyWuBMfxRoA6EbYG6g2H3DzCqKqz7HGJjkndNpVSmcP36dYKDg4mLsx51R0dHM2vWLABiYmIICgoiNjY2xeOOHz+e69evkytXLq5cucL48ePp3r07H330Efny5SN//geWnkqy9JQ4tgCljDHFjTEeQEdgfpq0xCM71H8L+m2EYrVh+Vsw9nE4tDJNmqOUSv8mT55Mu3btcHV1BWDlypVs27YNAA8PD+rXr38nkaSk7du389FHH/HOO+8wbdo0oqKiqFChApcuXcLHxyfF40HaTcedAWwAyhhjoo0xvUQkFggDlgF7gdkisjsp103xWlV5S0Ln2fD0bLgdC9+2g5md4YJDdcCUUhncoUOHyJ8/P35+flSqVIk8efJQsmRJLl269K/PTp8+ndatWwOwdu1aBg0axNy5c6lUqRJHjx6lTZs2TJ+evGVoibWjdevWdO/enVdeeYWlS5fSuHFjevXqxRtvvEHp0qUfqv/3JSJO8wJaAuP9/f0lucLDw+/9hZjrIpGfinxQUOT9AiLhH4rEXEt2nPTmvv12ctrv9GvPnj22XPfSpUtJ+nybNm1k9erVIiISHBwsO3bs+Ndnbt68KT4+Pv841rhxY9m5c+ed97GxsZIvX75ktNjxdtxt+vTpsnLlyjvvE/b7Xn+/QJQ48LM2PQ2OPzQRWQAsCAwM7J3iF3f3hKAh8FhH+PlNiPgQfpkOjT+Ess1B55srZZt3F+xmz8l//5afHHFxcbi6uhJQOAdvt3z0gZ/fvXs35cuXB2Dfvn2UKVPmX585d+4cuXLl+sex/fv3/+Ozrq6ueHh4/GuQukGDBpw+ffpf1xw6dOidOxhH23G3p59++oF9Sy6nShx378dhm5y+8OQUCOwBS16FWZ2hZD1o+gnkK2VfXKVUqrt+/To3btwgd+7cHD9+nLx58+Lh4cHVq1fp168fHh4ehISE0KxZs3+sxD5//jw5c+b8Vy2omzdv4unp+Y9jK1asSHY70opTJQ5b7zgSKh4Ez62BLRMg/EP4shbUfB6CX4EszjPVT6n0wJE7A0clZTrunj17KFeuHGBND/77z/PmzaN9+/a0bNmSDh060LlzZ+Li4rhx4waenp4cPXqUwoUL/+Na58+fJ3/+/MkqLHi/diRMYJ07d07ytZMjPc2qemgPPR03qVzdrGQxYCs81gHWfwGjAmHHbLDGXJRSGdjdj4eyZs3Ktm3b2LdvH9HR0RQtaq0e+HsWVaNGjVi7di0AZcuW5dy5c5QvX57169cDEB4eTrNmzVK0HX8nsAkTJjB/fupNQtU7jpTglR9aj4GqPWHxEJjXG6ImW4+vClVM1aYopVJOt27d7vy5bt26HDlyBABfX1+io6OpVKkSt2/fBiAsLIwRI0bQoEEDvLy82Lx58z+u9d133/Hhhx+maDt++OEHKlSoAPx/AksNTnXHkeZ8q8KzK6HVKDh3EMYHw6LBcC1hWS6lVEbWrl07vv/+e55//nlatmwJQOXKlQkNDb2zAPBuMTExtGnT5oED2kn1dwID7iSw1OBUdxypMjj+IC4uUKUblGtlzbzaPAF2zYP6b0KVZ8Al9X4rUErZI3v27EyZMuVfx3v27HnPz3t4ePzjriGltGvXjrCwMBYtWnQngaUGp0ocafao6l6y5oKmH1tJZPErsPAl2DoVmg2DotXTunVKKSdwvwRmN31UZTefR6H7QnhiElz5AyY1hB/6wuUzad0ypZRKFk0cqcEYqNDeKp5Y5yXYOdcqnrh+NMTdSuvWKaVUkjhV4kj16bhJlcULGrxjFU98pCb8/F8YWxuORKRxw5RSynFOlTgkpYsc2iWfP3SeA51mQtxNmNYaZnWFv35P65YppdQDOVXiyFCMgTJNod8mqPcGHFwOo6tD5CdwK/mbyCullN00caQ1d08Ietka/yjdGMKHwpjqsG+xrj5XSqVLmjjSi1xF4amvodt8cM8KMzvB9PZw7lBat0wppf7BqRJHuh8cd0SJYOi71irXfnwzfFkTlr8NN6+kdcuUUgpwssSRYQbHH8TVHWr1s4onVnwK1o2E0YGwY44+vlLKSVy4cCGtm5BsTpU4nI5XAWjzJfRaAV4+MO9ZmNocTu9K65YppR7SSy+9lNZNSDZNHBlB0WrQexW0GAln98JXdWHxy3A94/7GolRGsXTpUsqUKYO/vz8fffTRv76+f/9+KlWqdOeVI0cORo4cCYCfnx8VKlSgUqVKBAYG/uOa+/btY9iwYanWj5TkVLWqnJqLq7XrYEBrCP8fbJkIu76H+m9B5a5aPFEpG8TFxdG/f3+WL1+Or68v1apVo1WrVgQEBNz5TJkyZfjll1/ufL5IkSK0bdv2ztfDw8PJly/fP66bL18+unTpQlhYWOp0JIWl+zsOY0wJY8wkY8zctG5LupAtDzQfBn0iIV9pWDAQJtaH6Ki0bplSTmfz5s34+/tTokQJPDw86NixIz/99NN9P79y5UpKlixJsWLFEr3ujh07eOyxx1K6uanG1jsOY8xkoAVwVkTK33W8CfA54ApMFJF/3//FE5EjQK/USBzh+8+y9OgtDrgcdvgcd1cXmpYvRMGcng/+cEoqVBF6LIGdc+DnN63kUamzVdLEq0DqtkUpJ3XixIk7O/2Btf/Fpk2b7vv5mTNn0qlTpzvvjTE0atQIYwzPPfccffr0Aaw7jokTJ5IvX74728BmJHY/qpoKjAam/X3AGOMKjAEaAtHAFmPMfKwkknB7rJ4ictbmNt6xaMcp5u6Pgf37knTeh4v30T7Ql+eDS1I0TzabWncPxlizrso0tVacbxwLexdAyOtQvbc1O0spZ7DkNTi9M0UulTUu1tr2uWAFaHrf31kBkHvMYjTG3POzMTExzJ8//x+7/K1bt47ChQtz9uxZGjZsSNmyZQkKCqJVq1a0atXq4TqShmxNHCKy2hjjl+BwdeBQ/J0ExpiZQGsR+RDr7iTNfNCmPA1y/0ndunUdPuePyzcZv+YIc6KOM3vLcdpWLkL/UH/88mW3saUJZPGGRu9be38seQWWvQ7bpln7gZQITr12KOVkfH19OX78+J330dHRFC5c+J6fXbJkCVWqVMHHx+fOsb8/W6BAAdq2bcvmzZsJCgqyt9GpwNwro6ZoACtxLPz7UZUxpj3QRESejX/fFaghIvccJTLG5AWGYt2hTIxPMPf6XB+gD4CPj0/VmTNnJqu9V65cwcvLK8nn/XnjNkuO3iLieCyxt6FmIVeaFHenWI5UHrQWId+5TZQ8PImsN85yNn9tDpfswU3P/Imeltx+Z3Ta7/QrZ86c2LGbZ1xcnMP7c8fGxlKlShXmz59P4cKFCQkJYdKkSfd8vNS9e3caNGhAly5dALh69Sq3b9/G29ubq1ev0rp1a1599VUaNmyYov1xVMJ+Hzp0iISLpUNDQ7eKSGDCc/9FRGx9AX7ArrveP4mVAP5+3xUYlUKxWgLj/f39JbnCw8OTfa6IyJlL12Xooj1S7s0lUuzVhdJq1BqZtfl3uXrz1kNdN8liromEfyTyfgGRDwqKRH4iEnP9vh9/2H5nVNrv9GvPnj22XPfSpUtJ+vyiRYukVKlSUqJECfnggw/uHG/atKmcOHFCRESuXr0qefLkkb/++uvO1w8fPiwVK1aUihUrSkBAwD/OTQsJ+32vv18gShz4WZsW03GjgaJ3vfcFTqZBO2xRwNuT/zQrR/8Qf+Ztj+a7Tb/zyvc7eH/hHtpWKcLTNR6hbMEc9jfEPSuEvAqVOsGy/8KqD2D7dOvxVenG9sdXykk0a9aMZs2a/ev44sWL7/w5W7ZsnD9//h9fL1GiBL/++qvt7UsLaTEddwtQyhhT3BjjAXQE5qfEhSUdlRzJmc2dHrWL8/NLQczpW4sGAT7M3HKcJiPX0HniRnadSKV6WrkegQ7fQNcfwdUDvnsKpj8F5x2fOaaUUnezNXEYY2YAG4AyxphoY0wvEYkFwoBlwF5gtojsTqF46a7IoTGGan55+KxDJTa9Xp//NCvL3lOXaTl6LYNm/8LJv66nTkNKhlrFExt9AL+tt4onrnwPYq6mTnyllNOwNXGISCcRKSQi7iLiKyKT4o8vFpHSIlJSRIamYLx0c8dxL7mze9AnqCQRL4fQJ6gEC3ecInRYBMOW7efKzVj7G+DmAY8PgAFRUP4JWDMcRlezVqBr8USllIPS/crxpEiPdxz3ksPTndeblmPloGAaP1qQ0eGHCPk0nG83/satuNv2N8C7ILQdBz1/hmx5YW5PHvv1DTizx/7YSqkMz6kSR3q/40ioaJ5sfNGpMj/2r03xfNl548dd1B8eyfdbo4m7nQp3AI/UgD4R0OIzvK78BuPqWAutrv9lf2ylHCB6J2yLh/17darEkVHuOBKqVDQXs5+rxZTu1fD2dGPwnF9p9FkkC3ec5LbdCcTFFQJ7sqnGl1C1O2z+CkZVhW3fwO1UuPtR6j48PT05f/68Jo8UJiKcP38eT8/kl0lyquq4IrIAWBAYGNg7rduSVMYYQssWIKRMfpbtPs3wnw8Q9t12yhY8xJBGZahfrsB9Sx2khFj3HNBwBFR9xirZPj8Mtk6Bpp+Cb1Xb4ip1P76+vkRHR/PHH3+k6HVv3LjxUD80M6q7++3p6Ymvr2+yr+VUicMZGGNoUr4QDQMKsnDHST5bfoBnp0XxeMm8vNUywP41IIUeg57LYMcsWP4WTKxnlW2v/zZ4Jb76XKmU5O7uTvHixVP8uhEREVSuXDnFr5vepWS/9VFVOuXqYmhdqQgrBgXzXutH2XPqEs0+X8MbP+7kz6sx9gY3Bh7rCGFRUCsMfp1hPb7aOA7iUmH2l1IqXXOqxJHRBscd4ebqQrdafkQMCaFbLT9mbD5OyKfhTFl31P4ZWJ45oPFQeH49FKkCS1+Fr4Lg2Fp74yql0jWnShzOLFc2D95p9ShLB9blsaK5eHfBHpp+vobw/WftHzzMXwa6/gAdvoWbl619z+f0gIsn7I2rlEqXNHFkMKV8vJnWszoTuwUSd1voMWUL3SZvZu+pS/YGNgbKtYT+myD4Ndi/GEYHwpoREHvT3thKqXTFqRKHM41xJMYYQ4MAH5a9GMRbLQLYeeIizb9Yw6tzd3D20g17g3tkg9DXrQRSsh6sfBe+rAUHfrY3rlIq3XCqxOGMYxyJ8XBzoWed4kQOCaVn7eLM2x5NyLAIPl9xkOsxcfYGz+0HHadDl++tu5HvnoTvOsKfR+yNq5RKc06VODKrnNnceaNFACsGBRNcOj+frThA/eERLPj1pP3jH/4N4PkN0PA9OLYGxtSEle9r8USlnJgmDidSLG92xnapypy+tcid3YMBM7bTYfxG9py0efzDzQNqD7Sm7wa0hjXDYHR12P2DFk9Uyglp4nBC1fzyMD+sDv9rW4GDZy7TYtQa3vxxFxfsXv+RoxA8MQF6LIWsuWFOd5jWCs7utTeuUipVOVXiyCyD445wdTE8XeMRIoaE0q2WH99t/p3Q4RFM23CMWLvXfxSrBc9FQrNhcGoHjK0NS1+HG/p9UcoZOFXiyGyD447Imc2dd1o9yuIX6lKuYA7e+mk3zb9Yy7pD5+wN7OIK1XvDgG1QpRtsHGutPt8+XYsnKpXBOVXiUPdXpqA33/WuwbguVbgaE0vniZt47psofj9/zd7A2fNCy5HQJxxyF4ef+sHkRnBim71xlVK20cSRifxdQHHFoGBeblyGNQfP0WBEJB8v3cf1WJsHsQtXtoonthkHF36DCfVg/gtw9by9cZVSKe6+icMY42KM6WmMWWSM+dUYs9UYM9MYE5KK7VM28HR3pX+oP+FDQmjxWCHGRhzmtTXXmbct2t79P1xcoFIna+vamv1g+7cwqjJsnqDFE5XKQBK745gEPAJ8CIQDi+KPvWGMGZAKbbvDGNPGGDPBGPOTMaZRasZ2Zj45PBnxVCV+6Pc4eTwNg2b/Svtx69l1wuZBbM+c0OR/VvHEQo/B4iEwPhh+W29vXKVUikgscVQVkXdEZK2IvAg0EpHlQHOgn6MBjDGTjTFnjTG7EhxvYozZb4w5ZIx5LbFriMiPItIb6A50cDS2ckzlR3LzZk1PPnmiIr+dv0bL0Wt5fV4qlG8vUBa6zYcnv7a2q53SFL5/Fi6dtDeuUuqhJJY4bhljSgIYY6oAMQAichNIyvOMqUCTuw8YY1yBMUBTIADoZIwJMMZUMMYsTPAqcNepb8Sfp1KYizE8Va0oq4aE0OPx4syOOk7osFSYvmsMPNoGwrZA0CuwZz6MCoS1IyHW5sSllEqWxBLHy0C4MeYg8H38e4wx+YGFjgYQkdXAnwkOVwcOicgREYkBZgKtRWSniLRI8DprLB8DS0REp+PYKGdWd95qGcCSgXV5tLA1fbfFqLVsOGzzILZHNqj3X6t4YolgWPE2jK0FB1fYG1cplWTmfrWMjDGFgNNAXhF5qEn/xhg/YKGIlI9/3x5oIiLPxr/vCtQQkbD7nP8C8AywBfhFRMbd4zN9gD4APj4+VWfOnJmstl65cgUvL69knZuR3avfIkLUmThm7ovh/A2hWkFXOpbxIG9W+yfj5Tm/Ff9DE8l2/STn8tbgkH8vbmT1SfE4+v3OfDJr3x3pd2ho6FYRCXzgxUTkni9gCbAR+AgIAdzu99kHvQA/YNdd758EJt71viswKrnXv+s6LYHx/v7+klzh4eHJPjcjS6zf12NiZeTyA1L6v4ulzBuLZeTyA3I9Jtb+Rt26IbJmhMgHhUTeyy+yaqjIzaspGkK/35lPZu27I/0GosSBn7X3/dVRRJrGJ4wIoC2w0RgzzxjTxxjzyAMzUuKigaJ3vfcFHnpEVHTluC083V0Z2KAUKwcHU7+sT3z13UiW7jplb/VdtyxQ5yVr+m65lhD5MYypbo2DaPFEpdJMos8cROSGiCwVkYFi3b4MBtyA0caYzQ8RdwtQyhhT3BjjAXQE5j/E9QCtVWU339zZGNO5CjN618Qrixt9v91G10mbOXT2sr2BcxSG9pOg+yLIkgNmd4Vv2sAf++2Nq5S6pyQ9rBaRoyLypYi0Auo4co4xZgawAShjjIk2xvQSkVggDFgG7AVmi8juJLb9Xu3TO45UUKtkXha9UId3WgawI/ovmoxcw9BFe7h845a9gf3qwHOroemncHI7jH0clv0XbthcNl4p9Q/JGuU0xuwUazbUA4lIJxEpJCLuIuIrIpPijy8WkdIiUlJEhianHfdol95xpBI3Vxe61y5O+JAQ2lf1ZeLao9QbHsn3W21efe7qBjX6WMUTKz0NG8ZYe5//MkOLJyqVShIrOdLuPq8ngIKp2EaH6R1H6svrlYWPnqjIj/1qUzhXVgbPsVaf74y2OXlnzwetRkHvlZCzKPzYF6Y0gZO/2BtXKZXoHccsoBXWTKW7Xy0AT/ublnR6x5F2Hiuaix+ef5xP2lfk9z+v0WrMWl6ft4PzV27aG7hIVei1HFqPsfY7Hx8CC16EawmXDimlUopbIl/bAQwTkV0Jv2CMaWBfk5JPRBYACwIDA3undVsyIxcXw1OBRWlSviCfrzjI1PXHWLTjFIMblaFzjUdwc7Vp/YeLC1TuYs28ivgINn0Fe36Eem9A1R7W3iBKqRST2L/kF4H7jTq2taEtD03vONKHHJ7uvNkigKUD61LBNydvz7dWn288YvPqc8+c0ORD6LsWfMrDosFW8cTfN9obV6lMJrF1HGtE5Pf7fC3KviYln45xpC+lfLz5tpe1edTlG7F0HL+RsO+2ceridXsD+wTAMwug/RTrkdXkxjCvD1w+bW9cpTIJ3chJ2eruzaMG1i/F8j1nqDcskjHhh7gZG2dnYCjfziqeWHcw7P7B2rp23RdaPFGph+RUiUMfVaVfWT1cealhaVYMCqZuqXx8umw/jT5bzap9Z+wN7JEd6r8F/TZCsdqw/E1r/cehlfbGVcqJOVXi0EdV6V/RPNkY3y2QaT2r4+pi6Dk1ip5Tt3D03FV7A+ctCZ1nw9Oz4XYsfNsOZnbG87rNiUspJ5TYrCoAjDG5gG5YhQrvfF5EXrCvWcrZBZXOz9KBQUxdf5TPVxyk8WerebZuccLq+ZPN44H/WyZf6cZQPBg2jIY1w6kW9zN4HoE6L4J7VvviKuVEHLnjWIyVNHYCW+96pTv6qCpj8XBzoU9QSWvv84qF+DLiMPWGRTL/15P2Fk9094SgIRC2hfN5q0PkR1bxxL0LtHiiUg5wJHF4isggEZkiIl///bK9Zcmgj6oypgI5PBnRoRJz+9Yir5cHL8zYTsfxG9l32uYaVDl92fPoy9YMLA8vmNXFeoT1xwF74yqVwTmSOL4xxvQ2xhQyxuT5+2V7y1SmE+iXh/lhdRjatjwHzlym+RdreWf+bi5et7l4YvEgeG4NNPkYordaOw/+/CbctLnqr1IZlCOJIwb4FKvC7d+PqdLlOg6V8bm6GDrXKEb4kBCerv4I0zYco96wCGZt+d3+4ok1+8KArfBYR1j/hbX3+Y7Z+vhKqQQcSRyDAH8R8ROR4vGvEnY3TGVuubJ58H6b8iwYUIcS+bPz6vc7afvlOn45/pe9gb3yW3Wvnl1p7QMyrzdMaQqndtgbV6kMxJHEsRu4ZndDlLqXRwvnZPZztRjZoRKnLt6gzZh1vDL3V87ZXTzRN9BKHq1GwbkDVumSRYO1eKJSODAdF4gDfjHGhAN3/rWmx+m4xpiWQEt/f/+0bopKQcYY2lQuQoMAH0atPMjkdUdZsus0LzUoTbdaxewtnlilm1U8MfxD2DIBds2D+m9ClWe0eKLKtBz5F/cjMBRYTzqfjquzqpybVxY3Xm9WjqUvBlGpaC7eW7iHZl+sYf3hc/YGzpobmn1iDaAXKAcLX4IJoXD8YXZPVirjemDiiJ96O4P/TxjfpdfpuCpzKJnfi2k9q/NV16pci4nj6Qmb6P/dNk7+ZXPxxILlrX3Pn5gEV87CpIbww/NwWVefq8zlgYnDGBMCHATGAF8CB4wxQTa3S6lEGWNo/GhBVgwK5qUGpVmx5wz1h0cyetVBbtyyuXhihfYQFgV1XoKdc6ziietHQ5zN04aVSicceVQ1HGgkIsEiEgQ0Bj6zt1lKOcbT3ZWBDUqxYlAwwaXzM+znAzQeuZqVe22+C8jiBQ3esYonPlITfv4vjK0Nh8PtjatUOuBI4nAXkf1/vxGRA4C7fU36J2NMOWPMOGPMXGPM86kVV2UsRfNkY1zXqnzTqzpuLoZeX0fRY8pm+4sn5vOHznOg00yIuwnftIFZXeGve25lo5RTcCRxRBljJhljQuJfE3BwcNwYM9kYc9YYsyvB8SbGmP3GmEPGmNcSu4aI7BWRvsBTQKAjcVXmVbdUfpYMDOK/zcqx5dgFGn+2mk+W7uNaTKx9QY2BMk2h3yZru9qDy2F0dYj8FG7dsC+uUmnEkcTxPNZajheAgcAeoK+D158KNLn7gDHGFWu8pCkQAHQyxgQYYyoYYxYmeBWIP6cVsBbQTRTUA3m4udA7qASrBgenQfHEl63No0o3hvAP4MsasG+xrj5XTsUk9g8p/of81yLSJdkBjPEDFopI+fj3tYB3RKRx/PvXAUTkQweutUhEmt/na32APgA+Pj5VZ86cmaz2XrlyBS8vr2Sdm5E5c78PXojj270x/HbpNmVyu9AlIAtFva3fmezsd64LOyh1cDzZrx3nfJ6qHPLvxfVsRWyJlVTO/P1+kMzad0f6HRoaulVEHvxkR0QSfQHLAI8HfS6R8/2AXXe9bw9MvOt9V2B0IueHAF8AXwH9HYlZtWpVSa7w8PBkn5uROXu/Y+Nuy7cbj8lj7y6TEq8vkrd/2iV/XYuxv9+xMSLrx4j8z1fk3bwiP78lcuOyvTEd4Ozf78Rk1r470m8gShz4GevIyvHoGCJ9AAAgAElEQVRjwDpjzHzgzkijiIxw4Nx7Mfc4dt/bHhGJACIcurCuHFf38XfxxGblCzF8+X6mbTjG/F9P0ro4BN0WXFzu9b9lSgR2h1r9rCm8K96BdSOtwomN3ofyT1jjI0plMI6McZwEFsZ/1vuuV3JFA0Xveu8bH0Mp2+XO7sEHbSpYxRPzZWfKrphUKp5YANp8Cb1WWH/+vhdMbQ6ndz34XKXSmfsmDmPMN/F//EtE3k34eoiYW4BSxpjixhgPoCMw/yGud4doyRHloEcL52RO31r0qZgldYsnFq0GvVdBy8/h7F74qi4sfhmuX7A3rlIpKLE7jqrGmGJAT2NM7rs3cXJ0IydjzAysfTzKGGOijTG9RCQWCMMaO9kLzBaR3Q/bkfh4unWscpgxhscLu7FqSAjPBZVg3rYThA6LYMq6o8TG3bYvsIsrVO1u7f0R2Au2TLRWn2/9Gm7bGFepFJJY4hgHLAXK8s/ihg5v5CQinUSkkIi4i4iviEyKP75YREqLSEkRGfpwXfhHPL3jUEmWsHjiuwv20PyLtWw4fN7ewNnyQPNh8NxqyFcaFrwAE+tBtO6TptK3+yYOEflCRMoBk0WkhPz/Jk7pdiMnveNQD8O/gFU8cVyXqly5GUunCRsJ+24bpy7aXTyxAvRYAu0mwqVTMLE+/NjfKqSoVDrkSHXcDFPmQ+841MMyxtCkfEFWDg5mYP1SLN9zhnrDIhkTfoibsTYXT6z4JAyIgsdfgB2zrMdXG8dq8USV7ti0A07a0DsOlVI83V15qWFpVgwKJqh0Pj5dtp9Gn61m1T67iyd6W1N1+22wdiFc+hqMqwtHV9sbV6kkcKrEoXccKqUVzZONr7oGMq1ndVxdDD2nRtFz6haO2V48sRR0mQcdpsOtq/B1S5j9DFyMtjeuUg5wqsShlF2CSudn6cAgXm9alk1HztPos9V8uiwViieWawH9N0PIf+DAUhhdDVZr8USVthzZyKmdMeagMeaiMeaSMeayMeZSajQuqfRRlbKTh5sLzwWXZNWQEJpXLMSY8MPUHx7JAtuLJ2aFkFetBOJfH1Z9AF/WhP1L7YupVCIcueP4BGglIjlFJIeIeItIDrsblhz6qEqlBp8cnnzWoRJz+tYidzYPBszYTqcJG9l/+rK9gXMXgw7fQtcfrFImMzrA9Cfh/GF74yqVgCOJ44yI7LW9JUplMNX88rBgQB3eb1Oefacv0+yLNby7YDcXr9s8C6pkPei7Dhp9AL9tsO4+Vr4HMTaPuygVz9GNnGYZYzrFP7ZqZ4xpZ3vLlMoAXF0MXWsWI3xwCB2rFWXq+mPUHx7B7Kjj3L5t4+MrNw94fIA1fffRdrBmuDX+set73ftD2c6RxJEDuAY0AlrGv1rY2ajk0jEOlVZyZ/dgaNsKLAirQ7G82Xll7g7ajV3Pjmibiyd6F4R2X0HPZdZK9Lk9rRlYZ/bYG1dlao4sAOxxj1fP1GhcUukYh0pr5YvkZG7fWox46jGiL1yn9Zh1vPb9Ds7bXTzxkZrQJxKaj4Azu2BcHVjyGly3OXGpTMmRWVW+xpgf4vcOP2OM+d4Y45sajVMqIzLG0K6KL+FDgnm2TnHmbo0mdFgEX68/Zn/xxGq9YMA2qPoMbBpnrT7f9o0WT1QpypFHVVOwyp4XBooAC+KPKaUS4e3pzn+bB7D0xbpU9M3F2/N302LUWjYdSYXiiS0+g+ciIW9JmB8GkxrAia32xlWZhiOJI7+ITBGR2PjXVCC/ze1Symn4F/Dmm17VGdelCpdvxNJh/EZemLGd0xdtXsRX6DFr7KPtV9aK8wn14acwuPKHvXGV03MkcZwzxnQxxrjGv7oANv/KlDw6OK7SK6t4YiFWDArmhfqlWLr7NPWGRzA24rD9xRMf6whhUfB4GPw6A0ZVpUj0AoizcdW7cmqOJI6ewFPAaeAU0D7+WLqjg+Mqvcvq4cqghqVZOSiY2v75+HjpPpqMXEPEfptLqHvmsNZ9PL8eilSh1KGJ8FUQHFtrb1zllByZVfW7iLQSkfwiUkBE2ojIb6nROKWcVdE82ZjQLZCve1bHAN2nbOHZr6P4/fw1ewPnLwNdf2DXo6/BzcvWvudzesDFE/bGVU4lsT3HX4n/7yhjzBcJX6nXRKWcV3Dp/Cx9MYjXmpZl/eFzNPgskuE/7+d6jL2Pr87lrwX9N0Hwa7B/MYwOtBYRxto8bVg5hcTuOP4uMxLFv7eO1ekZSqUQDzcX+gaXZNXgEJqWL8ioVYdoMCKSxTtP2Vs80SMbhL5uJZCS9ayyJV/WhAM/2xdTOYXEto5dEP/HayLy9d0vrJXkqcYYk90Ys9UYky5XrCuVEgrm9OTzjpWZ1acm3p5u9Ju+jS6TNnHwjN3FE/2g43Rr/w/jCt89Cd91hD+P2BtXZViODI6/7uCxfzHGTI5fOLgrwfEmxpj9xphDxpjXHLjUq8BsR2IqldHVKJGXhQPq8F7rR9kZfZGmn6/h/YV7uHTD5uKJ/vWtwfOG78GxNTCmplXCPSZVf09UGYDb/b5gjGkKNAOKJBjTyAE4Oo9vKjAamHbXdV2BMUBDIBrYYoyZD7gCHyY4vydQEdgDeDoYU6kMz83VhW61/GheoRDDft7P5HVH+emXk7zWtCztKhfBxcXYFNgDag+ECk/B8resTaN+mQGNh0JAa2t6r8r0ErvjOIk1vnGDf45tzAcaO3JxEVkN/JngcHXgkIgcEZEYYCbQWkR2ikiLBK+zQChQE3ga6G2M0V0LVaaR1ysLH7aryE/9a+ObOytD5vzKE+PWszPa5rVKOQrBExOgx1LImhvmPAPTWsHZffbGVRmCedDgmzHGXUSSfY9sjPEDFopI+fj37YEmIvJs/PuuQA0RCXvAdboD50Rk4X2+3gfoA+Dj41N15syZyWrvlStX8PLySta5GZn2O/27LcK6E7HMORDD5RgI9nXjidIeeHsk/S4gSf2WOAqfXEbxo9NxjbvOiSLNOebXkTi37EmOmx5kpO95SnKk36GhoVtFJPBB17rvo6q7+BljPgQCuOtxkYiUcODce7nX/+UPnDoSX+oksa+PN8acAlp6e3tXDQkJSVbjIiIiSO65GZn2O2OoBwy8cYuRyw/y9YZjbD9/iyGNSvN0jWK4JuHxVdL7XR+uvgqr3qPo1q8pemEjNHwXKnYEl4z1ECCjfc9TSkr229Eih2OxxjVCscYrvnmImNFA0bve+2I9FntounJcZQY5PN15q2UASwbWJaBQDt78ySqeuOVYwqfCKSx7Xmj5OfQJt7ax/fF5mNwITm63N65KdxxJHFlFZCXWY63fROQdrF98kmsLUMoYU9wY4wF0xBo3eWhaq0plJqV9vPmudw3GPF2Fi9dieHLcBl6cuZ0zl2wunli4MvT8GdqMhQvHYHwoLBgIV9NlCTtlA0cSx434AemDxpgwY0xboIAjFzfGzAA2AGWMMdHGmF4iEguEAcuwFhnOFpHdyWy/UpmaMYbmFQuxYnAwYaH+LN55mnrDIvgq8jAxsXbu/eEClZ6GAVuhZj9rz49RVWDzBC2emAk4kjheBLIBLwBVga7AM45cXEQ6iUghEXEXEV8RmRR/fLGIlBaRkiIyNLmNv0c8fVSlMqVsHm4MaVyG5YOCqFUyLx8u2UeTz1cTecDmEuqeOaHJ/+D5dVCoIiweAuOD4bf19sZVacqRIodbROSKiETHbxvbTkQ2pkbjkkofVanMrlje7Ex8phpTelRDBJ6ZvJk+06I4/qfNi/gKlINu8+HJr63taqc0he+fhUspMnyp0pnEihyOjP/vAmPM/ISv1Gui4/SOQylLaJkCLH2xLq80KcOag+doMCKSEcsP2F48kUfbQNhmCHoZ9syHUYGwdiTExtgXV6W6xKbj/j1zalhqNEQplbKyuLnSL8SftpWL8L/F+/hi5UG+3xrNmy3KkcXW4onZod4b1hjI0v/Airdh+zfQ9GPwb2BfXJVqEity+HcF3DzARhGJvPuVOs1LGn1UpdS/FcqZlVGdKjOjd028srjR99ttDIu6waGzNhdPzFMCnp4JneeCCHz7BMx42pqJpTI0RwbHWwEHjDHfGGOaG2McWTSYJvRRlVL3V6tkXha9UId3WgZw5OJtmoxcw9BFe7hsd/HEUg2h3wao/zYciYDR1SH8f1o8MQNzZHC8B+APzMGqF3XYGDPR7oYppVKem6sL3WsX5+O62Xiiii8T1x6l3vBI5m2LtnfvD7csUHcQhG2Bci0g8mMYU90aB7EzrrKFQ7UC4mtVLcEqSLgVaG1no5JLH1Up5ZgcWQwft6/ID/1qUzinJ4Nm/0r7cRvYdcLmfzs5i0D7ydB9EWTJAbO7wjdt4I/99sZVKeqBiSN+74ypwCGgPTARKGRzu5JFH1UplTSViubih361+eSJihw7d5WWo9fy3x92cuGqzbOg/OrAc6uh6SdWyZKxj8Oy/8KNS/bGVSnCkTuO7sCPQGkReSZ+8Z4uDVXKSbi4GJ6qVpRVQ0J4ppYfM7ccJ3R4BN9u/I242zY+RnJ1gxrPwYBt1gysDWOsvc9/mQG3bVz1rh6aI2McHYHtQF0AY0xWY4y33Q1LDn1UpVTy5czqzjutHmXRC3UoW9CbN37cRavRa9n6m93FE/NBq1HQeyXk9IUf+8KUJnDyF3vjqmRz5FFVb2Au8FX8IV+sO5B0Rx9VKfXwyhbMwYzeNRnVqTJ/Xo3hibEbGDTrF85etrl4YpGq0GsFtB4D5w/D+BBY8CJcszlxqSRz5FFVf6A2cAlARA7iYJFDpVTGZIyh5WOFWTk4mP6hJVm44xT1hkUycc0RbsXZXDyxchereGKNvrBtmlU8cctEuG3jqneVJI4kjpvxW7wCEL+OQ+fPKZUJZPNw4+XGZfn5pSCqF8/DB4v20vTzNaw9eM7ewFlzQdOPoO9a8CkPiwZbxRN/T5dl8jIdRxJHpDHmP0BWY0xDrPUcC+xtllIqPfHLl53J3asx6ZlAbsXdpsukTTz/7VaiL9i8iM8nAJ5ZAO2nWI+sJjeGeX3g8ml746pEOZI4XgP+AHYCzwGLgTfsbFRy6eC4UvaqX86HZS8GMaRRacL3n6XBiEg+X3GQG7dsLp5Yvp21eLDuYNj9A4yqCuu+0OKJacSRWVW3sQbD+4lIexGZILYuMU0+HRxXyn6e7q6E1SvFysEh1C/rw2crDtDws0h+3n3a3tXnHtmh/lvQbyMUqw3L34RxteHwKvtiqntKrKy6Mca8Y4w5B+wD9htj/jDGvJV6zVNKpVdFcmVlTOcqfPdsDTzdXOnzzVaembKFw39csTdw3pLQeTY8PRvibsE3bWFWF7jwm71x1R2J3XG8iDWbqpqI5BWRPEANoLYx5qVUaZ1SKt173D8fiwfW5c0WAWz/7QJNRq7mwyV7uXLT5nXCpRtbdx/13oRDK63aVxEfwa3r9sZViSaObkAnETn69wEROQJ0if+aUkoB4O7qQq86xVk1JIQ2lYrwVeQR6g2L4MftJ+x9fOXuCUFDrPGPMs0g4kMrgexdqMUTbZRY4nAXkX/NuRORPwB3+5r0T8aYEGPMGmPMOGNMSGrFVUolXX7vLHz65GPM6/c4BXN68uKsX+jw1Ub2nLS5BlVOX3hyijUDy8MLZnW29v84d9DeuJlUYokjsekKDk1lMMZMNsacNcbsSnC8iTFmvzHmkDHmtQdcRoArgCcQ7UhcpVTaqvJIbn7sV5uP2lXg0B9XaDFqDW/9tIu/rtk8C6p4EDy3Bpp8DNFR8GUtWP4W3LR506pMJrHE8Zgx5tI9XpeBCg5efyrQ5O4DxhhXYAzQFAgAOhljAowxFYwxCxO8CgBrRKQp8CrwblI7qJRKGy4uho7VHyF8cAhdaxbj242/ETosghmbf7e/eGLNvtbq88c6wLrPrb3Pd8zWx1cpJLGtY11FJMc9Xt4i4tCjKhFZDSQsNFMdOCQiR+JXpM8EWovIThFpkeB1Nn46MMAFIEsy+qiUSkM5s7nzbuvyLBxQl1IFvHl93k7ajFnHtt8v2BvYK79V9+rZlZCjMMzrDVOa4nX5iL1xMwFj95IMY4wfsFBEyse/bw80EZFn4993BWqISNh9zm8HNAZyAWNFJOI+n+sD9AHw8fGpOnPmzGS198qVK3h5eSXr3IxM+525pFW/RYRNp+KYuT+Gv24KdYq48WRpD3JmMTYHvk3B0yspcWQa7rcuc7JwE44W70yse7os9G0LR77noaGhW0Uk8EHXSovE8STQOEHiqC4iA1IqZmBgoERFRSXr3IiICEJCQlKqKRmG9jtzSet+X70Zy6hVh5i09giebq4MbFCKZx73w93VoU1Jk+/6BaK/7Y/vySXgmQvqvwlVngEXV3vjpgOOfM+NMQ4lDpu/S/cUDRS9670vcDIlLqwlR5TKGLJnceO1pmVZ9mIQVYrl5oNFe2n2+RrWH7K7eGJuDpXqYw2gFygHC1+CCaFwfLO9cZ1MWiSOLUApY0xxY4wH0BGYnwbtUEqlsRL5vZjaoxoTugVyIzaOpyduot/0rZz4y+ZFfAXLW/uePzEJrpyFSQ3hh75w+Yy9cZ2ErYnDGDMD2ACUMcZEG2N6xW87GwYsA/YCs0Vkd0rE01pVSmU8xhgaBviw/KVgBjUszap9Z6k/PIJRK1OheGKF9hAWBbVfhJ1zreKJ60dbpUzUfdmaOESkk4gUEhF3EfEVkUnxxxeLSGkRKSkiQ1Mqnj6qUirj8nR35YX6pVgxKJjQMgUYvvwAjT5bzYo9Z+xdfZ7FCxq+a5UveaQG/PxfGFsbjkTYFzODS4tHVbbROw6lMj7f3NkY26Uq3/aqgYebC89Oi6LH1C0cPXfV3sD5/KHzXOg4A+JuwrTWMKsr/PW7vXEzIKdKHHrHoZTzqFMqH0sG1uWN5uWIOnaBxp+t5uOl+7hqZ/FEY6BsM+i3CULfgIPLYXR1iPwEbtm853oG4lSJQ+84lHIu7q4uPFu3BKuGBNPisUKMjThM/eGRzP/1pP3FE4Nftoonlm4E4UPhyxqwb7GuPsfJEofecSjlnAp4ezLiqUp8/3wt8nl78MKM7XQcv5G9p2wunpirKDw1Dbr9BG6eMLMTTH8Szh2yN24651SJQ+84lHJuVYvl4af+dRjatjz7z1ym+RdreGf+bi5es3kWVIkQ6LsWGv8Pjm+CL2vC8rfhps2bVqVTTpU4lFLOz9XF0LlGMSKGhNC5RjGmbThG6PAIZm35ndu2Fk90h1r9rem7FZ6EdSNhdDVrGm8me3zlVIlDH1UplXnkyubB+23Ks2BAHUrky86r3++k7Zfr+OX4X/YG9vaBtmOh13KrkOL3vWBqczi968HnOgmnShz6qEqpzOfRwjmZ07cWIztU4tTFG7QZs45X5v7KuSs37Q1ctDr0DocWI+HsXviqLix+Ga7bXPU3HXCqxKGUypyMMbSpXIRVQ0J4LqgE87adIHRYBFPWHSU27vaDL5BcLq4Q2MPa+yOwJ2yZaK0+3/o13LYxbhrTxKGUchpeWdx4vVk5lr4YRKWiuXh3wR5ajFrLxiPn7Q2cLQ80Hw59IiFfaVjwAkysZ+1C6IScKnHoGIdSCsC/gBfTelbnq65VuXIzlo7jNxL23TZOXbS5eGKhitBjCbSbAJdOwcT68GN/q5CiE3GqxKFjHEqpvxljaPxoQVYMCubFBqVYvucM9YZFsvBwDDdjbS6eWPEpGBAFj78AO2ZZj682fOk0xROdKnEopVRCnu6uvNigNCsGBRNUOh9zD96i8WerWbXP5hLqWbyh0fvQbwP4BsKy12FcXTi62t64qUATh1IqUyiaJxtfdQ1kSKAnri6GnlOj6DV1C8dsL55YCrrMgw7T4dZV+LolzOkOF6PtjWsjTRxKqUylfD5XlgwM4j/NyrLxyHkafbaaT5ft41qMzcUTy7WA/psh5D+wf4m1eHD1MIi1edqwDZwqcejguFLKER5uLvQJKsmqISE0r1iIMeFW8cQFthdPzAohr1oJxL8+rHofxtSAA8vsi2kDp0ocOjiulEoKnxyefNahEnP61iJ3Ng8GzNhOpwkb2X/6sr2BcxeDDt9C1x+sUibfPQXfdYDzh+2Nm0KcKnEopVRyVPPLw4IBdXi/TXn2nrpMsy/W8O6C3Vy8bvMsqJL1oO86aPQBHFtnFU9c+R7E2Dzu8pA0cSilFFbxxK41ixE+JIQO1Yoydf0x6g2LYHbUcXuLJ7p5wOMDrOm7j7aDNcOt8Y9d89Jt8URNHEopdZc82T34X9sKLAirQ7G82Xhl7g7ajV3Pr7YXTywI7b6Cnsuslehze1gzsM7ssTduMqT7xGGMcTHGDDXGjDLGPJPW7VFKZQ7li+Rkbt/HGf7kY0RfuE6bL9fx2vc7OG938cRHalqlS5qPgDO7YFwdWPIaXLc5cSWBrYnDGDPZGHPWGLMrwfEmxpj9xphDxpjXHnCZ1kAR4BaQcSc+K6UyHBcXwxNVfQkfEkyv2sWZuzWa0GERfL3+mP3FE6v1ggHboOozsGmctfp82zfponii3XccU4Emdx8wxrgCY4CmQADQyRgTYIypYIxZmOBVACgDbBCRQcDzNrdXKaX+xdvTnTdaBLD0xbpU9M3F2/N302LUWjalRvHEFp9BnwjIWxLmh8GkBnBiq71xH8DWxCEiq4E/ExyuDhwSkSMiEgPMBFqLyE4RaZHgdRbrLuPvAvc2FphRSqnE+Rfw5pte1RnbuQqXb8TSYfxGXpixndMXb9gbuHAla+yj7Vfw13GYUB9+CoOr5+yNex/G1sUugDHGD1goIuXj37cHmojIs/HvuwI1RCTsPudnA0YB14B9IjLmPp/rA/QB8PHxqTpz5sxktffKlSt4eXkl69yMTPuduWTWfkPK9f1mnLD4yC0WHb2Fq4FWJd1p5OeOu4tJgVben2vsNfyOzaLIiQXEuXpyzO9pThZuiri4JnqeI/0ODQ3dKiKBD2qDW9KanCLu9bd63+wlIteAXg+6qIiMN8acAlp6e3tXDQkJSVbjIiIiSO65GZn2O3PJrP2GlO174/ow+Pw13l+0hzl7zhD1pwdvtQwgtEyBFLn+/TWDP/bjsuQVSh2aQKlL66HZJ+BX575npGS/02JWVTRQ9K73vsDJNGiHUko9tEfyZmNCt0Cm9qgGQI8pW3j26yh+P3/N3sD5y0DXH+Gpb+DmZWvf87k94ZL9P07TInFsAUoZY4obYzyAjsD8lLiwlhxRSqWVkDIFWPZiEK82Kcv6w+do8FkkI37ez/UYm/f+CGgF/TdB8KuwbzFcPmVfvHh2T8edAWwAyhhjoo0xvUQkFggDlgF7gdkisjuF4mmRQ6VUmvFwc+H5kJKsGhxC0/IF+WLVIRqMiGTxzlP2Fk/0yAah/4FBe6BIVfvixLN7VlUnESkkIu4i4isik+KPLxaR0iJSUkSGpmA8veNQSqW5gjk9+bxjZWb1qYm3pxv9pm+jy6RNHDxjc/HEbHnsvX68dL9yPCn0jkMplZ7UKJGXhQPq8F7rR9kZfZGmn6/h/YV7uHQjY28h61SJQ+84lFLpjZurC91q+RE+JIQnA32ZvO4o9YZF8v3WaHuLJ9rIqRKH3nEopdKrvF5Z+LBdRX7qXxvf3FkZPOdX2o9bz64TGe/nlVMlDr3jUEqldxV9czHv+cf5tH1Ffv/zGi1Hr+U/P+zkwtWYtG6aw5wqcSilVEbg4mJ4MrAoq4aE0OPx4szacpyQYRF8s/E34jLA4yunShz6qEoplZHk8HTnrZYBLBlYl4BCOXjzx120HLWWLccSlvhLX5wqceijKqVURlTax5vvetdgzNNV+OtaDE+O28BLs37h7CWbiycmk1MlDqWUyqiMMTSvWIgVg4MJC/Vn0Y5ThA6L4KvIw8TEpv0eHHfTxKGUUulINg83hjQuw/JBQdQqmZcPl+yjyeerWX3gj7Ru2h1OlTh0jEMp5SyK5c3OxGeqMaV7NW7fFrpN3kyfaVEc/9Pm4okOcKrEoWMcSilnE1q2AMteCuLlxmVYc/AcDUZEMmL5AXuLJz6AUyUOpZRyRlncXOkf6s/KwcE0DPDhi5UHaTAikqW7TttbPPE+NHEopVQGUThXVkY/XYUZvWvilcWNvt9updvkzRw6eyVV26GJQymlMphaJfOy6IU6vN0ygF+O/0WTkav53+K9XE6l4olOlTh0cFwplVm4ubrQo3ZxwoeE8EQVX8avPkK94ZFs/c3+xYNOlTh0cFwpldnk88rCx+0r8mP/2pQt6E2xvNltj+lmewSllFK2q1Q0F9/0qpEqsZzqjkMppZT9NHEopZRKknT/qMoYUxfojNXWABF5PI2bpJRSmZqtdxzGmMnGmLPGmF0Jjjcxxuw3xhwyxryW2DVEZI2I9AUWAl/b2V6llFIPZvcdx1RgNDDt7wPGGFdgDNAQiAa2GGPmA67AhwnO7ykiZ+P//DTwrM3tVUop9QC2Jg4RWW2M8UtwuDpwSESOABhjZgKtReRDoMW9rmOMeQS4KCKXbGyuUkopB6TFGEcR4Phd76OBB80h6wVMSewDxpg+QB8AHx8fIiIiktW4K1euJPvcjEz7nblk1n5D5u17SvY7LRKHucexRKt0icjbD7qoiIw3xpwCWnp7e1cNCQlJVuMiIiJI7rkZmfY7c8ms/YbM2/eU7HdaJI5ooOhd732BkylxYRFZACwwxrQ1xvyWzMvkA86lRHsyGO135pJZ+w2Zt++O9LuYIxdKi8SxBShljCkOnAA6Yg18pxgRyZ/cc40xUSISmJLtyQi035lLZu03ZN6+p2S/7Z6OOwPYAJQxxkQbY3qJSCwQBiwD9gKzRWS3ne1QSimVcuyeVdXpPscXA4vtjK2UUsoeWnLk38andQPSiPY7c8ms/YbM2/cU67dJi20HlVJKZVx6x6GUUipJNHHcJSk1tDKye9UQM8bkMcYsN//X3v28WFXGcTbzrp4AAAO9SURBVBx/f5iMooQpyQgt+kELJWrahGALGyKsJFsUFAXua2FQRLWJAqEIqj+gohZFSWWJq8SMooWRaWgYRCAtFAciyTZG9WnxPOMdhok4zb3ndM/9vGC45zxzuHy/zDP3e8557v1e6Yf6eEmXMY6CpCsl7Zd0TNJ3krbX8V7nLukCSV9J+rbm/Vwdv0bSgZr3e5LO7zrWUZA0JemQpD11v/d5Szou6Yikw5K+rmNDm+cpHNWCHlp3AuuBByWt7zaqkXkT2Lxo7Clgn+3rgX11v2/+AB63vQ7YADxa/8Z9z/0sMGv7JmAG2CxpA/Ai8ErN+xdKh4Y+2k55B+e8Scn7NtszC96CO7R5nsIxcK6Hlu3fgXeBrR3HNBK2PwcWfzHxVgbdh98C7m01qBbYPmn7m7p9hvJisoae5+7it7q7ov4YmAXer+O9yxtA0lrgbuC1ui8mIO9/MLR5nsIxsFQPrTUdxdKFy22fhPICC6zuOJ6Rqs03bwYOMAG519s1h4E5YC/wI3C6fq4K+jvfXwWeBP6q+6uYjLwNfCLpYO3jB0Oc5//7L3JqUeMeWjGeJF0MfAA8ZvvXchLab7b/BGYkTQO7gHVLHdZuVKMlaQswZ/ugpE3zw0sc2qu8q422T0haDeyV9P0wnzxXHAMj66E1Jk5JugKgPs79y/FjSdIKStF42/aHdXgicgewfRr4jLLGMy1p/uSxj/N9I3CPpOOUW8+zlCuQvueN7RP1cY5yonALQ5znKRwD53po1XdZPADs7jimNu0GttXtbcDHHcYyEvX+9uvAMdsvL/hVr3OXdFm90kDShcDtlPWd/cB99bDe5W37adtrbV9N+X/+1PZD9DxvSRdJWjm/DdwBHGWI8zwfAFxA0l2UM5Ip4A3bOzoOaSRqD7FNlG6Zp4BngY+AncBVwE/A/bYXL6CPNUm3Al8ARxjc836Gss7R29wl3UhZDJ2inCzutP28pGspZ+KXAoeAh22f7S7S0am3qp6wvaXvedf8dtXd84B3bO+QtIohzfMUjoiIaCS3qiIiopEUjoiIaCSFIyIiGknhiIiIRlI4IiKikRSOiJZImpb0SNdxRCxXCkdEe6aBFI4YeykcEe15AbiufkfCS10HE/Ff5QOAES2pHXn32L6h41AiliVXHBER0UgKR0RENJLCEdGeM8DKroOIWK4UjoiW2P4Z+FLS0SyOxzjL4nhERDSSK46IiGgkhSMiIhpJ4YiIiEZSOCIiopEUjoiIaCSFIyIiGknhiIiIRlI4IiKikb8B3HhEC0vR8JUAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "%matplotlib inline\n", "import matplotlib.pyplot as plt\n", "import numpy\n", "\n", "neighbor = [[1, 3, 0, 0], [2, 4, 0, 1], [2, 5, 1, 2],\n", " [4, 6, 3, 0], [5, 7, 3, 1], [5, 8, 4, 2],\n", " [7, 6, 6, 3], [8, 7, 6, 4], [8, 8, 7, 5]]\n", "transfer = numpy.zeros((9, 9))\n", "\n", "site_0=numpy.zeros(50)\n", "fit=numpy.zeros(50)\n", "for k in range(9):\n", " for neigh in range(4): transfer[neighbor[k][neigh], k] += 0.25\n", "position = numpy.zeros(9)\n", "position[8] = 1.0\n", "for t in range(50):\n", " print t,' ',[\"%0.5f\" % abs(i- 1.0 / 9.0) for i in position]\n", " position = numpy.dot(transfer, position)\n", " fit[t]=0.75**t\n", " site_0[t]=abs(position[0]- 1.0 / 9.0)\n", "\n", "plt.clf()\n", "plt.semilogy(site_0, label=\"$\\pi_0(t)-\\pi_0^{eq}$\")\n", "plt.semilogy(fit, label=\"$0.75^t$\")\n", "plt.ylabel('Deviation from 1/9')\n", "plt.xlabel('t')\n", "plt.grid(True)\n", "plt.legend(shadow=False, fancybox=True)\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Dual pebble game with a finite probability of switching between the games to ensure irreducability of the transfer matrix." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Populating the interactive namespace from numpy and matplotlib\n", "transition red->blue at time = 176\n", "transition blue->red at time = 197\n", "transition red->blue at time = 198\n", "transition blue->red at time = 199\n", "transition red->blue at time = 211\n", "transition blue->red at time = 213\n" ] } ], "source": [ "%pylab qt\n", "import random, pylab\n", "\n", "random.seed('1234')\n", "sigma = 0.4\n", "epsilon = 0.4 # probability to switch from red to blue pebble, and vice versa\n", "pylab.figure()\n", "s_map_red = [(1.0, 1.0), (2.0, 1.0), (3.0, 1.0), \n", " (1.0, 2.0), (2.0, 2.0), (3.0, 2.0), \n", " (1.0, 3.0), (2.0, 3.0), (3.0, 3.0)] \n", "offset = 3.0\n", "s_map_blue = [(x+offset,y-offset) for (x,y) in s_map_red]\n", "neighbor = [[1, 3, 0, 0], [2, 4, 0, 1], [2, 5, 1, 2],\n", " [4, 6, 3, 0], [5, 7, 3, 1], [5, 8, 4, 2],\n", " [7, 6, 6, 3], [8, 7, 6, 4], [8, 8, 7, 5]]\n", "color = 'red' #chose 'red' or 'blue'\n", "site = 8\n", "tmax = 240\n", "for iter in range(tmax):\n", " period = 4\n", " if (iter%period) == 0:\n", "\t# Begin of graphical output\n", " pylab.clf()\n", " maxlength = len(str(tmax-1))\n", " number_string = str(iter).zfill(maxlength)\n", " if color == 'red': cir = pylab.Circle(s_map_red[site], radius=sigma, fc='r')\n", " if color == 'blue': cir = pylab.Circle(s_map_blue[site], radius=sigma, fc='b')\n", "\tpylab.figure()\n", " pylab.gca().add_patch(cir)\n", " pylab.plot([0.5, 3.5], [0.5, 0.5], 'r')\n", " pylab.plot([0.5, 3.5], [1.5, 1.5], 'r')\n", " pylab.plot([0.5, 3.5], [2.5, 2.5], 'r')\n", " pylab.plot([1.5, 1.5], [0.5, 3.5], 'r')\n", " pylab.plot([2.5, 2.5], [0.5, 3.5], 'r')\n", " pylab.plot([3.5, 3.5], [0.5, 3.5], 'r')\n", " pylab.plot([0.5+offset, 3.5+offset], [1.5-offset, 1.5-offset], 'b')\n", " pylab.plot([0.5+offset, 3.5+offset], [2.5-offset, 2.5-offset], 'b')\n", " pylab.plot([0.5+offset, 3.5+offset], [3.5-offset, 3.5-offset], 'b')\n", " pylab.plot([0.5+offset, 0.5+offset], [0.5-offset, 3.5-offset], 'b')\n", " pylab.plot([1.5+offset, 1.5+offset], [0.5-offset, 3.5-offset], 'b')\n", " pylab.plot([2.5+offset, 2.5+offset], [0.5-offset, 3.5-offset], 'b')\n", " pylab.title('t = '+ number_string)\n", " pylab.axis('scaled')\n", " pylab.axis([0.5, 6.5, -2.5, 3.5])\n", " pylab.xticks([])\n", " pylab.yticks([])\n", " number_string_filename = str(iter/period).zfill(3)\n", " pylab.pause(0.05)\n", " pylab.close()\n", "\t# End of graphical output\n", " newsite = neighbor[site][ random.randint(0, 3)]\n", " newcolor = color\n", " if (color == 'red') and (site == 2) and (newsite == 2):\n", " if random.random() < epsilon:\n", " newcolor = 'blue'\n", " newsite = 6\n", " print \"transition red->blue at time = \", iter\n", " if (color == 'blue') and (site == 6) and (newsite == 6):\n", " if random.random() < epsilon:\n", " newcolor = 'red'\n", " newsite = 2\n", " print \"transition blue->red at time = \", iter\n", " site = newsite\n", " color = newcolor" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Recurrent pebble game with $2\\times 2$ grid with small aperiodicity." ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Populating the interactive namespace from numpy and matplotlib\n" ] } ], "source": [ "%pylab qt\n", "import math, random, pylab\n", "\n", "sigma = 0.4\n", "epsilon = 0.1\n", "pylab.figure()\n", "s_map = [(1.0, 1.0), (2.0, 1.0)] \n", "neighbor = [[1], [0]]\n", "pos = 0\n", "tmax = 20\n", "for iter in range(tmax):\n", " # Begin of the graphics output\n", " pylab.figure()\n", " number_string = str(iter).zfill(len(str(tmax)))\n", " cir = pylab.Circle(s_map[pos], radius=sigma, fc='r')\n", " pylab.gca().add_patch(cir)\n", " pylab.plot([1.5, 1.5], [0.5, 1.5], 'b')\n", " pylab.title('t = '+ number_string)\n", " pylab.axis('scaled')\n", " pylab.axis([0.5, 2.5, 0.5, 1.5])\n", " pylab.xticks([])\n", " pylab.yticks([])\n", " #pylab.savefig('2x1pebble_epsilon'+number_string+'.png', transparent=True)\n", " pylab.pause(0.05)\n", " pylab.clf()\n", " pylab.close()\n", " # End of the graphics output\n", " newpos = neighbor[pos][0]\n", " if random.random() < epsilon:\n", " newpos = pos\n", " pos = newpos" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Basic pebble throw game with inhomogeneous steady state probabilities." ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "comparison: weight, histogram\n", "site: 0 weight: 3.0 histo: 3.00848\n", "site: 1 weight: 0.5 histo: 0.49825\n", "site: 2 weight: 1.0 histo: 0.99388\n", "site: 3 weight: 0.5 histo: 0.50321\n", "site: 4 weight: 1.0 histo: 1.00448\n", "site: 5 weight: 0.5 histo: 0.49447\n", "site: 6 weight: 2.0 histo: 2.01049\n", "site: 7 weight: 0.5 histo: 0.50085\n", "site: 8 weight: 1.0 histo: 0.98589\n" ] } ], "source": [ "import random\n", "\n", "#define the grid in terms of transitions:[right,up,left,down]\n", "histo = [0, 0, 0, 0, 0, 0, 0, 0, 0]\n", "neighbor = [[1, 3, 0, 0], [2, 4, 0, 1], [2, 5, 1, 2],\n", " [4, 6, 3, 0], [5, 7, 3, 1], [5, 8, 4, 2],\n", " [7, 6, 6, 3], [8, 7, 6, 4], [8, 8, 7, 5]]\n", "weight = [3.0, 0.5, 1.0, 0.5, 1.0, 0.5, 2.0, 0.5, 1.0]\n", "pos = 8\n", "n_iter = 1000000\n", "for iter in range(n_iter):\n", " new_pos = neighbor[pos][random.randint(0, 3)]\n", " if random.random() < weight[new_pos] / weight[pos]:\n", " pos = new_pos\n", " histo[pos] += 1 \n", "\n", "norm = sum(weight)\n", "print 'comparison: weight, histogram'\n", "for k in range(9): \n", " print 'site: ', k,' weight: ', weight[k], ' histo: ', norm * histo[k] / float(n_iter)" ] } ], "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 }