courses/IMTx-Queue-Theory/Week4_Lab_Continuous_Time_Markov_Chains (1).ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<p><font size=\"4\">  MOOC: Understanding queues</font></p>\n",
    "<p><font size=\"4\">  Python simulations</p>\n",
    "<p><font size=\"4\">  Week IV: Continuous time Markov chains </p>\n",
    "\n",
    "\n",
    "----------------------------"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this lab, we focus on the simulation of continuous time Markov chains. In the lab of week 2 we have simulated a M/M/1 queue. This week, we are going to study a $M/M/1/K$ queue to illustrate the effects of a finite buffer length. We will compute the loss probability in this model and observe the influence of the load $\\rho$ when the capacity $K$ is large. \n",
    "\n",
    "-----------------"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**1)** Complete the code of the function **MM1K** below. This function generates one trajectory of a $M/M/1/K$ queue. The function returns the instants of events (arrivals or departures), the number of customers in the system at these instants, as well as the number of arrivals and of lost customers. Customers are lost if the buffer is full when they arrive. Default parameters will be set as follows: MMM1K($K=3$, $\\lambda = 4$, $\\mu = 5$, $N0 = 2$, $Tmax=100$). $\\lambda$ and $\\mu$ are the arrival and departure rates, $K$ is the maximum number of customers in the system, $N0$ is the initial number of customers, and the evolution of the number of customers in the system is simulated over $[0,Tmax]$. Plot a trajectory of the number of customers in the system against time, obtained after running function **MM1K** with the default parameters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline              \n",
    "from pylab import *          "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def MM1K(K=3,lambda_ = 4.,mu = 5.,N0 = 2,Tmax=100):\n",
    "    N0       = min(N0,K)# enforcing buffer length constraint\n",
    "    p        = lambda_/(mu+lambda_) # probability that the next event is an arrival when N(t) > 0\n",
    "    T        = [0]  # list of instants of events (arrivals/departures)\n",
    "    N        = [N0] # initial number of customers in the system, list of number of customers at arrivals/departures\n",
    "    losses   = 0    # number of lost clients\n",
    "    arrivals = 0    \n",
    "    while T[-1]<Tmax:\n",
    "        if N[-1]==0:\n",
    "            # inter-event when N(t)=0:\n",
    "            tau      = -1./lambda_*log(rand()) \n",
    "            event    = 1\n",
    "        else:\n",
    "            tau   = -1./(lambda_+mu)*log(rand()) # inter-event when N(t)>0\n",
    "            event = 2*(rand()<p)-1 # +1 for an arrival, -1 for a departure\n",
    "        # Unlike in function MM1, when N[-1]==K, if a new client arrives this client is lost \n",
    "        # and the number of lost clients is incremented by 1\n",
    "        if event==1:\n",
    "            arrivals +=1    # event==1 corresponds to an arrival\n",
    "            if N[-1]==K:\n",
    "                ###########################\n",
    "                # supply value of events and update losses\n",
    "                # when a customer arrives while N[-1]==K\n",
    "                event   = 0 #loss cliente does not increase N \n",
    "                losses  += 1 \n",
    "                ###########################\n",
    "        N = N + [N[-1]+event]\n",
    "        T = T + [T[-1]+tau]\n",
    "\n",
    "    T    = T[:-1] # event after Tmax is discarded\n",
    "    N    = N[:-1]\n",
    "    return array(T),array(N),arrivals,losses\n",
    "\n",
    "#------------------\n",
    "T,N,arrivals,losses = MM1K(K=3,lambda_ = 4.,mu = 5.,N0 = 2,Tmax=10**3)\n",
    "V1 = losses/arrivals"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA3IAAADiCAYAAAALWfU/AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xm8JGV56PHf4wHXgRkJLoCMqEGjcRl1NEYTZ0S9AiomxCQaY4xxuYkbJho1fkiMW1yiJhqNN1w3NMYlaIga3KIOCffeKIMOCG7BBUFQMBFhUEHG5/5RVU5NTVV1dfepc04zv+/ncz7ndPVb7/vUW2+93c+p6urITCRJkiRJi+N6qx2AJEmSJGk6JnKSJEmStGBM5CRJkiRpwZjISZIkSdKCMZGTJEmSpAVjIidJkiRJC8ZETpIkSZIWjImcJEmSJC0YEzlJ+6SI+EZEPGiV2r5DRHwuIq6MiGdMsd4eMUfEeRGxdZQgB1rOGFZzn2h2EfGyiHjmasdxXRARn4mIn1/tOCQtBhM5SWtC+Sb+OxFxk9qyJ0bEtlUMayzPAbZl5gGZ+bpZK8nMn8/MbfMGM08CNU8MYyZuZd3XRMTBjeU7IiIj4ojG8kMj4qJp122s97SI2B4RV0fE2zriarZTT8wfFRHfi4gtc2z6ioqImwG/A/zdCrR1UET8U0RcFREXRMRv9ZQ9IiJOK/vz2xHx+ojYb4q2blru61uXjyMiXlTus7vNsQ1/HxGXRMQVEfGViHhio8irgBfNWr+kfYuJnKS1ZD/ghNUOYhrTvDmsuTVw3nLHor18HXh09SAi7gLcqKPsscBHZli3vt7FwEuAt/TE1Gynqv9xwBuAh2bm6T3rrzW/C5yWmT9cgbbeAFwD3AJ4DPDGnrNXfwtcChwCbAK2AE+Zoq1NwPcy84Lyn0unAEcB987Ms2eMH+BlwBGZeSBwHPCSiLhn7fkPAA+IiEPmaEPSPsJETtJa8pfAsyNiQ/OJ8r/jP1t7/LaIeEnt8Tci4o8j4pzyP/ZvjohbRMSHy0sY/zUibtqo9l4R8YXyv/ZvjYgblnUdGhHvi4jLIuLrzcsfy7aeGxHnAFe1JXMRcceI2BYRl5eXHx5XLv8k8ADg9RGxMyJu37Lu4RHx/rL9/4qI17d1Vv2MzsCYn132z/cj4j0RccOIeAewEfhgGc9zyvLPjYhvlX335Yh44IAYWtvoWK+1XWBT1/qTtrHFOyjOFlUeB7y9o+yxwGkzrPvT9TLz/Zl5KvBfPTE12yEingy8GnhIZv7ftpUi4u4R8dlyf7wnIt7dGP+Tjo/OvptnXeAYYK/EMyIOiIiTymPr0oj4w54+mahMpn4N+NPM3JmZZ1AkPY/tWOU2wHsz80eZ+W2K5HmaSxY3ATsiYiNwBvB94KjMvHTmjQAy87zMvLp6WP7crvb8j4CzgP8xTzuS9g0mcpLWku3ANuDZM67/a8CDgdsDDwc+DDwfOJhivmu+8X8M8BCKN1K3B06MiOsBHwTOBg4DHgg8MyIe0lj30cBDgQ2ZeW39iYjYv6zjY8DNgacD74yIO2TmUcC/A0/LzHWZ+ZXGukvAh4ALgCPKGN7dt9FTxPwbwNEUb3LvCvxuZj4W+Cbw8DKeV0bEHYCnAffKzAPKPvpGXwx9bbQVamu3b/0ptrHuP4ADy6R6CfhN4O+bhcr9dX/g49Os27Fep47yfwC8GHhgZm7vWO/6wKkUyeVBwD9SjPVBZuy7oeveBfhyy6qnAl8Fbgk8CnhVRNyyUfeHyn90tP18qFHf7YFdjePlbLqTs9cCj4qIG0fEYRQJ515nQnvcHdifYhy8IzN/LzOvaRaachuqdf42In4AfAm4hEZiD3wRmPnyTUn7DhM5SWvNnwFPj+KzN9P6m8z8TmZ+iyJZ+nRmfq78D/g/Ubw5q3t9Zl6Ymf8NvJQiObsXcLPMfFFmXpOZXwP+N8Wb0brXleu2XVJ2H2Ad8PKyjk9SJGePbinbdG/gUOCPM/Oq8ozCGRPWmSbmi8vt/SDFWYc2u4AbAHeKiP0z8xuZ+dUBsU/TxrTrD93GpurM2oMp3jh/q6XM/YGzM/PKKdftWq9LW/kHUyQLn+9Z7z4UScVfZ+aPM/MU4MyBbcLsfTdk3Q3AHtsfEQ8DyMxXZObV5fj/FkUy9lOZ+bDM3NDx87BGHOsozorVfR84oCPu0ymSvCuAiyj+SXTqgO2tbCrXPyczX9NVaMptqNZ5Shn3LwPvB65uFLmSol8lqZeJnKQ1JTPPpUh6njfD6t+p/f3DlsfrGuUvrP19AUUCdWvg0Pp/1inO6t2iZ92mQ4ELM/MnjfoPm7wJHA5c0DzLN8HQmL9d+/sH7N0fAGTm+cAzgT8HLi0v4zt0YCyD2phh/aHb2PQO4LcozuwNvaxy6Lpd63VpK//7FAnOmyIiOtY7FPhWZmZt2QVTtDtr3w1Z93vsnUwdB/xz9aA8q7eePY/Hae0EDmwsO5BGEllr76MUSdJNKM7I3xR4xZCGIuIGwB0pznreMUa4I2dm7ir/QXMrirOydQcAly93m5Kue0zkJK1FLwCexJ6Jzw+AG9ce73GZ1owOr/29keJmFRcCX2/8Z/2AzDy2sW7S7WLg8PINZb3+trNBTRcCG2O6m6gMjbnLXtuSmf+Qmb9E8UY+GfgmeEp9fdg00zZm5gUUNy45luKNfZtjgX+ZYd3W9Xq0lb+U4nLFX6a4QUebS4DDGonexkaZvuNjUt/Ns+45NM60Ab/Anp8TPAr4bmbucQlmFJ9f3dnx8+FGnV8B9ouII2vL7kb7TYMOoji2X1+eEfwv4K0U/T/EnSnOSv878KvAiyPiqLaCU25Dm/2ofUaudEeKy0YlqZeJnKQ1pzwj9B72/EzbDuC3ImIpIo6muAvdvJ4aEbeKiIMozjK8B/gMcEUUN/u4UdnenSPiXlPU+2ngKuA5EbF/FN+z9nAmfNat9BmKN+4vj4ibRHFDkvsNWGeemL8D3LZ6EMX33B1Vnpn4EcXZzF0D65rGHu1OMM82PoHiRhVXNZ+IiNsAN8jML02zbtt6EbFfFDdnWQKWyn2336R2MvNiimTn6Ij4q5YY/h9wLfCMso3jKS7Bres7Pib13TzrnlYvX34O8EjgkeX2/zxFgvrclu0+pvx8ZNvPMY2yV1Ek0y8qj4v7AY+gOGvarPe7FAn4H5T9tYHiZjU/TY6iuKHL25rrlu4OnJuZ12bmZ4GnAu8t9+HM2xARN4/iKybWlf34EIrLrT9ZK3MD4J4M/NylpH2biZyktepFFJdFVU6gSIYup7hJyTSfd+nyDxQ3JPla+fOSzNxVtrOJ4s3gd4E3UVwaNkgWN0U4juIGC9+leCP7Oz3JQn3dqv2fpbgZyEUUN9oYss6sMb+M4kYvl0fEsyk+H/fysp5vU9yw5fkD65pGs91O82xjZn41O24kQnHDms7LI3vWbVvvRIqk93nAb5d/nziwnQspkrlHRsTLGs9dAxxPcYnn9yjGQ/MMYefxMaDv5ln37cCxEVF9NcMdKW6Mcy5Fon4q8NLyc33zegrFV0BcCrwL+IPM/OkZufLsWDVOj6e4ac5lwPkUiXD9zpmHA/+no51NFMktAJn5doq54tSofc/lDJLiMsqLKPbjq4BnZuY/18ocR/EdkxfP0Y6kfUTsecm9JEn7jog4jeISvGk+6zb1erO201Pf24CLMvPESWXHFhF/AVyamX8dEb8N/GpmDr6r5kqL4i6gZwN3zcwfr3Y8dRHxaeAJ5WeFJanXLF9kK0nSdcU24FMrsN6s7ax5mVk/W3s3itvnr1nlGc47rnYcbTLzF1Y7BkmLw0srJUn7rMx8ZbZ/hcSyrjdrOwvorhRf1SBJGpmXVkqSJEnSgvGMnCRJkiQtGBM5SZIkSVowa+pmJwcffHAeccQRqx2GJEmSJK2Ks84667uZebNJ5dZUInfEEUewfXvXV/1IkiRJ0nVbRFwwpJyXVkqSJEnSgjGRkyRJkqQFYyInSZIkSQvGRE6SJEmSFoyJnCRJkiQtGBM5SZIkSVowJnKSJEmStGBM5CRJkiRpwZjISZIkSdKCMZGTJEmSpAVjIidJkiRJC8ZETpIkSZIWjImcJEmSJC0YEzlJkiRJWjAmcpIkSZK0YEzkJEmSJGnBmMhJkiRJ0oIxkZMkSZKkBWMiJ0mSJEkLxkROkiRJkhaMiZwkSZIkLRgTOUmSJElaMCZykiRJkrRgTOQkSZIkacGYyEmSJEnSgjGRkyRJkqQFYyInSZIkSQtmtEQuIm4YEZ+JiLMj4ryIeOFYbUmSJEnSvmS/Eeu+GjgqM3dGxP7AGRHx4cz8jxHblCRJkqTrvNESucxMYGf5cP/yJ8dqb2xbt8K2bcPKVYaUnzWWSSa13RZnfRuHbu8k89SzXDEstzH2cdc+Xc7tb7YxzRipys9bR72e5rpdY3JoW33trqVxVG173/heibE/dhurefwuZ9tdY3hMfeN+6LEztJ2V2kf1cV9Z7ra76p60ncsxX1R1jLVNfXPFGO3OY4xxNaTOafbz2GNxiGn33UocO83X3mn6c95yzXab+3MtjvU1ITNH+wGWgB0UCd0rOso8GdgObN+4cWOuRVu2ZC4tFb8zM9evL36alpYyofi9fv3u8ssdS1V318+ktutl6nVV29jc3nlinbWe5YphubX1Xdd4mKXOafbjPG1MM0bq2zhPHdXxsX793vu2a0w2l01S1V39HvNYnEXVD2190CwzZsxjt7Gax++WLbvn4SH6jt+2Y3Ps8dQ37tv6dZbjpKuuWU2aA+vH/iyxDtHVD13joS2WWeOqz5FjbFNXvWO1O48xjv0hdU4q0/aaPdZYHKJt361fv/s46SrffJ+5XLHA7nqruvv6ZdK80+zfKu4h7xmqn/r+XItjfWzA9hyQa415aSWZuQvYFBEbgH+KiDtn5rmNMicBJwFs3rx5zZ6xW7duurKbNo0Xy6ZN0/13sa+Oetn6Nk6zvX3mqWe5YlhubX23XHXWLWf9zTamGSNV+R075qtjwwbYWZ6jb9u3bf06a19v2lTEW/29llTb3je+V2Lsj93GWj1+p9U8Npf7uOxrs62tocfOEKuxj8aYP+ete9p5rauOMUyqd63NbzDOuBpS56Qy9fFRfz1biWO6K55py+/Ysft1dDktLe3d1pB4hvTfNHFX7e7Ysff+XItjfS1YkbtWZublwDbg6JVoT5IkSZKuy8a8a+XNyjNxRMSNgAcBXxqrPUmSJEnaV4x5aeUhwMkRsUSRML43Mz80YnuSJEmStE8Y866V5wB3H6t+SZIkSdpXTXVpZUTcNCLuOlYwkiRJkqTJJiZyEbEtIg6MiIOAs4G3RsRrxg9NkiRJktRmyBm59Zl5BXA88NbMvCfFjUskSZIkSatgSCK3X0QcAvwG4M1KJEmSJGmVDUnkXgh8FDg/M8+MiNsC/zluWJIkSZKkLkPuWnlJZv70BieZ+TU/IydJkiRJq2fIGbm/GbhMkiRJkrQCOs/IRcQvAvcFbhYRf1R76kBgaezAJEmSJEnt+i6tvD6wrixzQG35FcAjxwxKkiRJktStM5HLzNOB0yPibZl5wQrGJEmSJEnqMeRmJzeIiJOAI+rlM/OosYKSJEmSJHUbksj9I/C/gDcBu8YNR5IkSZI0yZBE7trMfOPokUiSJEmSBhny9QMfjIinRMQhEXFQ9TN6ZJIkSZKkVkPOyD2u/P3HtWUJ3Hb5w5EkSZIkTTIxkcvM26xEIJIkSZKkYSZeWhkRN46IE8s7VxIRR0bEw8YPTZIkSZLUZshn5N4KXAPct3x8EfCS0SKSJEmSJPUaksjdLjNfCfwYIDN/CMSoUUmSJEmSOg1J5K6JiBtR3OCEiLgdcPWoUUmSJEmSOg25a+ULgI8Ah0fEO4H7Ab87ZlCSJEmSpG5D7lr58Yj4LHAfiksqT8jM744emSRJkiSpVeellRHxc+XvewC3Bi4BLgY2lsskSZIkSaug74zcs4AnAa9ueS6Bo/oqjojDgbcDtwR+ApyUma+dMU5JkiRJUqkzkcvMJ5W/HzBj3dcCz8rMz0bEAcBZEfHxzPzCjPVJkiRJkuhJ5CLi+L4VM/P9E56/hOJyTDLzyoj4InAYYCInSZIkSXPo+/qBh/f8PGyaRiLiCODuwKdnCXK17dix5++dO/d8fuvW4mfXrvb1+mzdOlss85Rre66+jc3tbWrGXG1/17bs2DFsO+t1TIqhr46xbN26dzw7dhTjoTkmpjFpX7X176Q+n6aNumb/1w3ZxrZ92FVPfdva+rVt2aS2d+7c+zistzVp/Xl17ZPmdrb1ZVdf9NU7a4zTHl9929QX17TH73JqjoOubZjl+G1uV1fds+gb9/X91jdepm1v3rE1tA+bZZaj7T71/tm1a88xW8U8b/9Nsw/mmberv8fsr8q0bbSVX65jvyuWej8298GQPm4bi0PaHWLa/dwW05B1q/Lzjon6tu/a1f4a2rdO/XEzlmmOjXrZvvlkNV9X1qTMHPUHWAecBRzf8fyTge3A9o0bN+ZatGVL8Xv9+uJnaWn3si1bimVbtuxeVqmW99Vbr2tIHOvXDyvb1Xazjupxtay+HW11NGOub3+zfFVX9dzSUnf8W7ZkQvFTj2dSH/bFNkm9nbbHzbq7+qPZB0P3Ub3evuebfVj1Y1dMQ9poG69V/3VtQ3P/Neuob/+kvqyvX23TpPj6trXvOGw+31SNu2nGTl8MbcdBM7b6fq0e12NoPrd+/Z7Hx7wx9h3jbes0+6bvuK+eHzo+J5nlmGq23bUNVZ/3zZdty4fU3TaeJm1L27FaHzvNebrtGJxmvlyO/dQ17rvaqqvG9LxjpN5GV0zNPuuKqd4nbfui7fWz2UZXfJPGW7N88/Wp2V/V3DnNfp9k2tfRtvcEbWNq2uO4XnfbvmzOic391vV+pP64rj6f1ft1qHqbQ16f28Zh1W5b+23ju+qbeV67muOs7XVq0jqZe78/bnvPN83rc9vr4DSvW9cFwPYckGdN/PqBiPgL4JWZeXn5+KYUn307ccC6+wPvA96ZHZdiZuZJwEkAmzdvzsmp58rbtq34vWlT8Z+Adev2fH7Tpt1lmssnadY1yZA6J5Vre65aVt+Orjq6tr/rv0LVc2ec0R0TwNLS3m0P3d6u2JZT236uHs/zH7G+bay3V+/DSX0+TRt1Vf+1lR/St9V6O3Z0t1ltU/V7w4buMm1197Xdtl5V36S+Wo6x0xVjM7a2eJaWihjantu0qdj3Q/5TOiTG+pw2RFvfDBmD0x6/y6nZdts2VH0+bZxD6p51PDXrrj9um6fnNc080mVo/y1n3EPUx/pYbffNO03Tjrfmvj/99KnDm8m0Y7drfhgzlvqc2JyLhrTd3GdDjumhpn19rrdZnYka0n5VZt7Xrvq2Dx3LbfNU1/vjofVNeo2cts59Rd+llZVjqiQOIDO/Bxw7aaWICODNwBcz8zWzhyhJkiRJqhuSyC1FxA2qBxFxI+AGPeUr9wMeCxwVETvKn4kJoCRJkiSp38RLK4G/Bz4REW+l+P643wNOnrRSZp4BxHzhSZIkSZKaJiZymfnKiDgHeBBFYvbizPzo6JFJkiRJkloNOSNHZn4E+MjIsUiSJEmSBhjyGTlJkiRJ0hpiIidJkiRJC6YzkYuIT5S/X7Fy4UiSJEmSJun7jNwhEbEFOC4i3k3jDpSZ+dlRI5MkSZIktepL5P4MeB5wK6D5hd4JHDVWUJIkSZKkbp2JXGaeApwSEX+amS9ewZgkSZIkST2GfI/ciyPiOOD+5aJtmfmhccOSJEmSJHWZeNfKiHgZcALwhfLnhHKZJEmSJGkVDPlC8IcCmzLzJwARcTLwOeBPxgxMkiRJktRu6PfIbaj9vX6MQCRJkiRJwww5I/cy4HMR8SmKryC4P56NkyRJkqRVM+RmJ++KiG3AvSgSuedm5rfHDkySJEmS1G7IGTky8xLgAyPHIkmSJEkaYOhn5CRJkiRJa4SJnCRJkiQtmN5ELiKuFxHnrlQwkiRJkqTJehO58rvjzo6IjSsUjyRJkiRpgiE3OzkEOC8iPgNcVS3MzONGi0qSJEmS1GlIIvfC0aOQJEmSJA025HvkTo+IWwNHZua/RsSNgaXxQ5MkSZIktZl418qIeBJwCvB35aLDgFPHDEqSJEmS1G3I1w88FbgfcAVAZv4ncPMxg5IkSZIkdRuSyF2dmddUDyJiPyAnrRQRb4mIS/36AkmSJElaXkMSudMj4vnAjSLiwcA/Ah8csN7bgKPniE2SJEmS1GJIIvc84DLg88D/BE4DTpy0Umb+G/Dfc0UnSZIkSdrLxESu/FLwk4EXU3wVwcmZOfHSyuuqnTuL3zt2zF7H1q3FT1VHW11VmXm01du2rNqmrjrqcdRjbouvXn+zrV27irYmbVsznqF9Palcvc3q73pbk/phSPs7d+6up207m8umHUe7dhU/zXarOiftk6Z6PFW5vn4YajnqaNM27qY5Tqr1+/qrWt42XrrUj+d6O0PjGXqswp7H0JC4JtVZ75O2vm2bp7rmrL7nqzJdcQ+ZE4bs+6F92VzWPI669kuz/KQyfbEOiXM51fu/73Wn77iYpq3q96Tt6jteJh2HzfmrMs8cVJ/L+7Z5yLjq6u/mNjf/HjIW6q8FfcdcM5b68q52m/NC15zYtl1tMdRfG/u0jZvm+66+vqm3MWnOahqyLX3t7tw5uY62ftuxY8/H1X6tlrWNh2o7296Tzvvesa4t/iHH886de79XadY3zWtkVz3zvldezr5aNZnZ+wM8FLgQ2AacDnwTOGbSeuW6RwDnTijzZGA7sH3jxo251m3ZUvxev774u/rpKru0VJStL6uvW6+rWaZteVdbbbrqbYuzSz3epaVM6I8RinJtbS0t7a4Duretvnz9+j3LVm32xVmPrf58FVd9vzSXNfuiKtPWZlsM9bab/dPss6790aUeb3NZ1afN5/v2eTPO6mdpaVg8zfHYtQ+XU30/V4+r/TOkzbZ9VD+OqzHafK4+puv9WR8fzb6o749Jc0Rbn7fNMfXyzb7oG/PV466xUO+X+vbU66+2pz4P1OttO/ab+6et3qrNtuOvWUdze+pjv7mNze1q9mOzL7rGRZfmmOk6vtvGVld9k/bPpDHeLFMfE825vD7mqvFa7cO246lrG9riq9dX1Tkp7mabzXHRNkfXx0bb/pxV29ipYmr2Z1Nzfurqg7bxUR/TbX3WnKvqx2HfPqvWqWKpH+dVHzfbbe7Xtnjb5tDm8VfV27WP2+asttfOevnmOs1tbvZJ175qao65Ia8rzTmzWb7r9arrGKo/V/1d32/N7az3TX0/9x1zzXHct431OXvS2GyOhbZ628Z8X9m28VDfhua81baNXdrmlLUE2J4Dcq0hXwj+auABmXk+QETcDvgX4MPLlEieBJwEsHnz5jV/pm/btuL3pk2z17Fp0+56uuqap/5p663H0vZc/T8WS0t7Lm/+Z2ap9g2DzbbWrSt+f//73e21xbM08FsLJ/VZ1X7977Zls9QNu+Ou+mVI/0+7n9et27t/1q3b/d+vtm3oa6N6rop3x47+fhiqb0zNW2/zP2jTxNu2j+rLKl1jt0vzmN62DfYbMLtu2wYbNnTX2Ve+rS+amnFPOv63boUzztizfPV8/TicdEx2jbmqvra4+/q47XitHjfnk6HzXnNZfU5r7s821fMbNnSfaai3UR2nfZZj3p9Uf9+YgyLOrjiGbEPT0tLk46dvLE+ao6tlbftzVs1josuQcbW01D43N7e5bWy3xdX2etwVS7Oeevn6MVXt0/q+auu/tniHvKYtLQ2br5r1VjHU33c1X1ubcXa9txpybG3atPf8N1TXnNH1WtK2DVX/VM9V21uNobq2vhlD23gecjwPqW/Ia3c1HzXfa9b7Zh7L8X5ntQ35jNylVRJX+hpw6UjxSJIkSZIm6PyfcUQcX/55XkScBryX4msHfh04c1LFEfEuYCtwcERcBLwgM988d8SSJEmStI/ru/jn4bW/vwNsKf++DLjppIoz89FzxCVJkiRJ6tCZyGXm41cyEEmSJEnSMBM/jh8RtwGeTnEHyp+Wz8zjxgtLkiRJktRlyF0rTwXeDHwQ+Mm44UiSJEmSJhmSyP0oM183eiSSJEmSpEGGJHKvjYgXAB8Drq4WZuZnR4tKkiRJktRpSCJ3F+CxwFHsvrQyy8eSJEmSpBU2JJH7VeC2mXnN2MFIkiRJkia73oAyZwMbxg5EkiRJkjTMkDNytwC+FBFnsudn5Pz6AUmSJElaBUMSuReMHoUkSZIkabCJiVxmnr4SgUiSJEmShpmYyEXElRR3qQS4PrA/cFVmHjhmYJIkSZKkdkPOyB1QfxwRvwLce7SIJEmSJEm9hty1cg+ZeSp+h5wkSZIkrZohl1YeX3t4PWAzuy+1lCRJkiStsCF3rXx47e9rgW8AjxglGkmSJEnSREM+I/f4lQhEkiRJkjRMZyIXEX/Ws15m5otHiEeSJEmSNEHfGbmrWpbdBHgC8DOAiZwkSZIkrYLORC4zX139HREHACcAjwfeDby6az1JkiRJ0rh6PyMXEQcBfwQ8BjgZuEdmfm8lApMkSZIktev7jNxfAscDJwF3ycydKxaVJEmSJKlT3xeCPws4FDgRuDgirih/royIK1YmPEmSJElSU99n5PqSPEmSJEnSKhk1WYuIoyPiyxFxfkQ8b8y2JEmSJGlfMVoiFxFLwBuAY4A7AY+OiDuN1Z4kSZIk7St671o5p3sD52fm1wAi4t3AI4AvjNjmitqxo/i9aVN3mV27YOdO2Lp19zpt5XfsaC/TtXyaGOdZv17Prl2wtLTn8p0tt8DZtWtYW9PEVpWt2qvWqz9f77NmuWpfVb+r2Net23NZPaaqzK5d/dvRpt4vVZ3NGGH2/dFnaL/WY2zbj8vZ1nKr9ynMFn/X2G22U/+7Oa4mjY+hx0KXIeOk69ioj/nm+OtTzVlt5att7ZoH2ra3q7563FW9zeOvvh3147W+7ytt29jsvyF9Mc+xUK3bdqzD3q8HzfW7Yppmvqj3YX1MTDsOm/uiiqHvNa1r38zSZnNctJWZ5/ifJa6h/TlNPzS3bcjrTdfrcdeU79ogAAAIjklEQVQ+a8ZRHUvV767jelK8zdc52Pv4a8bZNT7rdfTtz2n39bSvtW2vAW3Haz2eSftsmter5nPTbG9zrumKe9p5oZrb2sbcLLE2xzxMjqGt7vqy5rzU9V6xK45FFpk5TsURjwSOzswnlo8fC/xCZj6tUe7JwJMBNm7ceM8LLrhglHjGUA2Qbdsml6lrK98sV5XpWj7UvOs362muv3Vr+7IhbU0TW9/B2Fy/q+y2bXvGVsXeXNZX91D1fqnX2Vw26/7oa3NovzZjnCWW5Rpf87Rb35fT1jFp7A4ZG/XybW30PT/JrHNMpRn3kDj6yk+aB9qe76qvbexMOnbbtqdvDHbtz2a5uuU+FoaOoUkx9T3fF8eQ+qv1+vZPtX7fa9q8x+WkNrv296z7bJ64qrYnlR/6Oth1/Exap6vdevttx0HXa+Akba+hXTFNmg/aTNqfs46poev0zRl9hs7R027fNNs7zVxTN817tWnmkEn1TTvm+95/TrO9dSv1nmUWEXFWZm6eWG7ERO7XgYc0Erl7Z+bTu9bZvHlzbt++fZR4JEmSJGmtG5rIXW/EGC4CDq89vhVw8YjtSZIkSdI+YcxE7kzgyIi4TURcH3gU8IER25MkSZKkfcJoNzvJzGsj4mnAR4El4C2Zed5Y7UmSJEnSvmLMu1aSmacBp43ZhiRJkiTta8a8tFKSJEmSNAITOUmSJElaMCZykiRJkrRgTOQkSZIkacGYyEmSJEnSgjGRkyRJkqQFYyInSZIkSQvGRE6SJEmSFoyJnCRJkiQtGBM5SZIkSVowJnKSJEmStGBM5CRJkiRpwZjISZIkSdKCMZGTJEmSpAVjIidJkiRJC8ZETpIkSZIWjImcJEmSJC0YEzlJkiRJWjAmcpIkSZK0YEzkJEmSJGnBmMhJkiRJ0oIxkZMkSZKkBWMiJ0mSJEkLxkROkiRJkhaMiZwkSZIkLRgTOUmSJElaMCZykiRJkrRgIjNXO4afiojLgAtWO44WBwPfXe0gdJ3l+NLYHGMak+NLY3OMaUxrcXzdOjNvNqnQmkrk1qqI2J6Zm1c7Dl03Ob40NseYxuT40tgcYxrTIo8vL62UJEmSpAVjIidJkiRJC8ZEbpiTVjsAXac5vjQ2x5jG5PjS2BxjGtPCji8/IydJkiRJC8YzcpIkSZK0YEzkekTE0RHx5Yg4PyKet9rxaPFFxOER8amI+GJEnBcRJ5TLD4qIj0fEf5a/b7rasWpxRcRSRHwuIj5UPr5NRHy6HF/viYjrr3aMWlwRsSEiTomIL5Vz2S86h2m5RMQflq+P50bEuyLihs5hmkdEvCUiLo2Ic2vLWuesKLyufO9/TkTcY/Uin8xErkNELAFvAI4B7gQ8OiLutLpR6TrgWuBZmXlH4D7AU8tx9TzgE5l5JPCJ8rE0qxOAL9YevwL4q3J8fQ94wqpEpeuK1wIfycyfA+5GMdacwzS3iDgMeAawOTPvDCwBj8I5TPN5G3B0Y1nXnHUMcGT582TgjSsU40xM5LrdGzg/M7+WmdcA7wYescoxacFl5iWZ+dny7ysp3gAdRjG2Ti6LnQz8yupEqEUXEbcCHgq8qXwcwFHAKWURx5dmFhEHAvcH3gyQmddk5uU4h2n57AfcKCL2A24MXIJzmOaQmf8G/Hdjcdec9Qjg7Vn4D2BDRByyMpFOz0Su22HAhbXHF5XLpGUREUcAdwc+DdwiMy+BItkDbr56kWnB/TXwHOAn5eOfAS7PzGvLx85lmsdtgcuAt5aX774pIm6Cc5iWQWZ+C3gV8E2KBO77wFk4h2n5dc1ZC/X+30SuW7Qs8xafWhYRsQ54H/DMzLxitePRdUNEPAy4NDPPqi9uKepcplntB9wDeGNm3h24Ci+j1DIpP6f0COA2wKHATSgudWtyDtNYFuo100Su20XA4bXHtwIuXqVYdB0SEftTJHHvzMz3l4u/U526L39fulrxaaHdDzguIr5BcTn4URRn6DaUlymBc5nmcxFwUWZ+unx8CkVi5xym5fAg4OuZeVlm/hh4P3BfnMO0/LrmrIV6/28i1+1M4MjyTknXp/iw7QdWOSYtuPLzSm8GvpiZr6k99QHgceXfjwP+eaVj0+LLzD/JzFtl5hEUc9YnM/MxwKeAR5bFHF+aWWZ+G7gwIu5QLnog8AWcw7Q8vgncJyJuXL5eVuPLOUzLrWvO+gDwO+XdK+8DfL+6BHMt8gvBe0TEsRT/zV4C3pKZL13lkLTgIuKXgH8HPs/uzzA9n+Jzcu8FNlK8kP16ZjY/mCsNFhFbgWdn5sMi4rYUZ+gOAj4H/HZmXr2a8WlxRcQmipvpXB/4GvB4in8MO4dpbhHxQuA3Ke7y/DngiRSfUXIO00wi4l3AVuBg4DvAC4BTaZmzyn8gvJ7iLpc/AB6fmdtXI+4hTOQkSZIkacF4aaUkSZIkLRgTOUmSJElaMCZykiRJkrRgTOQkSZIkacGYyEmSJEnSgtlvchFJkhZTRPwM8Iny4S2BXcBl5eMfZOZ9VyUwSZLm5NcPSJL2CRHx58DOzHzVasciSdK8vLRSkrRPioid5e+tEXF6RLw3Ir4SES+PiMdExGci4vMRcbuy3M0i4n0RcWb5c7/V3QJJ0r7MRE6SJLgbcAJwF+CxwO0z897Am4Cnl2VeC/xVZt4L+LXyOUmSVoWfkZMkCc7MzEsAIuKrwMfK5Z8HHlD+/SDgThFRrXNgRByQmVeuaKSSJGEiJ0kSwNW1v39Se/wTdr9WXg/4xcz84UoGJklSGy+tlCRpmI8BT6seRMSmVYxFkrSPM5GTJGmYZwCbI+KciPgC8PurHZAkad/l1w9IkiRJ0oLxjJwkSZIkLRgTOUmSJElaMCZykiRJkrRgTOQkSZIkacGYyEmSJEnSgjGRkyRJkqQFYyInSZIkSQvGRE6SJEmSFsz/B8IQn7eDl3EnAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7ff3b8fc85c0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Plotting the evolution of the number of clients in the system\n",
    "\n",
    "def step(x,y,Tmax=0,color='b'):\n",
    "    # plots a step function representing the number\n",
    "    # of clients in the system at each instant\n",
    "    if Tmax==0:\n",
    "        Tmax = max(x)\n",
    "    x = append(x,[Tmax])  # number of clients\n",
    "    y = append(y,[y[-1]]) # instants of events\n",
    "    for k in range(len(x)-1):\n",
    "        vlines(x[k+1],y[k],y[k+1],color=color)\n",
    "        hlines(y[k],x[k],x[k+1],color=color)\n",
    "\n",
    "K=3\n",
    "T,N,arrivals,losses = MM1K()\n",
    "rcParams['figure.figsize'] = [15,3]\n",
    "step(T,N)\n",
    "xlabel('Time')\n",
    "ylabel('Number of clients')\n",
    "lambda_ = 4.\n",
    "mu = 5.\n",
    "title('Number of clients in the M/M/1/K queue' \n",
    "      +r'($\\rho =%g$, $K=%g$)'%(lambda_/mu,K))\n",
    "axis(ymin=-1,ymax=4)\n",
    "yticks(range(4), range(4));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "**2)** Letting $K$ range from 1 to 11, plot the loss probability for $\\lambda = 4$ and for $\\lambda = 10$ (and $\\mu=5$). Remarks ? Compare it to the theoretical loss probability.\n",
    "\n",
    "Observe on the curves that when $\\rho<1$ ($\\lambda=4$), the blocking probability of the $M/M/1/K$ queue tends to 0 as $K$ increases since the system tends to behave as a stable $M/M/1$ queue. \n",
    "\n",
    "When $\\rho>1$ ($\\lambda=10$), the rate of arrivals exceeds that of departures and the corresponding $M/M/1$ queue is unstable. So, even if $K$ is large the loss probability of the $M/M/1/K$ queue does not tend to zero.\n",
    "\n",
    "The loss probability in a M/M/1/K queue is:\n",
    "$$\n",
    "\\pi_K=\\dfrac{1-\\rho}{1-\\rho^{K+1}}\\rho^K\n",
    "$$\n",
    "\n",
    "Clearly, when $\\rho<1$, $\\pi_K$ tends to 0 as $K$ tends to infinity, whereas $\\pi_K$ tends to $(\\rho-1)/\\rho$ when $\\rho>1$ and $K$ tends to infinity."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA20AAADgCAYAAACQNI5RAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAIABJREFUeJzs3Xl8leWd///X5yzZQxICYUlIwr4FJRDQzhStgkvFaqH9OtrpqIPUdqzf1nY61o4zHevPUWsdW9tav6OtaGuVjgoaNyy1tTLUskhQNgOoIAn7vmQ/5/r9cU4OJ8kJOYGEE8j72eZx7nPf133f1wlp4Z3PdV+XOecQERERERGRnsmT6A6IiIiIiIhI+xTaREREREREejCFNhERERERkR5MoU1ERERERKQHU2gTERERERHpwRTaREREREREejCFNhEREcDM7jKzp0/y3BvN7H9PcPx1M7shVlszO2pmw07mviIi0jv4Et0BERE5s5nZFmCuc+4Pie5LT+Wc++wJjmU0b5vZk0CVc+7fTke/RETkzKBKm4iInPUsRH/niYjIGUl/gYmISLcxs6+Y2WYz229m5WY2OLzfzOzHZrbbzA6Z2ftmVhI+doWZrTezI2ZWbWbfaefaN5rZUjP7WfgaH5jZ9Kjjb5nZf5rZUqAGGGZmg8P92B/u11daXTbFzH4XvvcqMzs36np3mNmH4WPrzWxW2y6dsC9z2/kczsxGmNnNwN8Dt4eHTL5sZv9iZi+0av8zM/tJR997ERE5eyi0iYhItzCzi4H7gGuAQcBWYH748KXABcAoIBv4O2Bf+NivgK865zKBEuCPJ7jNecBHQD/gP4AFZtY36vg/ADcDmeH7PwtUAYOBLwL3Rocr4GrgOaAv8Azwopn5w8c+BKYBWcAPgKfNbFAn+nJCzrnHgN8CDzjnMpxznwOeBi43s2wAM/MR+l79Jt7riojImU+hTUREusvfA08451Y55+qB7wGfMrNioJFQkBoDmHNug3NuR/i8RmCcmfVxzh1wzq06wT12Az9xzjU6534HVAIzo44/6Zxb55xrAgYCnwa+65yrc86tBn5JKNg1e9c597xzrhF4CEgBzgdwzj3nnNvunAuG77UJmNqJvnRa+HvyNvB/wrsuB/Y65949leuKiMiZRaFNRES6y2BC1S0AnHNHCVXT8p1zfwR+DjwC7DKzx8ysT7jpF4ArgK1m9mcz+9QJ7lHtnHNR77eG79tsW6v+7HfOHWnVPj9We+dckONVOczsejNbbWYHzewgoSpgv0705WQ9BXw5vP1lVGUTEel1FNpERKS7bAeKmt+YWTqQC1QDOOd+6pybDIwnNEzyX8L7VzjnrgbygBeB/znBPfLNzKLeF4bv2yw6RG0H+ppZZqv21VHvh0T11wMUANvNrAh4HLgVyHXOZQNrgeh7d9SXeLgY+14Ezgk/83cloSGUIiLSiyi0iYhIV/CbWUrUl4/QM2H/aGYTzSwZuBdY5pzbYmZTzOy88PNix4A6IGBmSWb292aWFR6ieBgInOC+ecA3zMxvZv8HGAu8Fquhc24b8BfgvnAfzwFuomUImmxms8P9vw2oB/4KpBMKVHsAzOwfCVXaTqovJ7ALaLFmm3OuDnie0PdzuXPuk05eU0REznAKbSIi0hVeA2qjvu5yzr0J/DvwArADGA5cG27fh1Dl6gChYYT7gAfDx/4B2GJmh4GvcXxoYCzLgJHAXuA/gS865/adoP11QDGhCthC4D+cc4ujjr9EaKKPA+F+zA4/o7Ye+C/gHULBagKw9BT7EsuvCD3Pd9DMXoza/1T4nhoaKSLSC1nL4fciIiJnBjO7kdCi3p9OdF+6m5kVAh8AA51zhxPdHxEROb1UaRMREenBws/WfRuYr8AmItI7+RLdAREREYktPHnLLkJDSC9PcHdERCRB4hoeaWaXAw8DXuCXzrn7Wx3/MXBR+G0akBeeWUtEREREREROQYehzcy8wEbgEkLr1awArgs/lB2r/f8FSp1zc7q4ryIiIiIiIr1OPM+0TQU2O+c+cs41APOBq0/Q/jrg2a7onIiIiIiISG8XzzNt+cC2qPdVwHmxGoYXHx0K/LGd4zcDNwOkp6dPHjNmTKc6KyIiIiIicrZ499139zrn+nfULp7QZjH2tTem8lrgeedczIVQnXOPAY8BlJWVuZUrV8ZxexERERERkbOPmW2Np108wyOrgCFR7wsILUoay7VoaKSIiIiIiEiXiSe0rQBGmtlQM0siFMzKWzcys9FADvBO13ZRRERERESk9+owtDnnmoBbgTeADcD/OOfWmdndZnZVVNPrCC382fEaAiIiIiIiIhKXuBbXds69BrzWat/3W72/q+u6JSIiIiIiIhDf8MhucXD7xyx6+WE27NtAQ6AhUd0QERERERHp0TpcXLu7jE9Ndb8dXszd13n5eIif4qxiRvcdzeic0NeovqPol9ovIX0TERERERHpbmb2rnOurKN2cQ2P7A7mINl5+I73claUDKFyfyUrdq7g1Y9ejbTJTcmNBLlRfUcxOmc0xVnF+D3+RHVbRERERETktEpYaMPA/ElMvOzv+ZvS0sjuA3UH2HhgI5X7K6k8UMnGAxt5esPTNAYbAfB7/IzIHsGonFHHK3N9R5OVnJWoTyIiIiIiItJtEjY8cmJRkfvLiy+SFhXY2tMYbOTjQx9Tub+yRaDbX7c/0mZA2oA2VbnCzEK8Hm93fgwREREREZGT0uOHR/r6948rsEGoujYqZxSjcka12L+3dm8kwDUHuqXVSwm4AAAp3hRG5oxsUZUblTOKjKSMLv88IiIiIiIi3SFhlbaysjK3cuXKLr9ufaCeDw9+GAlxGw9spPJAJYfqD0Xa5GfktxxemTOa/Mx8PJawyTRFRERERKSX6fGVtu6S7E1mXO44xuWOi+xzzrGrZleLoZWV+yt5a9tbOEKhNd2fzsjskYzuOzoS6EZmjyTNn5aojyIiIiIiInL2Vdo6o7apls0HNrcYXrnxwEaONh4FwDAK+xSGQlx4wpPROaMZmD4QM0to30VERERE5MzWayttnZHqS2VC/wlM6D8hss85R/XR6tDMlftDQys37NvA4q2LI236JPVp+Zxc31GMyB5Bsjc5ER9DRERERETOYr06tMViZhRkFlCQWcD0wumR/UcbjrLp4KYWSxEs2LSA2qZaALzmpbhPcWTmyuZhlv1T+6sqJyIiIiIiJ02hLU4ZSRmU5pVSmnd8xsugC7LtyLbjQW7/RlbvXs3rH78eaZOTnNMiyI3OGc2wrGH4vVogXEREREREOtarn2nrLofqDx2fuTIc6DYf2ExDsAEAn8fHsKxhLSpyo/uOpm9K3wT3XERERERETpd4n2lTaDtNmoJNbD289fjsleHK3J7aPZE2/VP7H6/KhQNdUZ8ifB4VREVEREREzjZdOhGJmV0OPAx4gV865+6P0eYa4C7AAe85577UqR6f5XweH8OzhzM8ezhXcEVk//66/ZGZK5sD3bIdy2gKNgGhJQyGZw9vUZUblTOKrOSsRH0UERERERE5jTqstJmZF9gIXAJUASuA65xz66PajAT+B7jYOXfAzPKcc7tPdN3eVmnrjMZAIx8d+qhFkNt4YCP76/ZH2gxKHxSZuXJ0TijMFfYp1ALhIiIiIiJniK6stE0FNjvnPgpfeD5wNbA+qs1XgEeccwcAOgpscmJ+rz80aUnf0Xxu+OeA0FIEe2v3RtaUax5euaR6CQEXAEJLGIzMHtlmBst0f3oiP46IiIiIiJyCeEJbPrAt6n0VcF6rNqMAzGwpoSGUdznnFrW+kJndDNwMUFhYeDL97bXMjP5p/emf1p9P5386sr8+UM/mg5sja8pV7q/kjS1v8PzG5yNtCjIKWqwpNzpnNPkZ+VqKQERERETkDBBPaIv1L/vWYyp9wEjgM0ABsMTMSpxzB1uc5NxjwGMQGh7Z6d5KG8neZMbnjmd87vjIPuccO4/tbFmVO7CRP37yR1z4jy7DnxF5Pq450I3IGUGqLzVRH0VERERERGKIJ7RVAUOi3hcA22O0+atzrhH42MwqCYW4FV3SS+kUM2NQxiAGZQziM0M+E9lf01gTWSC8+Xm58g/LqamsAcBjHgozCyMhrnl45YC0AarKiYiIiIgkSDyhbQUw0syGAtXAtUDrmSFfBK4DnjSzfoSGS37UlR2VU5fmT+Pc/udybv9zI/uCLkj10eoWwyvX7l3LG1veiLTJSs6KTHbSXJkbkT2CJG9SIj6GiIiIiEiv0mFoc841mdmtwBuEnld7wjm3zszuBlY658rDxy41s/VAAPgX59y+7uy4dA2PeRiSOYQhmUOYXjQ9sv9IwxE2HdgUCXIbD2zk+Y3PUxeoA8BrXoZmDW0xvHJ039H0S+2XqI8iIiIiInJW0uLaErdAMMAnRz6JzFzZHOh21eyKtOmb0rfF0MrRfUczNGsofo8/gT0XEREREel54p3yX6FNTtnBuoOhZ+SiqnKbD26mMdgIgN/jZ3j28FCICwe60TmjyU7JTnDPRUREREQSR6FNEqox2MjWQ1tDQS6qMre3dm+kTV5aXosQN6rvKIoyi/B6vAnsuYiIiIjI6dGVi2uLdJrf42dEzghG5IxgJjMj+/fV7ms5vPJAJe9sf4cm1wRAijeFEdkjGNV3VIvKXGZSZqI+ioiIiIhIQqnSJgnXEGjgo0MfHV9TLhzoDtYfX+ZvcPrgyMLgzZW5gswCPOZJYM9FRERERE6eKm1yxkjyJjGm7xjG9B0T2eecY3fN7sjC4M2B7u2qtwm6IABpvjRG5oxsMfHJqJxRpPnTEvVRRERERES6nCptckapbarlw4MfHl8gPFyZO9J4BADDGJI55PjsleFANyh9kBYIFxEREZEeRZU2OSul+lIp6VdCSb+SyD7nHDuO7Tg+vDJcmVu8dXGkTWZSZpvZK4dnDyfFl5KIjyEiIiIiEjeFNjnjmRmDMwYzOGMwFxVeFNlf01jDxgMbWwyvXLh5IbVNtUBoYfHiPsWRmStH5Yxi2LYmUt7fRNrUqaSVlibqI4mIiIiIRCi0yVkrzZ/GxLyJTMybGNkXdEGqjlRF1pSrPFDJe3ve4/UtrzOyyvH9ZwP4AhDwGfNuLmLf8H6k+lJJ86eR5kuLvKb6U1u8b3GsVXu/x6+hmSIiIiJy0hTapFfxmIfCPoUU9inkkqJLIvsPNxzm4589iD/4POYcFnBM3OZj6egUjjUdY0/tHmqbaqlprKGmqYb6QH3c9/SZr23Iiwp7zcdihcNYr6m+VFJ9qZo5U0RERKSXUGgTAfok9WHkxbP45Dcv4xob8fr9XHPd3dzYzhDJQDAQCnFNNZEgF3ltqqG2se2x1u131+xuEQRrmmoiM2PGozm8nSjcRV7bqQq2fvV7/V31LRURERGRLqLQJhKWVlpK4bwnqFm+grSpU074TJvX4yUjKYOMpIwuu79zjvpAfbshL7IvRkhsPna08WibMNipqqDH134QbBUG460OpvhSVBUUEREROQUKbSJR0kpLEzYBiZmR4kshxZdC35S+XXbdpmBTm4pedCiMVRVs/doVVcEThbwW+0/wfKCqgiIiItIbxRXazOxy4GHAC/zSOXd/q+M3Aj8CqsO7fu6c+2UX9lNETpLP4yMzKZPMpMwuu2brqmD0a3MIPFFVsKaphqMNoapg5Lym2k5XBdt7PjBmGGwvLJ5BVcGaioq4KsEiIiJydukwtJmZF3gEuASoAlaYWblzbn2rpr9zzt3aDX0UkR4mIVXB1uEwVigMVwVb73O4uPvQ0XOC7VUFTxQKu6IqWFNRwSf/OAfX0IAlJVE47wkFN+kx9AsFEZHuFU+lbSqw2Tn3EYCZzQeuBlqHNhGRU9JdVcG6QF3s4Nd6aGg7zxMeaTjCrppdLc5vCDZ06nO1u0RE1MyiJwyMf34N19AAwSCusYFd//tH0kYNxswwLLKsRPN25D/N+1vtM1rtjzoe+m+r/Vq2QtqhXyhIT6dfKsipCroggWCAxmAjARcgEAzQ5JpoCjYRcIHQa/S+YICAa9U+2ESTa4psB1wA3/qPGODzDYynD/GEtnxgW9T7KuC8GO2+YGYXABuBbznntsVoIyJyWplZpELWlVXBxmBjJAhGh7z2KoCxqoS7anbFXRUceczxfU8Qn4MmC3Ln4Xlsev7JLvs88epMwGu3baxQGd724ImrXbvhszNBNXp/jM/Rom0Hnznm9+Ekvz8e85z8Z46zbVd85uZ7Fry8nMKGeizoCDbU85dXH2d3yqfxmhePxxN6tfZfm786at+iXYz27Z0Tva1fPvQ++qVC9+muINPuvhjHYx1r3Y927xOrX1HXjL52Z0bsxKt5feBcry8/nvbxhLZY/w/XuucvA8865+rN7GvAU8DFbS5kdjNwM0BhYWE8/RMR6ZH8Hj/+JD99kvp02TWbq4Kxloqobaxl97nrSX5/M8dKivnSmCE4HM65yLmR/4S3W+yPta+Dtg5H6L+t9jffs722ra7fpm0HfYm7zzH61qnP0U6fY7UNBoOn9plP8GcV3Tbogqf2mbvw+xOvkV7H9z0u9AsFj+MX9jabli/pxE/+6XOicBczWHraBkgzaz9Yxgid7V671T08xD4/VnvDFIrjVLN8RdQohcZQxa0bQ1tXBpno470lyMTDa1685sXn8eH1ePGZL7LtNS9+jz/UxhNq4zNfZNvv8+P1ePGbP9Le5wmf33zNVvvaXMdOcN/oa0bdN9Y13VPP0xB8KmbQiiWe0FYFDIl6XwBsj27gnNsX9fZx4IexLuScewx4DKCsrCwxf9IiIj1UdFUwl9y2DYqmw5Wnv1/Se8Ub8OoufY/alSvxl5Uy75wSAi4Q+cdr0AWPv2/12vwVT3vnXNv97bRvc91W7Vu3jbUdcAGCwTj64AI0BBvatA8SbHN+zHsFw21bfZ5E/YP4VPTEUNwvYx9/6zOsyQh64dnU99m99N8VZHpokGnuQ7v3NF+P/OXAyaiZdgmfPPEsDhfXD425DtqZmY/QkMfphGaHXAF8yTm3LqrNIOfcjvD2LOC7zrnzT3TdsrIyt3Llynj6KCIiItKrtBcQW+zvJaE4VvvOhOKh2xoYuzXA5qEpfFKUGleQaR0iFGSkO9RUVDBs6tTqnY2NBR217bDS5pxrMrNbgTcITfn/hHNunZndDax0zpUD3zCzq4AmYD9w4yl9AhEREZFezMzwmZbTFTmbpZWWsqupaWc8bTustHUXVdpERERERKQ3M7N3nXNlHbXruavIioiIiIiIiEKbiIiIiIhIT6bQJiIiIiIi0oMptImIiIiIiPRgmpZIRERERKQHaGxspKqqirq6ukR3RWJISUmhoKAAv99/2u+duNB2dBdsWw5DpiasCyIiIiIiPUVVVRWZmZkUFxdr7bUexjnHvn37qKqqYujQoaf9/okbHnl4Bzx1VSi4iYiIiIj0cnV1deTm5iqw9UBmRm5ubsKqoAl8ps1BoAG2LElcF0REREREehAFtp4rkX82CQxtBt4kKJ6WuC6IiIiIiIj0cAkLbfs8ubxa+v/Yl3NuorogIiIiIiLS4yVsIpL9nhy+vsSPb+mbfGZ0f2aVFjB9bB4pfm+iuiQiIiIiItLjJKzSNjIvg0W3TeOmTw9lTfUhvv7MKqb85x/43oL3Wf7xfoJBl6iuiYiIiIj0WmvWrKGoqIhHH3200+cGAgFKS0u58soru6FnvVdCF9ceM7AP37tiLH+5Yzq/uWkql4wdwEurt3PNf7/DBT/6Ew/9vpKP9x5LZBdFRERERHqVCRMmMH/+fH796193+tyHH36YsWPHdkOvereEhrZmXo8xbWR/Hvq7iay4cwYPXXMuQ/ul87M/beaiB9/i848s5TfvbOHAsYZEd1VERERE5KyXl5fHunXrOnVOVVUVr776KnPnzu2mXvVecT3TZmaXAw8DXuCXzrn722n3ReA5YIpzbuXJdCg92cfsSQXMnlTAzkN1vLS6moUV1fz7S+u4+5X1fGZ0HrNL87l4bB7JPj3/JiIiIiJnnx+8vI712w936TXHDe7Df3xufFxt77jjDurr69m6dStFRUUATJs2jSNHjrRp++CDDzJjxgxuu+02HnjggZht5NR0GNrMzAs8AlwCVAErzKzcObe+VbtM4BvAsq7q3MCsFL564XC+euFw1m8/zIJVVbz03nYWr99FVqqfmecM4guT8plUmKM1LUREREREusCiRYs4duwYM2fOZN26dZHQtmRJ++srv/LKK+Tl5TF58mTeeuut09TT3iOeSttUYLNz7iMAM5sPXA2sb9Xu/wMeAL7TpT0MGze4D+MGj+OOz45h6Yf7WLiqigWrqnhm2ScU5abx+Yn5zJ6UT1FuenfcXkRERETktIm3ItbV6urquP322ykvL2fevHmsXbuWK664AjhxpW3p0qWUl5fz2muvUVdXx+HDh/nyl7/M008/fbo/wlkpntCWD2yLel8FnBfdwMxKgSHOuVfMrFtCWzOf18OFo/pz4aj+HK1vYtHanSysqOKnf9zEw29uYnJRDrNK87nynEFkpyV1Z1dERERERM4q99xzD9dffz3FxcVMmDCB8vLyyLETVdpmzJjBfffdB8Bbb73Fgw8+qMDWheIJbbHGHUbm4zczD/Bj4MYOL2R2M3AzQGFhYXw9PIGMZB9fnFzAFycXsONQLS9WbGfBqir+7cW13P3yei4ek8esSflcNDqPJF+PmHNFRERERKRHqqysZPHixSxduhQIzSJ57733JrhXAmDOnXg9NDP7FHCXc+6y8PvvATjn7gu/zwI+BI6GTxkI7AeuOtFkJGVlZW7lypOaq+SEnHOs236YBauqKX+vmr1HG8hO83PlOYOYPamA0iHZev5NRERERHqcDRs2aLr8Hq6r/4zM7F3nXFlH7eKptK0ARprZUKAauBb4UvNB59whoF/Ujd8CvnOys0eeKjOjJD+Lkvws/vWKMSzZvJcFq6p5bmUVT//1E4b2S+fzE/OZVZpPYW5aIrooIiIiIiIStw5Dm3OuycxuBd4gNOX/E865dWZ2N7DSOVd+4iskjs/r4aLReVw0Oo8jdY28vnYnC1ZV8eM/bOTHf9jIlOIcZpUWMHPCILLS/InuroiIiIiISBsdDo/sLt01PDIe1QdrebGimgWrqvhwzzGSfB5mjM1jVmkBF47qr+ffREREROS00/DInq8nD4886+Rnp/L1i0Zwy2eGs6b6UPj5t+28tmYnOWl+rjp3MLMmFXBuQZaefxMRERERkYTqlaGtmZlxTkE25xRkc+fMsby9cQ8LKqp5dsU2nnpnK8P6pTN7Uj6fL82nIEfPv4mIiIiIyOnXq0NbNL/Xw/SxA5g+dgCHaht5fc0OFlRU8+DvN/Lg7zcydWhfvjApn89OGESfFD3/JiIiIiIip4dCWwxZqX6unVrItVML2ba/hhcrqllYUc13X1jD919ax4xxA5hdms8Fo/rj9+r5NxERERER6T4KbR0Y0jeN/zt9JLdePIL3qg6xYFUVL7+3nVff30FuehKfO3cwsyflMyFfz7+JiIiIiEjXU5koTmbGxCHZ3H11Ccv+dQaPX1/G1KF9eWbZJ1z186Vc8uO3eeRPm6k+WJvoroqIiIiIdNq+ffuYOHEiEydOZODAgeTn50feb9y4kZKSkkR3MWLNmjUUFRXx6KOPdvrcQCBAaWkpV155ZTf0rHuo0nYSknweLhk3gEvGDeBQTSOvrtnBwooqfvRGJQ/+vpLzh+Yya1I+ny0ZSKaefxMRERGRM0Bubi6rV68G4K677iIjI4PvfOc7AGzZsqXL7+ecwzmHx9P5OtKECROYP38+3/72t/mnf/qnTp378MMPM3bsWA4fPtzp+yaKKm2nKCvNz5fOK+S5r/0Nb//LRdw2fRQ7DtVy+/PvM+U//8A3nq3gT5W7aQoEE91VERERETnbbFsOS/4r9NrNAoEAX/nKVxg/fjyXXnoptbWhEWZPP/00U6dOZeLEiXz1q18lEAgA8NBDD1FSUkJJSQk/+clPgFD4Gzt2LLfccguTJk3ipptu4uGHH47c48477+SnP/1pXP3Jy8tj3bp1nfoMVVVVvPrqq8ydO7dT5yWaKm1dqDA3jW/OGMk3po9g1ScHWVhRxSvv76D8ve30y0jmqvDzb+MH99HzbyIiIiJyarYth6eugkADeJPghnIYMrXbbrdp0yaeffZZHn/8ca655hpeeOEFJk+ezO9+9zuWLl2K3+/nlltu4be//S3jx49n3rx5LFu2DOcc5513HhdeeCE5OTlUVlYyb948fvGLX7BlyxZmz57NN7/5TYLBIPPnz2f58vgC6B133EF9fT1bt26lqKgIgGnTpnHkyJE2bR988EFmzJjBbbfdxgMPPBCzTU+m0NYNzIzJRTlMLsrh+1eO50+Vu1mwqorf/HULTyz9mFEDMphVWsDnSwczKCs10d0VERERkTPRliWhwOYCodctS7o1tA0dOpSJEycCMHnyZLZs2cLBgwd59913mTJlCgC1tbXk5eVx4MABZs2aRXp6OgCzZ89myZIlXHXVVRQVFXH++ecDUFxcTG5uLhUVFezatYvS0lJyc3M77MuiRYs4duwYM2fOZN26dZHQtmTJknbPeeWVV8jLy2Py5Mm89dZbp/KtOO0U2rpZks/DZeMHctn4gRysaeCV93ewYFUVP1z0AQ+88QF/MzyXWaUFXF4ykIxk/XGIiIiISJyKp4UqbM2VtuJp3Xq75OTkyLbX66W2thbnHDfccAP33Xdfi7bNwyFjaQ5yzebOncuTTz7Jzp07mTNnTof9qKur4/bbb6e8vJx58+axdu1arrjiCuDElbalS5dSXl7Oa6+9Rl1dHYcPH+bLX/4yTz/9dIf3TDQ903YaZacl8eXzi1hwy9/y1nc+wzcuHsm2/bV857n3mHLPH7htfgV/3riHQNAluqsiIiIi0tMNmRoaEnnxnd0+NLI906dP5/nnn2f37t0A7N+/n61bt3LBBRfw4osvUlNTw7Fjx1i4cCHTpsUOlbNmzWLRokWsWLGCyy67rMW1q6ur27S/5557uP766ykuLmbChAmsXbs2cmzJkiWsXr26zdeMGTPBubiYAAAgAElEQVS47777qKqqYsuWLcyfP5+LL774jAhsoEpbwhT3S+dbl4zithkjeXfrARZUVPPKe9t5cfV28jKTuXriYGaVFjBucJ9Ed1VEREREeqohUxMS1pqNGzeOe+65h0svvZRgMIjf7+eRRx7h/PPP58Ybb2Tq1FDf5s6dS2lpacxZKJOSkrjooovIzs7G6/UCEAwG2bx5M3379m3RtrKyksWLF7N06VIgNIvkvffe270fsgcw5xJT1SkrK3MrV65MyL17qrrGAH/6YDcLKqp5q3I3jQHHmIGZzJ6Uz9UT8xnQJyXRXRQRERGRbrJhwwbGjh2b6G6cdsFgkEmTJvHcc88xcuRIANauXcsTTzzBQw89lODetdTVf0Zm9q5zrqyjdnENjzSzy82s0sw2m9kdMY5/zczWmNlqM/tfMxt3Mp3u7VL8Xj47YRCPX1/Gsn+dwd1XjyfF7+Xe1z7gU/e9yT/8ahkLK6qoaWhKdFdFRERERE7Z+vXrGTFiBNOnT48ENoCSkpIeF9gSqcNKm5l5gY3AJUAVsAK4zjm3PqpNH+fc4fD2VcAtzrnLT3RdVdri99GeoyysqGZhRTVVB2pJS/Jy+fiBzJqUz98M74fXo+UDRERERM50vbXSdiZJVKUtnmfapgKbnXMfhS88H7gaiIS25sAWlg5oJo0uNKx/Bv986Wi+NWMUK7ceYMGqKl5ds4MFFdUM6JPM5yfmM2tSPmMG6vk3EREREZGzTTyhLR/YFvW+CjivdSMz+zrwbSAJuDjWhczsZuBmgMLCws72tdfzeIypQ/sydWhf7rpqPG9uCK3/9qv//Zj/fvsjxg3qw+xJ+Vw1cTB5mXr+TURERETkbBBPaIs19q5NJc059wjwiJl9Cfg34IYYbR4DHoPQ8MjOdVWipfi9zDxnEDPPGcS+o/W8/N52FlZUc8+rG7j3tQ1MG9mf2ZPyuXTcQFKTvInuroiIiIiInKR4QlsVMCTqfQGw/QTt5wOPnkqnpHNyM5K58W+HcuPfDmXz7qMsrKjixYrtfHP+atKTQpObzC7N5/xhuXj0/JuIiIiIyBklntC2AhhpZkOBauBa4EvRDcxspHNuU/jtTGATkhAj8jL4l8vG8M+XjGbZx/tZWFHFa2t28vy7VQzKSuHqifnMnpTPqAGZie6qiIiIiIjEocPQ5pxrMrNbgTcAL/CEc26dmd0NrHTOlQO3mtkMoBE4QIyhkXJ6eTzGp4bn8qnhufzgqhIWb9jFwlVVPL7kI/7fnz+kJL8Ps0oLuOrcwfTPTE50d0VEREREpB1aXLuX2XMk9Pzbgooq1lYfxusxLhjZj1mTCrh03ABS/Hr+TURERCQRNOV/z9eTp/yXs0j/zGTmfHoocz49lE27jrCgopoXK6r5xrMVZCT7uGLCQGaVFnDe0L56/k1EREREpAfwJLoDkjgjB2Ty3cvHsPS7F/PM3PO4vGQgr76/g+se/yvTHvgTP3rjAzbvPproboqIiIjIabRmzRqKiop49NHOzS04Z84c8vLyKCkpaXNs0aJFjB49mhEjRnD//fd3VVd7DQ2PlBZqGwL8fv1OFqyqZsmmPQQdnFOQxezSfD537mByM/T8m4iIiEh36EnDI9955x2+/e1v884778R9zttvv01GRgbXX389a9eujewPBAKMGjWKxYsXU1BQwJQpU3j22WcZN25cd3S9W2l4pPQIqUlerp6Yz9UT89l9uI7y97azYFU1d728nnte3cCFo/oze1IB08fm6fk3ERERkbNUXl4e69at69Q5F1xwAVu2bGmzf/ny5YwYMYJhw4YBcO211/LSSy+dkaEtURTapF15fVKYO20Yc6cN44Odh1m4qpoXV1fz5ge7yUzxMXPCIGaV5jOlWM+/iYiIiHSlHy7/IR/s/6BLrzmm7xi+O/W7cbW94447qK+vZ+vWrRQVFQEwbdo0jhw50qbtgw8+yIwZM9q9VnV1NUOGHF/2uaCggGXLlnWy972bQpvEZczAPnzvij7cfvkY3vlwHwtWVVH+3nbmr9hGQU4qs0rzmVWaz7D+GYnuqoiIiIicgkWLFnHs2DFmzpzJunXrIqFtyZIlJ3W9WI9jmekX/p2h0Cad4vUYnx7Zj0+P7Mc9DU28sS70/Nsjf9rMz/64mYlDspk9KZ8rzxlM3/SkRHdXRERE5IwUb0Wsq9XV1XH77bdTXl7OvHnzWLt2LVdccQVw8pW2goICtm3bFnlfVVXF4MGDu77zZzGFNjlpaUk+ZpUWMKu0gF2H63hpdTULVlXz/ZfWcffL67loTB6zS/O5eGweyT49/yYiIiLS091zzz1cf/31FBcXM2HCBMrLyyPHTrbSNmXKFDZt2sTHH39Mfn4+8+fP55lnnumqLvcKCm3SJQb0SeHmC4Zz8wXDWb/9MAsrqnhx9XYWr99FnxQfM88ZzBcm5TO5KEflcBEREZEeqLKyksWLF7N06VIAJkyYwL333hv3+ddddx1vvfUWe/fupaCggB/84AfcdNNN+Hw+fv7zn3PZZZcRCASYM2cO48eP766PcVbSlP/SbZoCQZZ+uI+Fq6p4Y90uahsDFPZN4/Ol+cwuzae4X3qiuygiIiLSY/SkKf8lNk35L2cdn9fDhaP6c+Go/hytb2LR2p0srKjiZ3/cxE/f3MSkwmxmTSrgc+cMIjtNz7+JiIiIiMSi0CanRUayjy9OLuCLkwvYcaiWl1ZvZ8GqKv79xbXc/fI6Lh6Tx6zSAi4a01/Pv4mIiIiIRFFok9NuUFYqX7twOF+9YBjrth9mYUU1L62u5o11u8hO83PlOYOYVVrApMJsPf8mIiIiIr2eQpskjJlRkp9FSX4W3/vsGJZs3svCVdU8t7KKp//6CcW5aeHZKfMpzE1LdHdFRERERBIirtBmZpcDDwNe4JfOuftbHf82MBdoAvYAc5xzW7u4r3IW83k9XDQ6j4tG53GkrpHX1+5k4apqfvyHjfz4DxspK8ph9qQCZk4YRFaaP9HdFRERERE5bToMbWbmBR4BLgGqgBVmVu6cWx/VrAIoc87VmNk/AQ8Af9cdHZazX2aKn2vKhnBN2RCqD9byYkU1Cyuq+deFa7irfB3Tx+Yxe1IBF47qT5LPk+juioiIiIh0q3gqbVOBzc65jwDMbD5wNRAJbc65P0W1/yvw5a7spPRe+dmpfP2iEdzymeGsrT7MC6uqePm97by+dic5aX4+d+5gZk8q4NyCLD3/JiIiIiJnpXhCWz6wLep9FXDeCdrfBLwe64CZ3QzcDFBYWBhnF0VCz79NKMhiQkEWd84cy5JNe3hhVTXzV2zj1+9sZVi/dGaV5vP50nyG9NXzbyIiIiJy9ohnbFms8kXMFbnN7MtAGfCjWMedc48558qcc2X9+/ePv5ciUfxeDxePGcAjX5rEyn+bwQ+/MIF+mcn81+KNTHvgT1zz3+8wf/knHK5rTHRXRURERM4Y+/btY+LEiUycOJGBAweSn58feb9x40ZKSkoS3cWINWvWUFRUxKOPPtqp8+bMmUNeXl6bz7Jo0SJGjx7NiBEjuP/++9s5O3HiCW1VwJCo9wXA9taNzGwGcCdwlXOuvmu6J3JifVL8/N2UQv7nq59iye0X8Z1LR7H3aD13LFhD2T1/4OvPrOLNDbtoDAQT3VURERGRHi03N5fVq1ezevVqvva1r/Gtb30r8j4pKanL7+ecIxg8uX+jTZgwgfnz5/PrX/+6U+fdeOONLFq0qMW+QCDA17/+dV5//XXWr1/Ps88+y/r169u5QmLEE9pWACPNbKiZJQHXAuXRDcysFPhvQoFtd9d3U6RjQ/qmcevFI3nz2xfy4tf/luumDOGdD/dx01MrOf/eN7mrfB3vVx3EuZiFYhEREZEzTk1FBXv/+zFqKiq6/V6BQICvfOUrjB8/nksvvZTa2loAnn76aaZOncrEiRP56le/SiAQAOChhx6ipKSEkpISfvKTnwCwZcsWxo4dyy233MKkSZO46aabePjhhyP3uPPOO/npT38aV3/y8vJYt25dpz7DBRdcQN++fVvsW758OSNGjGDYsGEkJSVx7bXX8tJLL3Xqut2tw2fanHNNZnYr8AahKf+fcM6tM7O7gZXOuXJCwyEzgOfCk0F84py7qhv7LdIuM2PikGwmDsnm364cx58r97Cgoopnln3Ck3/ZwvD+6cyeVMDnS/PJz05tefK25bBlCRRPgyFTE/MBREREROJQU1HBJ/84B9fQgCUlUTjvCdJKS7vtfps2beLZZ5/l8ccf55prruGFF15g8uTJ/O53v2Pp0qX4/X5uueUWfvvb3zJ+/HjmzZvHsmXLcM5x3nnnceGFF5KTk0NlZSXz5s3jF7/4BVu2bGH27Nl885vfJBgMMn/+fJYvXx5Xf+644w7q6+vZunUrRUVFAEybNo0jR460afvggw8yY8aMmNeprq5myJDjAwsLCgpYtmzZSXyHuk9c67Q5514DXmu17/tR27G/AyIJ5vd6mDFuADPGDeBQTSOvrtnBwooqfvRGJT96o5Lzh/VldmkBn50wkMw9FfDUVRBoAG8S3FCu4CYiIiI9Vs3yFbiGBggGcY2N1Cxf0a2hbejQoUycOBGAyZMns2XLFg4ePMi7777LlClTAKitrSUvL48DBw4wa9Ys0tPTAZg9ezZLlizhqquuoqioiPPPPx+A4uJicnNzqaioYNeuXZSWlpKbm9thXxYtWsSxY8eYOXMm69ati4S2JUuWdPpzxRqF1dNmJY8rtImcDbLS/HzpvEK+dF4hn+yrYWFFNQsrqrj9hff595fW8qOBb3JlUz0egqHgtmWJQpuIiIj0WGlTp2BJSbjGRszvJ23qlG69X3JycmTb6/VSW1uLc44bbriB++67r0Xb5uGQsTQHuWZz587lySefZOfOncyZM6fDftTV1XH77bdTXl7OvHnzWLt2LVdccQVwcpW2goICtm07Pll+VVUVgwcP7rAfp5NCm/RKhblpfHPGSL4xfQQV2w6ycFU1z79XxCXOh9+aCDgvP9+YR/2RDQzKSmFwdiqDs1IZnJ1C3/SkHvfbFxEREel90kpLKZz3RKjCNnVKt1bZ2jN9+nSuvvpqvvWtb5GXl8f+/fs5cuQIF1xwATfeeCN33HEHzjkWLlzIb37zm5jXmDVrFt///vdpbGzkmWeeaXHtX//61+Tn57dof88993D99ddTXFzMhAkTKC8/Pt3GyVTapkyZwqZNm/j444/Jz89n/vz5LfrREyi0Sa9mZkwqzGFSYQ4NV45j1V9GcGDdH/lrYCx/3l/A9o+30NDUclajZJ8nEuQGZaWSn53CoOxUBmWlkJ+dyqDsVDKS9T8tERER6X5ppaUJCWvNxo0bxz333MOll15KMBjE7/fzyCOPcP7553PjjTcydWpo1NLcuXMpLS1ly5Ytba6RlJTERRddRHZ2Nl6vF4BgMMjmzZvbTBpSWVnJ4sWLWbp0KRCaRfLee++Nu7/XXXcdb731Fnv37qWgoIAf/OAH3HTTTfz85z/nsssuIxAIMGfOHMaPH3+S35HuYYmaSa+srMytXLkyIfcWiZdzjn3HGthxsI7th2rZfrCWHYfqqD5Yy47w9q7DdQRb/c+oT4ovHOrCVbrsUJVuUFaoYjcwK4UkXzyTt4qIiEhvsWHDBsaOHZvobpx2wWCQSZMm8dxzzzFy5EgA1q5dyxNPPMFDDz2U4N611NV/Rmb2rnOurKN2KgeInICZ0S8jmX4ZyUwoyIrZpjEQZPeRerYfPB7qQtuh19XbDnKgprHVdaFfRnJ42GU4zGVHBbysFPplJOPxaBimiIiInL3Wr1/PlVdeyaxZsyKBDaCkpKTHBbZEUmgTOUV+r4f87NS2ywdEqW0IsP1Qbahid7D2+PahWjbuOsJblXuobQy0uq4xMBzo8sNVu0HZ4eGYWaFw1yfFp+frRERE5Iw1btw4Pvroo0R3o8dTaBM5DVKTvAzvn8Hw/hkxjzvnOFTbGKnO7ThUS/XBOnaEw93yj/ez63AdTa3GYaYneRkUVZ0b3OrZukFZKaT4vafjI4qIiIhIN1FoE+kBzIzstCSy05IYN7hPzDaBoGPPkfo2FbvmIZnrtx9m79H6NuflpicxKDslPPtl9HN2ode8zBS8GoYpIiIi0mMptImcIbye0HDJgVkpUBi7TV1jgF2HmydKqWtRsdu6r4Z3PtzHkfqmNtcdkBl6vm5Qc5jLajmJSk6aX8MwRURETgPnnP7O7aESNYEjKLSJnFVS/F6KctMpyk1vt83husbI83StK3bvVx3kjbV1NARaLnOQ4veEglxzoGvxbF3oNV3LHIiIiJySlJQU9u3bR25uroJbD+OcY9++faSkpCTk/vpXlkgv0yfFT5+BfkYPzIx5PBgML3MQDnLbw5W67eGg9/amPew+Uk/rXzZlpfqPz4aZHb0geahqNzArBb9XyxyIiIi0p6CggKqqKvbs2ZPorkgMKSkpFBQUJOTeCm0i0oLHY/TPTKZ/ZjLnFGTHbNMYCLLrcF0k0LUejvnuJwc4GGOZg7zM5BPOhpmbnqRlDkREpNfy+/0MHTo00d2QHkihTUQ6ze/1UJCTRkFOWrttahqaoqp0tVEzY9axYedh3vxgF3WNLYdhJnk9DMxKiTxXF3rOruV2nxR/d388ERERkR4lrtBmZpcDDwNe4JfOuftbHb8A+AlwDnCtc+75ru6oiJxZ0pJ8jMjLYERe+8scHKxpDFXpDrWs2G0/WMuyj/ez83AdgVbLHGQk+yLP0R0fjnl8ApWBWuZAREREzjIdhjYz8wKPAJcAVcAKMyt3zq2PavYJcCPwne7opIicfcyMnPQkctKTKMnPitkmEHTsPlLXYv266Ird2upD7DvW0Oa8fhlJLSZJyQ9X6Zq3+2cma5kDEREROWPEU2mbCmx2zn0EYGbzgauBSGhzzm0JHwvGuoCIyMnweoxBWakMykplclFOzDZ1jQF2HmqeBbNluPtozzH+d9NejjUEWpzj8xgD+qRE1qprDnjNM2TmZ6eSlaplDkRERKRniCe05QPbot5XAed1T3dERDonxe+luF86xf1iL3PgnONwXRM7wkscVIdDXfN2xScHee3QDhoDLYdhpvq9kQAXWbOuedmD8HZqkoZhioiISPeLJ7TF+lXzSa0sZ2Y3AzcDFBa2szqwiEgXMjOyUv1kpfoZM7BPzDbBoGPvsfqodetaVuwqd+5hz9G2yxzkpPmPV+miK3bhoDegTxcvc7BtOWxZAsXTYMjUrruuiIiI9GjxhLYqYEjU+wJg+8nczDn3GPAYQFlZWeKWFBcRieLxGHmZKeRlpnDukNjLHDQ0NS9zEHqerjoq1FUdqGXFlgMcqm25zIHHIC8zNAxzUHjSlOZw1/ycXW56UnzDMLcth6eugkADeJPghnIFNxERkV4intC2AhhpZkOBauBa4Evd2isRkR4myedhSN80hvRtf5mDY/VNkbXqdrSq2K3ffpg/rN9FfVOrZQ58ntAMmM2zYUZmxjxesctM8YcqbIEGcIHQ65YlCm0iIiK9RIehzTnXZGa3Am8QmvL/CefcOjO7G1jpnCs3synAQiAH+JyZ/cA5N75bey4i0sOkJ/sYkZfJiLzMmMedcxyoaQyvW1cbqdo1h7t3PtzLzsN1tFrlgMwUH9PTM/khPnw4Avj41ZZB7Ni/Fr/XE/6ydreTfB58nvB+n4ckrwef5/h2zPN9Hvzhc7we06QsIiIiCWSu9UMap0lZWZlbuXJlQu4tItJTNQWC7D5S37JiF67aZe1dxaja91hp46lwo2gMBGkKOBoCQRoDwTbP3HUVs9CC6v5w0PN7w8EvKuglhbeb9zWHQZ/XjgdDn+HzeEjytR8y21w3fE7ztt/rCV8j+vy211PIFBGRM4GZveucK+uoXVyLa4uIyOnh83rCwyRTmVzU+mjo/9NvjnGec45A0NEUDIe4piCNAUdjONDF2m4Ih77m/Q3hc5qCx7dDwTBIQ7vXavW+yXG0qSmy3RgI0hg8vh19z6bWJcWu/D56jDaVRp+1EzhD7XxR26HAeTyMhiqPUdvhY77oimb4nOZtv88T6UdS+ByfxyLb0UFT6waKiMiJKLSJiJwFzAyf1/B5Q8sgnAmCQRcKdAEXDobhANgUDAfH+IJiQ/j85v0N4fObt5vbRVclo8+vbQxwuK51+AzSGHTHt8PndhdPuJrZttJ4PPiFhrOGK4++qEDZuqLZuvLos8hQV3+MwBhdFY28b66Keo+HXX9U8PVpyKyIyGml0CYiIgnh8RjJHi/JZ8jfRM3VzOYA19RO1bJFpTMYX9UzduXThcPn8QDbFK5a1tY2RoJo5J4xrh3oxmpmyyGpoRCZkuQl1R/6Sgl/pSZ5SfV7QttR+1J8ntBrq/3Hz/WErpXkJcXnxaNqpIj0YmfIX5UiIiKJFV3NTOXMq2Y2Nlcdg8e324TNqKplQ3TVM+r8pmDLCmZzUKxvClLXGAh/hSqY+481RPbVhr/qGk+uYpnk80QCYWqSl+Rw6GsbEj0t3zcHP//x/ScKib6uXFtRRKSLKLSJiIicpVpUM5O78UadWPjdOUd9U5DahuggFw52DcHI+8hrw/EQGHnfFIicX9cY4EhdU+haDcevVdMYOKnJefxeaxn4/F5SoiqDzfuSI4Gw/ZDY4r3P26KyqAlzRKQzFNpERETk5HVy4Xez46Eopxu75VxoGGtdQ7BFyKttDFAXCX7ByL76xlZtGoORkNgcDg8ca2B79PGGUDg8mWGoHiMq4J04JLYMhJ6WobGDkJjs8ygcipwFFNpERETk5PXQhd/NjGSfl2Sflyz83XqvxkDweHWwVbUwVkisi6oaHg+NwUi18HBtI7ujr9EQoK4pNKvryUiNGv6ZEhUEYz1j2DIketqGykggbPmcYorfq1lQRbqRQpuIiIicvOJpoQpbc6WteFqie3TaNU/GkpnSveEwEHQtglx9q2ph6zDYZlhp9LDT8Ffzc4cth6We2nOHLZ4fjKr6xZxkJp4qos9LStLxZxr13KH0RgptIiIicvKGTA0NiYzzmTY5eV6PkZ7sI72bp1wNBl1kYpnaTobE1kNKm193H2mMPJ8Yfd2Tee7Q57FWVcPYIfF4JfD4WolJMZa3iLlUxonaavkLSQCFNhERETk1Q6YqrJ1FPB4LVcaSTt9zh7Wtq32tqoWRYaYtZiINH48aZnrgWENU0AxS3xhgdNMGprKet4NjWeVGdfnnSAoHPl/U8hfNayL6PBYJjC1DocVeJ9EXu13rNRRDbTsfOBUyz1wKbSIiIiJy2p2W5w63LYen7seFh+/WXreQ+oGT213yosVSFjGWxogsh9G8HmMwajt8rCnqGs1rKDaFt4/WN4XbNy/HEd5udZ+m07zGoi86SMYZOKNDps9z4sDZHDLjDZzHw6lCZjOFNhERERE5O4UnyrHwRDlp2/9C2vBPJbpXHWqzxmKw/SDZEBUyW6+r2NgqjDYHzqboUNnUtl1jVMhseZ9QyGwKBsPrNZ7ukHni7ejA6fd5IlXQjgJnqO2JA2fouq0DZ9vrdypkbltOfqYNjKepQpuIiIiInJ3O0IlyTtsai10kGHQ0BaPDYasg2dTyWJvKZRyBszlInihwHqtvimy3uU84ZDaFA3B3iQRBb3SlMjRcNSkq3I0LfMC/7/8eAzMsP57rKrSJiIiIyNlJE+WcFh6PkRQeTnkmcM61HQ4bCEYFvOYKYsvAGbNyGQ6cjQEXDoYdB87GQJDxB9/H5xrj7nNcoc3MLgceBrzAL51z97c6ngz8GpgM7AP+zjm3Je5eiIiIiIh0B02UI62YGUm+BIfMbQZP/Q4HcZX9OuypmXmBR4DPAuOA68xsXKtmNwEHnHMjgB8DP+xcr0VERERERHqJcBV411G3PZ7m8cTLqcBm59xHzrkGYD5wdas2VwNPhbefB6abpnoRERERERGJbchUqo+4nfE0jSe05QPbot5XhffFbOOcawIOAbnxdEBERERERETaF09oi1Uxaz32Mp42mNnNZrbSzFbu2bMnnv6JiIiIiIj0avGEtipgSNT7AqD12MtIGzPzAVnA/tYXcs495pwrc86V9e/f/+R6LCIiIiIi0ovEE9pWACPNbKiZJQHXAuWt2pQDN4S3vwj80TnXfQsgiIiIiIiI9BIWT7YysyuAnxCa8v8J59x/mtndwErnXLmZpQC/AUoJVdiudc591ME19wBbT/UDSEQ/YG+iOyESg342pSfTz6f0VPrZlJ5KP5tdq8g51+EQxLhCm/R8ZrbSOVeW6H6ItKafTenJ9PMpPZV+NqWn0s9mYpwZy5aLiIiIiIj0UgptIiIiIiIiPZhC29njsUR3QKQd+tmUnkw/n9JT6WdTeir9bCaAnmkTERERERHpwVRpExERERER6cEU2s5gZjbEzP5kZhvMbJ2ZfTPRfRKJZmZeM6sws1cS3ReRaGaWbWbPm9kH4f8P/VSi+yQCYGbfCv+dvtbMng0vqySSEGb2hJntNrO1Ufv6mtliM9sUfs1JZB97C4W2M1sT8M/OubHA+cDXzWxcgvskEu2bwIZEd0IkhoeBRc65McC56OdUegAzywe+AZQ550oIrY97bWJ7Jb3ck8DlrfbdAbzpnBsJvBl+L91Moe0M5pzb4ZxbFd4+QugfHfmJ7ZVIiJkVADOBXya6LyLRzKwPcAHwKwDnXINz7mBieyUS4QNSzcwHpAHbE9wf6cWcc28D+1vtvhp4Krz9FPD509qpXkqh7SxhZsVAKbAssT0RifgJcDsQTHRHRFoZBuwB5oWH7/7SzNIT3SkR51w18CDwCbADOOSc+31ieyXSxgDn3A4IFRCAvAT3p1dQaDsLmFkG8AJwm3PucKL7I2JmVwK7nXPvJrovIjH4gEn/f3v3FmNXXcVx/PuzlIiXlBD6UMAEbJoaJLECIi0GCS2RB1NqRHnBdO/Uz6kAAAWbSURBVIpJX0SNERNRsCZ4N7EPYDBeQo1WqbSjrahcFAsR1EjbgbYgGluVWi+F4g2kii4f9n/wpCkDlg77lPl+Xs6e/z577bXnPJyss/77v4Frq+rVwKM4vUdDoN0bdAFwEnAc8OIkF/eblaRhYNF2mEsyna5gW11Vo33nIzVnAYuT/Bq4Hjg3yVf7TUl60i5gV1WNz0xYS1fESX1bBOysqj1V9S9gFFjQc07S/v6YZBZAe/1Tz/lMCRZth7Ekobsn4/6q+kzf+Ujjquryqjqhqk6ku4n+tqry12INhar6A/BgkrltaCFwX48pSeN+C5yZ5EXtO34hLpKj4bMBWNq2lwLre8xlyjii7wT0rJwFvA3YmmSsjX2gqr7bY06SdDh4J7A6yZHADmBZz/lIVNVPk6wFNtOtEL0F+Hy/WWkqS/J14Bzg2CS7gBXAJ4BvJHk73Q8Nb+kvw6kjVdV3DpIkSZKkp+D0SEmSJEkaYhZtkiRJkjTELNokSZIkaYhZtEmSJEnSELNokyRJkqQhZtEmSZoUST6YZHuSe5OMJXntQcQYSXLcZOTXhySLk7y/bS9JcnLfOUmShp/PaZMkHXJJ5gNvBE6tqn1JjgWOPIhQI8A2YPchTK83VbWB7sG0AEuAG/HB3pKkp2GnTZI0GWYBD1XVPoCqeqiqdidZmOSb429Kcl6S0STTkqxKsi3J1iTvSXIhcDrdQ7DHkhyV5LQktyfZlOTmJLNanI1JVia5I8n9SV7T4v4yyUcOlGCS85NsTnJPkh+0sTOS3JVkS3ud28ZHkqxPclOSB5KsGIjzrZbP9iTLnyb+SJJrkiwAFgOfbtc2O8nmgWPnJNl0yD4NSdJhzU6bJGky3AJ8KMkvgO8Da6rqduA24LNJZlbVHmAZcB0wDzi+qk4BSHJ0Vf05yaXAZVV1d5LpwNXABVW1J8lFwEeBS9o5/1lVZyd5N7AeOA3YC/wqycqqeng8uSQzgS8AZ1fVziTHtF0/b2NPJFkEfAx4c9t3BnAK8BjwsyTfqaq7gUuqam+So9r4OrofRQ8UH4CquivJBuDGqlrbcvpLknlVNdb+L6uezQcgSXr+sNMmSTrkqurvdEXTcmAPsCbJSFUV8BXg4iRHA/OB7wE7gJcnuTrJ+cBfDxB2Ll3RdGuSMeAK4ISB/ePTDrcC26vq963TtwN42X6xzgTuqKqdLd+9bXwGcEOSbcBK4JUDx9xaVQ9X1T+AUeB1bfxdSe4BftLOM2eC+BP5IrAsyTTgIuBrz+AYSdIUYKdNkjQpqurfwEZgY5KtwFK67tF1wLeBx4EbquoJ4JEkrwLeALwDeCv/66CNC10xNv8pTrmvvf5nYHv87/2/7wLUAWJcBfywqt6U5MSW/5OXtP8lJjkHWATMr6rHkmwEXjhB/ImsA1bQdSM3DXYGJUlTm502SdIhl2RukjkDQ/OA3wBU1W66hUWuoE0BbAuVvKCq1gFXAqe24/4GvLRtPwDMbIuckGR6ksFO2P/jx8Drk5zUYo1PX5wB/K5tj+x3zHlJjmnTIJcAd7b3P9IKtlfQddgmij9o8NqoqseBm4Fr6QpbSZIAizZJ0uR4CfDlJPcluRc4GfjwwP7VwINVNb5y4vF0HbkxukLu8ja+CvhcG58GXAh8sk1HHAMWHExy7X665cBoi7Wm7foU8PEkd7bzDfoR3dTOMWBdu5/tJuCIdo1X0U2RnCj+oOuB97VFT2a3sdV0HbpbDua6JEnPT+luL5Ak6bmT5BpgS1V9qe9cnokkI8DpVXXpJJ/nMmBGVV05meeRJB1evKdNkvScakvZPwq8t+9chkl7FMJs4Ny+c5EkDRc7bZIkSZI0xLynTZIkSZKGmEWbJEmSJA0xizZJkiRJGmIWbZIkSZI0xCzaJEmSJGmIWbRJkiRJ0hD7L2Ksz0xcPXUzAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7ff3b5b96780>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "Ks        = arange(1,12,2) # system capacities\n",
    "Ploss_est = zeros(len(Ks)) # estimated loss probabilities\n",
    "\n",
    "########################################################\n",
    "# complete the value returned by function estimate_Ploss\n",
    "# that estimates the loss probability from the obbserved \n",
    "# number of arrivals and number of lost arrivals\n",
    "def estimate_Ploss(arrivals, losses):\n",
    "    return float(losses)/float(arrivals)\n",
    "########################################################\n",
    "mu = 5\n",
    "for lambda_ in [4,10]:\n",
    "    # estimated loss probabilities:\n",
    "    for index,K in enumerate(Ks):\n",
    "        T,N,arrivals,losses = \\\n",
    "        MM1K(lambda_=lambda_,K=K,Tmax=10**3)\n",
    "        Ploss_est[index] = estimate_Ploss(arrivals, losses)\n",
    "    plot(Ks,Ploss_est,label=\"$\\lambda$=%d\"%lambda_)\n",
    "    # loss probabilities:\n",
    "    rho = lambda_/mu\n",
    "    Ploss = (1-rho)/(1-rho**(Ks+1))*rho**Ks \n",
    "    plot(Ks,Ploss,'.',label=\"Theory, $\\lambda$=%d\"%lambda_)\n",
    "axis(xmin=1,xmax=11)\n",
    "xlabel(\"System capacity\")\n",
    "title(\"Loss probability\")\n",
    "legend(loc=(.85,.2))\n",
    "\n",
    "#--------------------------\n",
    "V2 = estimate_Ploss(2,1) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Your answers for the exercise"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "---------------------------\n",
      "RESULTS SUPPLIED FOR LAB 4:\n",
      "---------------------------\n",
      "V1 = 0.17\n",
      "V2 = 0.50\n"
     ]
    }
   ],
   "source": [
    "print(\"---------------------------\\n\"\n",
    "      +\"RESULTS SUPPLIED FOR LAB 4:\\n\"\n",
    "      +\"---------------------------\")\n",
    "results = (\"V\"+str(k) for k in range(1,3))\n",
    "for x in results:\n",
    "    try:\n",
    "        print(x+\" = {0:.2f}\".format(eval(x)))\n",
    "    except:\n",
    "        print(x+\": variable is undefined\")   "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}