HW_1.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "nbgrader": {
     "grade": false,
     "locked": false,
     "solution": false
    }
   },
   "source": [
    "<font size=\"4\" style=\"color:red;\"> **IMPORTANT: ** Only modify cells which have the following comment</font>\n",
    "\n",
    "```python\n",
    "# modify this cell\n",
    "```\n",
    "\n",
    "<font style=\"color:red;\"> Do not add any new cells when submitting homework. For Docker users, to test out new code, use the coding **scratchpad** by clicking the triangular icon in the bottom right corner of the screen. (**hotkey:** control-B)  </font>\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbgrader": {
     "grade": false,
     "locked": false,
     "solution": false
    }
   },
   "source": [
    "# Exercises:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbgrader": {
     "grade": false,
     "locked": false,
     "solution": false
    }
   },
   "source": [
    "**Note: ** Make sure you have read the *What is Probability?* notebook before attempting these exercises."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbgrader": {
     "grade": false,
     "locked": false,
     "solution": false
    }
   },
   "source": [
    "In this excercise you will write code to estimate the probability that $n$ flips of a fair coin will result in number of `\"heads\"` between $k_1$ and $k_2$.\n",
    "\n",
    "You should write the body of two functions:\n",
    "1. **`seq_sum`**: generates a random sequence of coin flips and counts the number of heads.\n",
    "2. **`estimate_prob`**: Using calls to `seq_sum`, estimate the probability of the number of heads being between $k_1$ and $k_2$. \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbgrader": {
     "grade": false,
     "locked": false,
     "solution": false
    }
   },
   "source": [
    "### Notebook Setup:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbgrader": {
     "grade": false,
     "locked": false,
     "solution": false
    }
   },
   "source": [
    "The folowing magic command downloads many python packages like *numpy* and allows the notebooks to plot graphs with *matplotlib*. \n",
    "\n",
    "<font color=\"red\">**DO NOT**</font> import other packages. You already have all the packages you need.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Populating the interactive namespace from numpy and matplotlib\n"
     ]
    }
   ],
   "source": [
    "%pylab inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Specifically, you can now use `random.rand(x)` which for some $x \\in N$ generates $x$ random numbers. You **will** use this command in your homework."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.2834288796977271"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "random.rand()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.4424793 , 0.1028654 , 0.27975932, 0.86693149])"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "random.rand(4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##  Exercise 1:\n",
    "\n",
    "Write a function, **seq_sum**, which generates $n$ random coin flips from a fair coin and then returns the number of heads. A fair coin is defined to be a coin where $P($heads$)=\\frac{1}{2}$ \n",
    "\n",
    "The output type should be a numpy integer, **hint:** use `random.rand()` \n",
    "\n",
    "<font  style=\"color:blue\"> * **Code:** *</font>\n",
    "```python\n",
    "x = seq_sum(100)\n",
    "print x\n",
    "print [seq_sum(2) for x in range(20)]\n",
    "```\n",
    "\n",
    "\n",
    "<font  style=\"color:magenta\"> * **Output:** *</font>\n",
    "```\n",
    "49\n",
    "[0, 1, 1, 1, 1, 2, 1, 2, 1, 1, 0, 0, 2, 1, 1, 1, 0, 0, 1, 1]\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* Write your code for seq_sum in the cell below"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# modify this cell\n",
    "import random\n",
    "\n",
    "def seq_sum(n):\n",
    "    \"\"\" input: n, generate a sequence of n random coin flips\n",
    "        output: return the number of heads \n",
    "        Hint: For simplicity, use 1,0 to represent head,tails\n",
    "    \"\"\"\n",
    "    total_heads = 0\n",
    "    total_tails = 0\n",
    "    count = 0\n",
    "    while count < n:\n",
    "        coin = random.randint(0, 1)\n",
    "        if coin == 1:\n",
    "            total_heads += 1\n",
    "            count += 1\n",
    "        \n",
    "        elif coin == 0:\n",
    "            total_tails += 1\n",
    "            count += 1\n",
    "    \n",
    "    return total_heads\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* if the following cell runs without error you receive some points."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "nbgrader": {
     "grade": true,
     "grade_id": "ex1",
     "locked": true,
     "points": "5",
     "solution": false
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "51\n"
     ]
    }
   ],
   "source": [
    "# checking function \n",
    "\n",
    "x = seq_sum(100)\n",
    "print(x)\n",
    "assert unique([seq_sum(2) for x in  range(0,200)]).tolist() == [0, 1, 2]\n",
    "\n",
    "#\n",
    "# AUTOGRADER TEST - DO NOT REMOVE\n",
    "#\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 2:\n",
    "\n",
    "Write a function, **estimate_prob**, that uses **seq_sum** to estimate the following probability:\n",
    "\n",
    "$$ P(\\; k_1 <= \\text{number of heads in $n$ flips} < k_2 ) $$\n",
    "\n",
    "The function should estimate the probability by running $m$ different trials of **`seq_sum(n)`**, probably using a *`for`* loop.\n",
    "\n",
    "In order to receive full credit **estimate_prob** <font color=\"red\">MUST</font> call **seq_sum** (aka: seq_sum is located inside the **estimate_prob** function)\n",
    "\n",
    "<font  style=\"color:blue\"> * **Code:** *</font>\n",
    "```python\n",
    "x = estimate_prob(100,45,55,1000)\n",
    "print x\n",
    "print type(x)\n",
    "```\n",
    "\n",
    "<font  style=\"color:magenta\"> * **Output:** *</font>\n",
    "```\n",
    "0.686\n",
    "<type 'float'>\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Modify this cell\n",
    "\n",
    "def estimate_prob(n,k1,k2,m):\n",
    "    \"\"\"Estimate the probability that n flips of a fair coin result in k1 to k2 heads\n",
    "         n: the number of coin flips (length of the sequence)\n",
    "         k1,k2: the trial is successful if the number of heads is \n",
    "                between k1 and k2-1\n",
    "         m: the number of trials (number of sequences of length n)\n",
    "         \n",
    "         output: the estimated probability \n",
    "         \"\"\"\n",
    "    successful_trial = 0\n",
    "    for i in range(m):\n",
    "        x = seq_sum(n)\n",
    "        if x>k1 and x<k2:\n",
    "            successful_trial+=1\n",
    "        elif x<k1 or x>k2:\n",
    "            successful_trial+=0\n",
    "    \n",
    "    #print(successful_trial)\n",
    "    estimated_probability = successful_trial/m\n",
    "    return estimated_probability\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.946\n"
     ]
    }
   ],
   "source": [
    "# this is a small sanity check\n",
    "# the true check for this function is further down\n",
    "\n",
    "x = estimate_prob(100,40,60,1000)\n",
    "print(x)\n",
    "assert 'float' in str(type(x))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Estimate vs. True Probability\n",
    "\n",
    "We can now check how to see how close these estimates are to the true probabilities."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Helper Functions\n",
    "\n",
    "These helper functions are used to calculate the actual probabilities. They are used to test your code.\n",
    "\n",
    "It is not required that you understand how they work."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": true,
    "nbgrader": {
     "grade": false,
     "locked": true,
     "solution": false
    }
   },
   "outputs": [],
   "source": [
    "def calc_prob(n,k1,k2):\n",
    "    \"\"\"Calculate the probability using a normal approximation\"\"\"\n",
    "    n=float(n);k1=float(k1);k2=float(k2)\n",
    "    z1=(k1-0.5*n)/(sqrt(n)/2)\n",
    "    z2=(k2-0.5*n)/(sqrt(n)/2)\n",
    "    return (erf(z2/sqrt(2))-erf(z1/sqrt(2)))/2\n",
    "\n",
    "from math import erf,sqrt\n",
    "def evaluate(n,q1,q2,m,r=100):\n",
    "    \"\"\"Run calc_range many times and test whether the estimates are consistent with calc_prob\"\"\"\n",
    "    k1=int(q1*n)\n",
    "    k2=int(q2*n)\n",
    "    p=calc_prob(n,k1,k2)\n",
    "    std=sqrt(p*(1-p)/m)\n",
    "    print(\"computed prob=%5.3f, std=%5.3f\" %(p,std))\n",
    "\n",
    "    L=[estimate_prob(n,k1,k2,m) for i in range(r)]\n",
    "    med=np.median(L)\n",
    "    print(\"ran estimator %d times, with parameters n=%d,k1=%d,k2=%d,m=%d\" %(r,n,k1,k2,m))\n",
    "    print('median of estimates=%5.3f, error of median estimator=%5.3f, std= %f5.3'%(med,med-p,std))\n",
    "    return L,med,p,std,abs((med-p)/std)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": true,
    "nbgrader": {
     "grade": false,
     "locked": true,
     "solution": false
    }
   },
   "outputs": [],
   "source": [
    "def test_report_assert(n,q1,q2,m,r=100):\n",
    "    k1=int(q1*n)\n",
    "    k2=int(q2*n)\n",
    "    L,med,p,std,norm_err=evaluate(n,q1,q2,m,r=100)\n",
    "    hist(L);\n",
    "    plot([p,p],plt.ylim(),'r',label='true prob')\n",
    "    plot([med,med],plt.ylim(),'k',label='median of %d estimates'%r)\n",
    "    mid_y=mean(plt.ylim())\n",
    "    plot([p-std,p+std],[mid_y,mid_y],'g',label='+-std')\n",
    "    legend();\n",
    "    print('normalized error of median=',norm_err,'should be <1.0')\n",
    "    title('r=%d,n=%d,k1=%d,k2=%d,m=%d,\\nnorm_err=%4.3f'%(r,n,k1,k2,m,norm_err))\n",
    "    assert norm_err<1.0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Testing your Functions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* We now test your functions. The graphs below show how close your estimated probability is to the true probability for various values of $k_1$ and $k_2$. You can see that your answer is never exactly the correct probability. \n",
    "* For full credit, the code below must run without error."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "nbgrader": {
     "grade": true,
     "grade_id": "ex2",
     "locked": true,
     "points": "5",
     "solution": false
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "#### test no. 1\n",
      "computed prob=0.954, std=0.021\n",
      "ran estimator 100 times, with parameters n=100,k1=40,k2=60,m=100\n",
      "median of estimates=0.940, error of median estimator=-0.014, std= 0.0208405.3\n",
      "normalized error of median= 0.6957692550136406 should be <1.0\n",
      "#### test no. 2\n",
      "computed prob=0.159, std=0.037\n",
      "ran estimator 100 times, with parameters n=100,k1=55,k2=100,m=100\n",
      "median of estimates=0.130, error of median estimator=-0.029, std= 0.0365355.3\n",
      "normalized error of median= 0.7843141288500418 should be <1.0\n",
      "#### test no. 3\n",
      "computed prob=0.146, std=0.035\n",
      "ran estimator 100 times, with parameters n=100,k1=47,k2=49,m=100\n",
      "median of estimates=0.070, error of median estimator=-0.076, std= 0.0353595.3\n",
      "normalized error of median= 2.163135793265594 should be <1.0\n"
     ]
    },
    {
     "ename": "AssertionError",
     "evalue": "",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mAssertionError\u001b[0m                            Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-10-479fcccfdbd9>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[1;32m      9\u001b[0m         \u001b[0mprint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'#### test no.'\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mi\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     10\u001b[0m         \u001b[0mi\u001b[0m\u001b[0;34m+=\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 11\u001b[0;31m         \u001b[0mtest_report_assert\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mn\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mq1\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mq2\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mm\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mr\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m100\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     12\u001b[0m \u001b[0mtight_layout\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m<ipython-input-9-2dd50d59ad97>\u001b[0m in \u001b[0;36mtest_report_assert\u001b[0;34m(n, q1, q2, m, r)\u001b[0m\n\u001b[1;32m     11\u001b[0m     \u001b[0mprint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'normalized error of median='\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mnorm_err\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m'should be <1.0'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     12\u001b[0m     \u001b[0mtitle\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'r=%d,n=%d,k1=%d,k2=%d,m=%d,\\nnorm_err=%4.3f'\u001b[0m\u001b[0;34m%\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mr\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mn\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mk1\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mk2\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mm\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mnorm_err\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 13\u001b[0;31m     \u001b[0;32massert\u001b[0m \u001b[0mnorm_err\u001b[0m\u001b[0;34m<\u001b[0m\u001b[0;36m1.0\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;31mAssertionError\u001b[0m: "
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAl0AAAHkCAYAAAAEr5oLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzs3XecVPW5x/HPo4BIkyqCiEAkBMFlUZohKlhRiRpjQWMhxhAx3sQbu7FwLRFzjXqNSRSjQY39GpRrbxDFKAiKSA2iKEgVBOm68Nw/fmeXw+zM7uxO212+79frvGbm1GfOnPmd5/xO+Zm7IyIiIiK5tUuhAxARERHZGSjpEhEREckDJV0iIiIieaCkS0RERCQPlHSJiIiI5IGSLhEREZE8UNIlaTGzUWb290LHkQ4z62Rmbmb1Ch2LiBSGyiypiZR0JTCzi8xsqpltMbOxSYYfYWZzzWyjmU0ws31jw3YzswfM7GszW2Zmv9kZ4jWznmb2spl9aWYZP/jNzLqa2ebEAtPMzjSzz8xsg5k9Y2YtqzHv75rZs2a20sxWR3F3q2acXczsOTNbF33338eGtTSzcVGsn5nZmdVZRiZqwrYhuVfbfueaEG+2yiwzmxiVVeujbl5s2CAz2xYbtt7Mzq3GMrJSZlX2nSsrs7JR/mbCzNqZ2XgzW2IhQe2UMLzCbaOi7Sqfdoqky6p29LAEuAl4IMl8WgP/AK4FWgJTgSdio4wCugL7AoOBy81sSPWirlXxfgs8CfysGtMm8yfgvXgPM+sB3AucDbQFNgJ/rsa8mwPjgW7RfKYAz1Z1JmbWAHgVeAPYC+gAxJPEPwHfRMv4CfCX6DvkU03YNqQaVGbVqjLrIndvEnWJydCS2LAm7v5gNeaflTKLyr9zyjIri+VvJrYBLwE/TjF8FCm2jTS2q/xx9zrZAQuBK4AZwBagXhWnvwkYm9BvBPCv2OfGwCbge9HnL4CjY8NvBB5PMf9RhD/AQ8A6YBbQJ4Pvm494/x69rw88BjwNNIiNs1/YpDL63YZF66VseVH/3wGPxj5/h1BANE0yj06Al/7mhD/pQqBnknFbRuO2qmKcI4C3UgxrHMX23Vi/h4HRKcYfSyjAXgTWA28TErk7ga+AuUDvmrptqMtOpzKr9pVZwETg/BTDBgGL05xPzsusir5zZWUW1St/fwosisqwC4C+0ba9Brg7g3VeL5p/p4T+KbeNyrarfHZ1vabrDOB4wpHCM2a2JkX3XJrz6wF8WPrB3TcAC4AeZtYCaB8fHr2vqHbjBOBxth/J3F06wMJpq5oWL2a2O/AMYadwmrt/U1kQ6X4XM2sG3ABcksZ3WUBUSFSy7J8CtwJHuvvMJKMcCixz91XR+GdWEOsaM+sYTTcAWGhmL0bV9RPN7IBo2HeBre7+79hyKlu3pwHXAK0J6/Yd4P3o8/8Ct8e+U43cNiQrVGbVojIrcktUBrxtZoMShu1pZsvN7FMzu8PMGqex7FyVWRWprMyqTvnbn1DzdDrhAPK3wJHRvE4zs8Oi+H9QSfw/qCz4NLaNlNtVZfPOtrp+0d5d7r4oej80C/NrAqxM6LcWaBoNK/2cOCyVSe7+AoCZPQxcXDrA3WtivM0I1bsfAr/26JChMlX4LjcC97v7IjNLHNYkIVaoPN6LgfOAQe6+OHGgmXUgVKmXnft390eBR9OItQOhCvsE4HXg18CzZva9asY6zt2nRXGNAy5094eiz08AF8VirInbhmSHyqzaVWZdAcwmJCDDgP8zs+IoKZkLFEev+wIPEg6eflHB/HJZZlWksjKrOmXaje6+GXjFzDYAj7n7CgAzewvoDfzT3ScRkvhM4y+NKVX8qbarvKrrNV2LKh+lStYT/sRxzQhV7etjnxOHpbIs9n4j0LCK13JUJtvxDgCKCFXOWW0p3cyKCUdBd6QYpaLvksplwJ9SFF5tgFeAP7v7Y1WPmE2EHdCL0ZHzbUAroHs1Y12eMO/Ez03IrmxvG5IdKrNqSZkF4O6T3X2du2/xcL3W28Bx0bBl7j7b3be5+6fA5cAplcwyl2VWRSors2p6mVbZtlGd+HOiriddZX+y6DTQ+hTdi2nObxbQKzbPxoRz27Pc/StgaXx49H5WdQKvofG+AtwCvG5mbbP8XQYRrgX43MyWAZcCPzaz91N8ly7AbkC8OjzR0cA1ZrbDhZdRVfQrwHh3vzlh2E8qiHV9rKp+BrHtK8G/gXpm1jXWr9rbQqIaum1IdqjMqj1lVjIOlKumT2NYqVyWWRWprMyqTvmbFjM7pJL4D6lsHmlsGym3q0zjr7J8X0SWr45wEeKR1ZiuHtCQ8Ed9OHpfenFjG0KV5I+j/rcC78amHQ38E2gBfI+wEQyJDXdCtTGUv1C8E7ELKWtyvIQ7QGYCraPPFs1//2iahsBuVfwejQgXj5d2txGuZWoTDe8BfA0cQrgI8u/ELqAlXIw+NnFdEv5oy4ATomHNCHf/VPtCzmg+3QhH+kcCuwL/SbhGoEE0/HHChbuNgYHR79AjIb5Osdhvis37fGBi7PN+QElN3ZbVZadDZVbO4iU3ZVZz4JjS+Al3/G0AukXDBwEdo2XtA0wA/habfiz5LbMq/M5UXGZVq/yNDV9MrLyIpr+mGt+hYbR8J5TBDdPZNtLYrkYRK3Nz+j/Px0IK0VH9AmxU9IPGu1Gx4UcSztFvIty50ik2bDfCbdBfE6pSfxMb1oFQldkqtpxsFGCFivcmYDrhTppOSWJYmOHvt8Pyon5nAp8TCrZngZaxYa8DP0+2LoE+0fc7Fjg3GraBUOVc2nWsRownAx9H628iUQEVDWtJuHh3QxTzmbFhh0TbZ/3o81hyk3TlZdtQl50OlVm5jjerZRZhR/5etMw1wLvAUbHhvyHcUbeRcNr4j8Tu9iPPZVZl35kKyqxoeLXK36hftpKuxPg9nW0jje3qfuDmfPzPLVqg5JiZnUXYKV9V6FjSUZvitfDMrA+BInf/ttDxVMbMrgFWuvu9hY6lOmrTtiHVV9t+59oUb20rsypSF76LmU0HjvDojtCcLktJl4iIiEju1fUL6UVERERqBCVdIiIiInmgpEtEREQkD5R0iYiIiOSBki4RERGRPFDSJXWGmbU0s3FmtsHMPjOzMysZ/0AzezN66vFyM/t1bNj3zWyKma0zsxmW0OiqmbUxs0ctNMj6lZk9kqvvJSJ1X1XKryRPzP/GzD6KDS82s7fMbK2ZLTaz61LM53ozczM7MhffScqr6w1eSxWZWT13L6kJy61GLH8iNDzbltDQ7PNm9qG7l2vqwcxaExrC/U/Ck+8bEB6uiJm1BMYDI4F/AGcQGrLt4qG5CaL+7xEast0I9KxCnCKSAztL+eXuxyYsayLwRqzXo8A4tjevNsnMprv7+Ng03yG0Bbm0CjFKhlTTVYeY2UIzuzSqmVlrZk+YWcNo2M/N7GMzW21m482sfWw6N7Nfmtl8YH6s34VmNj+q7bnRzL5jZu+Y2ddm9mT0ULzKYhpqZtOjGqF/mVlRQrxXmNkMYIOZ1UvWL83v3pjQxMO17r7eQ8v144GzU0zyG+Bld3/EQ2O169x9TjTs+8Byd3/K3be6+98JLdSfHC3raEKzHpe5+1p3/9bdP0gnThFJTuVXlcqv+LSdCK1cPBzr3Ql4JCq/FgCTCE35xN0NXEFI9CRPlHTVPacBQ4DOQBEw3MwOJ7RzdhrQDviM0M5W3ElAf0K7XKWGAAcBA4DLgTGE9sX2IdTsnFFRIGZ2IKFZhl8ArYB7gfFmtltstDOA44HmsaPCHfqZ2XNRoZesey6a5rvAVnePN8D6IeULmlIDgNVRQbrCzP7PtjcMa5RvmNbYXps1AJgHPGhmq8zsPTM7rKJ1ISJpUfm1XUXlV9w5wFvu/mms353AOWZW38y6AQcDr8W+26nAN+7+QhrzlyxS0lX33OXuS9x9NfB/hGrqnwAPuPv77r4FuAo4ODpCKnWLu692902xfre6+9dR9fZM4BV3/8Td1wIvAr0rieXnwL3uPjk64noQ2EIoBOPxLkpY7g793H2ouzdP0Q2NpmlCaNA0bi3QNEVsHQhtmv2a0Cjtp4TGXgH+BbQ3szOiQutcQov0jWLTHk1owHYv4A/As9EpSxGpPpVf21VUfsWdQ2i/Ne45wqnDTYT2Bu939/cAzKwJ8Dvg4jTmLVmmpKvuWRZ7v5HwZ25PODoEwN3XA6uAvWPjLkoyr+Wx95uSfG5SSSz7ApfEj+wIR5ntY+MkW26yfpVZDzRL6NeM0BhtMpuAce7+nrtvBv4L+L6Z7RG1v3Ui4RTkcsIR82uERltLp13o7vdHpxYfj2IeWI24RWQ7lV/bVVR+AWDhBp+9CNellvZrSbhe9QagISHmY8zswmiU/wIeTqgZkzxR0rVzWEIoQICy6wdaAV/ExslFI5yLCC23x4/sGrn7Y7Fxki13h35W/k6dePdiNNq/gXpm1jU2aS+g3EWokRkJyyl9bwDu/k937+vuLQnXVXQDpqSYVkRyR+VXaucC/4gS0VJdCKcqH3L3EndfTDgde1w0/AjgV2a2zMyWEZKyJ83sikqWJVmgpGvn8CjwUwu3Ee9GqFqe7O4Lc7zc+4ALzKy/BY3N7HgzS6fKvIy7H+vuTVJ0x0bjbCDcUXhDtJyBhNqqh1PM9m/Aj6J1Uh+4Fpjk7msAzKx3dGqxGXAbsNjdX46mHQe0MLNzzWxXMzuFcNT9dpXWjoikQ+VXEma2O3Aq5U8t/jsMtjPNbBcz2ws4nXCNGISkqyfh1G0xIan9BeHuSckxJV07AXd/nZBUPE24Pfg7wLA8LHcq4bqIu4GvgI+B4Tlc5IXA7sAKwvVZI0tvtzazQ8ys7GjQ3d8Argaej8bfD4g/F+dy4EvC0W474EexaVcDJwCXEq67uBI40d2/zNk3E9lJqfwqX35FTiKUPxMS4v6acKf1f0ZxTydc03ZzNHyVuy8r7YCtwFcJtWWSI+ausyQiIiIiuaaaLhEREZE8UNIlGTGzqyu5QFREpEZS+SX5ptOLIiIiInmgmi4RERGRPKiRDV63bt3aO3XqVOgwRCRPpk2b9qW7tyl0HNmg8ktk55NuGVYjk65OnToxderUQochInliZp9VPlbtoPJLZOeTbhmm04siIiIieaCkS0RERCQPlHSJiIiI5EGNvKYrmW+//ZbFixezefPmQociO5mGDRvSoUMH6tevX+hQRHZ62hdIIWW6P6g1SdfixYtp2rQpnTp1wswKHY7sJNydVatWsXjxYjp37lzocER2etoXSKFkY39Qa04vbt68mVatWulPJnllZrRq1UpH1SI1hPYFUijZ2B/UmqQL0J9MAPj888/5/PPP87Y8bXd1n5ntY2YTzGyOmc0ys19H/UeZ2RdmNj3qjit0rJL+fzLfZYXUfZnuD2pV0lVIa9as4c9//nOhw0hp4cKF9OzZs9Bh5MWmTZvYtGlTocOQuqUEuMTduwMDgF+a2f7RsDvcvTjqXihciFJVuSgrtC+QTCjpSlNFf7StW7fmJYZ8LUdkZ+PuS939/ej9OmAOsHdho5KaSPsCyYSSrjRdeeWVLFiwgOLiYi677DImTpzI4MGDOfPMMznggAPKHV3cdtttjBo1CoAFCxYwZMgQDjroIA455BDmzp1bbv6jRo3i7LPP5vDDD6dr167cd999AOWWA3D77bfTs2dPevbsyZ133lk2j5KSEs4991yKioo45ZRT2LhxYw7XiEjdZGadgN7A5KjXRWY2w8weMLMWKaYZYWZTzWzqypUr8xSpFIL2BZKJWnP34g4uvhimT8/uPIuLIbbRJho9ejQzZ85kerTciRMnMmXKFGbOnEnnzp1ZuHBhymlHjBjBPffcQ9euXZk8eTIXXnghb7zxRrnxZsyYwbvvvsuGDRvo3bs3xx9/PMAOy5k2bRp/+9vfmDx5Mu5O//79Oeyww2jRogXz5s3j/vvvZ+DAgZx33nn8+c9/5tJLL81svYjsRMysCfA0cLG7f21mfwFuBDx6/QNwXuJ07j4GGAPQp08fz1/EO7lK9gX7lCYbjRqlP0/tCySHVNOVgX79+lV62+j69ev517/+xamnnkpxcTG/+MUvWLp0adJxTzzxRHbffXdat27N4MGDmTJlSrnlTJo0iR/96Ec0btyYJk2acPLJJ/PWW28BsM8++zBw4EAAzjrrLCZNmpStrypS55lZfULC9Yi7/wPA3Ze7+1Z33wbcB/QrZIxSM2lfIOmqnTVdFRyF5FPjxo3L3terV49t27aVfS69pXTbtm00b9687KioIol3RZR+ji/HPfVBdKrpRaRiFv4s9wNz3P32WP927l66Z/wRMLMQ8UkKlewLFs2bB0C3bt1yGob2BZIu1XSlqWnTpqxbty7l8LZt27JixQpWrVrFli1beO655wBo1qwZnTt35qmnngLCH+XDDz9MOo9nn32WzZs3s2rVKiZOnEjfvn3LjXPooYfyzDPPsHHjRjZs2MC4ceM45JBDgHB79DvvvAPAY489xg9+8IOMvrPITmQgcDZweMLjIX5vZh+Z2QxgMPCfBY1SCk77AsmEkq40tWrVioEDB9KzZ08uu+yycsPr16/PddddR//+/Rk6dCjf+973yoY98sgj3H///fTq1YsePXrw7LPPJl1Gv379OP744xkwYADXXnst7du3LzfOgQceyPDhw+nXrx/9+/fn/PPPp3fv3gB0796dBx98kKKiIlavXs3IkSOz9O1F6jZ3n+Tu5u5F8cdDuPvZ7n5A1P+EWK2X7KS0L5BMWEVVlABm9gAwFFjh7j2jfk8ApfW1zYE17l6cZNqFwDpgK1Di7n3SCapPnz4+derUHfrNmTOH7t27pzN5rTRq1CiaNGmiix3TMC9Ppwzi5syZQ/d77w0fasjp7brEzKalWz7UdMnKL8mequwLClFWZCrlvqD0Ia8dO+Y/KNlBsm0w3TIsnWu6xgJ3Aw+V9nD302ML+gOwtoLpB7v7l2ksR6Rmy/YdsyIi6dIDoeuESpMud38zem5NOdHFp6cBh2c3rJ1P6XNcRERk56V9Qd2W6TVdhwDL3X1+iuEOvGJm08xsRIbLEhEREam1Mn1kxBnAYxUMH+juS8xsT+BVM5vr7m8mGzFKykYAdNQ5axEREaljql3TZWb1gJOBJ1KN4+5LotcVwDgqeLCgu49x9z7u3qdNmzbVDUtERESkRsrk9OKRwFx3X5xsoJk1NrOmpe+Bo9GDBUVERGQnVWnSZWaPAe8A3cxssZn9LBo0jIRTi2bW3sxeiD62BSaZ2YfAFOB5d38pe6HXboMGDaL0tvLjjjuONWvW5GxZW7Zs4cgjj6S4uJgnntixYvKpp56iR48e7LLLLiTe5n7LLbew33770a1bN15++eWy/i+99BLdunVjv/32Y/To0RnHN336dF544YWyz+PHj8/KfAHuvPNONfYqIjWW9gXb7RT7Anevcd1BBx3kiWbPnl2uX2122GGH+XvvvZeXZb3zzjt+6KGHJh02e/Zsnzt3brl4Zs2a5UVFRb5582b/5JNPvEuXLl5SUuIlJSXepUsXX7BggW/ZssWLiop81qxZGcX3t7/9zX/5y1+mPf7cuXN97ty5aY277777+sqVK6sbWpnZs2e7H3ZY6CTrgKleA8qebHTJyi/JnqrsC6pSVhRK2vuCuXNDl4G6ti+oimztC9yTb4PplmF6In2aFi5cyPe+9z3OP/98evbsyU9+8hNee+01Bg4cSNeuXcsaJN2wYQPnnXceffv2pXfv3mVPHN60aRPDhg2jqKiI008/nU2xZ6506tSJL78MjzI76aSTOOigg+jRowdjxowpG6dJkyb89re/pVevXgwYMIDly5eXi3H16tWcdNJJFBUVMWDAAGbMmMGKFSs466yzmD59OsXFxSxYsGCHabp37570wYHPPvssw4YNY7fddqNz587st99+TJkyhUefe522HTqxvkFL5q7YyKBjT+Tehx5nxuI1O3TPT/qAgYOOZP+iYg7sdzDPTpzCjMVruO2esezXrTvd9u/JQf2/z7RPVnDVb6/h0cceLzv6Gjt2LBdddBEAw4cPZ+TIkQwePJguXbrwz3/+k6uvvprjjjuO4cOHl8U7cuRI+vTpQ48ePbj++usBuOuuu1iyZAmDBw9m8ODBALzyyiscfPDBHHjggZx66qmsX78egCuvvJL999+foqIiPaBWpMA6Xfl8ym7xV5vKlTepuvVbSli/pWSHfpnSviDsC6ZMmcJ+++1Hly5daNCgAcOGDUv6hP0FCxYwZMgQDjroIA455BDmzp0LhJq1nj170qtXLw499FC++eYbrrvuOp544om09wXnnXce3bt3r1X7glrZ4PXFF1+cVqOhVVFcXMydlTxp/OOPP+app55izJgx9O3bl0cffZRJkyYxfvx4fve73/HMM89w8803c/jhh/PAAw+wZs0a+vXrx5FHHsm9995Lo0aNmDFjBjNmzODAAw9MuowHHniAli1bsmnTJvr27cuPf/xjWrVqxYYNGxgwYAA333wzl19+Offddx/XXHPNDtNef/319O7dm2eeeYY33niDc845h+nTp/PXv/6V2267rawNsHR88cUXDBgwoOxzhw4d+OKLL1ixagN7td+7rP+e7drz0QfTyk1/w5UXc80tt7Nv5+8w44Op3PzbS/nrE+O5987f85e/P03bdu35eu1a6jdowIWXXM2sGR/w+Nj7ABg7duwO8/rqq6944403GD9+PD/84Q/5+9//zk033bRDAXLzzTfTsmVLtm7dyhFHHMGMGTP41a9+xe23386ECRNo3bo1X375JTfddBOvvfYajRs35tZbb+X222/noosuYty4ccydOxczy2n1vohkz+9HXcW8WR+lHL6tZAsAu9Tbraxf490q3u3V2H1BURGtWrSoMfsCgH322WeH/pMnTy43/YgRI7jnnnvo2rUrkydP5sILL+SNN97ghhtu4OWXX2bvvfdmzZo1NGjQgBtuuIGpU6dy9913A5XvC95++23++te/0rdv31qzL6iVSVehdO7cmQMOOACAHj16cMQRR2BmHHDAASxcuBAI2fP48eO57bbbgNDC/Oeff86bb77Jr371KwCKioooKipKuoy77rqLcePGAbBo0SLmz59Pq1ataNCgAUOHDgXgoIMO4tVXXy037aRJk3j66acBOPzww1m1ahVr11bUWEBqnqR5KDNL2T9u44b1fDh1CpddMLys3zfffANAcd/+XHfJLzl66EkccewP04rlhz/8Ydl6btu2bdnRWI8ePVi4cCHFxcU8+eSTjBkzhpKSEpYuXcrs2bPLreN3332X2bNnM3DgwLKYDj74YJo1a0bDhg05//zzOf7448vWs4hIMgXZF3z2Ga1atKgx+4Jt27Yl7R+3fv16/vWvf3HqqaeW9duyJSTCAwcOZPjw4Zx22mmcfPLJacWSuC+I/wa1ZV9QK5Ouyo5CcmW33bYfLe2yyy5ln3fZZRdKSkqAsIE+/fTTSatpEzfIRBMnTuS1117jnXfeoVGjRgwaNIjNmzcDoRHV0ul33XXXsuXFpZMQpatDhw4sWrSo7PPixYtp3749G+utY9mSL8r6r1i6hD3b7rXDtNu2baPpHnvw5MtvlZvvtbfcwYwPpvLW669w2jGHJB0nUXw9J/4GJSUlfPrpp9x222289957tGjRguHDh5ettzh356ijjuKxx8o/Wm7KlCm8/vrrPP7449x999288cYblcYlIoV1+ahbKhz+7epQVtVvub12vqhD84yXW5B9QZSs1JR9AZCyf6lt27bRvHnzpGem7rnnHiZPnszzzz9PcXFxWmev6sK+QNd0ZdkxxxzDH//4x7KN/oMPPgDg0EMP5ZFHHgFg5syZzJgxo9y0a9eupUWLFjRq1Ii5c+fy7rvvVmnZ8WVMnDiR1q1b06xZs2p9jxNOOIHHH3+cLVu28OmnnzJ//nz69etHj14H8vnCBSz+/DO+/eYbXhr/Dw476tgdpm3StBl779ORV557Bggb+LzZ4RTAooWfUtS7D7+89Gqat2zFsiVf0LhJEzZuWF+tOAG+/vprGjduzB577MHy5ct58cUXy4Y1bdqUdevWATBgwADefvttPv74YwA2btzIv//9b9avX8/atWs57rjjuPPOO7N+6lpEdj51fV/Qt29f5s+fz6effso333zD448/zgknnLDDtM2aNaNz58489dRTQNgXfPjhh0C41qt///7ccMMNtG7dmkWLFu1QXldHbdgX1Mqarprs2muv5eKLL6aoqAh3p1OnTjz33HOMHDmSn/70pxQVFVFcXEy/fuWfEztkyBDuueceioqK6Nat2w7n0dMxatSosmU0atSIBx98sNJpxo0bx3/8x3+wcuVKjj/+eIqLi3n55Zfp0aMHp512Gvvvvz/16tXjT3/6E7vuuiv16tXjqht/z8izfsy2rVs56fSfsF+37uXm+7u77uPmqy/hvrtuo6SkhGNOOJlu+x/A7Tdfx+efLsDd6f+Dw+i2f0/ate/AA3+6k+LiYq666qoqfWeAXr160bt3b3r06EGXLl3KqowhXE9w7LHH0q5dOyZMmMDYsWM544wzyqq4b7rpJpo2bcqJJ57I5s2bcXfuuOOOKscgIhJX1/cFAHfffTfHHHMMW7du5bzzzqNHjx7l5vvII48wcuRIbrrpJr799luGDRtGr169uOyyy5g/fz7uzhFHHEGvXr3o2LEjo0ePrtP7AktWDVloffr08cTnhMyZM4fu3cvv3CW/snH3TyrpVvvPmzcPIGm1fa7MmTOH7iNHhg8TJ+ZtuTsLM5vm7n0KHUc2JCu/pGo6Xfl8ymH3ndCOth27pDWfXJ1eLIio3COP5Z4klywfSbcM0+lFERERkTxQ0iUiIiKSB0q6RERERPJAF9LLTiHTa9GWf7WJtZ+sAmBYkutNFo4+PqP5i4hI3aeaLhEREZE8UNIlIiIikgdKuvLozjvvZOPGjUmHxRv2FBGRukv7gp2Xkq5qmDhx4g6tmqeroj+aiIjULtoXSFXpQvoc2LBhA6eddhqLFy9m69atXHvttSxfvpwlS5YwePBgWrduzYQJE/jb3/7GLbfcQruvJ1a9AAAgAElEQVR27fjud7+7Q1tSIiJSu2lfIIkqTbrM7AFgKLDC3XtG/UYBPwdWRqNd7e4vJJl2CPA/wK7AX919dDaCvvili5m+LLttIhXvVcydQ7LTkPZLL71E+/btef75cJfb2rVr2WOPPbj99tuZMGECrVu3ZunSpVx//fVMmzaNPfbYg8GDB9O7d++sLF9EZGfw+3evYt6qj1IO31YSmnjZpd72JKbxbhXv9rQvkFxK5/TiWGBIkv53uHtx1CVLuHYF/gQcC+wPnGFm+2cSbKH179+f4uJizj//fMaPH09xcXFZ+1RxBxxwAK+99hpXXHEFb731FnvssUe5eU2ePJlBgwbRpk0bGjRowOmnn56vryEiIhnQvkCqq9KaLnd/08w6VWPe/YCP3f0TADN7HDgRmF2Nee0gW0chVTV58mQgnMcfO3YsY8eOBWDRokUUFxcDcMEFF3DBBRcwbdo0XnjhBa666iqOPvporrvuunLzM7O8xS4iqZnZPsBDwF7ANmCMu/+PmbUEngA6AQuB09z9q0LFKTu6fMAtFQ7PVduL2hdIdWVyTddFZnYOMBW4JElBtDewKPZ5MdA/1czMbAQwAqBjx44ZhJV/++yzD9Onbz/duWTJElq2bMlZZ51FkyZNyv6QTZs2Zd26dbRu3Zr+/fvz61//mlWrVtGsWTOeeuopevXqVaBvILLTKyGUY++bWVNgmpm9CgwHXnf30WZ2JXAlcEUB45QaTPsCqUx1k66/ADcCHr3+ATgvYZxkqbunmqG7jwHGAPTp0yfleLXBRx99xGWXXcYuu+xC/fr1+ctf/gLAiBEjOPbYY2nXrh0TJkxg1KhRHHzwwbRr144DDzyQrVu3FjhykZ2Tuy8Flkbv15nZHMKB44nAoGi0B4GJKOmSNGlfIInMvfL8Jjq9+FzphfTpDDOzg4FR7n5M9PkqAHevuD6YkHRNnTp1h35z5syhe/fulcYquZVpczoVSbfaf968eQB069Yt7Xln3AzQ55/QdNi5AAw7s/z9IGoGKDNmNs3d+xQ6Digr094EegKfu3vz2LCv3L1FkmniNfUHffbZZ/kJto7qlKSprVL3ndCOth27pDWfXJ1eLIio3KMK5Z7kRrJ8JN0yrFrP6TKzdrGPPwJmJhntPaCrmXU2swbAMGB8dZYnIpIPZtYEeBq42N2/Tnc6dx/j7n3cvU+bNm1yF6CI1GrpPDLiMUL1emszWwxcDwwys2LC6cKFwC+icdsTHg1xnLuXmNlFwMuER0Y84O6zcvItREQyZGb1CQnXI+7+j6j3cjNr5+5Lo4PNFYWLUERqu3TuXjwjSe/7U4y7BDgu9vkFoNzjJEREahILt4/dD8xx99tjg8YD5wKjo9dnCxCeiNQRteqJ9O6uW2sl79wdT30PiNQNA4GzgY/MrPT2s6sJydaTZvYz4HPg1ALFJxHHtS+QgknnOviK1Jqkq2HDhqxatYpWrVrpzyZ54+6UbPyaz9Z8ywGFDkZyxt0nkfyOa4Aj8hmLVOyzNd/SqtXX1GvUTPsCySt3Z9WqVTRs2LDa86g1SVeHDh1YvHgxK1eurHxkyZnlX23K2bznrNs9rfGWLVsGwLZt29Ked3XjdpzP1nzLHyd/FZ5nIiIF9cfJX/EfwL7Nv8RS5snBtg3h8ZG7rN9S1i/dcqbGico9qlDuSfY1bNiQDh06VHv6WpN01a9fn86dOxc6jJ3esRXcyp2pdB+7MHLkSCA8DTpduYxbRPLn6y3buPnNVWmNu+zRKwHYK/aYl1r7eJeo3KMK5Z7UPNV6ZISIiIiIVI2SLhEREZE8UNIlIiIikgdKukRERETyQEmXiIiISB4o6RIRERHJAyVdIiIiInmgpEtEREQkD5R0iYiIiOSBki4RERGRPKg1zQBJ3dcpzaZ6ln2yqkrj50MuY6m1zZaIiMgOKq3pMrMHzGyFmc2M9ftvM5trZjPMbJyZNU8x7UIz+8jMppvZ1GwGLiIiIlKbpHN6cSwwJKHfq0BPdy8C/g1cVcH0g9292N37VC9EERERkdqv0qTL3d8EVif0e8XdS6KP7wIdchCbiIiISJ2RjQvpzwNeTDHMgVfMbJqZjcjCskRERERqpYwupDez3wIlwCMpRhno7kvMbE/gVTObG9WcJZvXCGAEQMeOHTMJS0RERKTGqXZNl5mdCwwFfuLunmwcd18Sva4AxgH9Us3P3ce4ex9379OmTZvqhiUiIiJSI1Ur6TKzIcAVwAnuvjHFOI3NrGnpe+BoYGaycUVERETqunQeGfEY8A7QzcwWm9nPgLuBpoRThtPN7J5o3PZm9kI0aVtgkpl9CEwBnnf3l3LyLURERERquEqv6XL3M5L0vj/FuEuA46L3nwC9MopOREREpI5QM0AiIiIieaCkS0RERCQPlHSJiIiI5IGSLhEREZE8UNIlIiIikgdKukRkp2dmD5jZCjObGes3ysy+iB6LM93MjitkjCJS+ynpEhGBscCQJP3vcPfiqHshyXARkbQp6RKRnV7UJuzqQschInVbRg1ei4jUcReZ2TnAVOASd/8q2UhmNgIYAdCxY8c8hidV1enK53M274Wjj8/ZvKVuUNJVB+WyUBHZifwFuBHw6PUPwHnJRnT3McAYgD59+ni+AhSR2kWnF0VEknD35e6+1d23AfcB/Qodk4jUbkq6RESSMLN2sY8/AmamGldEJB06vSgiOz0zewwYBLQ2s8XA9cAgMysmnF5cCPyiYAGKSJ2gpEtEdnrufkaS3vfnPRARqdN0elFEREQkD5R0iYiIiORBWklXiiYyWprZq2Y2P3ptkWLac6Nx5pvZudkKXERERKQ2Sbemayzlm8i4Enjd3bsCr0efd2BmLQkXpPYn3G59farkTERERKQuSyvpStFExonAg9H7B4GTkkx6DPCqu6+OnuT8KsnbNxMRERGp0zK5pqutuy8FiF73TDLO3sCi2OfFUb9yzGyEmU01s6krV67MICwRERGRmifXj4ywJP2SNpGhZjR2Tqvrj+Eb+6RK03xzXBh/WYNyZ7Rz6sICLXfQ2P+u9rTFexVz55A7sxiNSP5Vp5yA/JcVmfxXK1U8PbyOHZS7ZdQCtb1My6Sma3npE5uj1xVJxlkM7BP73AFYksEyRURERGqlTGq6xgPnAqOj12eTjPMy8LvYxfNHA1dlsEypY1p+O6LK0yx7IRy17nXm6GyHU6E/R8sdluflThx+fF6XJ1LTVKecgPyXFTn9rw4aFF7vnJi7ZUjOpfvIiMeAd4BuZrbYzH5GSLaOMrP5wFHRZ8ysj5n9FcDdVwM3Au9F3Q1RPxEREZGdSlo1XSmayAA4Ism4U4HzY58fAB6oVnQiIiIidYSeSC8iIiKSB0q6RERERPJASZeIiIhIHijpEhEREcmDXD8cVUQy1OnK53M274Wj9TgKEZF8UdIlIiJSw737ySoAhuXgIEwHX/mj04siIiIieaCkS0RERCQPlHSJiIiI5IGSLhEREZE8UNIlIiIikgdKukRERETyQI+MEBGpo3L5jDcpL5fr+/GczVnySTVdIiIiInmgmq4C0RGo1AR62r2ISP5Uu6bLzLqZ2fRY97WZXZwwziAzWxsb57rMQxYRERGpfapd0+Xu84BiADPbFfgCGJdk1LfcfWh1lyMikmtm9gAwFFjh7j2jfi2BJ4BOwELgNHf/qlAxikjtl61ruo4AFrj7Z1man4hIPo0FhiT0uxJ43d27Aq9Hn0VEqi1bSdcw4LEUww42sw/N7EUz65FqBmY2wsymmtnUlStXZiksEZHKufubwOqE3icCD0bvHwROymtQIlLnZJx0mVkD4ATgqSSD3wf2dfdewB+BZ1LNx93HuHsfd+/Tpk2bTMMSEclUW3dfChC97lngeESklstGTdexwPvuvjxxgLt/7e7ro/cvAPXNrHUWlikiUmOopl5E0pGNpOsMUpxaNLO9zMyi9/2i5a3KwjJFRHJtuZm1A4heV6QaUTX1IpKOjJIuM2sEHAX8I9bvAjO7IPp4CjDTzD4E7gKGubtnskwRkTwZD5wbvT8XeLaAsYhIHZDRw1HdfSPQKqHfPbH3dwN3Z7IMEZFcM7PHgEFAazNbDFwPjAaeNLOfAZ8DpxYuQhGpC/REehHZ6bn7GSkGHZHXQESkTlPbiyIiIiJ5oKRLREREJA+UdImIiIjkgZIuERERkTxQ0iUiIiKSB0q6RERERPJASZeIiIhIHijpEhEREckDJV0iIiIieaCkS0RERCQPlHSJiIiI5IGSLhEREZE8UNIlIiIikgdKukRERETyIOOky8wWmtlHZjbdzKYmGW5mdpeZfWxmM8zswEyXKSIiIlLb1MvSfAa7+5cphh0LdI26/sBfolcRERGRnUY+Ti+eCDzkwbtAczNrl4flioiIiNQY2Ui6HHjFzKaZ2Ygkw/cGFsU+L476iYiIiOw0snF6caC7LzGzPYFXzWyuu78ZG25JpvHEHlHCNgKgY8eOWQhLREREpObIuKbL3ZdEryuAcUC/hFEWA/vEPncAliSZzxh37+Pufdq0aZNpWCIiIiI1SkZJl5k1NrOmpe+Bo4GZCaONB86J7mIcAKx196WZLFdERESktsn09GJbYJyZlc7rUXd/ycwuAHD3e4AXgOOAj4GNwE8zXKaIiIhIrZNR0uXunwC9kvS/J/begV9mshwRERGR2i5bz+kSEZEq6nTl84UOQUTySM0AiYiIiOSBki4RERGRPFDSJSIiIpIHuqZLRKQCZrYQWAdsBUrcvU9hIxKR2kpJl4hI5Qa7+5eFDkJEajedXhQRERHJAyVdIiIVc+AVM5sWtRErIlItOr0oIlKxge6+xMz2BF41s7nu/mZ8hCgZGwHQsWPHQsQoUiPl8ll0C0cfn7N554pqukREKuDuS6LXFcA4oF+Scca4ex9379OmTZt8hygitYSSLhGRFMyssZk1LX0PHA3MLGxUIlJb6fSiiEhqbYFxZgahvHzU3V8qbEgiUlsp6RIRScHdPwF6FToOEakbdHpRREREJA9qfU2X7owQERGR2kA1XSIiIiJ5UO2ky8z2MbMJZjbHzGaZ2a+TjDPIzNaa2fSouy6zcEVERERqp0xOL5YAl7j7+9Et1dPM7FV3n50w3lvuPjSD5YiIiIjUetWu6XL3pe7+fvR+HTAH2DtbgYmIiIjUJVm5psvMOgG9gclJBh9sZh+a2Ytm1qOCeYwws6lmNnXlypXZCEtERESkxsg46TKzJsDTwMXu/nXC4PeBfd29F/BH4JlU81EzGiIiIlKXZZR0mVl9QsL1iLv/I3G4u3/t7uuj9y8A9c2sdSbLFBEREamNqn0hvYV2Me4H5rj77SnG2QtY7u5uZv0ISd6q6i5TREREsiuXz7vMpdr4nM5M7l4cCJwNfGRm06N+VwMdAdz9HuAUYKSZlQCbgGHu7hksU0RERKRWqnbS5e6TAKtknLuBu6u7jEKrrdm/iIiI1Dx6Ir2IiIhIHijpEhEREckDJV0iIiIieaCkS0RERCQPlHSJiIiI5IGSLhEREZE8UNIlIiIikgdKukRERETyQEmXiIiISB4o6RIRERHJAyVdIiIiInmgpEtEREQkD5R0iYiIiOSBki4RERGRPFDSJSIiIpIHGSVdZjbEzOaZ2cdmdmWS4buZ2RPR8Mlm1imT5YmI5Ftl5ZyISLqqnXSZ2a7An4Bjgf2BM8xs/4TRfgZ85e77AXcAt1Z3eSIi+ZZmOScikpZMarr6AR+7+yfu/g3wOHBiwjgnAg9G7/8XOMLMLINliojkUzrlnIhIWjJJuvYGFsU+L476JR3H3UuAtUCrDJYpIpJP6ZRzIiJpqZfBtMlqrLwa44QRzUYAI6KP681sXppxtAa+THPcQlKc2dX6s1uH5jXOg0vf3Dq0KpPVmvVJluO0ql1MsG82l51FaZVhdaj8qknxZC2Wz6r2n02loOsmofypSb8T1Kx4shJLFcsvSLMMyyTpWgzsE/vcAViSYpzFZlYP2ANYnWxm7j4GGFPVIMxsqrv3qep0+aY4s0txZldtibMA0inn6kz5VZPiqUmxQM2KpybFAjUrnpoUSzKZnF58D+hqZp3NrAEwDBifMM544Nzo/SnAG+6etKZLRKQGSqecExFJS7Vruty9xMwuAl4GdgUecPdZZnYDMNXdxwP3Aw+b2ceEGq5h2QhaRCQfUpVzBQ5LRGqpTE4v4u4vAC8k9Lsu9n4zcGomy0hDlav0C0RxZpfizK7aEmfeJSvnsqimrfeaFE9NigVqVjw1KRaoWfHUpFjKMZ3tExEREck9NQMkIiIikgc1OulKo5mhjmY2wcw+MLMZZnZc1L++mT1oZh+Z2Rwzu6rAce5rZq9HMU40sw6xYeea2fyoOzdx2kLHaGbFZvaOma0ws6/M7PRcxZhJnLHhzczsCzNzM9uvJsYZbbevRNvm7Fw2j5VhnL83s1lRnHfpwcZVk8a6T9pMWq7KrzTiOdTM3jezEjM7JWFYVsup6sYSK49mRdtsxuVRJuslGl5a5tydaSyZxpPtsiXDWLJefqQRz2+i7z0jKtf2jQ3Ly762Uu5eIzvCRasLgC5AA+BDYP+EccYAI6P3+wMLo/dnAo9H7xsBC4FOaS73ImAqsAUYm2T4EcBcYCMwAegci7MJ4YaB9cAy4DfRNE8B50bvDwcejt63BD6JXltE71tUcT2lG+824B1gv9J1CewGPAB8DWwGHk0S43eBrsAo4GlgKdA8YRk9CRcafxk2qYx/88OjeL4q/c2Bq6P1WhIN2xR9pycT5vE/wKOEZyntl2QZewKPEW77Xwu8DfTPcNtcnbg84CVgdrSML4CJsWETgaOi902ARjn6D50eW1f/JOE/BBRH63gLMA34eex3/z7hcQmrou4LYFCu/u91rUuyjSQrvy4E7oneDwOeiN5Xu/zKMJ5OQBHwEHBKrH/G5VQWY/ku0DV6354k5VG+YokNLy1z7s7TdpMynmyWLRn+Tt8nlK27Rt07mZYfacYzuPQ7AyNj/6msbsOZdDW5piud5jccaBa934Ptz89xoLGFZ4PtDnxDSCzSsQS4iZCM7MDMWgP/AK4l/HhTgf8rjRP4LWEn+wfCj3+5mQ0hJDivR7OZEPsexwCvuvtqd/8KeBUYkmacVYn34Wj5k4C/s31djiIkVPsSnrp9ZBRvWYzu/m93nx/NchOwAmiTsKhvgScJbW1moh/wMXAl4Vb9hbE4fufuTYB5hMTxVkIBU7a+zOwgoC3wSgXLaBLN+yDCb/gg8LyZNalqnNFv3o/wm5eJtrvBhGS7JXA2cJiZfddCu3313P3V6Hutd/eNVVh2VexJSLhvIfwnyv5DFh5/8CwhKetOWA/XsH3bPJHwkMH+0XdsCQzKUZx1USbNpGVSflU7Hndf6O4zCAczcdkop7ISS7w8cvclJC+P8hILpF3m5CWeHJQtmawbBxoSkqPdgPrA8gxiSTeeCbHv/C7huXqQ/W242mpy0pVO8xujgLPMbDHh7qL/iPr/L6G2ZA2wklAblVah5e7/cPdnCEf3iU4GZrn7Ux7uzBxFSFpK530OIflo5e5zgPuA4YSM/MfROD8CmprZrcAVQD8zW2dmswgbbpWaGEknXkKi8lkUby9ga7Scc4Abo43wPeCDKN7SGBObbGpF+BN9ZmaPmdnTZtbA3ee5+/3RsjKxN+HPuYaQJG6k/PooXZdnAzNL4zSzXQjJ7mWJMzWzH5jZIjMbHP1hb3f3pe6+1cNDLRsA3aoY56Jop/hHwu8c9z3Cf2u1u28Fmkf9RxCO1NeY2T8snBb/bzP7LzN7ysz+Hm0LH0UJ2lUWTusuMrOjqxBfqaWEA4PSg5H4f2gQ4e7lCcAP3f0uQq1K6e9+CPBmNP37hN+jIIVULZVJM2n/C2wg/H6fA7e5e9KHSmc5nlxMm7P5mVk/wn93QSFiqajMKUQ8JC9bdi1ELO7+DqFsWRp1L0f7xExUNZ6fAS9Wc9qcqclJVzrNb5xBOKXWATiO8EywXQgZsRNqTQ4EPgVeMbM1Kbrn0oypB2GnH4Jx30DI3luYWQtCdffnsTg/jKa5lFDT8QFwGOFUzbZo2EeEnfJ4YGjptGb2XBbjtVi8CwgbW4Mo3tLvcymhtvDEWIwlsXk1IVQZXwCMI5ySOi064qhQFb5LI0IN1CWxfom/+aWEZLITIWkojfNC4AV3j/+xMLNjCKcTf+zuE5LEVhyti4+jzz+oINY1ZvYDtm+b/0lITBYlzpZQOxH/zbcQErt6hITmUqAvoaq8GPghoUayBSH5fZnw/9wbuAG4NxbznyuIb0ZCHIlK12cPYAY7bpvrCTv+EuAAoDHhSHFvoCPhNIKkJ5Nm0voRDozaEw4YLzGzLnmIJxfT5mR+ZtaO8H/5qbuXq4HKUyxJy5wMZRJPsrJleCFisXA9bXe2lx+Hm9mhGcRSpXjM7CygD/DfVZ0212py0pVO8xs/I9QslWbWDQmnRM4knAq7y92nE06r3ePuzVN06TbM1YSEU0mEWpm20TAIp2FK41wLNHX3Je5+srv3JpyChLATngtsimpDHgbalU7r7kOzGG98Xa4lFOarY59Lq+qvJSSRv436rYVwoSjwE+DfhARgAaGw25pOEFX4LscDy2OFWCMSfvMozo8J6+vyWJwHAxeZ2ULgtmj0vxKu+zvO3ackxhV9r4eB/yr9ru4+qYJYm7v7pGh9dgV+AVzH9irsUnOj9fguYQf6OqGKvUE07QdRjVsJ8Azhd3/L3V+O+j1FOGUy2t2/JVSjdzKz5lGMF1YQXzwxqug/1ARYm7BtTgF2idZFY8JDjte7+3rgDWD3bFwMu5OoSjNppaekS5tJOxN4yd2/dfcVhGtjMm3WJK3mjHIwbdbnF/1vnweucfd3M4gj01gSy5xzzGx0AeNJVrYcWKBYfgS8Gys/XgQGZBBL2vGY2ZGEfdgJ7r6lKtPmQ01OutJpfuNzwoXimFl3QtK1MurfkHAKqDHhx56bhZjWs/0aslJGuHamZfQ5HmczYJ2ZtY5q4ACuYvv1V3OAo6NasnqECwVfy0KcifGWrUtCwd6XcD1PaYyl13/tAayLxxit+3GEGrF9CbUdo909q0cJUY1TT6BeFOeuhNqs8QnjdSA8cPfBeJzu/hN37+junQhHehCOtJ5094+SLG93wvV477r7LVUMt/SasD8TkvsdWlqIEqXhhCRyGXAHoUZzcTRtCzMrvQ7lcMI2G7/eYRPwZSyp3RS9VuW6s9I4uxIORIwdt831QLOEbfNAQqII4Vqvg82snpnVJ9Rybsr2716HZdJM2ueEmgHLYvmVSXNGLxOVU1FZdXTUL++xxMqjh9z9qQxiyDiWJGXOQ+5e7o66fMVD8rJldoFi+ZxQg15afhxG2N9lotJ4zKw34azACdEBS6lsb8PV5wW4ej/djnDK8N+EmpXfRv1uIKxQCBeov01ICKYDR/v2uzY2EC7Enk045/4iYUeTrHsxybJvIuFuQMI1OW/HPjcmXHc0IoqzhNBMSGmcj0fdKcD8aJy/Emo9RhEuaj+PUHOzkFDdWS+aPqvxRutyPuG05h1Rv3XAtdH7UwjXhX1dGmPU/yzChfLLCEfhSwh3frRNEsN+JLl7MZ3vAlwc/Waro/W4LVru+wm/+f9E/cvWZZLlDY/W5aHRd744YfhuhD/co4SanfiwQyqIdT1wSDRe6Z2UJdF69Cj+22LrM/6bvwP8Ihp2FOHU3kfAWOBG4O+xGI4kuhM3+lwvmn+H6PM9FcQ3K8l/aBkhcYv/h64jJIHxONcRru8C+Bfhjsc5hP/Qy4QEteDlQm3pqLz8akio1fyYUMvYJerfJOo/K1r3l+Upnr7RNrGBUBbMik1bWk59TKjlLkgsbC+Ppse64kKtl9g8hpOFuxez8Dslli0NCvQ77UpIfkrLj9vztG5eIxzAlm4b43O1DVf7OxRqwTn/YiGJObIa09UjFIa3EE49NWR7ItSGcDrux1H/W4ntiIDRhB1VC8LF1EuBIbHhTnTbLFHSFRvWiVjSVZPjJZyGnAm0jj5bNP/92X7XSrlkqJLv0QjYK9bdRriguE3CeK8ANySZfhQ7PpbBCUlgR0KSeGHUvz6hhuuZqq7rhOXtmRCvE2okdo+GF0XroRHhKPjT+DqpZFuoMOmqQoy7RjFcQLj2rCFQPxrWgHBzxa8JSehF0ecG0fALCAXm3oTT0bOACxL+X8ML+R9Xp06dutrW1eTTi4VyDaFW4ErCUdWmqB/uvpKQwNxMeL5Rf3Y8tXQ9IQP/jJDM/Le7vwRlp8XWE45AanW87n4jIWl5zcxaEk47bmL73YubCHdMps3dN7r7stIuWvbm6DsQxbQ3ocr8oSSz2IdQ65k439JT0FeY2fmE02RDCdXLa8xsfdQdUsV4VyTEC+GUYOmpwLMJSeyKaPlHeXR9QQ63hURnE36LvxBq8DYR3Wnp4QaIkwh3sK4hHAWe5NtvjLiXkJx+REiwn4/6lZ7iacX2U5EiIpIGtb2YJ9HdFD3cPadPx8+WWhjvdOAId0/26Iwapbat20TRHZy/dPczCh2LiEhtoqRLREREJA90elFEREQkD5R0iYiIiOSBki4RERGRPFDSJSIiIpIHSrpERERE8kBJl9QJZrabmd1vZp+Z2Toz+8DMjq1g/J5m9rKZfWlmqRpNHWZmc8xsg5ktKH2Wl5ntb2ZTzeyrqHvNzPbP1XcTEZG6QUmX7CBqeLdGLLeKsdQDFhHa+NqD8NT8J82sU4rxvyU0lv6zFPEcRXiC/0+BpoQmhT6JBi8hNJ/TktCu4XhCc08iIiIpKemqQ8xsoZldamYzzGytmT1hZg2jYT83s4/NbLWZjTez9rHp3Mx+aWbzCe3wlfa70MzmRzVHN5rZd8zsHTP72syejJ5MXllMQ81supmtMbN/mVlRQrxXmNkMYEPUOGq5ful8d3ff4O6j3H2hu68qcXAAACAASURBVG9z9+cITe8clGL8ee5+P9ufop/ovwjNDb0bze8Ld/8imnZNtBwnNIG0ldDkkIiISEpKuuqe04AhQGdC+3/DzexwQtuMpwHtCM3+JNbMnERoJih+mmwIIWkZAFwOjAF+QmhypydQ4RPJzexA4AHgF4RmY+4FxpvZbrHRzgCOB5q7e0myfmb2XJS0JeueS7HstsB3SZ1UVRT3rkAfoE2UqC42s7vNbPeE8dYAm4E/Ar+r6nJERGTnoqSr7rnL3Ze4+2pC23nFhETpAXd/P2r/7yrg4IRTb7e4++pY24EAt7r71+4+i9D+3ivu/om7rwVeBHpXEsvPgXvdfbK7b3X3B4EthCQuHu+ihOXu0M/dh7p78xTd0MSFmll94BHgQXefW9kKS6ItoWHsUwhtFhZH3/Wa+Eju3pxwKvMi4INqLEdERHYiSrrqnmWx9xuBJkB7Qu0WAO6+HlgF7B0bd1GSeS2Pvd+U5HOTSmLZF7gkXjNFqCVrHxsn2XKT9UuLme0CPAx8Q0iGqqM0Afyjuy919y+B24HjEkd09w3APcBDZrZnNZcnIiI7ASVdO4clhAQIADNrTDjd90VsnFw0wrkIuDmhZqqRuz9WyXJ36GdmL5rZ/7d371FSlWe+x78PihJucpWg6DQqQQSaApuLh9EAXoJ3E+8RIzEeIsaTsGbGiSZRWUYTnOMg43giwaiYc7yFMShLjQYFkmAUAxERRVQMUQKCohjuCrznjyp6GrpbblW7u+H7WatWVe3L+z71VrP4rXfv2ntNLY/fVNkugHvIz1Sdm1L6bHeKTil9DCyppbaaNAKasm2IlSRpG4aufcODwDcjIlc4n+onwKyU0uIS93s3cGVE9I+8ZhFxekS02JVGUkqnppSa1/KoelmIu4BuwJnbHa6splBPE+CAwvsm251rdh/wvyLi4IhoDYwCnihse3JE9I6I/SKiJflZsI+BBbvyuSRJ+xZD1z4gpfQc+UsoPAosA44ELsqg39nkz+u6k3woeRsYXoq+IuIfyJ+wnwPerzITdklh/eGF94cXdvkH8ocRt55ovx5YWKXJHwN/At4kH6ZeBm4prGsFPAR8Aiwi/8vFoSmlDaX4bJKkvUPkf/UuSZKkUnKmS5IkKQOGLu2RiPjBjk5wlyRJHl6UJEnKhDNdkiRJGaiTmxvvSLt27VJZWVldlyEpI3PmzPkwpdS+ruuQpFKql6GrrKyM2bNn13UZkjISEX/d8VaS1LB5eFGSJCkDhi5JkqQMGLokSZIyUC/P6arJZ599xpIlS9iwwTutKFtNmjShU6dONG7cuK5LkSQ1YA0mdC1ZsoQWLVpQVlZGRNR1OdpHpJRYuXIlS5YsoXPnznVdjiSpAWswhxc3bNhA27ZtDVzKVETQtm1bZ1glSXuswYQuwMBVAu+++y7vvvtuXZdRr/l3J0kqhgYVuurSqlWr+NnPflbXZdRq8eLF9OjRY5f3W79+PevXry9BRZIkqSpD1076vNC1efPmTGrIqh9JklR8hq6ddO2117Jo0SJyuRzXXHMNM2bMYPDgwXz961+nZ8+e1WaabrvtNkaPHg3AokWLGDp0KMceeyzHH388b7zxRrX2R48ezaWXXsqQIUPo0qULd999N0C1fgDGjh1Ljx496NGjB+PGjatsY9OmTVx22WWUl5dz3nnnsW7duhKOiCRJ2hUN5teL2xg1CubOLW6buRxUCTDbGzNmDPPnz2duod8ZM2bw0ksvMX/+fDp37szixYtr3XfEiBGMHz+eLl26MGvWLK666iqmTZtWbbt58+bx4osvsnbtWnr37s3pp58OsE0/c+bM4b777mPWrFmklOjfvz9f/vKXad26NQsXLuSee+5h4MCBXH755fzsZz/jX/7lX/ZsXCRJUlE407UH+vXrt8PLCKxZs4Y//vGPnH/++eRyOb797W+zbNmyGrc9++yz+cIXvkC7du0YPHgwL730UrV+Zs6cyVe/+lWaNWtG8+bN+drXvsYf/vAHAA477DAGDhwIwLBhw5g5c2axPqokSdpDDXOm63NmpLLUrFmzytf7778/W7ZsqXy/9RIDW7ZsoVWrVpUzZJ9n+1/JbX1ftZ+U0i7vL0mS6p4zXTupRYsWrF69utb1HTp0YMWKFaxcuZKNGzfyxBNPANCyZUs6d+7MpEmTgHxoeuWVV2ps4/HHH2fDhg2sXLmSGTNm0Ldv32rbnHDCCTz22GOsW7eOtWvXMnnyZI4//nggf/mHF154AYCHHnqIf/zHf9yjzyxJkorH0LWT2rZty8CBA+nRowfXXHNNtfWNGzfmhhtuoH///pxxxhkcffTRleseeOAB7rnnHnr16kX37t15/PHHa+yjX79+nH766QwYMIDrr7+eQw45pNo2ffr0Yfjw4fTr14/+/ftzxRVX0Lt3bwC6devG/fffT3l5OR999BEjR44s0qeXJEl7Kj7vcNUuNRTRBPg9cCD5w5b/lVK6MSI6Aw8DbYA/A5emlD79vLYqKirS7Nmzt1m2YMECunXrVpRa66PRo0fTvHnzzE98X7hwIQBdu3bNtN+GJvO/v1Gj8s/15FB6qUXEnJRSRV3XIUmlVMyZro3AkJRSLyAHDI2IAcCtwO0ppS7Ax8C3itintHeaO7f4v9CVJNWpop1In/JTZmsKbxsXHgkYAny9sPx+YDRwV7H63VtsvaaXJEnaOxX1nK6I2C8i5gIrgKnAImBVSmlTYZMlwKG17DsiImZHxOwPPvigmGVJkiTVuaKGrpTS5pRSDugE9ANqOgmmxpPIUkoTUkoVKaWK9u3bF7MsSZKkOleSXy+mlFYBM4ABQKuI2HoYsxOwtBR9SpIk1WdFC10R0T4iWhVefwE4CVgATAfOK2x2GVDz9RIkSZL2YsWc6eoITI+IecCfgKkppSeA7wP/FBFvA22Be4rYZ4M1aNAgtl4W47TTTmPVqlUl62vjxo2cdNJJ5HI5HnnkkW3WPf3005xxxhk0atSI7S/T8dOf/pSjjjqKrl278swzz2yzT9euXTnqqKMYM2bMHtc3d+5cnnrqqcr3U6ZMKUq7AOPGjfPG35KkeqGYv16cB/SuYfk75M/vUi2qBo5SePnll/nss89qvBVRly5duOOOO7j11lu3Wf7666/z8MMP89prr7F06VJOOukk3nzzTQC+853vMHXqVDp16kTfvn0566yzOOaYY3a7vrlz5zJ79mxOO+00AM466yzOOuus3W6vqnHjxjFs2DCaNm1alPYkSdpdXpF+Jy1evJijjz6aK664gh49enDJJZfw7LPPMnDgQLp06VJ5c+q1a9dy+eWX07dvX3r37l159fn169dz0UUXUV5ezoUXXsj69esr2y4rK+PDDz8E4JxzzuHYY4+le/fuTJgwoXKb5s2b88Mf/pBevXoxYMAAli9fXq3Gjz76iHPOOYfy8nIGDBjAvHnzWLFiBcOGDWPu3LnkcjkWLVq0zT5HHnkkBx96OGs3buKt5auZt2QV85as4ue/fJhBp57Nwg/Ws7pxazp0KuPBJ57jwSeeo0OnMtYc0IY3Vqxj0Kln13iF/UWLFjF06FCOPfZYjj/+eN544w0AJk2aRI8ePejVqxcnnHACn376KTfccAOPPPJI5UzcxIkTufrqqwEYPnw4I0eOZPDgwRxxxBH87ne/4/LLL6dbt24MHz68sr+RI0dSUVFB9+7dufHGGwG44447WLp0KYMHD2bw4MEA/Pa3v+W4446jT58+nH/++axZk7/KybXXXssxxxxDeXl55heolSTtGxrkDa9HjRq1UzeQ3hW5XI5xO7j699tvv82kSZOYMGECffv25cEHH2TmzJlMmTKFn/zkJzz22GPccsstDBkyhHvvvZdVq1bRr18/TjrpJH7+85/TtGlT5s2bx7x58+jTp0+Nfdx77720adOG9evX07dvX84991zatm3L2rVrGTBgALfccgv/+q//yt13382PfvSjbfa98cYb6d27N4899hjTpk3jG9/4BnPnzuUXv/gFt912W+X9IHfG8veXUd7nvy8Q3qHjIax4fxkAXzzkv6/6cXDHQ/jbW/Or7T9ixAjGjx9Ply5dmDVrFldddRXTpk3jpptu4plnnuHQQw9l1apVHHDAAdx0003Mnj2bO++8E4CJEydu09bHH3/MtGnTmDJlCmeeeSbPP/88v/jFL+jbt29lmLzlllto06YNmzdv5sQTT2TevHl897vfZezYsUyfPp127drx4YcfcvPNN/Pss8/SrFkzbr31VsaOHcvVV1/N5MmTeeONN4iIkh7qlSTtuxpk6KornTt3pmfPngB0796dE088kYigZ8+eLF68GMjPpEyZMoXbbrsNgA0bNvDuu+/y+9//nu9+97sAlJeXU15eXmMfd9xxB5MnTwbgvffe46233qJt27YccMABnHHGGQAce+yxTJ06tdq+M2fO5NFHHwVgyJAhrFy5kk8++WT3PmwNt4eKCLZs2VLj8qrWrFnDH//4R84///zKZRs3bgRg4MCBDB8+nAsuuICvfe1rO1XKmWeeWTnOHTp02OY7WLx4Mblcjl/96ldMmDCBTZs2sWzZMl5//fVqY/ziiy/y+uuvM3DgQAA+/fRTjjvuOFq2bEmTJk244oorOP300yvHWZKkYmqQoWtHM1KlcuCBB1a+btSoUeX7Ro0asWlT/vqvKSUeffTRGu9luH042d6MGTN49tlneeGFF2jatCmDBg1iw4YNQP6G2lv332+//Sr7q6qm+2juqM/adOh4CMuX/q3y/fJlS2nf4YsAvF9l+YplS6vdmHvLli20atWqxtnI8ePHM2vWLJ588klyudxOzVhWHeftv4NNmzbxl7/8hdtuu40//elPtG7dmuHDh1eOW1UpJU4++WQeeuihauteeuklnnvuOR5++GHuvPNOpk2btsO6JEnaFZ7TVWRf+cpX+M///M/KAPTyyy8DcMIJJ/DAAw8AMH/+fObNm1dt308++YTWrVvTtGlT3njjDV588cVd6rtqHzNmzKBdu3a0bNlytz7Hl08+laen/JpPN25kybt/5d3Fi+iRO5buvfrw7uJFLHn3r3z26ac8PeXX1U56b9myJZ07d2bSpElAPuy88sorQP5cr/79+3PTTTfRrl073nvvPVq0aMHq1at3q06Av//97zRr1oyDDjqI5cuX85vf/KZyXdW2BwwYwPPPP8/bb78NwLp163jzzTdZs2YNn3zyCaeddhrjxo0r+qFrSZKggc501WfXX389o0aNory8nJQSZWVlPPHEE4wcOZJvfvOblJeXk8vl6Nev+g86hw4dyvjx4ykvL6dr164MGDBgl/oePXp0ZR9Nmzbl/vvv3+E+U6dO5cc/vpmPV33M1cMvpOsxPRn/wKMc1bUbp5xxDl8dMoD99t+fH9z8v9lvv/0AuO7H/8bIYeeyZfNmzrnwErp3716t3QceeICRI0dy880389lnn3HRRRfRq1cvrrnmGt566y1SSpx44on06tWLww8/nDFjxpDL5bjuuut26TMD9OrVi969e9O9e3eOOOKIysOHkD+37NRTT6Vjx45Mnz6diRMncvHFF1ce7rz55ptp0aIFZ599Nhs2bCClxO23377LNUiStCNR0yGpulZRUZG2v2bUggUL6NatprsKaU8sXLiQNRs30bhNjbfE3CnlnVoVsaL6KfO/v0GD8s8zZmTXZx2KiDkppYodbylJDZeHFyVJkjJg6JIkScqAoUuSJCkDhi5JkqQMGLokSZIyYOiSJEnKgKErQ+PGjWPdunU1rqt6k2dJkrT3MXTthhkzZjB8+PBd3u/zQpckSdq7eUX6Eli7di0XXHABS5YsYfPmzVx//fUsX76cpUuXMnjwYNq1a8f06dO57777+OlPf0rHjh350pe+tM19BSVJ0t6lQYauUU+PYu77xb0/Xu6LOcYNLc6NtJ9++mkOOeQQnnzySSB/T8WDDjqIsWPHMn36dNq1a8eyZcu48cYbmTNnDgcddBCDBw+md+/eRelfkiTVPx5e3AX9+/cnl8txxRVXMGXKFHK5HLlcjmeeeWab7Xr27Mmzzz7L97//ff7whz9w0EEHVWtr1qxZDBo0iPbt23PAAQdw4YUXZvUxJElSHWiQM13FmpHaVbNmzQLy53RNnDiRiRMnAvDee++Ry+UAuPLKK7nyyiuZM2cOTz31FNdddx2nnHIKN9xwQ7X2IiKz2iVJUt0qWuiKiMOAXwJfBLYAE1JK/xERo4H/CXxQ2PQHKaWnitVvfXDYYYcxd+5/H+5cunQpbdq0YdiwYTRv3rwynLVo0YLVq1fTrl07+vfvz/e+9z1WrlxJy5YtmTRpEr169aqjTyBJkkqtmDNdm4B/Tin9OSJaAHMiYmph3e0ppduK2Fe99uqrr3LNNdfQqFEjGjduzF133QXAiBEjOPXUU+nYsSPTp09n9OjRHHfccXTs2JE+ffqwefPmOq5ckiSVSqSUStNwxOPAncBAYM2uhK6Kioo0e/bsbZYtWLCAbt26FbdIsXDhQtZs3ETjNofudhvlnVoVsaL6KfO/v0GD8s8zZmTXZx2KiDkppYq6rkOSSqkkJ9JHRBnQG5hVWHR1RMyLiHsjonUp+pQkSarPih66IqI58CgwKqX0d+Au4EggBywD/r2W/UZExOyImP3BBx/UtIkkSVKDVdTQFRGNyQeuB1JKvwZIKS1PKW1OKW0B7gb61bRvSmlCSqkipVTRvn37YpYlSZJU54oWuiJ//YN7gAUppbFVlnesstlXgfm720epzj+TPo9/d5KkYijmrxcHApcCr0bE1usn/AC4OCJyQAIWA9/encabNGnCypUradu2rde3UmZSSqxcuZImTZrUdSmSpAauaKErpTQTqCkNFeWaXJ06dWLJkiV4vldxvf/++3y6aQuN1mzc7TYWrP5CESuqf5o0aUKnTp3qugxJUgPXYK5I37hxYzp37lzXZex1Ro4cyYvvrOSLXx+z220sHnN6ESuSJGnv5L0XJUmSMmDokiRJyoChS5IkKQOGLkmSpAwYuiRJkjJg6JIkScqAoUuSJCkDhi5JkqQMGLokSZIyYOiSJEnKgKFLkiQpA4YuSZKkDBi6JEmSMmDokiRJyoChS5IkKQOGLkmSpAwYuiRJkjJg6JIkScqAoUuSJCkDRQtdEXFYREyPiAUR8VpEfK+wvE1ETI2ItwrPrYvVpyRJUkNRzJmuTcA/p5S6AQOA70TEMcC1wHMppS7Ac4X3kiRJ+5Siha6U0rKU0p8Lr1cDC4BDgbOB+wub3Q+cU6w+JUmSGoqSnNMVEWVAb2AW0CGltAzywQw4uBR9SpIk1WdFD10R0Rx4FBiVUvr7Luw3IiJmR8TsDz74oNhlSZIk1amihq6IaEw+cD2QUvp1YfHyiOhYWN8RWFHTvimlCSmlipRSRfv27YtZliRJUp0r5q8XA7gHWJBSGltl1RTgssLry4DHi9WnJElSQ7F/EdsaCFwKvBoRcwvLfgCMAX4VEd8C3gXOL2KfkiRJDULRQldKaSYQtaw+sVj9SJIkNURekV6SJCkDhi5JkqQMGLokSZIyYOiSJEnKgKFLkiQpA4YuSZKkDBi6JEmSMmDokiRJyoChS5IkKQOGLkmSpAwYuiRJkjJg6JIkScpA0W54LdU3o54exdz359Z1GbsnV6h74qA6LWNP5L6YY9zQcXVdhiTVG850SZIkZcCZLu21GvQsy6BB+edxM+qyCklSETnTJUmSlAFDlyRJUgYMXZIkSRkwdEmSJGXA0CVJkpSBooWuiLg3IlZExPwqy0ZHxN8iYm7hcVqx+pMkSWpIijnTNREYWsPy21NKucLjqSL2J0mS1GAULXSllH4PfFSs9iRJkvYmWZzTdXVEzCscfmxd20YRMSIiZkfE7A8++CCDsiRJkrJT6tB1F3AkkAOWAf9e24YppQkppYqUUkX79u1LXJYkSVK2Shq6UkrLU0qbU0pbgLuBfqXsT5Ikqb4qaeiKiI5V3n4VmF/btpIkSXuzot3wOiIeAgYB7SJiCXAjMCgickACFgPfLlZ/kiRJDUnRQldK6eIaFt9TrPYlSZIasqKFLtWdsmuf3O19339nZRErkSRJtfE2QJIkSRkwdEmSJGXA0CVJkpQBQ5ckSVIGDF2SJEkZMHRJkiRlwNAlSZKUAUOXJElSBgxdkiRJGTB0SZIkZcDQJUmSlAFDlyRJUgYMXZIkSRkwdEmSJGXA0CVJkpQBQ5ckSVIGDF2SJEkZMHRJkiRloGihKyLujYgVETG/yrI2ETE1It4qPLcuVn+SJEkNSTFnuiYCQ7dbdi3wXEqpC/Bc4b0kSdI+p2ihK6X0e+Cj7RafDdxfeH0/cE6x+pMkSWpISn1OV4eU0jKAwvPBJe5PkiSpXtq/rgvYKiJGACMADj/88DquRlkqu/bJOu1/8ZjT67R/SdK+odQzXcsjoiNA4XlFbRumlCaklCpSShXt27cvcVmSJEnZKnXomgJcVnh9GfB4ifuTJEmql4p5yYiHgBeArhGxJCK+BYwBTo6It4CTC+8lSZL2OUU7pyuldHEtq04sVh+SJEkNVb05kX5fVdcnkUuSpGx4GyBJkqQMGLokSZIyYOiSJEnKgKFLkiQpA4YuSZKkDOzzv17c018PegsZf4EpSdLOcKZLkiQpA4YuSZKkDBi6JEmSMmDokiRJyoChS5IkKQOGLkmSpAwYuiRJkjJg6JIkScqAoUuSJCkDhi5JkqQMGLokSZIyYOiSJEnKgKFLkiQpA/tn0UlELAZWA5uBTSmliiz6lSRJqi8yCV0Fg1NKH2bYnyRJUr3h4UVJkqQMZDXTlYDfRkQCfp5SmrD9BhExAhgBcPjhh2dU1p4ru/bJui5Be2hPv8PFY04vUiWSpL1ZVjNdA1NKfYBTge9ExAnbb5BSmpBSqkgpVbRv3z6jsiRJkrKRSehKKS0tPK8AJgP9suhXkiSpvih56IqIZhHRYutr4BRgfqn7lSRJqk+yOKerAzA5Irb292BK6ekM+pUkSao3Sh66UkrvAL1K3Y8kSVJ95iUjJEmSMmDokiRJyoChS5IkKQOGLkmSpAwYuiRJkjJg6JIkScqAoUuSJCkDhi5JkqQMGLokSZIyYOiSJEnKgKFLkiQpA1nc8Lqkyq59sq5LkCRJ2iFnuiRJkjJg6JIkScqAoUuSJCkDhi5JkqQMGLokSZIyYOiSJEnKgKFLkiQpA5mErogYGhELI+LtiLg2iz4lSZLqk5KHrojYD/g/wKnAMcDFEXFMqfuVJEmqT7KY6eoHvJ1Seiel9CnwMHB2Bv1KkiTVG1mErkOB96q8X1JYJkmStM/I4t6LUcOyVG2jiBHAiMLbNRGxsKRV1a4d8GEd9V1VpnX89dYz6ryGz1Ef6qi1hri1hL1GtX8+9Xos9sA/FLk9Sap3sghdS4DDqrzvBCzdfqOU0gRgQgb1fK6ImJ1SqrCO+lFDfamjPtRQX+qoDzVIUkOUxeHFPwFdIqJzRBwAXARMyaBfSZKkeqPkM10ppU0RcTXwDLAfcG9K6bVS9ytJklSfZHF4kZTSU8BTWfRVBHV+iLOgPtRRH2qA+lFHfagB6kcd9aEGSWpwIqVq57RLkiSpyLwNkCRJUgb2+tC1o1sQRcSBEfFIYf2siCgrLD85IuZExKuF5yFV9plRaHNu4XFwiWooi4j1VfoZX2WfYwu1vR0Rd0RUv7ZAEeu4pEoNcyNiS0TkSjQWJ0TEnyNiU0Sct926yyLircLjshKPRY11REQuIl6IiNciYl5EXFhl3cSI+EuVsciVcCw2V+lnSpXlnQvf3VuF7/KAEo7F4O3+LjZExDm7MxaStE9IKe21D/In7i8CjgAOAF4Bjtlum6uA8YXXFwGPFF73Bg4pvO4B/K3KPjOAigxqKAPm19LuS8Bx5K+D9hvg1FLVsd02PYF3SjgWZUA58EvgvCrL2wDvFJ5bF163LuFY1FbHl4AuhdeHAMuAVoX3E6tuW6oaCuvW1NLur4CLCq/HAyNLWcd2389HQNNdHQsfPnz42Fcee/tM187cguhs4P7C6/8CToyISCm9nFLaej2x14AmEXFgljXU1mBEdARappReSCkl8v8ZnpNRHRcDD+2gr92uIaW0OKU0D9iy3b5fAaamlD5KKX0MTAWGlmosaqsjpfRmSumtwuulwAqg/c4OQDFqqE3huxpC/ruD/HdZsrHYznnAb1JK63amVknaF+3toWtnbkFUuU1KaRPwCdB2u23OBV5OKW2ssuy+wmGT63dwOGtPa+gcES9HxO8i4vgq2y/ZQZvFrmOrC6keuoo5Fru6b6nGYocioh/52aFFVRbfUjjsePsOQvqe1tAkImZHxItbD+mR/65WFb67nW2zWLfpuojqfxc7OxaStE/Y20PXztyC6HO3iYjuwK3At6usvySl1BM4vvC4tEQ1LAMOTyn1Bv4JeDAiWu5km8WsI78yoj+wLqU0v8r6Yo/Fru5bqrH4/AbyM2z/F/hmSmnrDNB1wNFAX/KH275fwhoOT/mrwn8dGBcRR+5mm8Uai57kr8W31a6MhSTtE/b20LUztyCq3CYi9gcOIn9uChHRCZgMfCOlVDmbkVL6W+F5NfAg+UM0Ra8hpbQxpbSy0Ncc8jMqXyps32kHbRatjirrq81mlGAsdnXfUo1FrQrB90ngRymlF7cuTyktS3kbgfso3VhsPbRJSukd8ufV9SZ/P8RWhe9uZ9vcozoKLgAmp5Q+q1LfroyFJO0T9vbQtTO3IJoCbP0l3HnAtJRSiohW5P9jvS6l9PzWjSNi/4hoV3jdGDgDmE/t9qSG9hGxX6GvI4Au5E9iXwasjogBhcN53wAeL9VYFPpvBJxP/pyfUo5FbZ4BTomI1hHRGjgFeKaEY1GjwvaTgV+mlCZtt65j4TnIn0tVkrEojMGBhdftgIHA64Xvajr57w7y32XJtx3u2wAAAQNJREFUxqKKauf57eJYSNK+oa7P5C/1AzgNeJP8LNEPC8tuAs4qvG4CTALeJv8ruCMKy38ErAXmVnkcDDQD5gDzyJ9g/x/AfiWq4dxCH68AfwbOrNJmBfn/yBYBd1K40G0p6iisGwS8uF17pRiLvuRnX9YCK4HXqux7eaG2t8kf1ivlWNRYBzAM+Gy7v4tcYd004NVCLf8PaF6iGv5HoZ9XCs/fqtLmEYXv7u3Cd3lgqcaisK4M+BvQaLs2d2ksfPjw4WNfeHhFekmSpAzs7YcXJUmS6gVDlyRJUgYMXZIkSRkwdEmSJGXA0CVJkpQBQ5ckSVIGDF2SJEkZMHRJkiRl4P8Dtz3ZO5+n8DIAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 720x864 with 3 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# checking functions\n",
    "\n",
    "m=100\n",
    "i=1\n",
    "figure(figsize=[10,12])\n",
    "for n in [100,1000]:\n",
    "    for q1,q2 in [(0.4,0.6),(0.55,1.00),(0.47,0.499)]:\n",
    "        fig=subplot(3,2,i)\n",
    "        print('#### test no.',i)\n",
    "        i+=1\n",
    "        test_report_assert(n,q1,q2,m,r=100)\n",
    "tight_layout()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "nbgrader": {
     "grade": true,
     "grade_id": "ex3",
     "locked": true,
     "points": "5",
     "solution": false
    }
   },
   "outputs": [
    {
     "ename": "AssertionError",
     "evalue": "estimate is incorrect. should be 0.002000, instead is 0.001000",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mAssertionError\u001b[0m                            Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-11-ad4c02648f8e>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[1;32m     13\u001b[0m     \u001b[0mb\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mfloat\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0msum\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0ms\u001b[0m\u001b[0;34m>=\u001b[0m\u001b[0mk1\u001b[0m \u001b[0;32mand\u001b[0m \u001b[0ms\u001b[0m\u001b[0;34m<\u001b[0m\u001b[0mk2\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0mn\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0ms\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mLog\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m/\u001b[0m\u001b[0mm\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     14\u001b[0m     \u001b[0mn_correct\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0msum\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mnn\u001b[0m\u001b[0;34m==\u001b[0m\u001b[0;36m100\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0mnn\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0ms\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mLog\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 15\u001b[0;31m     \u001b[0;32massert\u001b[0m \u001b[0ma\u001b[0m\u001b[0;34m==\u001b[0m\u001b[0mb\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m\"estimate is incorrect. should be %4f, instead is %4f\"\u001b[0m\u001b[0;34m%\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mb\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0ma\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     16\u001b[0m     \u001b[0;32massert\u001b[0m \u001b[0mm\u001b[0m\u001b[0;34m==\u001b[0m\u001b[0mlen\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mLog\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m'should call seq_sum %d times, called it %d times'\u001b[0m\u001b[0;34m%\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mm\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mlen\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mLog\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     17\u001b[0m     \u001b[0;32massert\u001b[0m \u001b[0mm\u001b[0m\u001b[0;34m==\u001b[0m\u001b[0mn_correct\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m'the parameter n should be %d but sometimes it was not.'\u001b[0m\u001b[0;34m%\u001b[0m\u001b[0mn\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mAssertionError\u001b[0m: estimate is incorrect. should be 0.002000, instead is 0.001000"
     ]
    }
   ],
   "source": [
    "# checking functions \n",
    "from numpy import random\n",
    "def seq_sum(n):\n",
    "    #Log.append(n)\n",
    "    s=sum(random.randn(n)>0.5)\n",
    "    Log.append((n,s))\n",
    "    return s\n",
    "\n",
    "n,k1,k2,m = 100,45,50,1000\n",
    "for r in range(10):\n",
    "    Log=[]\n",
    "    a=estimate_prob(n,k1,k2,m)\n",
    "    b=float(sum([(s>=k1 and s<k2) for n,s in Log]))/m\n",
    "    n_correct=sum(nn==100 for nn,s in Log)\n",
    "    assert a==b, \"estimate is incorrect. should be %4f, instead is %4f\"%(b,a)\n",
    "    assert m==len(Log), 'should call seq_sum %d times, called it %d times'%(m,len(Log))\n",
    "    assert m==n_correct, 'the parameter n should be %d but sometimes it was not.'%n\n",
    "    \n",
    "print(\"all good!\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "\n",
    "\n",
    "\n",
    "\n",
    "\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.5"
  },
  "toc": {
   "colors": {
    "hover_highlight": "#DAA520",
    "navigate_num": "#000000",
    "navigate_text": "#333333",
    "running_highlight": "#FF0000",
    "selected_highlight": "#FFD700",
    "sidebar_border": "#EEEEEE",
    "wrapper_background": "#FFFFFF"
   },
   "moveMenuLeft": true,
   "nav_menu": {
    "height": "138px",
    "width": "252px"
   },
   "navigate_menu": true,
   "number_sections": true,
   "sideBar": true,
   "threshold": 4,
   "toc_cell": false,
   "toc_section_display": "block",
   "toc_window_display": false,
   "widenNotebook": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}