{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "<script type=\"text/javascript\" src=\"../_static/linksdl.js\"></script>\n",
    "<div class='alert alert-info'>\n",
    "**This is a fixed-text formatted version of a Jupyter notebook.**\n",
    "\n",
    " You can contribute with your own notebooks in this\n",
    " [GitHub repository](https://github.com/gammapy/gammapy-extra/tree/master/notebooks).\n",
    "\n",
    "**Source files:**\n",
    "[using_numpy.ipynb](../_static/notebooks/using_numpy.ipynb) |\n",
    "[using_numpy.py](../_static/notebooks/using_numpy.py)\n",
    "</div>\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "# Rapid introduction on using numpy, scipy, matplotlib\n",
    "\n",
    "This is meant to be a very brief reminder. It is strongly suggested to refer to more detailed\n",
    "introductions and tutorials see for instance:\n",
    "- [A Whirlwind tour of Python](http://nbviewer.jupyter.org/github/jakevdp/WhirlwindTourOfPython/blob/master/Index.ipynb)\n",
    "- [Python data science handbook](http://nbviewer.jupyter.org/github/jakevdp/PythonDataScienceHandbook/blob/master/notebooks/Index.ipynb)\n",
    "- [Scipy lectures](http://www.scipy-lectures.org/)\n",
    "\n",
    "## Introduction\n",
    "\n",
    "Here we will look at :\n",
    "- basic features regarding array manipulation and indexing\n",
    "- do a bit of plotting with matplotlib\n",
    "- use a number of useful scipy features\n",
    "- see an example of vectorization with a simple Monte Carlo problem"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "## numpy: arrays, indexing etc\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([3, 4, 5])"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.array([3,4,5])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[1, 2],\n",
       "       [3, 4]])"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.array([[1, 2],[3,4]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([  1.,   2.,   3.,   4.,   5.,   6.,   7.,   8.,   9.,  10.])"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "### linearly spaced 1D array\n",
    "np.linspace(1.,10.,10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([  1.        ,   1.29154967,   1.66810054,   2.15443469,\n",
       "         2.7825594 ,   3.59381366,   4.64158883,   5.9948425 ,\n",
       "         7.74263683,  10.        ])"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "### log spaced 1D array\n",
    "np.logspace(0.,1.,10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0.,  0.,  0.,  0.,  0.])"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "### 1D array of zeros\n",
    "np.zeros(5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.,  0.,  0.],\n",
       "       [ 0.,  0.,  0.],\n",
       "       [ 0.,  0.,  0.]])"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "### 2D array of zeros\n",
    "np.zeros((3,3))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "### Types and casts\n",
    "\n",
    "See numpy [dtypes](https://docs.scipy.org/doc/numpy/reference/arrays.dtypes.html)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 1  1  1  2  2  3  4  5  7 10]\n"
     ]
    }
   ],
   "source": [
    "x_int = np.logspace(0.,1.,10).astype('int')   # cast array as int\n",
    "print(x_int)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2\n"
     ]
    }
   ],
   "source": [
    "x_int[1] = 2.34   # 2.34 is cast as int\n",
    "print(x_int[1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "dtype('<U1')"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "array_string = np.array(['a','b','c','d'])\n",
    "array_string.dtype    # 1 character string"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'b'"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "array_string[1]='bbbb'   # 'bbbb' is cast on 1 character string\n",
    "array_string[1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "b'bbbb'"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "array_string = np.array(['a','b','c','d'],dtype=np.dtype('S10'))\n",
    "array_string[1] = 'bbbb'   # 'bbbb' is cast on 10 character string\n",
    "array_string[1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "### array indexing & slicing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "x = np.arange(10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "9"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x[-1]   # last element"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([3, 4, 5])"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x[3:6]  # subarray"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([1, 3, 5, 7, 9])"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x[1::2] # stride"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([9, 8, 7, 6, 5, 4, 3, 2, 1, 0])"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x[::-1] # stride"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0,  1,  2,  3,  4],\n",
       "       [10, 11, 12, 13, 14],\n",
       "       [20, 21, 22, 23, 24],\n",
       "       [30, 31, 32, 33, 34],\n",
       "       [40, 41, 42, 43, 44]])"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x = np.array([np.arange(10*i,10*i+5) for i in range(5)])\n",
    "x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "first column :  [ 0 10 20 30 40]\n",
      "last row     :  [40 41 42 43 44]\n"
     ]
    }
   ],
   "source": [
    "print(\"first column : \", x[:,0])\n",
    "print(\"last row     : \", x[-1,:])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 0  1  2  3  4]\n",
      " [10 11 12 13 14]\n",
      " [20 21 22 23 24]\n",
      " [30 31 32 33 34]\n",
      " [41 42 43 44 45]]\n"
     ]
    }
   ],
   "source": [
    "b=x[-1,:]   # This is a view not a copy!\n",
    "b[:] += 1\n",
    "\n",
    "print(x) # the initial matrix is changed!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[False  True False  True False]\n",
      " [False  True False  True False]\n",
      " [False  True False  True False]\n",
      " [False  True False  True False]\n",
      " [ True False  True False  True]]\n"
     ]
    }
   ],
   "source": [
    "# Fancy indexing \n",
    "print(x % 2 == 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 0  0  2  0  4]\n",
      " [10  0 12  0 14]\n",
      " [20  0 22  0 24]\n",
      " [30  0 32  0 34]\n",
      " [ 0 42  0 44  0]]\n"
     ]
    }
   ],
   "source": [
    "x[x % 2 == 1] = 0\n",
    "print(x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "### Broadcasting"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 5.,  6.,  7.,  8.,  9.])"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x = np.linspace(1, 5, 5) + 4   # 4 is broadcast to 5 element array\n",
    "x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 5.,  6.,  7.,  8.,  9.],\n",
       "       [ 5.,  6.,  7.,  8.,  9.],\n",
       "       [ 5.,  6.,  7.,  8.,  9.]])"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "y = np.zeros((3, 5)) + x   # x is broadcast to (3,5) array\n",
    "y"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "## Plotting with matplotlib\n",
    "\n",
    "We will see some plotting:\n",
    "- Simple plots\n",
    "- Histograms with matplotlib"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "# This is for embedding figures in the notebook\n",
    "%matplotlib inline   \n",
    "import matplotlib.pyplot as plt\n",
    "plt.style.use('ggplot')         # Fancy style"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "## Vectorization or loops: A very simple MC\n",
    "\n",
    "We want to solve a simple statistical question. Assume a Poisson random process of mean mu. What is the density probability function pdf(n_val) of having at least one realization of the Poisson process out of N larger than n_val? \n",
    "\n",
    "See for instance [this paper](https://arxiv.org/pdf/0903.4373.pdf)\n",
    "\n",
    "While this problem has an analytical solution we would like to test it with a simple MC. \n",
    "\n",
    "We will first do it as one would do it with a C code and we will progressively vectorize the problem. We will use a timer to compare performance.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1 loop, best of 3: 12.1 s per loop\n"
     ]
    }
   ],
   "source": [
    "### Define the function\n",
    "def poisson_sample_maximum(mu, N, Ntrials):\n",
    "    \"\"\"\n",
    "    Generate a set of Ntrials random variables defined as the maximum of N \n",
    "    random Poisson R.V. of mean mu\n",
    "    \"\"\"\n",
    "    res = np.zeros(Ntrials)\n",
    "    ### Do a loop\n",
    "    for i in range(Ntrials):\n",
    "        ### Generate N random varslues \n",
    "        Y = np.random.poisson(mu, size=(N))\n",
    "        ### Take the maximum \n",
    "        res[i] = np.max(Y)\n",
    "\n",
    "    return res \n",
    "   \n",
    "mu = 5\n",
    "N = 10\n",
    "Ntrials = 1000000\n",
    "    \n",
    "%timeit values = poisson_sample_maximum(mu, N, Ntrials)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "It does work, but no so fast...\n",
    "\n",
    "To do it in a efficient and pythonic way we have to avoid loops as much as possible.\n",
    "\n",
    "The idea here will then be to do all trials at once requiring random.poisson to produce a 2D matrix of size Nxtrials"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1 loop, best of 3: 1.1 s per loop\n"
     ]
    }
   ],
   "source": [
    "### Define a better function\n",
    "def poisson_sample_maximum_better(mu, N, Ntrials):\n",
    "    \"\"\"\n",
    "    Generate a set of Ntrials random variables defined as the maximum of N \n",
    "    random Poisson R.V. of mean mu\n",
    "    \"\"\"\n",
    "    ### Generate N*Ntrials random values in N x Ntrials matrix\n",
    "    Y = np.random.poisson(mu,size=(N,Ntrials))\n",
    "    ### Return the maximum in each row\n",
    "    return np.max(Y,0)\n",
    "   \n",
    "mu = 5\n",
    "N = 10\n",
    "Ntrials = 1000000\n",
    "    \n",
    "%timeit values = poisson_sample_maximum_better(mu, N, Ntrials)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "We can now compare the distribution of MC simulated values to the actual analytical function.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1e-06, 1)"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAD+CAYAAAAzmNK6AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XtUW9edL/CvhHgKEEIyYDB+QOwmOI4dWcQ2sRMwJKa5\nvVkmbZmbTmfdLjdO7LbxeNLeG5zlm7bTespaLdesNRdP09Z12rROM3NXkzY3fapp3CQ0Do9gO+AH\n+BVjA0ISICEQSDrn/kFg7BpsISSdc3S+n7+is6Sj37YUvtp7n72PRhRFEUREpDpaqQsgIiJpMACI\niFSKAUBEpFIMACIilWIAEBGpFAOAiEilGABERCrFACAiUildrN7I5/PhRz/6EXQ6HVavXo0tW7bE\n6q2JiGgWCwqAQ4cOob29HQaDAQ0NDTPHOzo6cOTIEQiCgMrKSmzfvh3vv/8+Nm7cCKvVioMHDzIA\niIgktqAhoPLycjz33HM3HBMEAYcPH8Zzzz2HgwcP4t1330Vvby+cTifMZvPUm2o58kREJLUF9QBK\nSkpgt9tvONbT04O8vDzk5uYCAMrKytDS0gKTyQSn04nly5fjVtsP2Ww22Gw2AEB9ff1CyiMioluI\n+ByAy+WCyWSaeWwymdDd3Y1PfvKT+PGPf4z29nasX79+ztdXVVWhqqpq5vG1a9fCqsNsNsPhcIT1\nWiVju9VFre0G1Nv2UNqdn58f0rkiHgCz/brXaDRISUnBl770pZDO0draira2Njz11FORLo+IiD4W\n8QCYHuqZ5nQ6YTQa53UOq9UKq9Ua6dKIiOg6EZ+NLS4uRl9fH+x2OwKBAJqbm+f9x7y1tRUvvPBC\npEsjIqLrLKgH0NjYiK6uLng8HuzatQu1tbXYunUrduzYgQMHDkAQBFRUVKCwsHBe52UPgIgo+hYU\nAHv37p31uMVigcViCfu8nAMgIoq+mK0Eng/2AIiIok+WK7I4B0BEFH3sARARqZQsewBERBR9sgwA\nDgEREUUfh4CIiFRKlj0AIiKKPgYAEZFKyTIAOAdARBR9nAMgIlIpWfYAiIgo+hgAREQqJcsA4BwA\nEVH0cQ6AiEilZBkAdKOJgICm4/0Y9PqxNk+PtXlpWGlOhU6rkbo0IlIwBoDMeZ/6NP5lzQ50Za3A\nstE+/MK+GC+f0iJFp8Wa3FSszdNjw5IM5KQnSl0qESkMA0DGxvxBfPueL+KMYTn2nP4FHrB3wKNL\nxYdZxThVswcn+r1ouerFz0848C8PLUVRdorUJRORgjAAZGrMH8Q//7kXZw3LsLfrZWwePAEAyAiM\nY5PjQ2z64ZMAgN60RfjmPTvxzdeGUf9BE8wvvyFl2USkILK8CkjtvJNBfOPNKzjnGMczXUdn/vjP\nZsnYIP7XycPwa3X41j1fxMi4P4aVEpGSyTIA1HwZ6OhkEF9/8wp6nD78jy0FKBs8ddvXLB0bwL4P\nX4Q9xYj/+XoXJgJCDColIqWTZQBYrVZV3hA+KIj45ptXcHHIh2cfKMCmwoyQX1sycgl7u15GZ58H\nDe9eQ1AQo1gpEcUDzgHIRHDno+gwrsS5tTvx5TP/DuubrQjO8xybHB9ib3kRDr51AT9oHcCu0lxo\nNLxUlIhmxwCQkWO5FqQFxrHF3hH2OT6zNh+X7cP4ZZcLplQdateYI1ghEcUTBoBM+LSJOL7obmwZ\n6ECSEAj7PAM1Zfh7AM47/w4/x3rc9ZN/RsnIJST88NeRK5aI4oIs5wDU6Piiu+FLSMaDA20LPpcG\nwK5zv4TZN4QfrdyOIDgMREQ3YwDIxLFcC3LGXbhz5HJEzpcs+PGF8/8Pl9Lz8Yf8jRE5JxHFl5gF\nwMDAAP7t3/4NDQ0NsXpLxXCNB3DSuBIPDLRDi8hdvbNp8BTuHurByyu2we0Lf1iJiOJTSAFw6NAh\nPPHEE/jqV796w/GOjg784z/+I55++mm89tprtzxHbm4udu/eHX6lcewvl0YgaLQoH2iP6Hk1AJ7o\n/hXGdMn42QlHRM9NRMoX0iRweXk5qqur0dTUNHNMEAQcPnwY+/fvh8lkwr59+2C1WiEIAo4ePXrD\n63fv3g2DwRDZyuPIWxfdWOn+CPnjkf8jvXRsAI/0NuMNzRZsW5mFYu4XREQfCykASkpKYLfbbzjW\n09ODvLw85ObmAgDKysrQ0tKCmpoa1NXVRb7SOHVpyIeLQxN4IsK//q/3d5f+iLdXluMHLQOof3gp\n1wYQEYAFXAbqcrlgMplmHptMJnR3d8/5fI/Hg5dffhmXLl3Cq6++ipqamlmfZ7PZYLPZAAD19fUw\nm8O7jl2n04X92lj6xemLSNBqcL997v1+Fkof9OHLm4vwL7ZutDlEVN+1KGrvJRWlfN6RptZ2A+pt\neyTbHXYAiOLNk5W3+mWZkZGBJ5988rbnraqqQlZWFtrapi6HdDjCGxYxm81hvzZWgoKI358egGVx\nGgx+b1TfqzRHi1WmFPyfv5xHSZaItMSEqL5frCnh844GtbYbUG/bQ2l3fn5+SOcK+yogk8kEp9M5\n89jpdMJoNIZ7uhuoZS+gUwNjcI0HUL4i+vMjWo0GT5bmYtgXxCunnLd/ARHFvbADoLi4GH19fbDb\n7QgEAmhubo7YfXzVshvony+OIC1Ri9KC9Ji830pTKiqLDXj9jAu97omYvCcRyVdIQ0CNjY3o6uqC\nx+PBrl27UFtbi61bt2LHjh04cOAABEFARUUFCgsLI1KUGm4K7wsIeO+KB5uXZSJZp533xm/zFdz5\nKADgc4l6vL1xH/7jh/+BPWf+HQC4TQSRSoUUAHv37p31uMVigcViiWhBwFQPoK2tLa6Hgd674oEv\nIKIiBsM/18vye/HwtffwxpL78XeX/ohc31BM35+I5EOWW0GoYQ7gzxfdWJSmQ0lOaszf+9Erf4FW\nFPFaYXnM35uI5EOWARDvXOMBnOz34sEVBmgluCbfNOlGRX8r/rS4FK6kzJi/PxHJgyy3g47XIaDp\ncfhTOWshlPw9NvzkeQSbrkpSS81Hb+FPi0vxq8IH8IQkFRCR1GQZAPE+CXwhvQA6IYCl3n7Jasjz\nubBloAN/yN+IWl8AmSmy/CoQURRxCEgCF9LzsdTbj0Qx2tf+3NpjH/0ZEwlJeP0sJ4KJ1EiWARDP\n6wBEABcyCrBi9JrUpaBwzI6Ng6fwxrkhjPmlDSMiij1ZBkA8XwXkSM7CaKIeRR5pxv7/1qcvvwnv\npIDfnBuWuhQiijFZBkA8u5BRAAAoGpVHABSPXsW9i/X49WkXJgKC1OUQUQzJMgDieQjoYno+tKKA\nZaN9Upcy47N3mzAyEcQfetgLIFITWV76Ec9XAV1IL0D+2CBSBL/UpcxYnZOGkkWpePW0C9UrjUhM\n4P0CiNRAlj2AeHYhI182wz/X++zdJjjHAjh2aUTqUogoRmTZA4hXw4l6uJKzsMIj/RVA1wvufBT3\nAFhu3Ytf266hvL4RGnCTOKJ4xx5ADF1Ml9cE8PU0AB7pfReX0/Nx2rBc6nKIKAZkGQDxOgkstyuA\n/tYWewfS/WP4TcH9UpdCRDEgyyGgeJ0Evpiej9xxJ/QBn9SlzCpZ8KOyrwWvF26GMykTOVIXRERR\nJcseQLyaWgEsz1//07Zd+ytEaPCH/I1Sl0JEUcYAiBHvZBD9qWYUyWwC+G/l+VxY7zyDP+RvgD/I\nhWFE8YwBECMXh6buwSvX8f/rffJqM0aSMtD8kUfqUogoihgAMXLeNTXuL/chIABYO9SN/LFBvMH9\ngYjimiwDIB6vArow5INxwg3j5KjUpdyWFiI+ebUZZx3j6HHKc8KaiBZOlgEQj7uBXnD5FPHrf1p5\nfxtSdBq8cY73CiCKV7IMgHgzERDQ656UzRbQodAHfahYYcDbl9xw+wJSl0NEUcAAiIFLwxMQRGVM\nAF/vkVVG+AURfzzP/YGI4hEDIAYufDwBXCSDu4DNx9KsZKzJTcNvzw0hKIhSl0NEEcYAiIELQz6k\nJ2mxyKe88fT/ssqIwbEAWq/Kf/KaiOZHlltBxJsLrgkUGVOgtF32gzsfxXqNFuYNz+KN187BevJH\nALhLKFG8iGkAvP/++2hvb4fb7ca2bduwdu3aWL69JAKCiEvDE/jUJ4xSlxKWBFFAZV8LXlnxMOwp\nRuQosBdDRLMLOQAOHTqE9vZ2GAwGNDQ0zBzv6OjAkSNHIAgCKisrsX379jnPcd999+G+++7D6Ogo\nXnrpJVUEQO/IBAKCiCJjstSlhG1rfyv+fXkV3syz4r9d+qPU5RBRhIQcAOXl5aiurkZTU9PMMUEQ\ncPjwYezfvx8mkwn79u2D1WqFIAg4evToDa/fvXs3DAYDAOCXv/wltm3bFqEmyNv0CuDi7BSJKwnf\noolhrB3qxpt5Vnz2kg0JUhdERBERcgCUlJTAbrffcKynpwd5eXnIzc0FAJSVlaGlpQU1NTWoq6u7\n6RyiKOLnP/851q1bh6KiogWWrgwXhiaQnKDB4owkqUtZkMq+FjSs/jxOGe/AeqmLIaKIWNAcgMvl\ngslkmnlsMpnQ3d095/N/+9vf4tSpUxgbG0N/fz8efvjhm55js9lgs9kAAPX19TCbzWHVptPpwn5t\nJH3kuYZVOenIzVmEAamLWYD7HJ3I8HthW3wftsng3/VvyeXzjjW1thtQb9sj2e4FBYAo3nxtuEYz\n97UujzzyCB555JFbnrOqqgpZWVloa2sDADgcjrBqM5vNYb82UgRRxDn7KLYWZUpey0IlikE8MNCO\n3+dvwoWrA8hMltdAkBw+bymotd2AetseSrvz8/NDOteC1gGYTCY4nc6Zx06nE0bjwq92iZe9gPo8\nfvgCgqLH/69X2deCgFaHYxe5MpgoHiwoAIqLi9HX1we73Y5AIIDm5uaI3MoxXnYDnVkBbIyPAFju\n7Uex+wr+eH5k1t4fESlLyENAjY2N6Orqgsfjwa5du1BbW4utW7dix44dOHDgAARBQEVFBQoLCxdc\nlNLvCRzc+SgA4KPlD0GzrBL5dY8jKAYlrioyqvpb8EJmIXpcPqw0pUpdDhEtQMgBsHfv3lmPWywW\nWCyWiBUETPUA2traFD8MNJicBeOkG4lx8scfADYPdODIXZ/GH3tGGABECifLvYDiZQ5gMMWIRb74\nuquWPuhD2dIMvH3ZjYkA7xlMpGSyDIB44UjJirsAAICHirMw5hd4z2AihZNlAMTDJLAADRzJWTBP\nxN/eOatzUrE4IxG28/EXbkRqIsvdQJU+CQwAw0npCGh1cdkD0Gg0qCrKwksnBtHnmVT8KmcitZJl\nDyAeOJKzAEztoxOPKooyodUANt4tjEixZBkA8TAENJgytSDOHKfbJ5vSEmFZrMefLozwbmFECsUh\noCgZTPm4BxCHQ0DT6xzKF61B6+p/QEfd17BuqJs3iiFSGFn2AOKBIzkLaYFx6IM+qUuJGqvzNNIC\n4/hL7r1Sl0JEYZBlAMTLEJA5Dn/9Xy9JCKDMfhLvmdfAp02UuhwimidZBkA8LARzJGfF7QTw9R4c\n+AA+XTJazKulLoWI5kmWARAPpnoA8TkBfL27Ri7C7BvCMQ4DESkOAyAKxhOSMZqYFpcTwH9LCxEP\nDHyAjuxVGPYFpC6HiOaBARAFg3G+BuBvPTjQDkGTgLcvuaUuhYjmQZYBoPRJYMfHl4CqYQgIAArH\n7FjhuYpjDAAiReE6gChQWw8AmOoFvJhRgF73BJZkJktdDhGFQJY9AKUbTDEiQQjCOKGeX8Sb7Seg\n1QDHLqqnzURKxwCIAkdKFkwTI0iAerZIyJ50457cNBy75ObtIokUggEQBYPJRphVNPwzrXyFAQOj\nfpxxjEtdChGFgAEQBYMpWVikkgng620szEByggZvcRiISBFkGQBKvgooKIhwJWeqagJ4WmqiFhsK\nM/DOZTf8QQ4DEcmdLANAyVtBuMYDEDQJcb8P0FzKl2didFJA+7VRqUshotuQZQAo2aDXDwCqHAIC\ngHWL9TCkJOAtrgkgkj0GQITZpwNAhUNAAJCg1eCBZZlo6R3F6GRQ6nKI6BZkuRBMyRzeqf1w1LIK\n+HrTN4rZnFGI19c/jb9+69vY2t8GALxZDJEMsQcQYYNjfmT4vUgR/FKXIpk7PFeQO+7Eu4vWSl0K\nEd0CAyDCBr1+1U4AT9MAKBs8iZPGlfDo0qQuh4jmwACIsEGvX7UTwNe7334CQW0C3lt0t9SlENEc\nYjYH0Nvbi9/85jfweDxYs2YNHn744Vi9dcyIogi7N4A1Kp0Avt6K0WtYPDaId3PW4qG+96Uuh4hm\nEVIAHDp0CO3t7TAYDGhoaJg53tHRgSNHjkAQBFRWVmL79u1znmPJkiV48sknIQiCYhd53Y53UoAv\nIKh+CAiYGga6334Cv1y2FcOJ6TBJXRAR3SSkACgvL0d1dTWamppmjgmCgMOHD2P//v0wmUzYt28f\nrFYrBEHA0aNHb3j97t27YTAY0Nraitdeew3V1dWRbYVMDI5NXwLKISAAuH/wBP7v8ir8ddEafErq\nYojoJiEFQElJCex2+w3Henp6kJeXh9zcXABAWVkZWlpaUFNTg7q6ulnPM73P/3e+8x1s3rx51ufY\nbDbYbDYAQH19Pcxmc8iNuZ5Opwv7teE67XYCgCpuBRmKZd4BFHr78W7OWnwhyp+FFJ+3HKi13YB6\n2x7Jdoc9B+ByuWAy/WfH3mQyobu7e87nd3Z24vjx4wgEArj33rlvIF5VVYWqqqqZxw6HI6z6zGZz\n2K8N14W+qV/+nAT+T/fbT+KV5VU4+1EfTGmJUXsfKT5vOVBruwH1tj2Udufn54d0rrADYLY93zUa\nzZzPX716NVavXh3SuVtbW9HW1qa4/YDsXj8StRpk+r1SlyIbZYMn8IsVD6P5Iw/+653ZUpdDRNcJ\nOwBMJhOcTufMY6fTCaPRGJGilHpLyEGvH2a9DloV3QjmdpaMDWLZ6DW8fTmVAUAkM2GvAyguLkZf\nXx/sdjsCgQCam5sj9kdbqdtBD3r9WKSP3jCHUt1vP4GzjvGZjfKISB5C6gE0Njaiq6sLHo8Hu3bt\nQm1tLbZu3YodO3bgwIEDEAQBFRUVKCwsjEhRiu0BjAVgWayXugzZ2Ww/gaNFn8Q7l92oKeEFoURy\nEVIA7N27d9bjFosFFoslogUBypwD8AcFDI0HkMMewE3yfC4UZ6fgncseBgCRjMhyN1Al9gAcYx/v\nAqqX5T+p5DYvy8BPPhhEn2cSizOSpC6HiCDTvYCUOAcwcyMY9gBmtXlpJgDg3Y88EldCRNNk+XNV\niT2A6QDgENDsctIT8QlzCt657MZnVnMYiEgOZBkASjT48RCQKY3/pLMJ7nwUZUs248gdj+LKnh3I\nH59ayMIbxRBJh0NAETLo9cOYkoCkBFn+k8pCmf0kAODdnHskroSIAJkGgNVqVdQVQMD0IjAO/9yK\nadKNVSOX8Z55jdSlEBFkGgBKNOgNcAI4BBsdp3AxowADKZFZNU5E4WMARIAoinCM+TkBHIKNgx8C\nAHsBRDIgywBQ2hzAyEQQk0ERZk4A31aez4UVnqu8VSSRDMjyL5bSLgPlJaDzs9FxCi+vqIYrKROL\npC6GSMVk2QNQGi4Cm5/pYaDj5tC2Byei6GAARMCgd3obCAZAKArH7CjwDnAYiEhisgwApc0BDHr9\nSNFpkJEky39OWdro+BCdWUVw+wJSl0KkWrL8i6W0dQCOsQDMaYm3vCMa3Wjj4IcQNAl4/+qo1KUQ\nqZYsJ4GVIrjzUQDAyLpdyAQQ3Dn7ttl0s6LRq8gZd6H5Iz2qirOkLodIlWTZA1AaT2Ia7wM8TxoA\nGxwf4kS/F97JoNTlEKkSAyAC3Il6ZDAA5m3T4CkEBKCVw0BEkmAALJAADXsAYVrl/gjGVB3+eoX3\nCCCSgiwDQElXAY3pUiBoEpDpH5O6FMXRQsTGJelou+bFRECQuhwi1ZFlACjpKiB34tRN4DkEFJ5N\nSzMwGRTRfo3/fkSxJssAUBJPYhoAcAgoTHfnpCEjScthICIJMAAWyDPTA+AQUDgStBrctyQDLVdH\n4Q9yGIgolhgAC+T+uAfAIaDwlS3NwJhfwIl+hihRLDEAFmh6DoBDQOFbm5eGVB2HgYhijQGwQJ5E\nPXRCAKnBCalLUazEBC0s+Xq0XB2FIIpSl0OkGgyABZpaBDYG7gK0MBuWpGPEF8Q5h0/qUohUI6YB\n4PP58Oyzz6KtrS2WbxtVXAQWGesL0pGgAY73chiIKFZC2gzu0KFDaG9vh8FgQENDw8zxjo4OHDly\nBIIgoLKyEtu3b7/leX71q19h06ZNC6tYZrgNxMJMb6iXCmD12p043pKJzx96Agk//LW0hRGpQEgB\nUF5ejurqajQ1Nc0cEwQBhw8fxv79+2EymbBv3z5YrVYIgoCjR4/e8Prdu3fj8uXLWLJkCfx+f2Rb\nIDFPoh7LRvukLiMulDo6cXjldlxLNaNQ6mKIVCCkACgpKYHdbr/hWE9PD/Ly8pCbmwsAKCsrQ0tL\nC2pqalBXV3fTOTo7OzExMYHe3l4kJSXh3nvvhVZ78wiUzWaDzWYDANTX18NsNs+7UQCg0+nCfm2o\nBsAeQCSVOrpweOV2vG9ejXvn+dnF4vOWI7W2G1Bv2yPZ7rDvB+ByuWAymWYem0wmdHd3z/n8xx9/\nHADw1ltvISMjY9Y//gBQVVWFqqqqmccOhyOs+sxmc9ivDVUQGnh1qdwHKEJyJoaxwnMVLaYSfHqe\nn10sPm85Umu7AfW2PZR25+fnh3SusCeBxVku1wvljljl5eVYv379LZ+jlM3gxnSpEDRaTgJHUKmz\nC2cNyzDCW0USRV3YAWAymeB0OmceO51OGI3GiBSllM3gpheBpbMHEDGlji4IGi1aeI8AoqgLOwCK\ni4vR19cHu92OQCCA5uZmWK3WiBSllB6AmxvBRVzR6FWYfUN4v5cBQBRtIc0BNDY2oqurCx6PB7t2\n7UJtbS22bt2KHTt24MCBAxAEARUVFSgsjMy1G1arNWJhEk0ebgMRcRpMDQP9SZ+NiYCAZB3XKhJF\nS0gBsHfv7Dc7t1gssFgsES0ImOoBtLW1yX4YiPcCiI5SRxd+W3A/Ovq92LAkQ+pyiOJW2FcBRZNS\negDcCC46Vg9fQFqiFu/3jjIAiKJIlv1rpcwBeBL1SAr6kSzE1+I2qSWKQazP16OldxRBgZvDEUUL\newAL4E5MQ4bfy43gouC+JRl4+7IH55zjuGtRmtTlEMUlWfYAlMKTqOfwT5Ssz9dDpwWvBiKKIlkG\ngHKGgNI4ARwl+qQE3J2ThuMMAKKo4RDQArgT9Sj2DUtdRlwK7nwUpfmb8MNVNfhozxdRMD4IANwl\nlCiCZNkDUAoPN4KLKqvzNADgfXOJxJUQxScGQJiCgojRxDRkBLgNRLQsmt4cjgFAFBWyDAAlzAF4\nJoMAuAYg2kqdXTiXuQweHa8EIoo0WQaAEjaDc08wAGJhvfM0BI0W7dmfkLoUorgjywBQAo9vKgA4\nBxBdxZ6ryJr0oM10l9SlEMUdBkCY2AOIDS1EWJxn8EH2KgQ1/LoSRZIs/49SwhzAdABk8F4AUWd1\nnoY3MQ1nMpdJXQpRXOE6gDC5J6buWMUhoOhbO9QNnRBAm+ku3CN1MURxRJY9ACXwTASREpxAssBb\nF0ZbanACJcMX0Mp5AKKIYgCEyT0R5K//GLI6T6NXn4t+z6TUpRDFDQZAmNwTQWRy/D9mplcFt17j\n3kBEkcIACBN7ALGV53OhYMyOlqv8NyeKFAZAmDwTQV4CGmPrnafx4cAYxv2C1KUQxQVZBoBSLgPl\nJaCxZXWeRkAQcaKfwUsUCbwMNAz+oIgxv8AhoBi7c+QS9IlatFwdxcZC3iuYaKFk2QOQO24EJw2d\nKGDdYj3aro5CEHmvYKKFYgCEwe2buvafARB7pQXpGPIFccE1IXUpRIrHAAgDt4GQjiVfDw2A1qu8\nHJRooRgAYeAQkHQMKTqsMqeihQFAtGAMgDC4uRW0pEoL9Ohx+TA0zm04iBYiZgHQ2dmJ559/Hj/4\nwQ/Q2dkZq7eNCg+HgCRlLUgHALRxVTDRgoR0GeihQ4fQ3t4Og8GAhoaGmeMdHR04cuQIBEFAZWUl\ntm/fPuc5NBoNUlJS4Pf7YTKZFl65hNwTQaQlapEoBqUuRZWWZyXDlKZDy9VRVBVnSV0OkWKFFADl\n5eWorq5GU1PTzDFBEHD48GHs378fJpMJ+/btg9VqhSAIOHr06A2v3717N+68804899xzGB4exk9/\n+lPs2bMnsi2JIfdEEJnJCVKXoVoajQalBel46+II/EEBiQkcySQKR0gBUFJSArvdfsOxnp4e5OXl\nITc3FwBQVlaGlpYW1NTUoK6ubs5zpaenw+/3L6Bk6bkngshgAEjKkq/H77qH0TU4jrV5eqnLIVKk\nsFcCu1yuG4ZyTCYTuru753z+8ePHceLECXi9XlRXV8/5PJvNBpvNBgCor6+H2WwOqz6dThf2a29n\nLNgLU3pqVM5Ntzb9mVZkGvHdd67h9JCAyrvNUf285Uyt7QbU2/ZItjvsABBnWYmp0WjmfP6GDRuw\nYcOG2563qqoKWVlZaGtrAwA4HI6w6jObzWG/9naGvD4s1nPYQQoDNWUz/33X2p145+0+fPZ//z1y\nX22O2uctZ9H8nsudWtseSrvz8/NDOlfYf8VMJhOcTufMY6fTCaPRGO7pbmC1WvHUU09F5FzRwDkA\nebjXeRZX9HkYTOZEMFE4wg6A4uJi9PX1wW63IxAIoLm5OWIbuMl5N9DJoABfQGQAyIDFdQYA8EH2\nJySuhEiZQhoCamxsRFdXFzweD3bt2oXa2lps3boVO3bswIEDByAIAioqKlBYWBiRouS8G+j0GoDM\nZFlupKoqS8bsMPuG0J79CfyD1MUQKVBIf8X27t0763GLxQKLxRLRgoCpHkBbW5ssh4HcMwHAHoDU\nNAAsrrP4S8698Ad5kxii+ZLlTKac5wAYAPJicZ6BT5eMk9fcUpdCpDiyDAA5m9kHKIUBIAdrhs9D\nJwTw3uUhqUshUhxZBoCcJ4HZA5CX1OAE7hy5hPcuMQCI5kuWM5lKmATOSGIAyIXFdRY/dd6BQa8f\ni/SJUpenwkLTAAAI00lEQVRDpBjsAcyTeyKA9CQtErRzL3qj2Lp3+nLQPm7PTTQf7AHME/cBkp+l\n3gHkpCeh7dooHr6Di8KIQiXLHoCccRWw/GgAbFxuxIm+MQQE3iyeKFQMgHliAMjTxmVGjAcEnBkc\nl7oUIsWQZQDIeQ7AMxFEBlcBy461MAsJGt4ljGg+ZPmXTO5zAOwByI8+WYe7ctLQfs2L/36v1NUQ\nKYMsewByNREQMBnkRnBytX6xHpeGJ+AcU/YNh4hihQEwD1wEJm+W/Kk7g7Vf4+WgRKFgAMwDA0De\nlmUlw5SqQxsDgCgkspwDkOtuoAwA+Zq+U9i6T3wGf/WswcTPvgydOLVDaMIPfy1laUSyJcsegFx3\nA3X7AgC4EZycWZxnMKZLxbnMpVKXQiR7sgwAuXLzZjCyd89wD7RiEB3GVVKXQiR7DIB5cE8EodUA\n+kT+s8mVPuDDSvcVdGQzAIhuh3/J5sEzEYQ+KYEbwcncOtc5nM9YAo8uTepSiGSNATAPXASmDOuG\nzkHUaHHSeIfUpRDJmiwHs+V2FVBw56MAAPfancjQ6hDcOfs9kkke7vD0Qu8fQ0f2Ktw/eFLqcohk\nS5YBINetINyJeuT6XFKXQbeRIAq4Z6gHHcZV4N6gRHPjENA8eBLTkO7nbpNKsG7oHJwpWehNy5G6\nFCLZYgDMw3hCCtICPqnLoBCsdZ0DAF4NRHQLDIAQiQB8CUlIDU5IXQqFIGdiGAVjdpzgegCiOTEA\nQuRLSIKo0TIAFGSt6xw+zCrCZFCQuhQiWWIAhGg8IQUAkMohIMVY5zqHyYQknOZdwohmFbOrgARB\nwCuvvILx8XEUFRWhvLw8Vm8dEeMJyQDAHoCCrB65AJ0QQEefF2vz9FKXQyQ7IQXAoUOH0N7eDoPB\ngIaGhpnjHR0dOHLkCARBQGVlJbZv3z7nOVpbW+FyuZCeng6TybTwymNsXDcVAGkMAMVIDU7izpHL\n+KBPz7uEEc0ipAAoLy9HdXU1mpqaZo4JgoDDhw9j//79MJlM2LdvH6xWKwRBwNGjR294/e7du3Ht\n2jWsWrUKDz30EBoaGrBmzZrItiTKZnoAAQaAkqwbOoufDRVjaDwAY6osl70QSSak/yNKSkpgt9tv\nONbT04O8vDzk5uYCAMrKytDS0oKamhrU1dXddI7s7GzodFNvp9Uqb+phOgBS2ANQlLWubvysCOjo\n86KiyCB1OUSyEvZPIpfLdcNQjslkQnd395zP37BhA3784x/jzJkzuOuuu+Z8ns1mg81mAwDU19cj\nPz8/3BIX9NobvNGKQgCfBgBsjcw5KSYKAbRIXUSURex7rkBqbXuk2h32T3FRvHmRvUYz9y6ZycnJ\n2L17N3bs2IHq6uo5n1dVVYX6+nrU19eHWxoAzNoLUQO2W13U2m5AvW2PZLvDDgCTyQSn0znz2Ol0\nwmg0RqQoIiKKvrADoLi4GH19fbDb7QgEAmhubpblBm5ERDS7hG984xvfuN2TGhsb8corr8DpdMJm\nsyEtLQ3FxcXIy8vDv/7rv+J3v/sdtmzZgo0bN8ag5NAVFRVJXYIk2G51UWu7AfW2PVLt1oizDeYT\nEVHcU971mEREFBEMACIilYrLpZHz2aJCyWbbomN0dBQHDx7E4OAgFi1ahH/6p39Cenq6xJVGlsPh\nQFNTE4aHh6HRaFBVVYVHHnkk7ts+OTmJr3/96wgEAggGg9i4cSNqa2tht9vR2NiI0dFRrFixAk8/\n/fTMost4IggC6urqkJ2djbq6OlW0+8tf/jJSUlKg1WqRkJCA+vr6yH7PxTgTDAbFr3zlK2J/f7/o\n9/vFr33ta+KVK1ekLisqOjs7xfPnz4vPPPPMzLGXXnpJfPXVV0VRFMVXX31VfOmll6QqL2pcLpd4\n/vx5URRFcWxsTNyzZ4945cqVuG+7IAji+Pi4KIqi6Pf7xX379olnz54VGxoaxHfeeUcURVF84YUX\nxN///vdSlhk1r7/+utjY2Ch+5zvfEUVRVEW7v/SlL4kjIyM3HIvk9zzuhoCu36JCp9PNbFERj0pK\nSm5K/paWFjz44IMAgAcffDAu2240GmeugkhNTUVBQQFcLlfct12j0SAlZWpb8mAwiGAwCI1Gg87O\nzpkr8MrLy+Ou3cDUOqP29nZUVlYCmFqIqoZ2zyaS3/P46i9h/ltUxJuRkZGZBXlGoxFut1viiqLL\nbrfj4sWLuOOOO1TRdkEQ8Oyzz6K/vx/btm1Dbm4u0tLSkJCQAGBqzy2XyyVxlZH34osv4vOf/zzG\nx6fu7eDxeFTRbgA4cOAAAOChhx5CVVVVRL/ncRcA4jy3qCDl8vl8aGhowBe+8AWkpaVJXU5MaLVa\nfPe734XX68X3vvc9XL16VeqSoq6trQ0GgwFFRUXo7OyUupyY+ta3voXs7GyMjIzg29/+dsT3Poq7\nAFD7FhUGgwFDQ0MwGo0YGhpCZmam1CVFRSAQQENDA7Zs2YINGzYAUE/bAUCv16OkpATd3d0YGxtD\nMBhEQkICXC4XsrOzpS4vos6ePYvW1lZ88MEHmJycxPj4OF588cW4bzeAmTYZDAaUlpaip6cnot/z\nuJsDUPsWFVarFceOHQMAHDt2DKWlpRJXFHmiKOL73/8+CgoK8KlPfWrmeLy33e12w+v1Api6IujU\nqVMoKCjA6tWr8d577wEA3nrrrbj7vn/uc5/D97//fTQ1NWHv3r24++67sWfPnrhvt8/nmxny8vl8\nOHnyJJYuXRrR73lcrgRub2/HT37yEwiCgIqKCjz22GNSlxQVjY2N6OrqgsfjgcFgQG1tLUpLS3Hw\n4EE4HA6YzWY888wzcXUpJACcOXMGzz//PJYuXTozvPf4449j5cqVcd32y5cvo6mpCYIgQBRFbNq0\nCZ/5zGcwMDBw0+WQiYmJUpcbFZ2dnXj99ddRV1cX9+0eGBjA9773PQBTk/6bN2/GY489Bo/HE7Hv\neVwGABER3V7cDQEREVFoGABERCrFACAiUikGABGRSjEAiIhUigFARKRSDAAiIpX6/3UoptgE239/\nAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10f0a4f98>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "values = poisson_sample_maximum_better(mu,N,Ntrials)\n",
    "\n",
    "### Make and plot the normalized histogram\n",
    "### We define the binning ouselves to have bins for each integer\n",
    "bins = np.arange(0, 10 * mu)\n",
    "histo = plt.hist(values, bins=bins, normed=True, log=True)\n",
    "\n",
    "### Now compare to the analytical solution\n",
    "from scipy.special import gammaincc\n",
    "\n",
    "### Define a lambda function to compute analytical solution\n",
    "proba = lambda nv, Nr, mu_p : gammaincc(nv + 1, mu_p) ** Nr - gammaincc(nv, mu_p) ** Nr\n",
    "\n",
    "x = 0.5 * (bins[:-1] + bins[1:])\n",
    "y = proba(bins[:-1], N, mu)\n",
    "plt.plot(x, y)\n",
    "plt.ylim(1e-6,1)   # restrict y range"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "## Exercices\n",
    "\n",
    "- write a vectorized function that takes an array of int and returns an array where square integers are replaced by their square root and the others are left unchanged"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "7.0\n",
      "[  0.   1.   2.   3.   2.   5.   6.   7.   8.   3.  10.  11.  12.  13.  14.\n",
      "  15.   4.  17.  18.  19.  20.  21.  22.  23.  24.   5.]\n"
     ]
    }
   ],
   "source": [
    "### A solution\n",
    "def replace_square(n):\n",
    "    sqrt_n = np.sqrt(n)\n",
    "    return n + (sqrt_n == sqrt_n.astype(int))*(-n + sqrt_n)\n",
    "\n",
    "print(replace_square(7.0))\n",
    "print(replace_square(np.arange(26)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "7.0\n",
      "[  0.   1.   2.   3.   2.   5.   6.   7.   8.   3.  10.  11.  12.  13.  14.\n",
      "  15.   4.  17.  18.  19.  20.  21.  22.  23.  24.   5.]\n"
     ]
    }
   ],
   "source": [
    "### or using where\n",
    "def replace_square2(n):\n",
    "    sqrt_n = np.sqrt(n)\n",
    "    return np.where(sqrt_n == sqrt_n.astype(int), \n",
    "                    sqrt_n, n)\n",
    "        \n",
    "print(replace_square2(7.0))       \n",
    "print(replace_square2(np.arange(26)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.3"
  },
  "nbsphinx": {
   "orphan": true
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}