2016-08-17-notes.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# UBC Scientific Software Seminar\n",
    "\n",
    "##  August 17, 2016\n",
    "\n",
    "Today's Agenda:\n",
    "\n",
    "1. SciPy Example: Fitting Wind Speed Distributions to the Weibull Distribution\n",
    "2. A tour of scikit-learn (by examples)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import pandas as pd\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1. SciPy Example: Fitting Wind Speed Distributions to the Weibull Distribution\n",
    "\n",
    "It is well-known in the [wind energy industry](http://wind-data.ch/tools/weibull.php?lng=en) that wind speed distributions (measured in m/s) are modelled by the [Weibull distribution](https://en.wikipedia.org/wiki/Weibull_distribution):\n",
    "\n",
    "![weibull](https://wikimedia.org/api/rest_v1/media/math/render/svg/847f7f8870ffa9eabda2b09dc503ff4cb7753695)\n",
    "\n",
    "where, in the case of wind speed, the independent variable $x$ is wind speed measured in m/s, the coefficient $\\lambda$ is approximately equal to the average wind speed and $k$ is usually in the range [1,2].\n",
    "\n",
    "Let's plot the wind speed distribution at the Vancouver Airport using pandas and matplotlib, and fit the distribution to the Weibull distribution using SciPy."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Plot the wind speed distribution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "wind = pd.read_csv(\"2016-08-10-notes/wind.csv\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Year</th>\n",
       "      <th>Month</th>\n",
       "      <th>Day</th>\n",
       "      <th>Hour</th>\n",
       "      <th>Pressure</th>\n",
       "      <th>Wind Direction</th>\n",
       "      <th>Wind Speed</th>\n",
       "      <th>Temperature</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>2010</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>100.87</td>\n",
       "      <td>E</td>\n",
       "      <td>22.0</td>\n",
       "      <td>5.5</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>2010</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>100.81</td>\n",
       "      <td>E</td>\n",
       "      <td>19.0</td>\n",
       "      <td>5.4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>2010</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>2</td>\n",
       "      <td>100.77</td>\n",
       "      <td>E</td>\n",
       "      <td>26.0</td>\n",
       "      <td>5.7</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>2010</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>3</td>\n",
       "      <td>100.69</td>\n",
       "      <td>ESE</td>\n",
       "      <td>26.0</td>\n",
       "      <td>5.5</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>2010</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>4</td>\n",
       "      <td>100.58</td>\n",
       "      <td>ESE</td>\n",
       "      <td>26.0</td>\n",
       "      <td>5.8</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   Year  Month  Day  Hour  Pressure Wind Direction  Wind Speed  Temperature\n",
       "0  2010      1    1     0    100.87              E        22.0          5.5\n",
       "1  2010      1    1     1    100.81              E        19.0          5.4\n",
       "2  2010      1    1     2    100.77              E        26.0          5.7\n",
       "3  2010      1    1     3    100.69            ESE        26.0          5.5\n",
       "4  2010      1    1     4    100.58            ESE        26.0          5.8"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "wind.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Year</th>\n",
       "      <th>Month</th>\n",
       "      <th>Day</th>\n",
       "      <th>Hour</th>\n",
       "      <th>Pressure</th>\n",
       "      <th>Wind Direction</th>\n",
       "      <th>Wind Speed</th>\n",
       "      <th>Temperature</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>54774</th>\n",
       "      <td>2016</td>\n",
       "      <td>4</td>\n",
       "      <td>1</td>\n",
       "      <td>7</td>\n",
       "      <td>102.15</td>\n",
       "      <td>ENE</td>\n",
       "      <td>4.0</td>\n",
       "      <td>8.7</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>54775</th>\n",
       "      <td>2016</td>\n",
       "      <td>4</td>\n",
       "      <td>1</td>\n",
       "      <td>8</td>\n",
       "      <td>102.15</td>\n",
       "      <td>ESE</td>\n",
       "      <td>17.0</td>\n",
       "      <td>10.1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>54776</th>\n",
       "      <td>2016</td>\n",
       "      <td>4</td>\n",
       "      <td>1</td>\n",
       "      <td>9</td>\n",
       "      <td>102.15</td>\n",
       "      <td>ESE</td>\n",
       "      <td>11.0</td>\n",
       "      <td>11.7</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>54777</th>\n",
       "      <td>2016</td>\n",
       "      <td>4</td>\n",
       "      <td>1</td>\n",
       "      <td>10</td>\n",
       "      <td>102.15</td>\n",
       "      <td>SE</td>\n",
       "      <td>10.0</td>\n",
       "      <td>13.5</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>54778</th>\n",
       "      <td>2016</td>\n",
       "      <td>4</td>\n",
       "      <td>1</td>\n",
       "      <td>11</td>\n",
       "      <td>102.15</td>\n",
       "      <td>SSW</td>\n",
       "      <td>13.0</td>\n",
       "      <td>13.4</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "       Year  Month  Day  Hour  Pressure Wind Direction  Wind Speed  \\\n",
       "54774  2016      4    1     7    102.15            ENE         4.0   \n",
       "54775  2016      4    1     8    102.15            ESE        17.0   \n",
       "54776  2016      4    1     9    102.15            ESE        11.0   \n",
       "54777  2016      4    1    10    102.15             SE        10.0   \n",
       "54778  2016      4    1    11    102.15            SSW        13.0   \n",
       "\n",
       "       Temperature  \n",
       "54774          8.7  \n",
       "54775         10.1  \n",
       "54776         11.7  \n",
       "54777         13.5  \n",
       "54778         13.4  "
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "wind.tail()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<class 'pandas.core.frame.DataFrame'>\n",
      "RangeIndex: 54779 entries, 0 to 54778\n",
      "Data columns (total 8 columns):\n",
      "Year              54779 non-null int64\n",
      "Month             54779 non-null int64\n",
      "Day               54779 non-null int64\n",
      "Hour              54779 non-null int64\n",
      "Pressure          54776 non-null float64\n",
      "Wind Direction    53180 non-null object\n",
      "Wind Speed        54777 non-null float64\n",
      "Temperature       54777 non-null float64\n",
      "dtypes: float64(3), int64(4), object(1)\n",
      "memory usage: 3.3+ MB\n"
     ]
    }
   ],
   "source": [
    "wind.info()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since we're only interested in wind speed, we can select the \"Wind Speed\" column:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0    22.0\n",
       "1    19.0\n",
       "2    26.0\n",
       "3    26.0\n",
       "4    26.0\n",
       "Name: Wind Speed, dtype: float64"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "wind_speed = wind[\"Wind Speed\"]\n",
    "wind_speed.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice that wind speed is measured in km/h in this dataset. Let's check to see if there are any missing values:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "6166   NaN\n",
       "6361   NaN\n",
       "Name: Wind Speed, dtype: float64"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "wind_speed[wind_speed.isnull()]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are only 2 missing values and we can fill them in with the values from the previous measurement using the [fillna](http://pandas.pydata.org/pandas-docs/stable/generated/pandas.DataFrame.fillna.html) DataFrame method:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "wind_speed = wind_speed.fillna(method='ffill')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To plot the wind speed distribution, we can plot a [histogram](http://matplotlib.org/api/pyplot_api.html#matplotlib.pyplot.hist) of wind speeds with bins of size 4 km/h up to the maximum wind speed:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "78.0"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "wind_speed.max()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAEACAYAAABfxaZOAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAEwBJREFUeJzt3XGsnfV93/H3x3UggRY3W1d7hWBaIKWzuhimOu5Qldul\nHbaj1n8OqxutpW2WBku0ShUs08TNn/lja4OoStEoKlUyULKmcSq3MxE5UqtpjjMwMLCDWRrHOHCX\nqCURThUR/N0f54EdTux7z7333Hse+/d+SUc+z3N+z3M+59x7P/fx75zn3FQVkqQ2bJh1AEnS+rH0\nJakhlr4kNcTSl6SGWPqS1BBLX5IaMlHpJ9mV5ESSF5LcfYEx9yU5meRYku0j6zcl+XSS40meS/L+\naYWXJC3PkqWfZANwP3AbsA3Yl+SmsTG7geur6kbgAPDAyM2fAA5V1c8A7wOOTym7JGmZJjnS3wGc\nrKpTVfU68Ciwd2zMXuARgKo6AmxKsjnJVcAvVNXD3W3fr6rvTC++JGk5Jin9q4HTI8svdesWG3Om\nW/eTwLeSPJzkySQPJnnXagJLklZurV/I3QjcAvxuVd0CfBe4Z43vU5J0ARsnGHMGuHZk+Zpu3fiY\n91xgzOmq+nJ3/TPAhV4I9kOAJGmZqirLGT/Jkf5R4IYkW5NcBtwOHBwbcxC4AyDJTuDVqlqoqgXg\ndJL3duM+CDy/SPheX+69996ZZzCnOc1pzjcvK7HkkX5VvZHkLuAww18SD1XV8SQHhjfXg1V1KMme\nJC8CZ4H9I7v4MPDJJO8Avjp2myRpHU0yvUNV/Tnw02Prfn9s+a4LbPs08HMrDShJmh7PyF2Gubm5\nWUeYiDmny5zTZc7ZykrnhaYtSfUliyRdDJJQa/BCriTpEmHpS1JDLH1JaoilL0kNsfQlqSGWviQ1\nxNKXpIZY+pLUEEtfkhpi6UtSQyx9SWqIpS9JDbH0Jakhlr4kNcTSl6SGWPqS1BBLX5IaYulLUkMs\nfUlqiKUvSQ2x9CWpIZa+JDXE0pekhlj6I7ZsuY4kq7ps2XLdrB+GJF1QqmrWGQBIUrPOkgRYbYYw\n68chqQ1JqKosZxuP9CWpIROVfpJdSU4keSHJ3RcYc1+Sk0mOJbl5ZP3Xkjyd5KkkX5pWcEnS8m1c\nakCSDcD9wAeBbwBHk3yuqk6MjNkNXF9VNyZ5P/B7wM7u5nPAXFX9zdTTS5KWZZIj/R3Ayao6VVWv\nA48Ce8fG7AUeAaiqI8CmJJu72zLh/UiS1tgkZXw1cHpk+aVu3WJjzoyMKeDxJEeT/KuVBpUkrd6S\n0ztTcGtVvZzk7zEs/+NV9ZfrcL+SpDGTlP4Z4NqR5Wu6deNj3nO+MVX1cvfvN5N8luF00XlLf35+\n/q3rc3NzzM3NTRBPktowGAwYDAar2seS79NP8kPAVxi+kPsy8CVgX1UdHxmzB7izqj6UZCfwO1W1\nM8kVwIaqei3JlcBh4GNVdfg89+P79CVpGVbyPv0lj/Sr6o0kdzEs7A3AQ1V1PMmB4c31YFUdSrIn\nyYvAWWB/t/lm4LNJqruvT56v8CVJ68Mzct+egdUf6b8T+N6q9rB581ZeeeVrq8wh6VK3kiN9S//t\nGZjG9I5TRJLWQ9MfwzCND0uTpEvdJXOk36ejdI/0Ja2Hpo/0JUlLs/QlqSGWviQ1xNKXpIZY+pLU\nEEtfkhpi6UtSQyx9SWqIpS9JDbH0Jakhlr4kNcTSl6SGWPqS1BBLX5IaYulLUkMsfUlqiKUvSQ2x\n9CWpIZa+JDXE0pekhlj6ktQQS1+SGmLpS1JDLH1JaoilL0kNsfQlqSETlX6SXUlOJHkhyd0XGHNf\nkpNJjiXZPnbbhiRPJjk4jdCSpJVZsvSTbADuB24DtgH7ktw0NmY3cH1V3QgcAB4Y281HgOenkliS\ntGKTHOnvAE5W1amqeh14FNg7NmYv8AhAVR0BNiXZDJDkGmAP8F+mllqStCKTlP7VwOmR5Ze6dYuN\nOTMy5reB3wJqhRklSVOycS13nuRDwEJVHUsyB2Sx8fPz829dn5ubY25ubi3jSdJFZTAYMBgMVrWP\nVC1+AJ5kJzBfVbu65XuAqqqPj4x5APhiVT3WLZ8APsBwLv+fA98H3gX8CPDHVXXHee6nlsqyRE5W\n/5+J/uxjNc+FpDYkoaoWPZgeN8n0zlHghiRbk1wG3A6MvwvnIHBHF2In8GpVLVTVR6vq2qr6qW67\nJ85X+JKk9bHk9E5VvZHkLuAww18SD1XV8SQHhjfXg1V1KMmeJC8CZ4H9axtbkrQSS07vrBend96+\nj758XST111pN70iSLhGWviQ1xNKXpIZY+pLUEEtfkhpi6UtSQyx9SWqIpS9JDbH0Jakhlr4kNcTS\nl6SGWPqS1BBLX5IaYulLUkMs/V66nCQrvmzZct2sH4CknvLz9N++l0tkH34ev9QCP09fkrQoS1+S\nGmLpS1JDLH1JaoilL0kNsfQlqSGWviQ1xNKXpIZY+pLUEEtfkhpi6UtSQyx9SWrIRKWfZFeSE0le\nSHL3Bcbcl+RkkmNJtnfrLk9yJMlTSZ5Ncu80w0uSlmfJ0k+yAbgfuA3YBuxLctPYmN3A9VV1I3AA\neACgqr4H/GJV3QxsB3Yn2THdhyBJmtQkR/o7gJNVdaqqXgceBfaOjdkLPAJQVUeATUk2d8vf7cZc\nDmxk9Z87LElaoUlK/2rg9MjyS926xcaceXNMkg1JngJeAR6vqqMrjytJWo2Na30HVXUOuDnJVcCf\nJPkHVfX8+cbOz8+/dX1ubo65ubm1jidJF43BYMBgMFjVPpb8y1lJdgLzVbWrW74HqKr6+MiYB4Av\nVtVj3fIJ4ANVtTC2r/8InK2q/3ye+/EvZ01tH/7lLKkFa/WXs44CNyTZmuQy4Hbg4NiYg8AdXYid\nwKtVtZDkx5Js6ta/C/hl4MRyAkqSpmfJ6Z2qeiPJXcBhhr8kHqqq40kODG+uB6vqUJI9SV4EzgL7\nu83/PvCH3TuANgCPVdWhtXkokqSl+IfR376XS2QfTu9ILfAPo0uSFmXpS1JDLH1JaoilL0kNsfQl\nqSGWviQ1xNKXpIZY+pLUEEtfkhpi6UtSQyx9SWqIpS9JDbH0Jakhlr4kNcTSl6SGWPqS1BBLX5Ia\nYulLUkMsfUlqiKUvSQ2x9CWpIZa+JDXE0pekhlj6ktQQS1+SGmLpS1JDLH1JaoilL0kNmaj0k+xK\nciLJC0nuvsCY+5KcTHIsyfZu3TVJnkjyXJJnk3x4muElScuzZOkn2QDcD9wGbAP2JblpbMxu4Pqq\nuhE4ADzQ3fR94Derahvw88Cd49tKktbPJEf6O4CTVXWqql4HHgX2jo3ZCzwCUFVHgE1JNlfVK1V1\nrFv/GnAcuHpq6SVJyzJJ6V8NnB5ZfokfLO7xMWfGxyS5DtgOHFluSEnSdGxcjztJ8sPAZ4CPdEf8\n5zU/P//W9bm5Oebm5tY826XpcpKsag+bN2/llVe+Np04kqZiMBgwGAxWtY9U1eIDkp3AfFXt6pbv\nAaqqPj4y5gHgi1X1WLd8AvhAVS0k2Qj8KfBnVfWJRe6nlsqyRE5g5dt3e7lE9jGdDKv5ekhae0mo\nqmUd4U0yvXMUuCHJ1iSXAbcDB8fGHATu6ELsBF6tqoXutj8Anl+s8CVJ62PJ6Z2qeiPJXcBhhr8k\nHqqq40kODG+uB6vqUJI9SV4EzgK/AZDkVuDXgGeTPMXw8POjVfXna/R4JEmLWHJ6Z70kqYcffnjF\n2+/fv5/ZT6v0ZR9O70gtWMn0Tq9K/8orf31F254799f87d9+ntmXbV/2YelLLbjoS3/lRfU0w3eD\nzrps+7IPS19qwVq9kCtJukRY+pLUEEtfkhpi6UtSQyx9SWqIpS9JDbH0Jakhlr4kNcTSl6SGWPqS\n1BBLX5IaYulLUkMsfUlqiKUvSQ2x9CWpIZa+JDXE0pekhlj6ktQQS1+SGmLpS1JDLH1JaoilL0kN\nsfQlqSGWviQ1xNKXpIZY+pLUkIlKP8muJCeSvJDk7guMuS/JySTHktw8sv6hJAtJnplWaK2Hy0my\nqsuWLdfN+kFIGrNk6SfZANwP3AZsA/YluWlszG7g+qq6ETgA/N7IzQ932+qi8j2gVnVZWDi1/rEl\nLWqSI/0dwMmqOlVVrwOPAnvHxuwFHgGoqiPApiSbu+W/BP5mepElSSs1SelfDZweWX6pW7fYmDPn\nGSNJmjFfyJWkhmycYMwZ4NqR5Wu6deNj3rPEmAnMj1yf6y6SJIDBYMBgMFjVPlJViw9Ifgj4CvBB\n4GXgS8C+qjo+MmYPcGdVfSjJTuB3qmrnyO3XAZ+vqp9d5H5q+ALgSjwNbGfl27+V4hLZRx8yDPex\n1PeXpJVLQlVlOdssOb1TVW8AdwGHgeeAR6vqeJIDSf51N+YQ8FdJXgR+H/g3I6E+BfwP4L1Jvp5k\n/3ICSpKmZ8kj/fXikf4099GHDMN99OX7S7oUrcmRviTp0mHpS1JDLH1JaoilL0kNsfQlqSGWviQ1\nxNKXpIZY+lpDq/tMfj+PX5q+ST57R1qhNz+Tf2UWFpZ1zomkCXikL0kNsfQlqSGWviQ1xNKXpIZY\n+pLUEEtfkhpi6UtSQyx99djqTu7yBC/pB3lylnpsdSd3gSd4SeM80pekhlj6ktQQS1+SGmLpS1JD\nLH1JaoilL0kNsfQlqSGWviQ1xNLXJc6zeqVRnpGrS5xn9UqjJjrST7IryYkkLyS5+wJj7ktyMsmx\nJNuXs60kaX0sWfpJNgD3A7cB24B9SW4aG7MbuL6qbgQOAA9Muu3FZTDrAJeYwawDTOgdF8UU0WAw\nWPP7mAZzztYkR/o7gJNVdaqqXgceBfaOjdkLPAJQVUeATUk2T7jtRWQw6wCXmMGsA0zodYZTRCu/\nLCycWvOUF0tJmXO2Jin9q4HTI8svdesmGTPJtlIDVveCsi8ma1rW6oXcFb3yddVVv7KiOzt37tu8\n9tqKNpXWyepeUF5YeCfJ0j9WH/vYxy5424YNV3Du3HdXnGFa+7jyyk3Mz8+vah9auUlK/wxw7cjy\nNd268THvOc+YyybY9i3f+c6fThBnMdN4l8VS+7jwD9X65rgYMkyyj6Wez4vlcazXPlZutWU9rX2c\nPfvtiX6B9cFiv0QvVpOU/lHghiRbgZeB24F9Y2MOAncCjyXZCbxaVQtJvjXBtgBU1cXxXSBJF7El\nS7+q3khyF3CY4WsAD1XV8SQHhjfXg1V1KMmeJC8CZ4H9i227Zo9GkrSoVK3uxBVJ0sVj5h/D0NeT\nt5I8lGQhyTMj696d5HCSryT570k2zTJjl+maJE8keS7Js0k+3LesSS5PciTJU13Ge/uWcVSSDUme\nTHKwW+5dziRfS/J095x+qcc5NyX5dJLj3ffo+/uWM8l7u+fxye7fbyf5cN9ydln/XZL/neSZJJ9M\nctlyc8609Ht+8tbDDHONugf4QlX9NPAE8O/XPdUP+j7wm1W1Dfh54M7uOexN1qr6HvCLVXUzsB3Y\nnWRHnzKO+Qjw/MhyH3OeA+aq6uaq2tGt62POTwCHqupngPcBJ+hZzqp6oXsebwH+EcMp6s/Ss5xJ\nfgL4t8AtVfUPGU7P72O5OatqZhdgJ/BnI8v3AHfPMtNYvq3AMyPLJ4DN3fUtwIlZZzxP5j8Bfqmv\nWYErgC8DP9fHjAzfYfY4MAcc7OvXHfgr4O+OretVTuAq4P+cZ32vco5l+6fAX/QxJ/ATwCng3V3h\nH1zJz/qsp3cutpO3fryqFgCq6hXgx2ec522SXMfwSPp/Mvwm6E3WbsrkKeAV4PGqOtq3jJ3fBn6L\nt7+pvo85C3g8ydEk/7Jb17ecPwl8K8nD3dTJg0muoH85R/0z4FPd9V7lrKpvAP8J+DrDt75/u6q+\nwDJzzrr0L3a9eRU8yQ8DnwE+UlWv8YPZZpq1qs7VcHrnGmBHkm3nyTTTjEk+BCxU1TEWf1N9H77u\nt9ZwOmIPwym9X6BnzyfDo9FbgN/tsp5l+L/5vuUEIMk7gF8FPt2t6lXOJD/K8GNstjI86r8yya+d\nJ9eiOWdd+pOc+NUnC91nCpFkC/B/Z5wHgCQbGRb+H1XV57rVvcxaVd9h+KE7u+hfxluBX03yVeC/\nAv8kyR8Br/QsJ1X1cvfvNxlO6e2gf8/nS8Dpqvpyt/zfGP4S6FvON+0G/ldVfatb7lvOXwK+WlV/\nXVVvMHzd4R+zzJyzLv23TvxKchnDk7cOzjjTqPD2I76DwG90138d+Nz4BjPyB8DzVfWJkXW9yZrk\nx958R0GSdwG/DBynRxkBquqjVXVtVf0Uw+/FJ6rqXwCfp0c5k1zR/c+OJFcynId+lv49nwvA6STv\n7VZ9EHiOnuUcsY/hL/s39S3n14GdSd6ZJAyfz+dZbs4evHCyC/gKcBK4Z9Z5RnJ9CvgGww9N+TrD\nE87eDXyhy3sY+NEe5LwVeAM4BjwFPNk9p3+nL1mBn+1yHQOeAf5Dt743Gc+T+QP8/xdye5WT4Vz5\nm1/vZ9/8uelbzi7T+xge3B0D/hjY1NOcVwDfBH5kZF0fc97L8IDpGeAPgXcsN6cnZ0lSQ2Y9vSNJ\nWkeWviQ1xNKXpIZY+pLUEEtfkhpi6UtSQyx9SWqIpS9JDfl/LDmcWhcl4oUAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f3d18f0be48>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "data, bins, _ = plt.hist(wind_speed,bins=np.arange(0,82,4),normed=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The output of the matplotlib.pyplot function [hist](http://matplotlib.org/api/pyplot_api.html#matplotlib.pyplot.hist) is a tuple with three entries: the frequency for each bin, the edges of each bin, and (silent) list of patches (which we don't need). We've unpacked these into variables `data` and `bins` (and forgotten the patches `_`). Notice that we used the keyword argument `normed=True` to get the frequency as a fraction instead of a hard count of the wind speed data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([  1.87298052e-02,   4.87093594e-02,   5.19724712e-02,\n",
       "         4.90790266e-02,   3.48217383e-02,   2.00852516e-02,\n",
       "         1.11447818e-02,   7.01911316e-03,   4.44969788e-03,\n",
       "         1.75249639e-03,   9.03630953e-04,   5.79601672e-04,\n",
       "         3.14901696e-04,   2.41881013e-04,   1.09531025e-04,\n",
       "         4.56379269e-05,   2.73827562e-05,   9.12758539e-06,\n",
       "         0.00000000e+00,   4.56379269e-06])"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0,  4,  8, 12, 16, 20, 24, 28, 32, 36, 40, 44, 48, 52, 56, 60, 64,\n",
       "       68, 72, 76, 80])"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "bins"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Fit the wind speed distribution to the Weibull distribution\n",
    "\n",
    "The data that we want to fit is given by the NumPy array `data` however we want to prepend a data point to the array to represent 0 frequency at wind speed 0."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[  0.00000000e+00   1.87298052e-02   4.87093594e-02   5.19724712e-02\n",
      "   4.90790266e-02   3.48217383e-02   2.00852516e-02   1.11447818e-02\n",
      "   7.01911316e-03   4.44969788e-03   1.75249639e-03   9.03630953e-04\n",
      "   5.79601672e-04   3.14901696e-04   2.41881013e-04   1.09531025e-04\n",
      "   4.56379269e-05   2.73827562e-05   9.12758539e-06   0.00000000e+00\n",
      "   4.56379269e-06]\n"
     ]
    }
   ],
   "source": [
    "Y = np.append(0,data)\n",
    "print(Y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We want the corresponding wind speed values to be the midpoint of each bin and 0 as the first value:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 0  2  6 10 14 18 22 26 30 34 38 42 46 50 54 58 62 66 70 74 78]\n"
     ]
    }
   ],
   "source": [
    "X = np.append(0,bins[1:] - 2)\n",
    "print(X)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.text.Text at 0x7f3d184bdba8>"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAEKCAYAAAD+XoUoAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAHd5JREFUeJzt3X+cHXV97/HXe41QCRJUJHYXWOOioGgNWEMqV90UU5PQ\na/po6xVMikGvzeNRA9za64Xarrt57K1XbrECohdpkYBIof6ONbWJD1kVNTEKEYSEhGVdw26IgoAQ\nrUX2c/+Y72Ynh7O75+ye3XNO5v18PPaRMzPfmfnsbPZ95nznO7OKCMzMrBha6l2AmZnNHoe+mVmB\nOPTNzArEoW9mViAOfTOzAnHom5kViEP/MCbpCUkvnuK6t0l6Z20rqi1J75D0rWlu40eS3lCjet4u\n6au56RFJL6nFttP2pvzzrGDbH5R0UXr9Rkl7Z2g/3ZI+NcHybZJePhP7toxDv0lIulTSppJ5eyR9\npWTebkn/DSAinhsRP56BWuZJuk7SPkmPS9ol6X/Vej8VKnujiaT2FLq/SF/7JG2U9KZDVo54ZUR8\nc6Id5LY14e9LRNwcEcsmq60S5d50Z/DneRzwZ8An8rur8T6+mjv2E23774HeWu7bDuXQbx7fBH5P\nkgAkvQiYA5xeMq8jtZ1JHwHmAqdExDzgLcD9M7zPqQhgXkQcA7wa+BrwBUnnV7kdpW1p3AbSs8ZZ\nrxmsATZFxK9nYuOSjgJeA3yjguZfBpZIOn4majGHfjPZDhwBLEzTrwduA+4rmdcfEQ/Bod0Lkq6X\ndLWkf01nvt+VtGB045KWStop6VFJH2XiwHotcHNE/AIgInZHxOdz2xqRdKGkfkk/lfR/8ytLeqek\neyU9IunfJJ2UW3aqpM1p2U5Jb80te346W39c0layN7jJKNX404i4CugBDtYjaUDS76fXr5W0PW1/\nn6TLU7PRsHosHbszU9fS7ZL+QdLDQPc43U3nlDsOpd0c+U8Tkv432c/y6rS/q3LHdfTneYykG9N2\nByT9TW5b75D0LUl/L+nnaf/5TyClljNBIEu6KHWDtY52/Uh6n6T9koYkrZS0XNJ9kh6W9Nclmzgb\n+HZEPJWmj5R0Q/re7pZ0xmjD9MbzA+DNE9Rr0+DQbxLpF2YbMNr//AayM/rby8w7uFrJZt4GdAPH\nAv3A3wFIegHwOeD9wHFp2VkTlLMV+KCkNZJOHqfNHwFnpK+Vo10VklYCl6blLwS+BfxzWnYUsBm4\nKdVxLvBxSaembX4c+CUwH3gXMJVrDp8Hjpd0SpllVwJXpE8vHcC/pPmjx/eYiDgmIral6TPJPuEc\nTzqWPPOYlz0O47QNgIj4W7Ljsi7t76Iy7a8Gngu8GOgEzpd0QW75ImAn8AKyLpPryny/o15FdvLw\nDJI+AJwPvCEihtPsF5GdgLSS/X/6R2AVcDrZseqS1J7bzAog3w35X4GbgXlkZ/YfK9ntTrJPZjYD\nHPrN5RuMBdDryYLh9pJ5+TO20rP1L0TEDyJiBPg0Y58QVgA/iogvRMTTEXEF8NAEdawjC+b3APco\nu45Qeib5oYh4PCIeBK4Azkvz1wL/J306GAE+BCyUdCLwh8BARNwYmR+SvRm9NfWn/zHQFRH/ERH3\nADdMUON4RoPr+WWW/SdwsqQXRMQvI+J7JctLj+dQRHw8IkYm6BoZ7zhMxWg3XgvZG/ilqc5B4MNk\n/fKjBiPik5E9XOsG4EUTdJkcCzxRMq9F0oeBNwGdEfHz3LL/BD4YEU8Dt5C9QV+RarkXuJdDQ3sF\nkL8edXtE/Huq7VPA75Ts+4lUk80Ah35z+SbwXyQ9DzguIvqB7wCvS/NeycT9+fkg/yVwdHrdCpSO\n1hh39EZE/DoiPhQRryU7k/wM8BlJ+V/UB3OvB9M+ANqBK1O3w8+BR8jOYNvSssWjyyQ9Cryd7Mz+\nhWTXMEq3W6229O8jZZa9CzgF2KVsFMk5k2yrkhEu4x2H6TiO7Fj8pGTbbbnpgz/riPgV2RvG0ZT3\nKNmnhrxjgXeTvUE/WbLskRh7UuOv0r8/zS3/1ei+JL0SeCwihsrVRvb/8Ld06EXy5wKPjVOrTZND\nv7l8l7Ffxm8DRMQTZGev7yY785xKEO4DTiqZd2IlK6ZA+CDZhd0FuUX59dsZO8PeC6yNiOenr+dF\nxNERsTUt6ytZdkxErAN+BjxVst3Smivxx8D+iNhd5nvpj4i3R8QLyfr9PyvpOYw/2qSSES6l9Y4e\nhwPAUbllv13Fth8mOxb5LpR2YKh880ndBbysZN7PyT55bZD0uiluF555ll+JlwM/nMY+bQIO/SYS\nEf8BfB94L1nXzqhvp3lTHbXzFeAVkv5I0rMkXUx2dl2WpL+V9LuSni3pSOB/kJ0t5vuF3yfp2NRt\ncxFZNwDANcD7Jb0ibWuepD9Ny/4VeJmk1ZLmpO3/rqRTUlfQ54EeSc9J679jku9LjHWJHC9pHdBF\ndk2h3Pe1StnwRYDHyYJ3hOwNZ4TKLhyXyh+Hixk7DjuAN0g6UdK8MjXtB8qO8U/H4l+Av5N0dOo/\n/0uyrpKp2ER2XaB0P98k66v/nKTXTnHbpf355RzsNkv/n14DbJni/mwSDv3m8w2yro7bc/O+leaV\njsCoaKx1RDwCvBW4jOwssoP0SWK8VYDrycJwiGx0xjkR8ctcmy+RjcK4g+xi3SfTvr5I1o9/i6TH\nyM4yl6VlTwJ/QHYBdzh9fQg4Mm3zQrKP/vvS9j452bcGPCrpidx+/jQibihpM2oZ2TWKX5ANS31b\n6sr6FdmF2m+nbqdFk+w3v+3xjsPXgFtTXdvTsrwrya5lPCLpijK1XkTWNfIA2Zv9TRFx/SS1jOdG\nYHkK3ENXyup8F7BR0sJnrFl+2wHZGzrZWft3Jth36fpvAW4bHYFmtadK/ohKukh3BdmbxHURcVmZ\nNleRDf06AKyJiB1p/jzgn8j6m0eAd+ZGP9hhSNIIcHJEPFDvWqwyaZjo6LDWWm3zrcCfRMS5Vazz\nXeBd6YKwzYBJQz9dYNlNdjY3THZWcm5E7Mq1WU42vOwcSWcCV0bE4rRsA/CNiLhe0hzgqNHx3XZ4\ncugbgLI7cJ/wSV5jmVNBm0XAntELhJJuAVYCu3JtVpJ9RCQitqV+2vlkV/FfHxFr0rLfAA78w5//\nBqeNdg1Zg6mkT7+NQ4emPcihQ8PKtRlK8xYADyu7G/QOSdem0RB2GIuIZ/ks36wxzfSF3DlkdyJ+\nLCLOILvwVHbkhJmZzbxKuneGOHQ89Ak8czzwEIeOR8632RsR30+vPwtcUm4nktwlYGZWpYio6sF+\nlZzpbye7Nb1d0hFkw+k2lrTZSPZ8DiQtJrsDb39E7Af2Shq98eNsslu0xyu+ob+6u7vrXoPrdJ2u\n03WOfk3FpGf6EfF0uqllM2NDNndKWpstjmsjYpOkFZLuJxuymX/w00XApyU9m2xM8QWl+zAzs9lR\nSfcOEfFVsmeS5Od9omR63Tjr/pDsUbxmZlZnviO3Cp2dnfUuoSKus7ZcZ225zvqq6I7c2SApGqUW\nM7NmIImYgQu5ZmZ2mHDom5kViEPfzKxAHPpmZgXi0DczKxCHvplZgTj0zcwKxKFvZlYgDn0zswJx\n6JuZFYhD38ysQBz6ZmYF4tA3MysQh76ZWYE49M3MCsShb2ZWIA59M7MCceibmRWIQ9/MrEAc+mZm\nBTKn3gU0g4GBQbq6NjA0NEJbWwu9vWtYsKC93mWZmVVNEVHvGgCQFI1SS97AwCBLl36U/v71wFzg\nAB0d3WzZcqGD38zqShIRoWrWcffOJLq6NuQCH2Au/f3r6eraUMeqzMymxqE/iaGhEcYCf9RchodH\n6lGOmdm0OPQn0dbWAhwomXuA1lYfOjNrPk6uSfT2rqGjo5ux4M/69Ht719StJjOzqaroQq6kZcAV\nZG8S10XEZWXaXAUsJ0vHCyLizjT/x8DjwAjwVEQsGmcfDXkhF8ZG7wwPj9Da6tE7ZtYYpnIhd9LQ\nl9QC7AbOBoaB7cC5EbEr12Y5sC4izpF0JnBlRCxOyx4AXhMRj06yn4YNfTOzRjRTo3cWAXsiYjAi\nngJuAVaWtFkJ3AgQEduAeZLmj9ZV4X7MzGyGVRLGbcDe3PSDad5EbYZybQLYImm7pHdPtVAzM5u+\n2bgj96yI2CfphWThvzMibp+F/ZqZWYlKQn8IOCk3fUKaV9rmxHJtImJf+vdnkr5A1l1UNvR7enoO\nvu7s7KSzs7OC8pqDH+VgZtPV19dHX1/ftLZRyYXcZwH3kV3I3Qd8DzgvInbm2qwA3pMu5C4GroiI\nxZKOAloi4klJc4HNwPqI2FxmP4fthVw/ysHMZsKMXMiNiKeBdWSBfQ9wS0TslLRW0p+nNpuAAUn3\nA58A/iKtPh+4XdKdwFbgy+UCf6YNDAyyevV6lizpZvXq9QwMDM7q/v0oBzNrFBX16UfEV4FTSuZ9\nomR6XZn1BoCF0ylwusqdZW/dOrtn2X6Ug5k1isN+KGUjnGX7UQ5m1igO+9RphLNsP8rBzBrFYf9H\nVMbOsvPBP7tn2QsWtLNly4V0dV2ee5SDL+Ka2ew77P+IikfOmNnhakaevTNbZnLIph+YZmaHI4e+\nmVmB+M8lmpnZhBz6ZmYF4tA3MysQh76ZWYE49M3MCsShb2ZWIA59M7MCceibmRWIQ9/MrEAc+mZm\nBeLQNzMrEIe+mVmBOPTNzArEoW9mViAOfTOzAnHom5kViEPfzKxAHPpmZgXi0DczKxCHvplZgTj0\nzcwKxKFvZlYgDn0zswKpKPQlLZO0S9JuSZeM0+YqSXsk7ZC0sGRZi6Q7JG2sRdFmZjY1k4a+pBbg\nauDNwGnAeZJOLWmzHOiIiJcCa4FrSjZzMXBvTSo2M7Mpq+RMfxGwJyIGI+Ip4BZgZUmblcCNABGx\nDZgnaT6ApBOAFcA/1axqMzObkkpCvw3Ym5t+MM2bqM1Qrs1HgPcBMcUazcysRubM5MYlnQPsj4gd\nkjoBTdS+p6fn4OvOzk46Oztnsjwzs6bS19dHX1/ftLahiIlPwCUtBnoiYlmavhSIiLgs1+Ya4LaI\nuDVN7wLeSNaXvxr4DfAc4LnA5yPi/DL7iclqMTOzMZKIiAlPpktV0r2zHThZUrukI4BzgdJROBuB\n81MRi4HHImJ/RLw/Ik6KiJek9b5eLvDNzGx2TNq9ExFPS1oHbCZ7k7guInZKWpstjmsjYpOkFZLu\nBw4AF8xs2WZmNhWTdu/MFnfvmJlVZ6a6d8zM7DDh0DczKxCHvplZgTj0zcwKxKFvZlYgDn0zswJx\n6JuZFYhD38ysQBz6ZmYF4tA3MysQh76ZWYHM6PP0rbYGBgbp6trA0NAIbW0t9PauYcGC9nqXZWZN\nxA9caxIDA4MsXfpR+vvXA3OBA3R0dLNly4UOfrOC8gPXDmNdXRtygQ8wl/7+9XR1bahjVWbWbBz6\nTWJoaISxwB81l+HhkXqUY2ZNyqHfJNraWsj+Pk3eAVpb/SM0s8o5MZpEb+8aOjq6GQv+rE+/t3dN\n3Woys+bjC7lNZHT0zvDwCK2tHr1jVnRTuZDr0Dcza1IevWNmZhNy6JuZFYhD38ysQBz6ZmYF4tA3\nMysQh76ZWYE49M3MCsShb2ZWIA59M7MCqSj0JS2TtEvSbkmXjNPmKkl7JO2QtDDNO1LSNkl3Srpb\nUnctizczs+pMGvqSWoCrgTcDpwHnSTq1pM1yoCMiXgqsBa4BiIhfA0si4nRgIbBc0qLafgtmZlap\nSs70FwF7ImIwIp4CbgFWlrRZCdwIEBHbgHmS5qfpX6Y2R5L9eUY/YMfMrE4qCf02YG9u+sE0b6I2\nQ6NtJLVIuhN4CNgSEdunXq6ZmU3HjP9h9IgYAU6XdAzwRUmviIh7y7Xt6ek5+Lqzs5POzs6ZLs/M\nrGn09fXR19c3rW1M+mhlSYuBnohYlqYvBSIiLsu1uQa4LSJuTdO7gDdGxP6SbXUBByLiH8rsx49W\nNjOrwkw9Wnk7cLKkdklHAOcCG0vabATOT0UsBh6LiP2SjpM0L81/DrAU2FVNgWZmVjuTdu9ExNOS\n1gGbyd4krouInZLWZovj2ojYJGmFpPvJ/p7fBWn13wZuSCOAWoBbI2LTzHwrZmY2Gf/lLDOzJuW/\nnGVmZhNy6JuZFYhD38ysQBz6ZmYF4tA3MysQh76ZWYE49M3MCsShb2ZWIA59M7MCceibmRWIQ9/M\nrEAc+mZmBeLQNzMrEIe+mVmBOPTNzArEoW9mViAOfTOzAnHom5kViEPfzKxAHPpmZgXi0DczKxCH\nvplZgTj0zcwKxKFvZlYgDn0zswJx6JuZFYhD38ysQBz6ZmYFUlHoS1omaZek3ZIuGafNVZL2SNoh\naWGad4Kkr0u6R9Ldki6qZfFmZladOZM1kNQCXA2cDQwD2yV9KSJ25dosBzoi4qWSzgSuARYDvwHe\nGxE7JB0N/EDS5vy6NrsGBgbp6trA0NAIbW0t9PauYcGC9nqXZWazZNLQBxYBeyJiEEDSLcBKIB/c\nK4EbASJim6R5kuZHxEPAQ2n+k5J2Am0l69osGRgYZOnSj9Lfvx6YCxxg69Zutmy50MFvVhCVdO+0\nAXtz0w+meRO1GSptI+nFwEJgW7VFWm10dW3IBT7AXPr719PVtaGOVZnZbKrkTH/aUtfOZ4GLI+LJ\n8dr19PQcfN3Z2UlnZ+eM11YkQ0MjjAX+qLkMD4/Uoxwzq1JfXx99fX3T2kYloT8EnJSbPiHNK21z\nYrk2kuaQBf6nIuJLE+0oH/pWe21tLcABDg3+A7S2ehCXWTMoPRlev3591duo5Ld9O3CypHZJRwDn\nAhtL2mwEzgeQtBh4LCL2p2WfBO6NiCurro6sH3r16vUsWdLN6tXrGRgYnMpmDOjtXUNHRzdZ8AMc\noKOjm97eNXWrycxmlyJi8kbSMuBKsjeJ6yLiQ5LWAhER16Y2VwPLyBJlTUTcKeks4JvA3UCkr/dH\nxFfL7CNKayl34bGjwxcep2N09M7w8AitrR69Y9bMJBERqmqdSkJ/NpQL/dWr1/PpT/9PSrsjVq26\nnJtu6p7V+szMGs1UQr+hO3N94dHMrLYaOvTHLjzm+cKjmdlUNXR6+sKjmVltNXSfPvjCo5nZeA67\nC7lmZja+w+5CrpmZ1ZZD38ysQBz6ZmYF4tA3MysQh76ZWYE49M3MCsShb2ZWIA59M7MCceibmRWI\nQ9/MrEAc+mZmBeLQNzMrEIe+mVmBOPTNzArEoW9mViAOfTOzAnHom5kViEPfzKxAHPpmZgXi0Dcz\nK5A59S7AmsvAwCBdXRsYGhqhra2F3t41LFjQXu+yzKxCioh61wCApGiUWqy8gYFBli79KP3964G5\nwAE6OrrZsuVCB79ZHUgiIlTNOu7esYp1dW3IBT7AXPr719PVtaGOVZlZNSoKfUnLJO2StFvSJeO0\nuUrSHkk7JJ2em3+dpP2S7qpV0VYfQ0MjjAX+qLkMD4/Uoxwzm4JJQ19SC3A18GbgNOA8SaeWtFkO\ndETES4G1wP/LLb4+rWtNrq2tBThQMvcAra3+wGjWLCr5bV0E7ImIwYh4CrgFWFnSZiVwI0BEbAPm\nSZqfpm8HHq1dyVYvvb1r6OjoZiz4sz793t41davJzKpTyeidNmBvbvpBsjeCidoMpXn7p1WdNZQF\nC9rZsuVCurouZ3h4hNbWFnp7fRHXrJl4yKZVZcGCdm66qbveZZjZFFUS+kPASbnpE9K80jYnTtJm\nUj09PQdfd3Z20tnZWe0mzMwOW319ffT19U1rG5OO05f0LOA+4GxgH/A94LyI2JlrswJ4T0ScI2kx\ncEVELM4tfzHw5Yh41QT78Th9M7MqzMg4/Yh4GlgHbAbuAW6JiJ2S1kr689RmEzAg6X7gE8Bf5Iq6\nGfgO8DJJP5F0QTUFmplZ7fiOXDOzJuU7cs3MbEIOfTOzAnHom5kViEPfzKxAfHOWzTo/k9+sfjx6\nx2aVn8lvVjsevWMNz8/kN6svh77NKj+T36y+HPo2q/xMfrP68m+azSo/k9+svnwh12bd6OidsWfy\ne/SO2VRM5UKuQ9/MrEl59I6ZmU3IN2dZU/INXmZT4+4dazq+wcss4+4dKwTf4GU2dQ59azq+wcts\n6hz61nR8g5fZ1Pm3xJqOb/AymzpfyLWm5Bu8zHxzlplZoXj0jpmZTcihb2ZWIL4j1wrLd/VaEblP\n3wrJd/Xa4cB9+mYV8l29VlTu3rFCqsVdve4esmbk0LdCGrurNx/8ld/VW657aOtWdw9Z46vof7ik\nZZJ2Sdot6ZJx2lwlaY+kHZIWVrOu2Wyb7l29teoeGhgYZPXq9SxZ0s3q1esZGBisan2zqkXEhF9k\nbwz3A+3As4EdwKklbZYDX0mvzwS2VrpubhvR6G677bZ6l1AR11mZBx74caxa1RNLlnwgVq3qiQce\n+HHZduXq7Oz8QEA842vJkg9Utf+Ojr8KeDKt/2R0dPzVuHVM9n28+tXnT/h9VLKNzs6Jj0Ut1q/3\nz71SzVBnys1Jczz/VUnoLwb+LTd9KXBJSZtrgLflpncC8ytZN7dsyv9hZ0t3d3e9S6iI66ytcnWu\nWtWTC+s4GNqrVvVUvN1abOPQN47uKb1xTPfNp9r1x/u5T/eNp9bbaG9/w4y/AU53GzMV+n8CXJub\nXg1cVdLmy8DrctNbgDMqWTe3bMpnOrOlmUOqETVznbU4S6/Fp4VD3zi66/LmU+36M3U8a7+N6t9E\nZ/v7mEroz9SQzarGjY7xsDlrDgsWtLNly4WsWnU5S5Z0s2rV5VVfxK3FI6JrMQpputuoRQ21uEbS\nCNtohBomM+nNWZIWAz0RsSxNX0r27nJZrs01wG0RcWua3gW8EVgw2bq5bfjOLDOzKkWVN2dVMmRz\nO3CypHZgH3AucF5Jm43Ae4Bb05vEYxGxX9LDFaw7pcLNzKx6k4Z+RDwtaR2wmWw0znURsVPS2mxx\nXBsRmyStkHQ/2efVCyZad8a+GzMzm1DDPHvHzMxmXt2fvdOoN29Juk7Sfkl35eY9T9JmSfdJ+ndJ\n8+pZY6rpBElfl3SPpLslXdRotUo6UtI2SXemGrsbrcY8SS2S7pC0MU03XJ2Sfizph+mYfq+B65wn\n6TOSdqb/o2c2Wp2SXpaO4x3p38clXdRodaZa/1LSjyTdJenTko6ots66hr6kFuBq4M3AacB5kk6t\nZ00515PVlXcp8LWIOAX4OvDXs17VM/0GeG9EnAb8HvCedAwbptaI+DWwJCJOBxYCyyUtaqQaS1wM\n3JubbsQ6R4DOiDg9IhaleY1Y55XApoh4OfBqYBcNVmdE7E7H8QzgNWRd1F+gweqU1ApcCJwREb9D\n1j1/HtXWWe0Yz1p+UcXNW3Wqrx24Kze9C5ifXr8I2FXvGsvU/EXgTY1aK3AU8H3gtY1YI3AC2X0m\nncDGRv25AwPAC0rmNVSdwDFAf5n5DVVnSW1/AHyrEesEWoFB4Hkp8DdO5Xe93t07bcDe3PSDaV6j\nOj4i9gNExEPA8XWu5xCSXkx2Jr2V7D9Bw9SaukzuBB4CtkTE9karMfkI8D4gf7GrEesMYIuk7ZL+\ne5rXaHUuAB6WdH3qOrlW0lE0Xp15bwNuTq8bqs6IGAY+DPwEGAIej4ivUWWd9Q79ZtcwV8ElHQ18\nFrg4Ip7kmbXVtdaIGImse+cEYJGk08rUVNcaJZ0D7I+IHUx8g2Ej/NzPiqw7YgVZl97rabDjSXY2\negbwsVTrAbJP841WJwCSng28BfhMmtVQdUo6FlhJ1gPRCsyVtKpMXRPWWe/QHwJOyk2fkOY1qv2S\n5gNIehHw0zrXA4CkOWSB/6mI+FKa3ZC1RsQvgD5gGY1X41nAWyQ9APwz8PuSPgU81GB1EhH70r8/\nI+vSW0TjHc8Hgb0R8f00/TmyN4FGq3PUcuAHEfFwmm60Ot8EPBARP4+Ip8muO7yOKuusd+gfvPFL\n0hFkN29trHNNeeLQM76NwJr0+h3Al0pXqJNPAvdGxJW5eQ1Tq6TjRkcUSHoOsJTsoXwNUyNARLw/\nIk6KiJeQ/V/8ekT8GdmzpdakZnWvU9JR6ZMdkuaS9UPfTeMdz/3AXkkvS7POBu6hwerMOY/szX5U\no9X5E2CxpN+SJLLjeS/V1tkAF06WAfcBe4BL611Prq6bgWHg1+lgX0B2AeVrqd7NwLENUOdZwNNk\nj62+E7gjHdPnN0qtwKtSXTuAu4C/SfMbpsYyNb+RsQu5DVUnWV/56M/77tHfm0arM9X0arKTux3A\n54F5DVrnUcDPgOfm5jVind1kJ0x3ATeQPbK+qjp9c5aZWYHUu3vHzMxmkUPfzKxAHPpmZgXi0Dcz\nKxCHvplZgTj0zcwKxKFvZlYgDn0zswL5/81pS5ToJ2XGAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f3d180e6f28>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(X,Y,'bo')\n",
    "plt.title('Wind Speed Distribution (km/h)')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "SciPy has a subpackage called [scipy.optimize](http://docs.scipy.org/doc/scipy/reference/optimize.html) which contains many methods for optimization and root finding. We will use the [curve_fit](http://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.curve_fit.html#scipy.optimize.curve_fit) method to fit the wind speed distribution to the Weibull distribution and so let's import the method:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from scipy.optimize import curve_fit"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can display the documentation for `curve_fit` in the notebook using `?`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "curve_fit?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We have the arrays `X` and `Y` for the fit and so we just need to define a function for the Weibull distribution. Notice that the first variable is the independent variable $x$ and the following parameters are the coefficients in the Weibull distribution:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def f(x,k,ell):\n",
    "    '''Weibull distribution with parameters k and ell'''\n",
    "    return np.where(x < 0, 0., (k/ell)*(x/ell)**(k-1) * np.exp(-(x/ell)**k))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We've used the NumPy function [numpy.where](http://docs.scipy.org/doc/numpy/reference/generated/numpy.where.html) to define a step function. Let's plot our function to make sure it corresponds with the [definition](https://en.wikipedia.org/wiki/Weibull_distribution):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7f3d1811eb70>"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAe8AAAFCCAYAAADYC1X9AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3Xl4VNX9x/H3yT4hC4SEsIclEUQUlC2yI6iIGlBBBFRE\nRUuV0qKVFsSC7U9FcSlFbUtlUxZxRVlEEQOygwooCgQCIQQCJAGyMFnn/P64E5hskGQmsyTf1/Pc\nZ7Z7z/1O0Hxyzr33XKW1RgghhBCew8vVBQghhBCiaiS8hRBCCA8j4S2EEEJ4GAlvIYQQwsNIeAsh\nhBAeRsJbCCGE8DAS3kI4iVLqXaXUNOvzfkqpZDvasiil2lifL1BKvWhHWy2UUplKKVXdNkq157Dv\nWU7bvZVSvzmqPSE8lYS3EBVQSv1FKbWm1HsJSqnVpd47pJS6/2rtaa0naK3/z/YtO8qr1LZKqbFK\nqUJrOGcqpY4opeYrpWJs6krWWofoq0z6YG3r+6sW5sDvaftHirXtzVrra6vbnhC1hYS3EBXbBNxc\n3CNVSjUGfIAbS73X1rquM1Wll7xVax0ChAKDADPwg1KqQzX2ebWAd/TvFJlFSohySHgLUbFdgB/Q\n2fq6D/AdcLDUe0e01qkASqn2SqmvlVLpSqnflFIjihsrZ3hbKaX+qpQ6q5RKVEqNtvngO6XUozav\nK9XrvRJtOKq1fgrYCMywth1l7eF6WV8/Yu2hF/fURyml2gPvYvwxk6WUyrD5Tu8opVYrpbKA/o76\nnkqpjRh/MOyz1jKi9DC89ef9nVLqnFLqZ6XU3aV+3nOVUqus229TSrW252cohLuQ8BaiAlrrAmAH\n0Nf6Vl+MHvbmct5DKRUIfA18AIQDDwDvWIOvPI2BMKAp8AjwX9vh7PJKqu53KcenGH94lGjb+h3+\nCdxu7a33BPZorQ8AvwO2aa2DtdZhNtuOAv6utQ4GtpSzr2p9T611P+vr663D+h+VqtUH+BL4CogA\n/gAsKdX2SOBvQH3gCGA7nC+Ex5LwFuLKNnI5qPsA31MyvPtY1wG4CziqtV5s7eXuBT4BRlA+DUzX\nWhdorTcBq4GrHjt3kJMYgVqeIuB6pVSA1vq01vpqJ4it1FpvB9Ba55Xzub3fs6JDBDcD9bTWs7TW\nhVrr74BVGH9MFPtMa/2D1toCLOHyiIkQHk3CW4gr2wT0Vko1AMK11keArUBP63sduXy8OwqIVUpl\nWJdzwGggsoK2z2mtc21eJ2H0Tp2hGZBR+k2t9UWM3uoE4JRS6kulVLurtHW1s8lr6ns2KWffSRjf\nrViqzfOLQJAD9iuEy0l4C3Fl2zCGXMdjHRLWWmdh9FzHAyla6yTruslAvNY6zLo0sA73Pl1B2w2U\nUiab1y2t7QLkAIE2nzV2zNe55B6MUYQytNbfaK1vs+7zIPDf4o8qaOtqw/k19T1PAi1KvdcSSKlC\nG0J4JAlvIa7A2mPcDUymZNhtsb5ne5b5KuAapdSDSikfpZSvUqrrFXquCphpXa8PcCewwvrZHuBe\npZRJKRUNPGbH1yg+M95LKdVKKfUvoB/WE9ZKrdNIKRVnPfZdAGQDFus6p4HmSinfauy/ut8zFWhD\n+XYAF5VSz1l/3v0xDl0sq2J9QngcCW8hrm4jxglRm23e+976XvHxbrTW2cBtGCeqnbQurwD+FbR7\nCjhnXe994EmtdYL1szcxwjMVWIBxEpytqpy8FquUygQuYJwtHwR001r/Wk57Xhh/lKQAaRjH9idY\nP9sA7AdSlVJnqrB/e77nDGCx9TDEcNsPrCcU3g0MsdY6F3jIpm25zEzUWuoq8zJUrhGlBgNvYfyP\n/57Welapz0Mw/qdsCXgDr2utF9q9YyGEEKIOsju8rdeGHgIGYvxlvQt4wHppSfE6fwVCtNZ/VUqF\nYxxHi9RaF9q1cyGEEKIOcsSweXcgQWudZB3GWg4MLbWOBoKtz4OBdAluIYQQonocEd7NKHm5xglK\nXqoBxrGoDkqpk8BeYJID9iuEEELUSc46Ye124CetdVPgRuBtpZRcbymEEEJUg48D2kjBOBGtWHPK\nXmc5DngZQGt9RCl1FGiPcQlOCUopOUNUCCFEnaK1rtIteR3R894FRFtvbuCHcZnMF6XWScK4mxFK\nqUjgGiCxoga11rJcZfnb3/7m8ho8YZGfk/ys5OckPyt3X6rD7p631rpIKfU0xg0Zii8V+00p9aTx\nsf4v8A9goVJqn3Wz57TWZaZmFEIIIcTVOWLYHK31V0C7Uu/9x+b5KYzj3kIIIYSwk8yw5qH69+/v\n6hI8gvycKk9+VpUjP6fKk59VzXHIDGuOpJTS7laTEEIIUVOUUugqnrDmkGFzIYQQjteqVSuSkpKu\nvqLwCFFRURw7dswhbUnPWwgh3JS1R+bqMoSDVPTvWZ2etxzzFkIIITyMhLcQQgjhYSS8hRBCCA8j\n4S2EEEJ4GAnvWiYvD154Ae66Cz780NXVCCGEqAkS3rXM738PP/4IDz8M06fDzJmurkgIISpn7969\nPPvss5Vad+nSpbz++uuMHDmS5cuX13BlV1aVuh1FrvOuRbZtg6++ggMHIDgY+veHbt3gppvg7rtd\nXZ0QorZJTEzkiSee4L777mPChAl2tfXGG2+wefNm6tevf9V1jxw5Qnp6Os888wxpaWnExMQQGxtL\nq1atrrrtnDlzSE1N5aWXXrKr3mJVqduRpOddi7z1FkydagQ3QKNGsHgxPPUUZGe7tjYhRO3Tpk0b\nQkNDGTRokN1tTZ48maFDh1Zq3f379/Paa68BEB4eTnR0NLt3l7nDdLkmTpzIihUrOHPmTLVrtVWV\nuh1Jet61RFaW0et+552S7/frB336wJtvGsPoQgjhKBaLhaNHjxITE1Pms8TERObNm1diYpLi50op\nYmNjiYuLq9Z+hwwZwpo1ay69PnXqFNHR0ZXaVinFmDFjWLx4cZmh7pqs2dEkvGuJL780Qrphw7Kf\nvfAC9O4Nf/zj5V65EMLzqSrNyVWx6k7itnv3brp160ZSUhJ79+5l4sSJHDhwAJPJRJs2bXj55Zcd\nU2ApPj4+dOzYEYBVq1bRtWtXOnfuXOntx44dS1xcXJnwrsmaHU2GzWuJTZvg1lvL/6xdOxg0qGyv\nXAjh2bR2zFJd69evJy8vj8zMTOLi4i4Ft7NkZmayaNEiPvjggyptl5aWhtlsZteuXTVUWc2Tnnct\nsWULPP54xZ9PmWKctPbMM+Aj/+pCCAfYsGEDkydPZvr06cyePbvE0LXtELQtRw5Bv/rqq8ybN4+g\noCCSkpKIioq66jbr1q0jISGB559/nvnz59OtWzen1uwo8mu8Fjh/Ho4dg06dKl6nc2do0QJWrYJh\nw5xWmhCiljKbzWRmZjJkyBA2b97M/v372blzJ6NHjwaqPwRd+sYdiYmJtG7dukygzp07l3vuuYe8\nvDx27dqF2WwmKiqKcePGoZRi/vz5ZdpetmwZe/bsYdasWWRnZ/PCCy/w1ltv4e/vb1fN5dVd02TY\nvBbYtQu6dAFf3yuvN2ECvPuuc2oSQtRu+/btY+DAgQD07NmTQ4cO0bRp02q3N3fuXN577z3i4+OZ\nOXMmWVlZAAwfPpw9e/aUWHfLli1MmjSJ7t2706RJE2JjY2nbti0AycnJ9O7du0z727dvZ/369cya\nNQuAoKAghg0bZvc14hXVXdPklqC1wD//CQkJMHfuldfLzYXmzeGHH6ASo0tCCBeTW4JWTUFBAZ07\nd2bfvn14e3u7upwy5JagooSDB42T0q4mIACGD4dly2q+JiGEcDZfX1/279/vlsHtaBLetcCBA9C+\nfeXWffBBeP99+84wFUII4VoS3rVAVcK7Z0+4eBH27q3ZmoQQQtQcCW8Pl5lpLM2aVW59Ly8YMwaW\nLKnZuoQQQtQcCW8Pd/gwtG1rhHJljR4Ny5fL0LkQQngqCW8Pl5JinEFeFR06QFAQVHIefyGEEG5G\nwtvDpaRUfsjc1r33wqefOr4eIYQQNU/C28PZE96ffCJD50II4YkkvD1cSgpUZ1Kjm26CvDzYv9/x\nNQkhhKhZEt4e7uTJ6vW8lZKhcyGE8FQS3h6uusPmAPfdJ+EthHAfe/fuLXOP7av56aefeOWVV2qo\noqtbunQpr7/+OiNHjrR7nvSqkPD2cPaE9803w6lTxh3JhBCiqhITExk0aBDvOuCOR2+88QYzZ84k\nIyOj0ttorZk+fTr5+fmV3mbOnDlMnTq1OiWWceTIEdLT03nmmWd4++23mTBhAsec9AtVwtuDmc3G\nbGkNG1Zve29vuOMOWLPGsXUJIeqGNm3aEBoayqBBg+xua/LkyQwdOrRK23zyyScMGDCgSttMnDiR\nFStWcObMmSptV579+/fz2muvARAeHk50dDS7nXQNrtzP24OlpUF4uHH8urruvBMWLIDf/95xdQkh\n6gaLxcLRo0eJiYkp81liYiLz5s0rcSet4udKKWJjY4mLi6v2vtPT0/Hy8iI8PJycnJxKb6eUYsyY\nMSxevLjMEH1Vax4yZAhrbHo/p06dIjo6utrfqSokvD3YuXPQoIF9bdx2Gzz2mNGDDwx0TF1CCOdQ\nM+34y92G/lv1rhndvXs33bp1Iykpib179zJx4kQOHDiAyWSiTZs2vPzyyw6przyffvop48ePZ9Gi\nRVXeduzYscTFxZUJ76rW7OPjQ8eOHQFYtWoVXbt2pXPnzlWupzokvD2YI8I7NBS6dIENG+CuuxxT\nlxDCOaobuo6yfv168vLyyMzMJC4ujltvvRWTyVTj+925cyc9evSo9vZpaWmYzWZ27dpFt27d7K4n\nMzOTRYsW8cEHH9jdVmVJeHswR4Q3GEPnq1ZJeAshqmbDhg1MnjyZ6dOnM3v27BJDxrZD0LYcMWy+\nY8cOzGYza9euZcuWLeTm5vLFF19Uqr1169aRkJDA888/z/z580uEd3VrfvXVV5k3bx5BQUEkJSUR\nFRVVre9VFUq72RRbSintbjW5qwULYONGWLjQvnYOHIBbb4Xjx+07fi6EcCzbY6/uxmw2069fP3bu\n3MnUqVPp0aMHOTk5jB49utptLlq0iPj4eBYsWHDpvcTERFq3bl0mUIvNnDkTpRQvvPACAOPGjUMp\nxfz588usu2zZMvbs2cOsWbPIzs7m2muv5fDhw/j7+1e75rlz53LzzTfTvHlzjh8/jtlspm/fvuWu\nW9G/p/X9Kv32lbPNPVhGhmN63u3agZ8f/Pyz/W0JIeqGffv2MXDgQAB69uzJoUOHaFqd6R6t5s6d\ny3vvvUd8fDwzZ84kKysLgOHDh7Nnz55yt/noo49YuXIlK1eu5OOPPwYgOTmZ3r17l1l3+/btrF+/\nnlmzZgEQFBTEsGHD7Lo2e8uWLUyaNInu3bvTpEkTYmNjadu2bbXbqwrpeXuw5583Qtf6B6dd/vAH\naNwYHHT5oxDCAdy55+2OCgoK6Ny5M/v27cPb29vV5ZQhPW8BOO6YNxjHvVevdkxbQgjhCr6+vuzf\nv98tg9vRJLw9mCPDu18/+OUXSE93THtCCCFqjoS3B3NkeAcEGAH+zTeOaU8IIUTNkfD2YI4Mb4Db\nb4d16xzXnhBCiJoh4e3Baiq85fwYIYRwbxLeHszR4R0dDSaTXDImhBDuTsLbg124YExv6kgydC6E\nEO5PwttDFRRAUZFxopkjSXgLIYT7c0h4K6UGK6UOKKUOKaWmVLBOf6XUT0qpX5RS3zliv3VZTg7U\nq+f46UwHDIAdO4z2hRBCuCe7w1sp5QXMBW4HrgNGKaXal1onFHgbuEtr3REYYe9+67qcHAgKcny7\nISFw003GnOlCCCHckyN63t2BBK11kta6AFgODC21zmjgE611CoDWOs0B+63TsrONnndNGDxYhs6F\nEKIyNm3aRG5uLnl5eXz//fdO268jwrsZkGzz+oT1PVvXAGFKqe+UUruUUg85YL91WvGweU2Q495C\nCFfYu3cvzz77bKXWXbp0Ka+//jojR4606+Yi9ho7diyBgYFERUWRkZHhtP06637ePsBNwC1APWCb\nUmqb1vqwk/Zf69TUsDlA587GZWjHjkGrVjWzDyGE50tMTOSJJ57gvvvuY8KECXa19cYbb7B582bq\n169/1XWPHDlCeno6zzzzDGlpacTExBAbG0urSvzCmjNnDqmpqbz00kt21Vts2rRpDB48mCZNmjh1\nTnVHhHcK0NLmdXPre7ZOAGla61wgVym1CegElBveM2bMuPS8f//+9O/f3wFl1i41OWzu5WXc33vd\nOnjyyZrZhxDC87Vp04bQ0FAGDRpkd1uTJ0+mYcOGbKzECTf79+/ntddeY+LEiYSHhxMdHc3u3bsr\nFd4TJ04kJiaGP/7xjzRq1Mjuun19fWnevHmVtomPjyc+Pt6u/ToivHcB0UqpKOAU8AAwqtQ6K4F/\nKaW8AX+gB/BGRQ3ahrcoX032vMEYOv/sMwlvIUTFLBYLR48eJSYmpsxniYmJzJs3r8RtMIufK6WI\njY0lLi6uWvsdMmQIa9asufT61KlTREdHV2pbpRRjxoxh8eLFZYboq1Pzzp070VqTnp5OTExMpb5T\n6U7pzJkzK1V7CVpruxdgMHAQSAD+Yn3vSeAJm3WeBfYD+4CJV2hLi6tbuFDrhx+uufZTU7UODdU6\nP7/m9iGEuLKr/j40ZjO2f6mmHTt26CeeeEIfO3ZMr1y5Urds2VJfvHix2u0tXLhQjxs3rkrbfPnl\nl3ro0KFV2ubIkSP6uuuuq9I2Ffnxxx8vPe/UqZM+f/58hetW9O9pfb9KueuQ67y11l9prdtprWO0\n1q9Y3/uP1vq/NuvM1lpfp7W+QWv9L0fsty67eBECA2uu/chIaN3auOZbCOGmHBXf1bR+/Xry8vLI\nzMwkLi6OAwcOYDKZHPgFrywzM5NFixbxwQcfVGm7tLQ0zGYzu3btsruGTp06XXreoEEDu4fDK8tZ\nJ6wJB8vNNeYhr0nFZ5337l2z+xFCeKYNGzYwefJkpk+fzuzZs0sMXdsOQdvSDhg2L/bqq68yb948\ngoKCSEpKIioq6qrbrFu3joSEBJ5//nnmz59Pt27dql3zkiVLWLNmDUuWLAEgOzvbaSetKe1mt5BS\nSml3q8kdvfwyZGYajzXlu+9gyhTYubPm9iGEqJjtsVd3Yzab6devHzt37mTq1Kn06NGDnJwcRo8e\nXe02Fy1aRHx8PAsWLLj0XmJiIq1bty4TqHPnzuXmm2+mefPmHD9+HLPZTN++fRk3bhxKKebPn1+m\n/WXLlrFnzx5mzZpFdnY21157LYcPH8bf379a9W7ZsoWioiL69u1LTk4OHTt2ZP/+/QRWMCxa0b+n\n9f0qzZcpc5t7qNxcx89rXlqvXnDwIKTJlDpCiFL27dvHwIEDAejZsyeHDh2iadOm1W5v7ty5vPfe\ne8THxzNz5kyysrIAGD58OHv27Cmx7pYtW5g0aRLdu3enSZMmxMbG0rZtWwCSk5PpXc5w4fbt21m/\nfj2zZs0CICgoiGHDhtl1jXivXr1ITk7mrbfeYtq0aSxfvrzC4HY06Xl7qClTICzMeKxJcXEwapSx\nCCGcy5173u6ooKCAzp07s2/fPqdec11Z0vMW5OVBNUd6quT22+Grr2p+P0IIYS9fX1/279/vlsHt\naBLeHsoZw+ZgzHP+9dd2nZAqhBDCwSS8PZSzwrttW2MymH37an5fQgghKkfC20M5a9gcZOhcCCHc\njYS3h3JWzxuMoXMJbyGEcB8S3h7KmeHdvz/s3g3WKzeEEEK4mIS3h8rNdd6weVAQdO9uTNoihBDC\n9SS8PVRenvN63iBD50II4U4kvD2UM4fN4XJ4yyVjQgjhehLeHsqZw+YAHTsavf3Dh523TyGEEOWT\n8PZQ+fng5+e8/Skll4wJIURpmzZtIjc3l7y8PL7//nun7VfC20MVFoKvr3P3OXiwcYtQIYSoCXv3\n7uXZZ5+t0jY//fQTr7zySg1VdHVjx44lMDCQqKgoMjIynLZfCW8PVVgIPk6+G/ugQbBpkzFkL4QQ\niYmJDBo0iHfffdfutt544w1mzpxZpQDUWjN9+nTy8/Mrvc2cOXOYOnVqdUos17Rp0zh+/DgpKSkM\nHTrUYe1ejYS3h3JFeIeFwXXXwebNzt2vEMI9tWnThtDQUAYNGmR3W5MnT65y+H3yyScMGDCgSttM\nnDiRFStWcObMmSptVxFfX1+aN2/u9JuhOPnXv3AUV4Q3XB46d8D/q0IID2exWDh69CgxMTFlPktM\nTGTevHklboNZ/FwpRWxsLHFxcdXed3p6Ol5eXoSHh5OTk1Pp7ZRSjBkzhsWLF5cZoq9OzTt37kRr\nTXp6OjExMXZ9p6qQ8PZQrgzvxx+H115z/r6FECWp+HiHtKP796/Wdrt376Zbt24kJSWxd+9eJk6c\nyIEDBzCZTLRp04aXX37ZIfWV59NPP2X8+PEsWrSoytuOHTuWuLi4MuFdnZoff/xxbrzxRgA6d+5M\nv379CA0NrXJNVSXh7aFcFd5du8LJk3DiBDRv7vz9CyEuq27oOsr69evJy8sjMzOTuLg4br31Vkwm\nU43vd+fOnfTo0aPa26elpWE2m9m1axfdunWzq5ZOnTpdet6gQQPi4+OdcuxbwttDuSq8vb3h1luN\ne3w/+qjz9y+EcB8bNmxg8uTJTJ8+ndmzZxMdHX3pM9shaFuOGDbfsWMHZrOZtWvXsmXLFnJzc/ni\niy8q1d66detISEjg+eefZ/78+SXCu6o1L1myhDVr1rBkyRIAsrOznXbsW2k3mzJLKaXdrSZ35ONj\nnPXtigBfuBDWrIEVK5y/byHqEttjr+7GbDbTr18/du7cydSpU+nRowc5OTmMHj262m0uWrSI+Ph4\nFixYcOm9xMREWrduXSZQi82cOROlFC+88AIA48aNQynF/Pnzy6y7bNky9uzZw6xZs8jOzubaa6/l\n8OHD+FdzxqstW7ZQVFRE3759ycnJoWPHjuzfv5/AwMBy16/o39P6fvlfsAJytrkH0hqKioxesCvc\ndhusX2/0/oUQddO+ffsYOHAgAD179uTQoUM0bdq02u3NnTuX9957j/j4eGbOnEmW9TaGw4cPZ8+e\nPeVu89FHH7Fy5UpWrlzJxx9/DEBycjK9e/cus+727dtZv349s2bNAiAoKIhhw4axfPnyatfcq1cv\nkpOTeeutt5g2bRrLly+vMLgdTXreHqiw0JjX3JXh2akT/PvfcPPNrqtBiNrOnXve7qigoIDOnTuz\nb98+p1+6VRnS867jXHW825bcZUwI4W58fX3Zv3+/Wwa3o0l4eyAJbyGEqNskvD2QO4R3r15w4ACk\npbm2DiGEqIskvD2QO4S3nx/0729cMiaEEMK5JLw9kDuEN8Bdd8GqVa6uQggh6h4Jbw/kLuE9ZIhx\n3FsuGRNCCOdygwgQVeUu4d2sGbRqBVu3Qt++rq5GiNonKiqqwslJhOeJiopyWFtuEAGiqtwlvMEY\nOl+9WsJbiJpw7NgxV5cg3JQMm3sgdwtvOe4thBDOJeHtgdwpvLt2NS4XS0x0dSVCCFF3SHh7IHcK\nby8v48S11atdXYkQQtQdEt4eyJ3CGy4f9xZCCOEcEt4eyN3C+9ZbYcsWyM52dSVCCFE3SHh7IHcL\n75AQiI01bhMqhBCi5kl4eyB3C2+Qs86FEMKZJLw9kDuG9513wpo1YLG4uhIhhKj9JLw9kDuGd3S0\nMXz+44+urkQIIWo/CW8PVFgI7niv+aFD4fPPXV2FEELUfhLeHqioyD3De9gwCW8hhHAGCW8PpDW4\n470KevSA9HRISHB1JUIIUbtJeHsodwxvLy8ZOhdCCGeQ8PZAWru6gorJ0LkQQtQ8h4S3UmqwUuqA\nUuqQUmrKFdbrppQqUErd64j91lXuOmwOMGAA/PorpKa6uhIhhKi97A5vpZQXMBe4HbgOGKWUal/B\neq8A6+zdZ13nzuHt7w+DB8MXX7i6EiGEqL0c0fPuDiRorZO01gXAcmBoOetNBD4Gzjhgn3WaO4c3\nyNC5EELUNEeEdzMg2eb1Cet7lyilmgLDtNbvAm4cO57B3cP7jjtg82bIzHR1JUIIUTs5a56utwDb\nY+FXjJ4ZM2Zcet6/f3/69+9fI0V5KncP75AQ6N0b1q6FkSNdXY0QQriX+Ph44uPj7WpDaTtPXVZK\nxQIztNaDra//Amit9SybdRKLnwLhQA7whNa6zJFRpZS2t6ba7sMP4ZNPYMUKV1dSsf/+FzZsgOXL\nXV2JEEK4N6UUWusqdckcMWy+C4hWSkUppfyAB4ASoay1bmNdWmMc9/59ecEtKsfde95gXO/91Vdg\nNru6EiGEqH3sDm+tdRHwNPA1sB9YrrX+TSn1pFLqifI2sXefdZ0nhHdkJNx4I6yTawuEEMLhHHLM\nW2v9FdCu1Hv/qWDdRx2xz7rME8Ib4P77jaH9YcNcXYkQQtQuMsOaB/KU8L73XuMe3zJ0LoQQjiXh\n7YE8JbwjI6FLF+PYtxBCCMeR8PZAnhLeACNGuPdZ8UII4YkkvD2QK8Jba02GOYMzOWe4kHuByl7O\nd++9xvXeMnQuhBCO46xJWoQDOSO8C4oKWHt4LasPrWZL8hYOZxwmwCcAX29f8grzsGgLbcPaEtss\nll4tezEkZgjhgeFl2mnUCLp2NQL8XrkdjRBCOISEtweqyfDOyc9h7s65vLH9DWLCYhjRYQTju4yn\nQ0QHAn0DL62XmZfJwbSDbDuxjZUHVzJx7UR6NOvBozc+yn3X3oevt++ldYuHziW8hRDCMeyeYc3R\nZIa1q1uwADZuhIULHdvu2oS1jP9yPD1b9GRG/xl0iOhQ6W1z8nNYdWgV7+5+l4SMBCb1mMRT3Z6i\nnl89zp6F6Gg4dQoCA6/elhBC1CWummFNOJmje96FlkL+sPYPTFg9gQ/u/YAVI1ZUKbgB6vnVY2TH\nkcQ/Es/aMWvZfXI30f+KZs6OOYSG5dO9u3HZmBBCCPtJeHsgR4Z3Tn4O93x4DwfTD7Lnd3vo36q/\n3W3eEHkDK0asYO2YtaxJWMON/7mRm+7ZyJIl9tcrhBBCwtsjOSq8cwtzuXvZ3TQIaMCqUauoH1Df\n/kZtdG7pL9G2AAAgAElEQVTcmbVj1vKPAf9gae5DrPZ7mEMpZx26DyGEqIskvD2UveFdZCnioc8e\nomFgQxYMXVDiBDNHUkpxz7X38NvTv9I6shHd53di9aHVNbIvIYSoKyS8PZAjet4z4mdwNucs79/z\nPt5e3o4p7AqC/IKYfetsWu5cztNrn+Z3q35HTn5Oje9XCCFqIwlvD2RveK9PXM/8PfP5cPiHBPgE\nOK6wq7j9dji1vS8rB+/FXGim+/+6cyDtgNP2L4QQtYWEtwey50q6tItpPPzZwywetpjIoEjHFVUJ\nfn7Gnca++CiERcMWMTl2Mn0W9GHFfpk/VQghqkLC2wPZ0/N+7pvnuP+6+xnYZqBji6qkhx6CDz4w\nvsNjNz3G1w9+zV+//SuT1k6ioKjAJTUJIYSnkfD2QNUN783HN/P1ka/5+4C/O76oSurRAwoL4Ycf\njNc3NrmR3eN3k5CRwJClQzhnPuey2oQQwlNIeHug6oR3kaWIp9Y8xZu3v0mwf3DNFFYJSsGDDxq9\n72INTA34YtQXXBdxHTe/dzMJ6Qkuq08IITyBhLcHqk54r9i/gnq+9RjeYXjNFFUFY8bAsmVQYDNK\n7uPlw1uD3+JPsX+i94LexB+Ld1l9Qgjh7iS8PVBVw7vQUsiMjTN4ccCLKDe4EXhMDFxzDawu53Lv\nJ7s+ydJ7lzLy45Es3LPQ6bUJIYQnkPD2QFUN76U/L6VxUGMGtnbNSWrleewx+N//yv9sYJuBbHxk\nIzPiZ/Dy9y9X+t7hQghRV0h4e6CqhLfWmlc2v8ILfV9wi153sREjYOtWSEkp//P24e3Z+thWlv2y\njElfTcKiLc4tUAgh3JiEtweqSnhvOLoBby9vbml9S80WVUX16hkBvmhRxes0DW7KpnGb2Ht6L6M+\nGUVeYZ7zChRCCDcm4e2BqhLec3fN5eluT7tVr7vYY4/Be++B5Qqd6voB9Vn34DoKLYUMWTqEzLxM\n5xUohBBuSsLbA1U2vJPOJ7EpaRNjbhhT80VVQ7duRg9848YrrxfgE8CK4Sto17Ad/Rf252yO3JlM\nCFG3SXh7oMqG97wf5/HQDQ8R5BdU80VVg1JXPnHNlreXN28PeZs7Y+6k78K+pGRWcLBcCCHqAAlv\nD1SZ8LZoCx/s+4BHOj/ilJqq68EHjUvGMjKuvq5Sir/f8nfGdR5HnwV9SDyXWPMFCiGEG5Lw9kCV\nCe8tx7cQ5BdEp8hOzimqmho2hLvuggULKr/Nc72e4889/0y/hf347exvNVecEEK4KQlvD1SZ8P5g\n3wc8eMODbnmiWmlPPQXvvnvlE9dKm9BtAi/d8hK3LL6Fn079VHPFCSGEG5Lw9kBXC+9CSyGfHviU\nBzo+4Lyi7BAbCyEhsG5d1bZ7qNNDvD3kbQYvGczW5K01U5wQQrghCW8PdLXw3pq8lRYhLWhVv5XT\narKHUkbv++23q77tvdfey+Jhixm6fCjrE9c7vjghhHBDEt4e6Grh/fmBzxnabqjzCnKAUaNg+3Y4\nerTq294efTuf3v8poz8ZzZcHv3R8cUII4WYkvD1UReGttWblwZUMbe9Z4R0YCGPHGse+q6NPVB/W\njFnD+C/H8+EvHzq2OCGEcDMS3h7oSj3v39J+o9BS6PZnmZdnwgTjrHOzuXrbd23alW8e+obJX09m\nwU9VOH1dCCE8jIS3B7pSeH9z5Btua3ObR5xlXlp0NHTvDkuWVL+N6yOv57ux3zFj4wzm7pzruOKE\nEMKNSHh7oCuF97dHv2VQm0HOLciBnnkG3nijapeNlXZNw2vY+MhG3tr+FrM2z3JccUII4SYkvD1Q\nReFdaClkU9Imt7uDWFUMGAABAbB2rX3ttKrfik3jNrFo7yKmb5gu9wQXQtQqEt4eqKLw3pWyi1b1\nWxFRL8L5RTmIUkbve/Zs+9tqGtyUjY9sZFXCKiavmywBLoSoNSS8PVBFGfTt0W8Z2Hqgc4upAfff\nD0eOwA8/2N9WRL0Ivhv7HdtTtvPkqicpshTZ36gQQriYhLcHqqjnve3ENnq37O38ghzM1xcmTYLX\nX3dMe/UD6vP1g1+TkJHA2M/HUmgpdEzDQgjhIhLeHqi88NZas/3Edno07+Gaohxs/HhjutSkJMe0\nF+wfzJrRa8gwZ3D/R/eTV5jnmIaFEMIFJLw9UHnhnZCRQJBfEE2Dm7qmKAcLCYHHH3fMse9iJl8T\nnz/wOV7Ki2EfDuNiwUXHNS6EEE4k4e2Bygvv7Se2c3Pzm11TUA2ZPNm45vvUKce16eftx/LhywkP\nDGfIkiFk5WU5rnEhhHASCW8PVF54b0veRmzzWNcUVEMiI40pU197zbHt+nj5sGjYItqHt2fQ+4M4\nZz7n2B0IIUQNk/D2QOWF946UHbUuvAH+/GdYuBDOnHFsu17Ki3fvfJfeLXozYNEAzuQ4eAdCCFGD\nJLw9UOnwLigq4Le037gh8gbXFVVDmjaF0aMdd+a5LaUUs2+bzdB2Q+m3sB8pmSmO34kQQtQACW8P\nVDq8D6UfomVoSwJ9A11XVA2aMgXmzYO0NMe3rZRi5oCZjOs8jj4L+nD0XDXuSSqEEE4m4e2BSof3\nvtP7uL7R9a4rqIa1aAEjRjj2zPPSnuv1HM/c/Ax9F/blQNqBmtuREEI4gEPCWyk1WCl1QCl1SCk1\npZzPRyul9lqXzUqp2ps0TlA6vH8+83OtHDK39fzzRu/75Mma28dT3Z/i/275P/ov7M/2E9trbkdC\nCGEnu8NbKeUFzAVuB64DRiml2pdaLRHoq7XuBPwDmGfvfuuyutbzBqP3/eij8OKLNbufhzs9zPyh\n84lbFsfqQ6trdmdCCFFNjuh5dwcStNZJWusCYDkw1HYFrfV2rfUF68vtQDMH7LfOqos9b4C//AU+\n/hgOHarZ/QyJGcIXo77gsS8eY+GehTW7MyGEqAZHhHczINnm9QmuHM6PA3be8LFusw3v87nnSb+Y\nTusGrV1blBM0bGhM3DJ9es3vK7Z5LPGPxDMjfgavbH5F7kgmhHArPs7cmVJqADAOuOLdM2bMmHHp\nef/+/enfv3+N1uVpbMP7t7O/cW3EtXipunHu4aRJEBNj3HGsS5ea3Vf78PZsfWwrgz8YzKmsU7w5\n+M0683MWQtSc+Ph44uPj7WpD2dujUErFAjO01oOtr/8CaK31rFLr3QB8AgzWWh+5QntaejlX9swz\n0KQJPPssvL/3fdYeXsvS+5a6uiyn+fe/4cMPYcOG8u+u5mjnc88zdPlQmgQ1YdGwRfj7+Nf8ToUQ\ndYZSCq11lX6bOaIbsQuIVkpFKaX8gAeAL0oV1hIjuB+6UnCLyisOrYSMBKLDol1bjJM9/jikp8Mn\nnzhnf/UD6rPuwXUUWAoYvGQwGeYM5+xYCCEqYHd4a62LgKeBr4H9wHKt9W9KqSeVUk9YV5sOhAHv\nKKV+UkrttHe/dZntsPnhjMPEhMW4tiAn8/GBf/7TGHkwm52zzwCfAFYMX8FNjW+i53s9OZIhf4MK\nIVzHIce8tdZfAe1Kvfcfm+fjgfGO2JcoG951recNMGAAdOtm3LTkhRecs09vL29ev/11YhrG0Gt+\nLz6+/2N6t7zi6RtCCFEj5OwbD1Qc3lprEjISiGlYt3rexWbPhjlz4Phx5+73d11/x6Jhi7j3w3tZ\n+nPdOddACOE+JLw9UHF4n8s9h0VbaGhq6OqSXCIqCp56yhg+d7bbo29nw9gNTP12Ki9ufFEuJRNC\nOJWEtwcqDu/kC8m0DG2JcsYp127qL3+BH3+EVaucv++OjTqy/fHtrDq0ioc+ewhzgZMOwAsh6jwJ\nbw9UHN4nMk/QIqSFq8txKZMJ/vMf+P3vISvL+ftvHNSY+EfisWgLfRb0IflC8tU3EkIIO0l4e6Di\nEdrkzGSahzR3bTFuYOBAY5k2zTX7D/QNZMm9S3ig4wP0+F8PNiVtck0hQog6Q8LbA0nPu6zZs+Gj\nj2C7i24GppTi2Z7PsnDYQkZ8NIK3d74tx8GFEDVGwtsDXTrmLT3vSxo2hDffNCZwyc11XR23tb2N\nrY9u5d8//JvHv3icvMI81xUjhKi1JLw9UImed6j0vIuNHAnt2xv3/naltmFt2fbYNi7kXaDPgj4c\nO3/MtQUJIWodp96YRDiG7dnmTh02Lyw05iU9cwbOnoULF+DiRcjJubzk5RkFWizGY/Hi7Q3+/sYS\nEFDyeXAwhIYaS0jI5Uefqv3nqZQx73mnTnDnncZELq4S5BfERyM+4s3tb9J9Xnfm3T2Poe2HXn1D\nIYSoBAlvD1Qc3ilZKTQLcfCt0XNy4OefjSUx8fJy7BicOwcNGkCjRhARYYRsvXolF39/8PIyCrRd\nioqMYD93zhjXzsszltxc4zTxCxdKLllZxqnkoaEQHk5ho0ZkNW1KZuPGZEZGkhkeTmaDBlwICSEn\nJIT8kBDy/f0psFi4daHmns81jzW24OOn8VIKLyj30U8pTF5emLy9jcfixfo62NubUB8f6vv4YPLy\nqtJleUopJt88mZ4tevLAxw+wMWkjrwx6BT9vP8f+mwkh6hy77yrmaHJXsat74gm4/qaLPJfeEPM0\nO64t1hqOHIHvvoONG2H3bmO6sg4d4PrrIToa2rQxlqgoCA+vcm/46iVozhYUcDQ3l6NmM8fz8kjN\nz+d0fj6pZjOnc3NJLSzkvMVCsNaEFBYSkpdHSG4uITk5hGRmUi8zE7/MTPyKivALCMDXZOLs+UBy\ndRAduwajQ0KwhIZiCQ7GEhSExcsLC1CkNfkWC2bbpajo0vOLRUVkFxVxoaiI84WFFGpNfWuQh3p7\nX3oe5utLI19fIvz8jEeb5+G+vvh6eZFhzmDs52M5m3OWD4d/SFT9KIf+HIUQnqs6dxWTnrcH0hpy\ndHr1ZlYrKoItW4xbcn3+uTEUPmAA9O8PU6YYB419fR1ec57FwoGLF/klJ4dfcnL4NSeHRGtg+3t5\n0ToggDYmEy39/Wnm789NQUE09vMj0ro09PXF+0q9Xq0hMxNSUiAlhbzEFN6dfpLbftpLhxDjPU6e\nNIb7IyONP0ZatSr5GBUFrVsbPf4KvsOFwkLOFxZeejxfWEhGYSFn8vM5Yjaz7cIFzhYUcKaggLP5\n+aQXFhLs7U2Ery9N2r9IflYiHb76Jw/E3MrgFl1oZv2+Tfz88POSU1CEEJUj4e2BtIaLOp3wwPDK\nb3T6NPzvf8aMJg0awPDhsGaN0ct28Axt+RYLe7Kz2Z6ZyY7MTH7MzuZYbi5tAgLoWK8eHevV4+HG\njWkbEEBrk4lQR/Tmlbp83LxDB/yBwf2gTx/4+mu48UbregUFRpAfOwZJScaybRssX248T042fj62\ngd6qFbRujX/btjRq1YpGgYGVLsuiNecKC42RhPx8UvKasLV+I5YeWseGjLOEN7iW0wWFpObnU9/H\nh2b+/jT186NlQABR/v5EBQTQKiCAqIAAGvv54VWHZ9MTQlwmw+Ye6LHHIPTG9ewNeZlvH/72yiuf\nPAl//7sRTsOHw4QJcNNNDq3HXFTElgsXWH/uHBsvXGBfdjbRJhOxISH0CAmhS3Aw7QMD8XdBz/LD\nD2HqVPjhB6hfvxIbFBVBaqoR5MUBf+wYHD1qHGI4cQKaNIG2bY3DCW3bllxCQytVV1ZeFpO+msT3\nx7/n/Xvep3uzHpzJzyclP5+UvDyS8/JIys3lWG4uSdblfGEhLayhXhzotuHe3N//yqMTQgi3VJ1h\ncwlvD/Too+B744eca/wJK0asKH+l/Hx46SX417+MDaZMMY5ZO8jBixdZmZbG1xkZ7MjKolO9egxs\n0IAB9evTNTiYIAcfG7fHH/5gZPBnnxnn0tmloMA4L+DIEeNEviNHSi7+/mUDvTjomzYtU8Cnv33K\n71f/ngldJzCt7zR8vCr+uZmLijiel3cp0G2D/VhuLmkFBbQMCKCtyUS0yURb6/O2JhNtAgII8Pa2\n88sLIWqChHcdMW4cFHR+m5Do/bxz5ztlV/jlFxg1yjh++8470Nz+iVy01vyQlcWnaWl8lpZGZmEh\nQ8PDGRIWRt/69Qlxo7AuLT/fOKw/cCC8+GIN7khr45i6bZjbBvyFC0aIx8SUWFKbBDN25185n3eB\n+XHzua7RddXafW5REUdzczlsNnPEbOZIbi5HzGYOm80cz80lws+PtgEBRrDbLgEB1K+B8xyEEJUj\n4V1HPPIIXOj8Ih1vKODvt/y95IcrVxrTjL32Gowda/fx7JS8PD44fZpFqankWSyMiIjgnogIugUH\ne9Tx1zNnoEcP+Mc/YMwYFxWRnQ2HD0NCwuVH66Kzssho1oCtAWcJu74HPQY8hE+79kbAR0ba/e9Y\npDXJxcFuE+pHrEuAlxfXBAZyjclEu8DAS8+jTSZM0mMXokZJeNcRY8fCyU5/4M6ebfhj7B8vf/D+\n+8Y9Mj/7DLp3r3b7Fq35OiODOSkpbM/M5L6ICMZGRtIrNNSjbz/6yy9wyy3Gj6dXL1dXU0pmJhw+\nTNrebaxZM4f6yWe5pbAFQUknjWvho6PL9NiJiTGut7fz30Rrzen8fA6ZzRy6ePHS40GzmaNmM439\n/MoN9pYBAXKMXQgHkPCuI8aOhUPXj+Gp2+7gwRseNN784gv43e/g22/h2mur1e7FoiIWpaYyJyWF\nAC8vJjVrxv2NGhFYi3peX31lHHbYvNk4FO2OtNYs+XkJz3z9DOM6j+OFTn8g8FhKiZ76paWoqPxQ\nj442Jny3M1wLLRaS8vKMMLcJ9kNmM2cLCmgbEHApzK8JDDTC3WQi3NfXo//QE8KZJLzriIcfhj0d\nBzPrnkncEXOHcTz15pvhyy+NseEqMhcV8e+TJ3k1OZnuwcFMbtGCvh7ey76S//4XXnkFvv8emjl4\ngjpHOp19mklfTWJHyg7+OfifxLWLK7tSenr5oZ6QYJwcV16wx8QYl8PZKaeoiMM2YX7w4sVLPXbg\nck/d+tguMJAYGYYXogwJ7zrioYdgc4eufPjwO3Rv0tUYAx41yjitugqKtGbeyZP8PSmJ7iEhzGjV\nik5BQTVUtXuZNQsWL4ZNm4wOqjtbn7iep9c8TUzDGP45+J+0adDm6htpDWlpFQe7n9/lofjo6JKL\nnT12rTXpBQUcLB5+twb6oYsXSczNpZGv76Uwb2fTY2/h7+9R51EI4SgS3nXEgw/CN9e2ZfPvviLm\nyy3GxCtbtlTpOqhN58/zh4QE6vv48Hp0NF2Cg2uwYvf0l7/Ahg3GkQZ3//r5Rfm8se0NZm+dzcTu\nE5nSewoBPgHVa0xr4wy+w4fLLgkJxjrlhXp0tDGvvR0BWzwMf9BmGL74+bnCQmKKw9ymt36NySRn\nw4taTcK7jhgzBlZdG8mBp3+gSZd+sGgR9O5dqW0zCgqYdPgwG8+fZ3bbtoyIiKi1w+NXozU89RTs\n2QNr11Z6fhWXOn7hOH9a9yd+OPkDLw98mZEdR+KlHDj5jdaQkVHyrHjbJT+//FCPjjYmr7Hjv6Ws\nwkISbML8oM2QfKCXV4kwL37eJiAAX5lWVng4Ce86YvRo+Kx9MOnt3iXwrbeN6T0rYW16OuMPHuS+\niAheatOGenLsEa2Now3bthnTqIaFubqiytl4bCPPfP0MXsqL1297nT5RfZyz43PnjHMsygv37OyK\ng71Zs2rPkKO15mR+folj6sUBn5KXR8uAgDJD8O1MJiL9/OrsH6bCs0h41xEPjNKsaO9D0Z441LBh\nxunnV1BgsfDnI0f4PC2N+e3bc4sDTlaqTbSGP/8ZvvnGWBo1cnVFlWPRFpb9vIypG6bSpUkXXhn0\nCtc0vMZ1BWVmXg720kPx588bE9RER1+eda5168t3rAuo3iGAPIuFIzZhbjsMX6B1uUPw1wQG1qor\nKITnk/CuI4aPMrMuuj5ZcwKMX5ZXmPY0NS+P+3/9lWBvbz649loayLHDcmkNf/ubMQX8mjVGxngK\nc4GZOTvm8NrW1xjWfhjP932eVvVbubqsknJySgb70aOX7xWfnGxcr14c5m3alHzeuHG1huPTCwpK\nnDBX3HM/kptLhPWkOdsh+HYmEy3k2nXhAhLedcSw0Wl4e7fmk6QbjdOlK3Dw4kVu37uXRxo35oVW\nreRM3kr4z39gxgxjIpfYWFdXUzXnzOd4Y9sbvLP7HUZ0GMG0PtNoEdrC1WVdXVGRcae34jAvXooD\nPivLuLNbcZjbhnvr1lDFKySKtCYpN7dMT/3gxYukFxYSbTJd6q3b9tzlD19RUyS864g7xyQxMOV6\nJvf7E8ycWe46uzIzifvlF15q3ZpxTZo4uULPtnq1MQXtv/8N993n6mqqLu1iGrO3zmbej/MY3XE0\nz/V6zjNCvCLZ2UaQ2/bWi8P96FHjUoHSvfbiW7q2aGHcLKayu7I9aa7UpW4mL68SPfXi521NJrkX\nu7CLhHcdcfuDvzJlV3dueW0pxJWduOOHrCzu2LeP/7VrR5wD7yRWl/zwA9x7r3FZ3osvgiceIj2T\nc4ZXt7zK/J/mE9cujud6PUeHiA6uLsuxtDZu4WrbU09MvHyv9pQU47p123uzFz8vXirRc9dacyo/\nv8wJc4fMZpJzc2kREECM9UYvtnd0ay13cxOVIOFdRwx8aDfvf9mTpj8eMHoaNn7LyeGWvXt5NyaG\nYRERLqqwdjhzBh54AHx9YelS95/MpSLnzOd4Z9c7/Gvnv+jRvAdTek2hZ4ueri7LOYqKjHvaF4d5\n8VJ8r/bjx8FkKj/Ui5ewsCsec8+3njSXYHOjl+IbwBzPzSXSz69EqNve1c2d78YnnEfCu47o99B3\nrF82EN8cc4khwXMFBXT94QemR0XxiAyVO0RhoTGZy8cfGzOy9e3r6oqqz1xgZsGeBczeOpvGQY2Z\n2H0i93W4Dz9vP1eX5jrFt3GtKNyTkow/AIqDvGVL4xa7LVqUfDSZym2+0GLheF5emVAvfl3P2/vS\nbVltQz3aZCJC5oevMyS864i77v+A99c8SoPs/EvvWbTm7p9/5prAQN70pFOlPcTq1TB+vHFV3syZ\nxuyinqrQUsiXB79k7q65/Hr2V8bfNJ4nuzxJsxA3nujdlS5cuBzkycnGcuLE5ccTJ4yh9/JCvfix\nnIDXWpOan18m1Itv1ZqvNW0DAmhjHX5vFRBA64AAWptMtAoIkHkaahEJ7zpi7N2zmfbj37km5cKl\n9149fpxV6el826mTzDhVQ86cgcceMw6jLl4MHTu6uiL7/Xr2V97Z9Q5Lf17KwDYDebLLk9zS+hbH\nztpW2xXPI1861G0fU1KMgK8o3Js1g6ZNSxx/P1dQwBGzmaO5uRzNzeWYzeOx3FyCvb0vB3pxuFuD\nPcrfX461exAJ7zrij4OmMDJlMTf/dgqAQxcv0vPHH9ndpQutKhi+E46hNcybB9OmweOPw/TpEBjo\n6qrsl5mXyft73+d/P/2PDHMGYzuN5ZHOj1TuJiji6oqH5ysK95MnjcXb25hmtmnTkkup9ywmE6fz\n80sE+qVHs5nkvDzCfX0vhbttsLcKCKC5vz/+8ke+25DwriP+PGAid2WspN/e4wAM2rOHu8PDmdS8\nuYsrqztOnYI//Ql27YK334bBg11dkePsSd3Dgp8WsPSXpXSI6MC4zuO499p7CfEPcXVptZvWxix1\np05dDvPSS/Fn/v7lhnrxe0VNmnCyYUOOQplgT8rL42ReHg18fGgZEEALf//LjzbPI/38ZG4IJ5Hw\nriOm9h3PoIsbuGX3ETacO8eThw7xW7du+Mhf0k731Vfw9NPGTbhefRWuv97VFTlOflE+qw6tYuGe\nhWxM2sgtrW9h5HUjueuauwjyqxu3jnVLWhvTzZYX6qXfCwgwZqiLjLy8NG5MUWQkpxs35nh4OMmh\noRwPCCC5qIjjeXkk5+VxPDeX84WFNCsV6C0DAmjp73/puZwt7xgS3nXE9D4P07dgB4O2HaDXTz/x\ndLNmjI6MdHVZdVZ+Prz7Lvzf/8HddxsntNW2QZBz5nN8fuBzVvy6gq3JW7mt7W3c3+F+BkcPJtjf\nze+nWlcVh3xqKpw+fXkp7/WZM0bQW8OdyEhymzThRIsWHG/ShOSGDTkeHGyEvLc3xwsKOJ6bi7dS\nNPP3NxY/P5panzf187v0fqSvr3QsrkLCu46Y2et+eqjfCFm9mXEHDvBr9+4yH7MbOH8eXn7ZOCY+\nYgRMmVLmMvxaIf1iOp8d+IyPfv2Irclb6dmiJ3dfczd3X3M3UfWjXF2eqI7ioC8v3Mt7z88PHR7O\n+ZYtSWnVipSmTTkZGUlKw4akhIZy0mQixc+PFKVIs1iIsIZ501KPts/r+/jU2UvjJLzriP+LHUpn\nUzIfv7uU6wIDebZlS1eXJGykpcFbbxnTq955Jzz3HFx3naurqhmZeZl8feRrvjz0JWsS1tAkqAl3\nX3M3g6MH06N5j7p9DXltpbUx3/zZs8aSlnbF54Xp6aQGBnKybduyQR8SwsnAQFJ8fSnw8qKplxeR\nfn40DgykcUAAjf38aOznZ7xn87y2nWwn4V1HvNp9MG3Dsnl02ssk9OhBI0++6LgWO3/eOJnt7bfh\nmmvgqadg2DBjxrbaqMhSxPYT2/ny0Jd8k/gNCekJ9Inqw6DWgxjYZiDXN7q+zvas6rzcXEhPv2LI\nZ1+4wMnCQk5rTaq3N6kREZxu0oTURo1IbdiQ1Pr1OV2vHqdNJupZLDS2WGisFI19fYkMCKBxvXo0\nDg2lcb16l8I+wkOG7CW864g3ut5CdteWrH3iD2y76SZXlyOuIj/fuEvZ228bd8V89FF46CEj0Guz\n9IvpfHfsO75N/Jb1R9eTmZdJv6h+9GrRi14te9G5cWd8vOSEJ1EOreHiRSPwMzJKLJaMDM5lZnI6\nN5fUggJStSbVy4vTPj6kmkykhoWRGh5OaoMGpAcFEZKfT0R+PhGFhUQAEV5eRPj4EO7vT4TJRERQ\nEBGhoUQ0aEBEcLBLro+X8K4j5tzUi/hRQ+g4fDQvtm7t6nJEFezbBwsXwrJlxkybDz4I999vnCdU\n29XloFUAAAzXSURBVCWdT2JT0ia2JG9h8/HNHL9wnG7Nuhlh3qIX3Zp1I8wU5uoyhSfTGszmS0Ff\nlJ7OuXPnOJuZydmcHNLy8zlbWMhZrTnr5cVZHx/O+vtz1mTibFAQaSEh+BUVEZGdTYTZTEReHuFF\nRURYLER4eRHu40OYvz8NTSbC6tUjLDiYsPr18W/QAOrXNyZ9qMbokoR3HfFu5+68/tc/seD2O+hT\nv76ryxHVUFgI69fD++8bU6927AhDhxrD6jExrq7OOc6Zz7HtxDa2HN/CluQt/HjqRxoGNqRLky7c\n1OQmujTpQpemXQgPlDvjCefQZjNZ585xtjjws7M5e/EiZ62hn2axcA7I8PIiw9eXDD8/0gMD8S0s\nJCwri7DMTMLMZsJycwkrKKBhYSFhFgthQJi3N2HW8A8LCCCsXj1MwcGokBDUoEES3nXBfzvdyMTX\nXyW9/wCC5DpLj5eXBxs2wMqVxlK/Ptx6KwwcCP36Ga/rAou2cDjjMD+c/IEfThnLj6d+pH5AfW6I\nvIHrIq7juojr6NioI+3D22PyldkEhetprckpKiKjsJCMixfJyMwkIzubjJwc0nNzycjLI6OggAyL\nhQytyVCKDB8f0q0nv4SZzZyKi5Pwrgtm9+3LK1OeI+3Ou1xdinAwi8W4l/i33xrL9u3QoQMMGAA3\n3ww9ehiX4dYVFm3hSMYRfj7zM/vP7Gf/WWM5nHGY5iHNLwV6TMMYosOiiQ6LJrJepJwYJ9ye1hqz\nxUJGQQEtTCYJ77rgL/cOY/WIkfw8apSrSxE1LC8Ptm2D776DHTtg504ICTFCvEcPuPFGY1a38Do2\nslxQVMDhjMNGmJ/Zz5FzRziccZjDGYe5WHDxUpBHh0XTtkFbWoa2pGVoS1qEtpDZ4YTbcdkxb6XU\nYOAtwAt4T2s9q5x15gB3ADnAI1rrPRW0JeF9FU+NfZCfBwxk0yPjXF2KcDKtISHBCPIdO2DvXvj5\nZ2NyrOuvN46dd+xoHDePjjamvq5rndALuRdKhPnhjMMcv3Cc5Mxkjl84jsnHRIvQFkaYhxiPzUOa\n0ySoCZFBkUTWi6RhYEO5s5pwGpeEt1LKCzgEDAROAruAB7TWB2zWuQN4Wmt9p1KqB/BPrXVsBe1J\neF/F2KeeJKNzN74c/7irSxFuQGvjjpM//2wsv/wChw8bl6VlZRmzvLVta4R5y5aX70DZrJkR7nXp\ntAmtNenmdCPMLyRfCvXkzGRSs1M5nX2a1OxUsvKziAiMoHFQYyKDIo3HepE0qteIhqaGhJnCSiwN\nTA1kQhpRbdUJb0f8b9sdSNBaJ1mLWA4MBQ7YrDMUWAygtd6hlApVSkVqrU87YP91TnZQECGFFleX\nIdyEUkYgN28Od9xR8rOsLCPEjxwxAj0hAeLjjTtRpqQY82OEh18O8vDwsktEhPHYoIExZO/JcwIp\npQgPDCc8MJybmlQ8R0JeYR5ncs5wOuf0pUA/nXOapPNJ7EndQ4Y549KSbk7nnPkcJl9TiWAPDQgl\nxD+EYL/gko/+l1+Xfh7oGyg9flEpjgjvZkCyzesTGIF+pXVSrO9JeFdDVnAwjSwyOiGuLjgYOnc2\nlvIUFhrTVp84YUxZnZZ2efKr3367/Dotzbh0NjPT6KmHhFS8BAYai8lUucXf35h1rrzFVZNj+fv4\n0yK0BS1CW1Rq/f9v7+5i5KrLOI5/fzvbhSJaXxAUChRCfAkJEkwUJSbEl0iqIlcE1ChyYyII0cSo\naCKXamKUqDcoEjQQI9VEYkSBYGO8QDRCqFAQKfKqNUJpV2iX7e7jxZntvnR2Z9ttmTns95OczP+c\n85/ZZ0+7+9tz9syz0zXN+MT4vEDfuWcn4y+Os2tiF+MT4zyz+xkefe7Redt2Texi/MXxfeM9e/cw\n1hlj7Zq1rB1dy1Frjlp8vGDbWGds3rKms2a/bf2W0ZFRRkdG6aRDZ6Qzb+wPFcNlKC+YXXD11YMu\nYag9ctqpvH3HnkGXoZeB0dHZs/blqGo6Xe7a1XvZuROef77pk7FjR/OXKXfvXnqZmIDJyd5Lp7N4\nsI+NzQZ8p9M8HsrxzA8OM/cMJIuPYYRkHck64JQe+5vHIwJHAK9f5LVIMcUEe/MCUyO72cvuZpzd\n7O0uk7zAruzm2e64mbObacaZ4sWey3R6b99/3yTFFNPs7T52x5mCCiN0CB1G6DDC6L5xGJ2zb/54\n35zqEEZIQhgBmseQZlzNuOc++u9jzpz+zwOYvUq9cFvP9fSft+Tzl3zegTsU4f0UMPcvY6zvbls4\n58Q+c/Z5cPPmfeNjNmzgmA0bVlrjy8oZD2/jogs/MegytAols2fMh7srXFVzZWCxYJ+cbFrPTk/P\nLlNTh25c1SwztSw2Xun++XND1ZFUHbms1zocFnvt6ZreP9S745nAnxlPLTavimKaonmE2fW544Pa\nt8Rr177xzK8bZz/J2jeuJdeXO285688++gg7/rltwfYDcyhuWOsAD9HcsPYv4G7g4qraOmfORuCy\n7g1rZwPf9YY1SZIGdMNaVU0luRy4jdm3im1N8plmd11bVb9JsjHJP2jeKuZ7nCRJOkg2aZEkaYAO\n5szb2wclSWoZw1uSpJYxvCVJahnDW5KkljG8JUlqGcNbkqSWMbwlSWoZw1uSpJYxvCVJahnDW5Kk\nljG8JUlqGcNbkqSWMbwlSWoZw1uSpJYxvCVJahnDW5KkljG8JUlqGcNbkqSWMbwlSWoZw1uSpJYx\nvCVJahnDW5KkljG8JUlqGcNbkqSWMbwlSWoZw1uSpJYxvCVJahnDW5KkljG8JUlqGcNbkqSWMbwl\nSWoZw1uSpJYxvCVJahnDW5KkljG8JUlqGcNbkqSWMbwlSWoZw1uSpJYxvCVJahnDW5KkljG8JUlq\nGcNbkqSWMbwlSWoZw1uSpJYxvCVJahnDW5KkljG8JUlqmRWFd5LXJLktyUNJfpdkXY8565PcmeT+\nJFuSXLGSjylJ0mq30jPvLwN3VNWbgTuBr/SYsxf4QlWdDrwLuCzJW1b4cVe9zZs3D7qEVvA4LZ/H\nank8TsvnsTp8VhreHwVu6I5vAC5YOKGq/l1V93bH/wO2Aies8OOuen5RLI/Hafk8VsvjcVo+j9Xh\ns9LwPraqtkMT0sCxS01OsgE4E/jTCj+uJEmr1mi/CUluB46buwko4Gs9ptcSr3M0sAm4snsGLkmS\nDkKqFs3b/k9OtgLnVtX2JG8Afl9Vb+0xbxT4NXBrVV3T5zUPviBJklqoqnIg8/ueefdxC3AJ8E3g\nU8CvFpn3Y+CBfsENB/4JSJK02qz0zPu1wM+BE4HHgAur6rkkbwR+WFUfTnIO8AdgC81l9QKuqqrf\nrrh6SZJWoRWFtyRJeukNXYe1JN9KsjXJvUl+keRVg65pmCQ5L8mDSf6e5EuDrmdY2RzowCQZSfLX\nJLcMupZhlmRdkpu736PuT/LOQdc0jJJ8PsnfktyX5MYkY4OuaVgkuS7J9iT3zdnWt+HZQkMX3sBt\nwOlVdSbwML0bv6xKSUaA7wMfBE4HLrbhzaJsDnRgrgQeGHQRLXAN8Jvujblvo+lboTmSHA98Djir\nqs6gubfqosFWNVSup/kePtdyGp7NM3ThXVV3VNV0d/UuYP0g6xky7wAerqrHqmoS+BlNoxwtYHOg\n5UuyHtgI/GjQtQyz7lXA91TV9QBVtbeqdg24rGHVAV7RfafRUcDTA65naFTVH4EdCzb3bXi20NCF\n9wKXArcOuoghcgLwxJz1JzGQ+rI5UF/fAb7IEn0aBMApwH+TXN/9FcO1SdYOuqhhU1VPA98GHgee\nAp6rqjsGW9XQO6CGZzCg8E5ye/d3ITPLlu7jR+bM+SowWVU3DaJGvTzYHGhpST4EbO9epUh3UW+j\nwFnAD6rqLOAFmsudmiPJq2nOJE8GjgeOTvKxwVbVOn1/kF7p+7wPSlV9YKn9SS6huYz33pekoPZ4\nCjhpzvr67jb10L1ktwn4aVUt1oNgtTsHOD/JRmAt8MokP6mqTw64rmH0JPBEVf2lu74J8KbR/b0f\n2FZVzwIk+SXwbsATscVtT3LcnIZn/+n3hKG7bJ7kPJpLeOdX1cSg6xkyfwZOS3Jy9+7Ni2ga5ai3\nZTcHWq2q6qqqOqmqTqX5/3Snwd1b97LmE0ne1N30PrzJr5fHgbOTHJkkNMfJG/vmW3iVa6bhGSzd\n8GyfgZx59/E9YAy4vfl3566q+uxgSxoOVTWV5HKaO/JHgOuqyi+KHrrNgT4ObElyDzYH0qFxBXBj\nkjXANuDTA65n6FTV3Uk2AfcAk93Hawdb1fBIchNwLvC6JI8DXwe+Adyc5FK6Dc/6vo5NWiRJapeh\nu2wuSZKWZnhLktQyhrckSS1jeEuS1DKGtyRJLWN4S5LUMoa3JEktY3hLktQy/wcdwMndqBK/HwAA\nAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f3d1811eba8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(8,5))\n",
    "x = np.linspace(-2,10,1000)\n",
    "labels = []\n",
    "for k,ell in [(1.2,1),(1.4,2),(1.2,5),(1.4,5)]:\n",
    "    y = f(x,k,ell)\n",
    "    plt.plot(x,y)\n",
    "    labels.append('$k = {0}$, $\\lambda = {1}$'.format(k,ell))\n",
    "\n",
    "plt.title('Weibull Distribution')\n",
    "plt.ylim([-0.2,0.8])\n",
    "plt.legend(labels)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we are ready to fit the data. According to [industry standards](http://wind-data.ch/tools/weibull.php?lng=en) we know the parameter $k$ should be between 1 and 2 and $\\lambda$ should be around the mean wind speed (~15 km/h in our case) therefore we can specify an intial guess for the curve fitting algorithm:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "params, _ = curve_fit(f,X,Y,p0=(1.5,15))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([  1.85706203,  15.02485969])"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "params"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.text.Text at 0x7f3d184fc048>"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAEKCAYAAAD+XoUoAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XecVPW9//HXhyrIgi02VFQsyGIEC67i4qokAibBeG3Y\nS5R7o9GbYtT4M5J69Sa5lmhivJaIiQGxhXjt0RVEg1gA2aUJiDSx66qolM/vj+9ZGYbZ3ZndmT1n\nZt7Px2MeO3PO95zz2Vn4zHe+7Zi7IyIi5aFD3AGIiEj7UdIXESkjSvoiImVESV9EpIwo6YuIlBEl\nfRGRMqKkX8LMrMHMdm3lsU+b2Tn5jSi/zOxMM5vSxnPMNrOheYrnFDN7NOX1ejPbPR/njs7X6r9n\nFuf+tZldFD0/3MyWFug6V5nZXc3sn2Zm+xTi2hIo6RcJM7vMzB5O27bAzP4vbdt8MzsRwN0r3P31\nAsTSy8xuM7OVZvahmc01sx/n+zpZyjjRxMz6REn3o+ix0swmmdmwjQ52H+Duk5u7QMq5mv3/4u53\nu/vwlmLLRqYP3QL+PbcBTgf+lHq5PF/j0ZT3vrlz/wb4RT6vLRtT0i8ek4FDzMwAzGx7oBMwKG1b\n36hsIV0LbA7s7e69gG8BrxX4mq3hQC937wnsBzwJPGBmZ+R4HovOZU0WMOvYxHHF4CzgYXf/vBAn\nN7PuwAHAM1kU/wdwhJltW4hYREm/mEwHugADo9fVwNPAvLRtC939Tdi4ecHM7jCzG83soajm+7yZ\n7dZ4cjP7mpnNMbP3zez3NJ+wDgLudvePANx9vrvfn3Ku9Wb2PTNbaGZvmdl/px5sZueYWb2ZvWtm\nj5jZLin7+pnZ49G+OWZ2Qsq+raLa+odm9i/CB1xLLIrxLXe/ARgLfBmPmS02syOj5weZ2fTo/CvN\n7LdRscZk9UH03h0cNS09a2b/Y2bvAFc10dx0TKb3Ib2ZI/XbhJn9kvC3vDG63g0p72vj37OnmY2L\nzrvYzK5IOdeZZjbFzH5jZu9F10/9BpJuBM0kZDO7KGoG27Gx6cfMLjGzVWa23MxGmdkIM5tnZu+Y\n2eVppzgKmOrua6LXXc3szuh3e9XM9m8sGH3wvAQc3Uy80gZK+kUi+g8zDWhsfx5KqNE/m2Hbl4el\nneYk4CpgC2Ah8CsAM9sauA/4CbBNtG9IM+H8C/i1mZ1lZns0UeZYYP/oMaqxqcLMRgGXRfu/AkwB\n/hbt6w48DvwliuNk4A9m1i865x+AT4HtgHOB1vQ53A9sa2Z7Z9h3PXBd9O2lL3BPtL3x/e3p7j3d\nfVr0+mDCN5xtid5LNn3PM74PTZR1AHf/f4T35cLoehdlKH8jUAHsCtQAZ5jZ2Sn7BwNzgK0JTSa3\nZfh9G+1LqDxswsx+CpwBDHX3FdHm7QkVkB0J/57+FzgVGER4r640sz4ppxkJpDZDfhO4G+hFqNnf\nlHbZOYRvZlIASvrF5Rk2JKBqQmJ4Nm1bao0tvbb+gLu/5O7rgb+y4RvCSGC2uz/g7uvc/TrgzWbi\nuJCQmC8A6iz0I6TXJK929w/dfRlwHTA62j4G+K/o28F64GpgoJntDHwDWOzu4zyYSfgwOiFqTz8O\nuNLdP3P3OuDOZmJsSmPi2irDvi+APcxsa3f/1N1fSNuf/n4ud/c/uPv6ZppGmnofWqOxGa8D4QP8\nsijOJcDvCO3yjZa4++0eFte6E9i+mSaTLYCGtG0dzOx3wDCgxt3fS9n3BfBrd18HjCd8QF8XxVIP\n1LNx0h4JpPZHPevuj0Wx3QV8Ne3aDVFMUgBK+sVlMnCYmW0JbOPuC4HngEOjbQNovj0/NZF/CvSI\nnu8IpI/WaHL0hrt/7u5Xu/tBhJrkRGCimaX+R12W8nxJdA2APsD1UbPDe8C7hBps72hfVeM+M3sf\nOIVQs/8KoQ8j/by56h39fDfDvnOBvYG5FkaRHNPCubIZ4dLU+9AW2xDeizfSzt075fWXf2t3X034\nwOhBZu8TvjWk2gI4j/AB/XHavnd9w0qNq6Ofb6XsX914LTMbAHzg7sszxUb4d7iZbdxJXgF80ESs\n0kZK+sXleTb8Z5wK4O4NhNrreYSaZ2sS4Upgl7RtO2dzYJQQfk3o2N0tZVfq8X3YUMNeCoxx962i\nx5bu3sPd/xXtq03b19PdLwTeBtaknTc95mwcB6xy9/kZfpeF7n6Ku3+F0O5/r5l1o+nRJtmMcEmP\nt/F9+ATonrJvhxzO/Q7hvUhtQukDLM9cvEWzgL3Str1H+Ob1ZzM7tJXnhU1r+dnYB5jZhmtKM5T0\ni4i7fwa8CPyA0LTTaGq0rbWjdv4P6G9mx5pZRzO7mFC7zsjM/p+ZHWhmnc2sK/CfhNpiarvwJWa2\nRdRscxGhGQDgZuAnZtY/OlcvMzs+2vcQsJeZnWZmnaLzH2hme0dNQfcDY82sW3T8mS38XsaGJpFt\nzexC4EpCn0Km3+tUC8MXAT4kJN71hA+c9WTXcZwu9X24mA3vwwxgqJntbGa9MsS0Csg4xj96L+4B\nfmVmPaL28+8Tmkpa42FCv0D6dSYT2urvM7ODWnnu9Pb8TL5sNov+PR0APNHK60kLlPSLzzOEpo5n\nU7ZNibalj8DIaqy1u78LnABcQ6hF9iX6JtHUIcAdhGS4nDA64xh3/zSlzN8JozBeJnTW3R5d60FC\nO/54M/uAUMscHu37GPg6oQN3RfS4GuganfN7hK/+K6Pz3d7Srwa8b2YNKdc53t3vTCvTaDihj+Ij\nwrDUk6KmrNWEjtqpUbPT4Baum3rupt6HJ4EJUVzTo32prif0ZbxrZtdliPUiQtPIIsKH/V/c/Y4W\nYmnKOGBElHA3PijEeS4wycwGbnJk5nM7hA90Qq39uWaunX78t4CnG0egSf5ZNjdRiTrpriN8SNzm\n7tdkKHMDYejXJ8BZ7j4j2t4LuJXQ3rweOCdl9IOUIDNbD+zh7ovijkWyEw0TbRzWmq9zngD8m7uf\nnMMxzwPnRh3CUgAtJv2og2U+oTa3glArOdnd56aUGUEYXnaMmR0MXO/uVdG+PwPPuPsdZtYJ6N44\nvltKk5K+AFiYgdugSl6ydMqizGBgQWMHoZmNB0YBc1PKjCJ8RcTdp0XttNsRevGr3f2saN9aQAm/\n9OkenNLYNCQJk02bfm82Hpq2jI2HhmUqszzathvwjoXZoC+b2S3RaAgpYe7eUbV8kWQqdEduJ8JM\nxJvcfX9Cx1PGkRMiIlJ42TTvLGfj8dA7sel44OVsPB45tcxSd38xen4vcGmmi5iZmgRERHLk7jkt\n7JdNTX86YWp6HzPrQhhONymtzCTC+hyYWRVhBt4qd18FLDWzxokfRxGmaDcVfKIfV111VewxKE7F\nqTgVZ+OjNVqs6bv7umhSy+NsGLI5x8zGhN1+i7s/bGYjzew1wpDN1IWfLgL+amadCWOKz06/hoiI\ntI9smndw90cJa5KkbvtT2usLmzh2JmEpXhERiZlm5OagpqYm7hCyojjzS3Hml+KMV1YzctuDmXlS\nYhERKQZmhhegI1dEREqEkr6ISBlR0hcRKSNK+iIiZURJX0SkjCjpi4iUESV9EZEyoqQvIlJGlPRF\nRMqIkr6ISBlR0hcRKSNK+iIiZURJX0SkjGS1nn45++wzeOwxmDkT3GGffWD4cOjZM+7IRERyp5p+\nE9zh9tuhTx+4/vqQ/NesgTvvhF13hZ/9DFavjjtKEZHcqKafwRdfwFlnwbx58OijMGjQxvsXL4Yf\n/hAOPRTuuw923z2WMEVEcqaafpp16+D44+HTT2Hq1E0TPsBuu4Vkf/bZUF0Nc+a0f5wiIq2hmn6a\nSy4JzTYPPwydOzddzgwuugi23BKGDYPnngtNQSIiSaakn+LBB+Hvf4cXX2w+4ac6/XR4+2049lh4\n9lnYfPPCxigi0ha6R27k3Xdh333hnnvgsMNyO9Y99AGsWwd/+UtBwhMR2URr7pGrpB8ZMwa6doUb\nbmjd8atXh/b/n/8cTjwxv7GJiGSipN9Ks2fDkUeG0TpbbtlEoYaGUHDAAKioyFjkhRfgm9+EGTNg\nhx0KF6+ICLQu6Wv0DnD55XDFFS0k/OpqGDo0/GxoyFhs8GA499wwnFNEJInKPunPmgUvvwz//u/N\nFJo9G+rqYO1aqK8Pz5twxRVhqOczz+Q/VhGRtir7pH/11fCf/xna85s0YABUVoYhPf37h+dN2Hxz\n+N3v4MILw2eEiEiSlHWb/qJFoUlm0aIs1tJpaAg1/MrKJtv0G7lDTU0Y0XP22fmKVkRkYwXryDWz\n4cB1hG8Gt7n7NRnK3ACMAD4Bznb3V6LtrwMfAuuBNe4+uIlrtHvS/8EPoEuXUNvPt6lT4ZRTYP78\nFr5FiIi0UkGSvpl1AOYDRwErgOnAye4+N6XMCOBCdz/GzA4Grnf3qmjfIuAAd3+/heu0a9L/7DPY\neWeYNq1wa+d84xvw9a+HmbsiIvlWqNE7g4EF7r7E3dcA44FRaWVGAeMA3H0a0MvMtmuMK8vrtKsH\nHoCBAwu7WNqvfgX/9V/wySeFu4aISC6ySca9gaUpr5dF25orszyljANPmNl0MzuvtYHm2y23wHkF\njma//eCQQ8ISzSIiSdAeNfAh7r4/MBK4wMxyXOQg/xYuDH2yo9K/rxTApZfCb38b1uIXEYlbNguu\nLQd2SXm9U7QtvczOmcq4+8ro59tm9gChuejZTBcaO3bsl89ramqoqanJIrzcTZgQlk9ujw7Wgw8O\nTUgTJsBppxX+eiJSumpra6mtrW3TObLpyO0IzCN05K4EXgBGu/uclDIjgQuijtwq4Dp3rzKz7kAH\nd//YzDYHHgd+5u6PZ7hOu3Xk7rcf/P73YYJte3jssTBLd9bUBjrUN7+Ug4hItlrTkdtiTd/d15nZ\nhYSE3Thkc46ZjQm7/RZ3f9jMRprZa0RDNqPDtwMeMDOPrvXXTAm/Pc2ZA++8k/tKmm3x9a9Drw4N\nNAyqptfSaKz/lClK/CLS7spncla0YNqvJw3g7c8quPbawl0qk0d++jzDfjmUzr42zOydPBmqqto3\nCBEpKVplsymNC6bV1TG3YyUfPjSFg4e1by3783caWLB9Nf2tng6V/VXTF5E20yqbTUlZMG23z+s5\nqHvTC6YVStdtKrj/+1O47tuTlfBFJDZlVdNfN7ue5T37s8uSeJLusmXw1a/CkiXK+SLSdqrpN6Wi\nAqZM4fv7T2bG7+OrZe+0Exx1FIwbF8vlRUTKpKZPqOzvuCOsWBFvLbu2Fi64ILQ4WU6fzyIiG1NN\nvxlPPBGWRIi7WeXww8Ps3GnT4o1DRMpT2ST9hx4K96+Nm1m4peKtt8YdiYiUo7Jo3nGH3r3D0Pg9\n9ijIJXLy5puwzz7wxhvxf/MQkeKl5p0mzJ0b5kP17Rt3JMH224c7a02YEHckIlJuyiLpP/VUGDWT\npI7T73xHTTwi0v7KJukfeWTcUWzs6KPDuP1XX407EhEpJyWf9NetC8Mkk5b0O3UKN07/85/jjkRE\nyknJJ/2ZM2G77cIY/aQ5/XS4+25YuzbuSESkXJR80v/nP5NXy2+0997h5uz//GfckYhIuSj5pN/Y\niZtUp58Od90VdxQiUi5Kepz+unWw1VbhnrjbbJPXU+fNW2/BXnuFTt0ePeKORkSKicbpp5k9G3bY\nIbkJH2DbbcNdvO6/P+5IRKQclHTSnzoVhgyJO4qWqYlHRNqLkn4CfOtb8OKLsHx53JGISKkr6aT/\n7LPFkfS7dYPjjgvDN0VECqlkk/6yZfDpp6GTtBicfjr85S9xRyEipa5kk/7UqXDooclab6c51dXw\nzjswZ07ckYhIKSvppF8MTTuNOnaEE07QypsiUlglm/Sfe664kj7AySeHpJ+QqRMiUoJKMul/9hnU\n18OgQXFHkpuDD4bVq2HWrLgjEZFSVZJJf9as0IHbvXvckeTGDE46CcaPjzsSESlVJZn0X3wRDjww\n7ihaR008IlJIJZv0Dzoo7ihaZ+DAsNb+9OlxRyIipagkk/706cVb0zfbUNsXEcm3rJK+mQ03s7lm\nNt/MLm2izA1mtsDMZpjZwLR9HczsZTOblI+gm/PJJ7BoEey7b6GvVDiNSX/9+rgjEZFS02LSN7MO\nwI3A0UAlMNrM+qWVGQH0dfc9gTHAzWmnuRioz0vELXjlFRgwALp0aY+rFUb//rDllmGugYhIPmVT\n0x8MLHD3Je6+BhgPjEorMwoYB+Du04BeZrYdgJntBIwEbs1b1M0o5k7cVGriEZFCyCbp9waWprxe\nFm1rrszylDLXApcA7TIeZfr04u3ETXXiiXDvveFGMCIi+dKpkCc3s2OAVe4+w8xqgGZXwhk7duyX\nz2tqaqipqcn5mi++CJdfnvNhibPnnrD99qGJZ+jQuKMRkSSora2ltra2Tedo8XaJZlYFjHX34dHr\nywB392tSytwMPO3uE6LXc4HDCW35pwFrgW5ABXC/u5+R4Tptvl3ihx/CTjvBBx+EtWyK3S9/GW6n\neMMNcUciIklUqNslTgf2MLM+ZtYFOBlIH4UzCTgjCqIK+MDdV7n7T9x9F3ffPTruqUwJP19eeQW+\n+tXSSPgAxx8P992nUTwikj8tJn13XwdcCDwO1AHj3X2OmY0xs/OjMg8Di83sNeBPwHcLGHOTZs6E\n/faL48qF0a9fuLH788/HHYmIlIoWm3faSz6ad845JyxaNmZMnoJKgJ/9LDRXXXtt3JGISNIUqnmn\naJRaTR9CE8+996qJR0Tyo2SS/tq14a5TxTwTN5P+/aFHD63FIyL5UTJJf968MHJn883jjiS/zDbU\n9kVE2qpkkn4pNu00OuEEmDhRyy2LSNsp6ReBffcNawm99FLckYhIsSuZpD9rVukmfTXxiEi+lEzS\nnzkzTMwqVY1JX008ItIWJZH033473FB8l13ijqRwBg0KwzZnzow7EhEpZiWR9Btr+ZbTFIXiYhY6\ndNXEIyJtUTJJv1Tb81Mdf7xG8YhI25RE0p89u/QmZWVy4IHw2Wfh9xURaY2SSPr19WHmaqnTKB4R\naauiT/ru5ZP0AU4c0cCCcc9DQ0PcoYhIESr6pL90KVRUhBuJl7yGBgb/sJo/vz6Uzw6qVuIXkZwV\nfdKvq4PKyrijaCezZ2P1dXRhLZ0W1IdfXkQkB0Wf9MupaYcBA6CyknWdOvNal/5l9GknIvlS9Em/\nrGr6FRUwZQr2zGS+2WsK81ZUxB2RiBSZok/6ZVXTB6iooMOhVQw/oYKJE+MORkSKTVEn/XIbuZOq\ncbllEZFcFHXSX7Ys3DRlq63ijqT9DRkCb70F8+fHHYmIFJOiTvpl1Z6fpmNH+Ld/U21fRHJT1Em/\nXJt2GqmJR0RyVdRJv5xr+gCHHQarVsGCBXFHIiLFoqiTfrnX9Dt2hOOOU21fRLJXtEnfHebMgX32\niTuSeKmJR0RyUbRJ/+23Q013m23ijiRe1dWwciW89lrckYhIMSjapD9vHuy1V9xRxK+xiUfLLYtI\nNoo66e+9d9xRJIOaeEQkW1klfTMbbmZzzWy+mV3aRJkbzGyBmc0ws4HRtq5mNs3MXjGzV83sqnwF\nrqS/wdChYaLaokVxRyIiSddi0jezDsCNwNFAJTDazPqllRkB9HX3PYExwM0A7v45cIS7DwIGAiPM\nbHA+AlfS36BjR/j2t1XbF5GWZVPTHwwscPcl7r4GGA+MSiszChgH4O7TgF5mtl30+tOoTFegE5CX\n23rPn682/VRq4hGRbGST9HsDS1NeL4u2NVdmeWMZM+tgZq8AbwJPuPv01ocbrFkDr78Oe+zR1jOV\njsMPhzfeUBOPiDSvU6Ev4O7rgUFm1hN40Mz6u3t9prJjx4798nlNTQ01NTUZz7l4Mey4I2y2Wf7j\nLVadOoUmnnvvhR//OO5oRKQQamtrqa2tbdM5zL351hYzqwLGuvvw6PVlgLv7NSllbgaedvcJ0eu5\nwOHuvirtXFcCn7j7/2S4jrcUS6N//AP+8Ad45JGsipeNJ5+Eyy+H6W3+LiUixcDMcHfL5Zhsmnem\nA3uYWR8z6wKcDExKKzMJOCMKogr4wN1Xmdk2ZtYr2t4N+BowN5cAM1F7fmY1NaHZa/HiuCMRkaRq\nMem7+zrgQuBxoA4Y7+5zzGyMmZ0flXkYWGxmrwF/Ar4bHb4D8LSZzQCmAY9FZdtEI3cyS23iERHJ\npMXmnfaSS/PO4YfDT38KRx1V4KCK0BNPwBVXwAsvxB2JiBRaoZp3Ekc1/aYdcURo3lETj4hkUnRJ\n/8MPoaEhjN6RTXXqFNbimTAh7khEJImKLunPnw977gkdii7y9nPKKfC3v8UdhYgkUdGlzoULQ9KX\nplVXw7vvhjuLiYikKrqk/9pr0Ldv3FEkW4cOcPLJqu2LyKaKLukvXKjlF7IxenRI+gkZnCUiCVF0\nSV81/ezsv39YfVNDN0UkVdEl/YULlfSzYaYOXRHZVFFNzvrkk3BP3E8+0eidbMybF5ZmWLYs1PpF\npLSU/OSsRYtg112V8LO1995hPkMbF+UTkRJSVOlTnbi5a+zQFRGBIkv66sTN3UknwQMPwOefxx2J\niCRBUSV9deLmbuedYcAAePTRuCMRkSQouqSv5p3cqYlHRBoVVdJX807rnHBCqOl/9FHckYhI3Iom\n6a9ZA8uXh9E7kputtw5LLk+cGHckIhK3okn6S5aE4YddusQdSXE680y48864oxCRuBVN0lfTTtuM\nHAlz5oS5DiJSvoom6asTt226dAkrb951V9yRiEiciibpq6bfdmeeCePGaeVNkXJWNElfNf22O+AA\n2GwzePbZuCMRkbgUVdJXTb9tzNShK1LuimKVTXfYfHNYtQoqKto5sBKzYkWYobt8OXTrFnc0ItIW\nJbvK5ltvQffuSvj5sOOOMHgwPPhg3JGISByKIum//jrstlvcUZSOM85QE49IuSqapK+ZuPlz7LEw\nfTq88UbckYhIeyuKpL94sZJ+PnXvHhZhu/32uCMRkfZWFElfzTv5d955IemvWxd3JCLSnrJK+mY2\n3Mzmmtl8M7u0iTI3mNkCM5thZgOjbTuZ2VNmVmdmr5rZRa0JUs07+bfffrD99vDYY3FHIiLtqcWk\nb2YdgBuBo4FKYLSZ9UsrMwLo6+57AmOAm6Nda4EfuHslcAhwQfqx2VDzTmGcdx787//GHYWItKds\navqDgQXuvsTd1wDjgVFpZUYB4wDcfRrQy8y2c/c33X1GtP1jYA7QO5cA168PHY59+uRylDSroQGe\nf56Tj2mgthZWrow7IBFpL9kk/d7A0pTXy9g0caeXWZ5exsx2BQYC03IJcNUq6NkzTM6SPGhogOpq\nGDqUipHVnDaqgTvuiDsoEWkvndrjImbWA7gXuDiq8Wc0duzYL5/X1NRQU1Ojpp18mz0b6upg7Vqo\nr+c/LqrjG7+s4rLLoENRdOuLlK/a2lpqa2vbdI4Wl2EwsypgrLsPj15fBri7X5NS5mbgaXefEL2e\nCxzu7qvMrBPwEPCIu1/fzHUyLsNw993w97/DhAm5/3KSQWNNv74e+vfHJ09h0NAKfvtbGDYs7uBE\nJBeFWoZhOrCHmfUxsy7AycCktDKTgDOiIKqAD9x9VbTvdqC+uYTfHA3XzLOKCpgyBSZPhilTsJ4V\nnH8+3Hxzy4eKSPFrMem7+zrgQuBxoA4Y7+5zzGyMmZ0flXkYWGxmrwF/Av4DwMyGAKcCR5rZK2b2\nspkNzyXA11+HPbcPHY80NORyqDSlogKqqr5czOj00+Gpp2DZspjjEpGCS/wqm986ooG/LKmm59I6\nqKwMtVStvJZ33/se9OoFv/xl3JGISLZKcpXNrgtm02Ppho5H6uriDqkkXXBBGLP/+edxRyIihZTo\npL9+PTz99gB8n0ro3Bn69w+1fcm7fv3CLN177ok7EhEppEQn/RUroPNWFXScuqHjUU07hfO978GN\nN8YdhYgUUqKT/pcjd9I6HqUwRo4MN6x54YW4IxGRQkl80tfErPbTsWNo21dtX6R0tcuM3NbSbNz2\nd845sMceYT2eHXaIOxoRybdE1/SXLNFCa+1tq63glFPg97+POxIRKYREJ/2lS2GXXeKOovz84Adw\nyy2aCydSihKd9N94Q0k/DrvvDkceCbfdFnckIpJviZ2R6x4G66xYEZZWlvb1wgtwwgnw2mthioSI\nJE9Jzch9//0wmkQJPx6DB4dO9IkT445ERPIpsUlf7fnxu+QS+M1vwrcuESkNiU36b7wBO+8cdxTl\nbeTIsBbPk0/GHYmI5Etik75q+vHr0AEuvxx+/nPV9kVKRWKTvmr6yTB6NLz5JjzzTNyRiEg+JDbp\nq6afDJ06wRVXhNq+iBS/xCZ91fST49RTwzpIU6bEHYmItFVik75q+snRuTP85Ceq7YuUgkROzlq3\nDrp1C8sAdO0ac2ACwBdfwF57wV//CkOGxB2NiEAJTc56882w8JcSfnJ06QI//Wmo8SekniAirZDI\npK+mnWQ644xwk5VHHok7EhFprUQmfXXiJlOnTvDrX4ex++vXxx2NiLRGIpO+avrJdeyxob/lb3+L\nOxIRaY1EJn3V9JPLDK6+Gq68MnTuikhxSWTSV00/2WpqoF8/uOmmuCMRkVwlMumrpp98v/tdaN9/\n++24IxGRXCQy6S9dqqSfaA0N7PPB85x7YgNXXBF3MCKSi8RNzvriC+jRA1avDjdRkYRpaIDqaqir\nY22/SvZeNYWJj1aw//5xByZSfgo2OcvMhpvZXDObb2aXNlHmBjNbYGYzzGxQyvbbzGyVmc3K5lor\nV8K22yrhJ9bs2VBXB2vX0mlePf9zbh0XXaQJWyLFosWkb2YdgBuBo4FKYLSZ9UsrMwLo6+57AmOA\nP6bsviM6NisrVkDv3tmWlnY3YABUVoYFefr35xuXVrJ6dVieQUSSL5ua/mBggbsvcfc1wHhgVFqZ\nUcA4AHefBvQys+2i188C72cb0PLlSvqJVlERltucPBmmTKHjFhX88Y/wox/BO+/EHZyItCSbpN8b\nWJryelm0rbkyyzOUycqKFbDjjq05UtpNRQVUVYWfhJuojx4NP/xhzHGJSIsSN3pHNf3i9ItfhLtr\nPfFE3JEKHLwWAAALhUlEQVSISHM6ZVFmOZA6VWqnaFt6mZ1bKNOisWPH8sgj0Lcv1NbWUFNTk+sp\nJCY9esAf/wj//u/w6qvQvXvcEYmUntraWmpra9t0jhaHbJpZR2AecBSwEngBGO3uc1LKjAQucPdj\nzKwKuM7dq1L27wr8w933beY67u4ceWRYvnfYsDb8VhKbU0+FrbeGG26IOxKR0leQIZvuvg64EHgc\nqAPGu/scMxtjZudHZR4GFpvZa8CfgO+mBHU38Bywl5m9YWZnN3c9Ne8UtxtvhAcfhMceizsSEckk\nUZOz1q93KipCZ27PnnFHJK311FNh7f2ZM0OtX0QKo+jvnNXQEH5Gg0KkSB15JJx0Epx/viZtiSRN\nopJ+Y9OO5fS5JUn0q1/BggVw661xRyIiqbIZvdNu1J5fOjbbDO65JyzTs//+cMABcUckIpCwmr6W\nYCgt/frBH/4Axx8P770XdzQiAglL+suXazZuqTnhBPj2t+G003RfXZEkSFTSV02/NF1zDXz8cbjF\n4pcaGuD55zf03otIu0hU0lebfmnq3Bnuuw/Gj4fbb2fDmvxDh4afSvwi7SZxHblq3ilNX/kKPPxw\nyPP7fjybg6I1+amvD+vzV1W1fBIRabNE1fTVvFPa9t47jOg58ecDWL37hjX5qayMOzSRspGoGbmd\nOzsNDdC1a9zRSCGNHw8//X4DT15fxy4jKjUbT6SVWjMjN1HNOz16KOGXg5NPhk8/reCwH1VReyDs\nrpwv0m4SlfS33z7uCKS9nHMOfPYZHHVUWId/l11aPkZE2i5RSX+77eKOQNrTd78La9aEATyPPRYm\nc4lIYSUq6aumX34uvhh69YIjjoB//AMOPDDuiERKW6JG76imX57OOgtuvhlGjszxdoua4CWSMyV9\nSYRRo8IErjPOgOuvz2JJZk3wEmmVRCV9Ne+Ut+rqUHG//XY491z4/PNmCs+eHSZ1pU7wEpEWJSrp\nq6Yvu+4KU6fCRx/BoYfC/PlNFBwwIEzq0gQvkZwkKumrpi8Q5mtMnAjf+Q4MGQJ33JGhuaeiAqZM\ngcmTw09N8BLJSqJm5C5b5lqGQTZSVwejR8Puu8NNN2mZDpFURX+P3G23jTsCSZrKSpg+HQYODI+b\nboJ16+KOSqR4Jaqmn5RYJJnq62HMGPj0U/jNb8IN2EXKWWtq+kr6UlTcQ3v/5ZeHVTuvuQb23Tfu\nqETiUfTNOyItMYMTTwy1/qOPhmHDwu0Yp02LOzKR4qCkL0Wpa9ewhMPixWHRthNPDM09DzwQ1vPJ\nmmb1SplR846UhDVrwg1a/vjH8EFwzjlw9tlh1E+TGmf11tWFHmMN/ZQio+YdKVudO8Opp8Kzz8Lj\nj4fJXVVVcNBBodN3yZIMB2lWr5Qh1fSlZK1dC7W14RvA/feH+y8ffXR4HHYYbLYmqunX14dZva2p\n6Tc0hA+PAQP0LUHaXcFG75jZcOA6wjeD29z9mgxlbgBGAJ8AZ7n7jGyPjcop6UvBrF0bxvs/9lh4\nzJ4dlnEeOqiBYTvUsc/xlWyzWysSvpqHJEYFad4xsw7AjcDRQCUw2sz6pZUZAfR19z2BMcDN2R5b\nTGpra+MOISuKc1OdOsEhh8DYsaHfdulSuPRSoKKCXzxRRd+BFey6K3zzm2E46N13w6xZsHp1M3Hm\nq3koT53J+rvnV7HEmats2vQHAwvcfYm7rwHGA6PSyowCxgG4+zSgl5ltl+WxRaNY/hEozpZtsQUM\nHw4/+1noA3jvvfDz7LOhWzd48EE46STYckv4xjdqGTIETjsNrrwSbr0VHnoIXv5iAF/sWYm3ZdG3\nfC0R3dBA7V13te2Do60fPlker3+f8crmzlm9gaUpr5cRknlLZXpneewGDQ36eiyx6NgR9torPI47\nbsP29evhRz+CY48No4IWLw6rgK5aBW++WcHH709h2/V1LFtaSbfBFWy5ZbgT2BZbbPzo1Qs23zx8\noDQ+uneHrefPpm9dHR3WrsXr61k9vY4Oh1bRqVOIybL54t74wTFrFrz0Uuv7JtrSVJWvpq589JHk\n6xxLl7Y+JyXl98igULdLzKmN6UvV1WoXlUTp0AF69gwV8aFDM5WoYP36Kt59F956Cz78ED74YOPH\nu+/CwoVh+YjVqzc8Pv0U7OMB3NGhkt2pZ/66/nztm5V8sC4MQV2/PjRLde4cfqY+T9026LPZjFtS\nBzhfzKzn3P3rmN2jig4dwoeGGRmfp24b0DCba2fV0dnXsmZWPT8aVsfcLaqyeo/MYJ8PZvPfKcdf\n8rXMx5vBggWZJ9N1W9vAb6ZVs8vHdbzRo5IfV01hdafsckHjh2O3tQ389/PV7PxxHUt7VPLjQzY9\nR0sfpN3WNnDNc9Wsa5jFontf4rJDs4+j8firn6tm54Y6llZU5nx8vs7RlBY7cs2sChjr7sOj15cB\nntoha2Y3A0+7+4To9VzgcGC3lo5NOYd6cUVEcpRrR242Nf3pwB5m1gdYCZwMjE4rMwm4AJgQfUh8\n4O6rzOydLI5tVeAiIpK7FpO+u68zswuBx9kw7HKOmY0Ju/0Wd3/YzEaa2WuEIZtnN3dswX4bERFp\nVmImZ4mISOHFvgyDmQ03s7lmNt/MLo07nkZmdpuZrTKzWSnbtjSzx81snpk9Zma94owximknM3vK\nzOrM7FUzuyhpsZpZVzObZmavRDFelbQYU5lZBzN72cwmRa8TF6eZvW5mM6P39IUEx9nLzCaa2Zzo\n3+jBSYvTzPaK3seXo58fmtlFSYszivX7ZjbbzGaZ2V/NrEuuccaa9BM+eesOQlypLgOedPe9gaeA\ny9s9qk2tBX7g7pXAIcAF0XuYmFjd/XPgCHcfBAwERpjZ4CTFmOZioD7ldRLjXA/UuPsgd28cBp3E\nOK8HHnb3fYD9gLkkLE53nx+9j/sDBxCaqB8gYXGa2Y7A94D93f2rhOb50eQap7vH9gCqgEdSXl8G\nXBpnTGnx9QFmpbyeC2wXPd8emBt3jBlifhAYltRYge7Ai8BBSYwR2Al4AqgBJiX17w4sBrZO25ao\nOIGewMIM2xMVZ1psXwemJDFOYEdgCbBllPAnteb/etzNO01N6kqqbd19FYC7vwkk6q6+ZrYroSb9\nL8I/gsTEGjWZvAK8CTzh7tOTFmPkWuASILWzK4lxOvCEmU03s+9E25IW527AO2Z2R9R0couZdSd5\ncaY6Cbg7ep6oON19BfA74A1gOfChuz9JjnHGnfSLXWJ6wc2sB3AvcLG7f8ymscUaq7uv99C8sxMw\n2MwqM8QUa4xmdgywysNigc0NIU7C332Ih+aIkYQmvWoS9n4SaqP7AzdFsX5C+DaftDgBMLPOwLeA\nidGmRMVpZlsQlrHpQ6j1b25mp2aIq9k44076y4FdUl7vFG1LqlXRmkKY2fbAWzHHA4CZdSIk/Lvc\n/e/R5kTG6u4fAbXAcJIX4xDgW2a2CPgbcKSZ3QW8mbA4cfeV0c+3CU16g0ne+7kMWOruL0av7yN8\nCCQtzkYjgJfc/Z3oddLiHAYscvf33H0dod/hUHKMM+6k/+XELzPrQpi8NSnmmFIZG9f4JgFnRc/P\nBP6efkBMbgfq3f36lG2JidXMtmkcUWBm3YCvAXNIUIwA7v4Td9/F3Xcn/Ft8yt1PB/5BguI0s+7R\nNzvMbHNCO/SrJO/9XAUsNbO9ok1HAXUkLM4Uowkf9o2SFucbQJWZbWZmRng/68k1zgR0nAwH5gEL\ngMvijiclrruBFcDn0Zt9NqED5cko3seBLRIQ5xBgHTADeAV4OXpPt0pKrMC+UVwzgFnAFdH2xMSY\nIebD2dCRm6g4CW3ljX/vVxv/3yQtziim/QiVuxnA/UCvhMbZHXgbqEjZlsQ4ryJUmGYBdwKdc41T\nk7NERMpI3M07IiLSjpT0RUTKiJK+iEgZUdIXESkjSvoiImVESV9EpIwo6YuIlBElfRGRMvL/AXoO\nIRuWh62wAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f3d17eb7ba8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "x = np.linspace(0,80,1000)\n",
    "y = f(x,params[0],params[1])\n",
    "plt.plot(x,y)\n",
    "plt.plot(X,Y,'r.')\n",
    "plt.title('Wind Speed Distribution (km/h)')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**EXERCISE**\n",
    "\n",
    "Plot the squared error\n",
    "\n",
    "$$\n",
    "\\sum_{i=0}^{N} \\left( y_i - f(x_i,k,\\lambda) \\right)^2\n",
    "$$\n",
    "\n",
    "for a grid of values for $k$ and $\\lambda$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x7f3d17dd7470>"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEACAYAAAC+gnFaAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztvW/MdctVH/abc5733us/YKc2xQ4uThFxEkod6oCDAm1f\n6oZY1AIJ9YOJFFoSpaiq2ygiCIpU2URIoY3itCLJB4zj1i0uSSlJTD9Qx6KOZEVpndrGdmJjGsD/\nsWsg1Pb1vfd9zp5+OGeeM3vt9Xf27HP28zxnSa/es2fWWjN7P2ev9VtrZtZJOWdc6EIXutCF7i9t\nzj2BC13oQhe60Hnp4ggudKELXeie08URXOhCF7rQPaeLI7jQhS50oXtOF0dwoQtd6EL3nC6O4EIX\nutCF7jmZjiCl9OaU0mdTSh8g7f95SunDKaUPppR+UpB9dUrpIymlj6aUfqTXpC90oQtd6EL9KFnn\nCFJK3w7giwDemnN++aHtIYAfA/BdOefrlNILc86fJ3IbAB8F8CoAnwbwHgCvzTl/pPtdXOhCF7rQ\nhZrJjAhyzu8G8Luk+T8F8JM55+sDz+cngsArAfxazvljOedHAH4OwPfMnO+FLnShC12oM7WuEbwM\nwL+TUvonKaX/I6X0zQzP1wD4RHX9yUPbhS50oQtdaEV0NUPu9+WcvzWl9C0A/i6Ar+s3rQtd6EIX\nutCpqNURfALALwBAzvk9KaUhpfSCnPNvVzyfAvC11fVLDm0spZQuRY8udKELXShIOec0V4fXEaTD\nv0J/H8C/B+AfpZReBuABcQLAfnH461NKLwXwGQCvBfB92iB/9/D/AzLBVH2OtCcA2+oz5Su0xTFH\nVsvQudR8AHBVPZHNBthuxteb6voBedJX2zFvrfgNTwI//pVj/tGk6A1w1wnTxN8W478ivVmtjeri\n+HDg4xKOHK+k46DnDf8CeMMfdMzFo0+bm7cfGH9RNCpvTCT5uiH/V/SG9wJveEVAx2zTsF56w/8N\nvOGPBYVqmDnMGHwg//eWydU/S+cApP81MA+FTEeQUnobgIcAXpBS+jiA1wP42wDeklL6IICnAXz/\ngffFAN6Uc35NznmXUnodgHdg//V8c875w9pYz8fUCRTq1U6N+YiPGOfasF9VSjYJ2G55OSQg1df0\nCTsMehoAfIWTnxp3z5gcD9cmGTKPbFSe63scx+egyWjG3/qGe6CQpj+qq5W/fhYtY0Up6sSWpl31\n+XEAFCj1oOsTynr5dzZLDzK/TjnnPy10/RmG9zMAXlNd/xKAP+SdzAsP/28wRvLbQPsGwLYyjLUB\nv9oC6dC3rT5vEkHzlUyiKLE2CrRPM8CUV7v+AoDnYWzgKT+Hfuk1h3C5qMCD9iMRwAY8IpUQt9T+\nOPbogCMNvc9F/h5UHUX8nmhDm8cTkJ9FRJd0fZuIexYUbc9B/YVy9X8kcd0SAVi8Zfwe98XQ0rgi\nRM/HNOVTPm83R8MtGfTtdm/UC9VIPUl5Ifoyt/BxvLXRDKZ1Hr4YY8RjpYE4w8vx0Jef++u3Rg9A\nPAIw2h/+QUyR35Lo34P8I29M9O1Sxn/4RzCOCDQddzgtBAAPvwG+Z0FpDuIvNCDmFKLjeng9jiNI\n5oGyU1FKKf/Lx44plzrHfrU9frc31Rd9hNZrY0gXBySjvBVkgPFLTA0t7aPIXuK19Fqof8koQEL2\nnrbSTg1QNAKQIglNxkLmms66fy5P4av/t2hOPp9+55aic0YOC6Hfie6WcVrWHLyRgjdKGID0l067\nWHwSet7zDh9oUh9Mu8SjGXRJhhrKViegzcNzPdcJcDl+ToZ+bbxpIO7bwhlaSV5rl4y8hJKthWFL\n1jLsZU6eRWHvGsINmnHy1zIIjBPRWchzr2sgKWfeA8/WaD+qrxhtr1y5D4vfyzeTVuUI8BUYG9fa\nQFD0z33W0je1AZSiB41P0ydde3mtyMC69ujgxqXzBHgDKBlFrk1KTUTbNbSspT8s9O5B4dtOPHRM\nLz8dIyonzYHTcRuMP0eSYZy7uEqRfsQAR9YUat5eEcIMWp8jmBMNSDya4ZcQvuZULDRv8faOBDzb\nQjnj7t06Ki0Sz4kENHTcGgloKSIPwvege+8iMXXKHuqF/i1nf5eJcwRz0HStz6undgQemUh0sFBk\nsC5H8HzwaNyD7Gk7FzFoMtQISHwROUu2XNdkbQ9tQf6eCIKTk9okJxCJECSjrCFtLRKwELrlBDwI\nv3YA3vUCDy+dQ0SGG5f729wX4gxla5TQGh0sERkUvoW2k67LEUj7pHtFA97F4QjCX3IdwJKnc+Nk\n6BwlHqmNjifxeXklvRq/p09D3RbS90YLMPTUfDDmJMm0oHePo76P1Cs6oE7Aq8O7blA7AY13wfWC\ndTqC2mDTz2DatQhAkqFOgcpIDkMz1BH0zxkfK1rwIH0pJaSNo7VR4vi0aIAjCX1b0YDVJ5GF9iNr\nAh6KrB8QmRyVKyQY/yU33qyVWP/bgqQzeaReHbWhjshcIoIDcY4givQ1A9+yHkCNrRY5RKIMC/1D\nuK6/mR4dHFqm9+SV42QleU1v0ePh94wh6Sv9GkI+tGcrYgAwpP0/lQ79O6/jqPTttpZyScdYbtgk\n5NSo647QZjdM/+QNW+U3Q0ZqiAa2gbWFK080kLGvOHAv1gi4sgq90kKaQddOBGvnAKz8/lxHEC1P\n4ZHxykltXoMvyUs6Kn72u6707bZyHwAMml4Aw9VG7QeAvEnIG5tvP5/jDfqzCOOH1fq+D8LDXcdp\noWWIdXlXwGYGfB4H3r4TYaOge2ePfWOarnXeI99y8d26HUESPgN+R8ChfqpLiw6kFBSns762tnzS\ncTwonUv7WDISgufmQqMNbVGYjJHBt+8YpJ4T3170DBvm1d6kPdLlXvsD/24zVbpPtSQM2CjGcC+3\nMxcaSgTv3zKk66znmarPnjHssQvp9377aSvmVNrumupLTj0jua3xy4/Ix7+Kwjvie5y7zyddc7No\nXY7guYf/l1go1vh6pYasUhCRCIK75gyzJcONe7inTNoysT1D2rfRr+lwxbXxxma42goonm8vyJjv\n2z+AaB9QDLL+Shc0rUfoG2T4IoNap5d/rpxEO2ORY62OwkpwbTskzadr7X6d48jBh+z3vHqUccSe\nss49z110BFylSY/xp4vF3EIvRbiSTBlHigIk/RKvJNuD/9CWicyOfLPzBhgqnnyIXmrknRMwbDaj\n8XJKQErYkXxz3mwmiLXgForYrXYmk4uMjdC+n7iElouB1mSlfipvOYOik41QGN795+lz8MgVkubt\nJU4nN8a5nQIPEaaUZi6Fc9DFq5PK2nJHZL/nfVzk3I74fOPPpXU6AmnxV1ogBnzOQ0P3Xj5NzuIV\n+LNyvduS7czbgyGv2wgaH7abyWIhh8AHXDFtHN+0bSegZ45XGqvwc3o0XRa616IKa0zKY/Htgsh9\n7loAtwawtOHOhwRFb5qTw/cS5zqj404jBnvNYIwrvZHCHgJZ+pZ4but0BNLaAM3JS45BQv5LRQVk\nbUBD6BN0viHI/JALv7lOaZL7poicIugjcqZfsbHx5ND6vi1N2gCYY9S8O6GdQ/NyBKAhfH4cz3iW\nbspj8x2jDS/qp4ZVug+vfKH2tYW7QbIhjUcPFIYcdchofkPeMC2+qnl96F/XN4fW5QieC3tL6JJR\nAZPzz1VfvcVwgtQpkr+qjAjJk1N07kHm9HrnQO5zET8XQQCRKEBC9C3rA3rkIPX5+o+IV48Sjl8W\n63WkyN2/e6gv4rfWBubqX4q8bnHOGsGcaCG6ruBdH/BEEvMShDytyxE8B/Kum57nCer0yna8N3zY\nHK8z6cvbDYaScjnk1EsKZthQVL5lEbl0PRj9+XCjEZ2S3qkeGd1z+qaoN7G8FqKPRgGDqm8r9tV6\n5/LU87D46nl5+aX5tMjLeo8RTCt552GvonioDxKesybAySeH8a/nfYxWphGFJzrwRhAttC5HUCKC\nOpXjSRFV71quHMBuUxn8VB3ySeTwTkrYlZ8oSwnXdSompQpVJbIdcHp9PUrZ1Aub+6/EUOkabxOk\n11N9+/GmJ7T0OXF6vLrktoF9xflxAAmZcvdcy3Bpjn3fdBzPeMd+fv4RniOfPFeZH9DuwZatyYP6\n22gJ7FnTcrGIHCn4x9xUC7xe+TSR0fmPW05lvvF21vuQGiprBATl3xjzKr9eb23Mm/E2xxq515+H\ntD9xWV7uGrVr6Jqi4DHvFA0Phzx74a1lpwh9O9G954tFD1y0wOf/7TYO9UuIW9sVNOXVIgEZ0XvW\nCuydQABn2IYKHevrAEdja0cBYwPvRc+9IwDvekGPKKMneXfESOsBrWMstWNoum5gI/7jmyXz9YwK\nVuUI8rOB6wrpswheQu9pjHBrYzxGtWPUOuXz6Jjq2f952yKHNqTvQf4cOubQPN8GcKiVn+t+HC+y\n1xF/Vvr2cnKtCRmZ+yMBwELrdSTgQePj8byHzCRE7pdvH+O81A/19jhsNtXhiQp8KP7KyeeJHObQ\nqhzBl567xUAQ/JDKIt5xN8txV8sUZdefc2XIj0Y9TZAph9qPn6eIfUBCwhThZtBo4Lgj+Kh3ei1F\nD/VcjvsG6HzKPY1l9v9PUT5GbdvJeJnlq/fNY8JLT8VmMudj33aiW5pvPYbUR+dFOepFYmr8E8aI\nWctCe9cO6jnVesscKBVN0gIxN2eObzqez7iXe/Ho1eYg8bfz+tcVaoMbmTP31y5wiBJtoUjdikw0\npP+UwMfpzKDrBP22ka7KEXzxseccPkkIftzeHgFI+iTUnkDz/X60ThGutVYg65suVdlI3875x5B/\nxnRRWRpH0guAPINxnx4N6Ih/7toB4I8a9Pnw/PyaiE+2UGxtoW2M81E/xHslGuilowHfGoEv95/v\n3xrBU3j8gFLHaJz7PIzap0hekzmi+Ok4Y4M/7auDvmGC+PUtnwP5k9Zz5/n3uX9dZ1kf4Hhiums+\nOiZAkf/0Kz/mldunGdYjsqfUY9uoLH9cjJrOi9cjzVPm9VHvYnH1vXlJewZLUJs7kw9dUXqaXMuu\n3Xeeegp/dDnaJ20xpc9B4qOrhj1pVY7gSzfFho4IM4/Q5kZA8DU/NXqSjK4D4ri6npJiGCpZfXeQ\nHRFkZjxuDlZEML0vSZesDwCrk7sHvp3XU7fLqFmOCup58ih9Y/Qf9QN9owIuHRaLDHS+vmsGlHpF\nDMu6GHktoIztG5/XY+0di+0u8iD8sU5+/vo9x2hVjuApPCFGAbtR+xXLo0UEuwmfjvo5uVJwTOKN\nHgKjB744hMvpmPJQPXwkQeUkXbRtDvqvSSpLsXc6bcjfQv2eEhDeg2JRtN8aHXCyrXokKmU61kU0\nDdJGutH2V1OKRwCUfCUjAPmQWc3Xo8CeRKtyBF/Cfo2ALuRKn3c3nxPGYb4kMy5HvO+rF6Oni8K4\naUugi8vjMWqUrMkWBE4DvbE8RvrHOmqe3YRnHElIcnZbkZ0uJu95p+NIOmg7N/8jHZG4VI6ikLa4\n3NpP583Pkdd31GtFCHx/dBunt/7P2raHRsi7lTS6lVLS69EzNe7yVoM0GUvifRzjheDTHipblSOg\nEcFuEh3YawFy5DBGyd4ogOqxogAP6p+P+O3IgY8Kprt5aJTD3aM2rjTnvY5YxODt4+cwzvdz5EH+\nERQfRfwlbcm1R4gWrpPGukvk30Xk+xEZTacHeU8RvVUu2uadQqUp31IRwqocwTEioAuIqUK+YxSv\nIf9C4+jhuFAsIf06KpAQtxZBUFk6b4l/rI9H82OdXBTgRfxTRM6jfX77phZFaPxSMTo67rRvHvrn\nUD1F/nKBujGfNB437rHNyjTLsp4xveSNIs5Fc1FuwmPqU/Lol1aHNFn+nL2E6ukP30z5aCTB8fQs\nRb0qR/BlPAFAWwPg1wNoWWJJnl5r+X9NjqaYSo67lt3z+a7pXMo9jGXsvD9Xnnkvt1Hlyj15owC5\nBLVclO44jq2fjqONJcl60LoX0UeRP7dF1fvKamh/7mvfspPotDR/N4x1d9ES0jVpCDyynuCJEFp4\n5tCqHMFTeNbMSOBYaplD8FrOvwX5+1C7HQnEc/++stNAaxRQIiZpfYCPArylp/PNs52iZE8fnQPX\nPy8CGLdH0X+0BEVNvcpRSzSQZ7AOksspSNSyJrC/68eadEVLUnjLUHiigzFE1I4+ttOqHMGTePbh\nk3zQ67o6iAXU6Mt/OIwexPJuLx1GxjZ+yIyWfxjnjKfbSvf3R7dkTts8ZSW8uvZt0rbRqU5pi+n4\nXrQxj+2AfABtOj7XL285lfvHPPwcZF67xARvdOeWmWjTdfvJv6vIZyxbS1Dw0QEvMz3YxvN5Dq8t\ndbhsVY7gqZvUUI26+cJwFAVTmYIoOfRff97Lcjr2WTpt/aC+HheCmxais9cEqD4O5Y8Ls3kjAa6E\nA1fkLRIJSO1Z5E+mHr6Poniaq/dECbI8Nwedb4qorfy/HHn41g32vL4Ko8Pg03dbaLPxGztvyQXt\nhEisXETRp+/wecbg9ZS3jpa1iNIKHYFWXsJT1I1+rv/sYx3DpE9G/1oUAcC8vp70S/o4HfRZePVw\n9ym3jecpjSHNz26Pof25SH/c78nWespP3PDu/FHDjcy1JqPLxvV46dSOox+K3V55jKE93mYzILGP\n1IgKtvzhM5Z3Ytw5h2Ij/p6HyGpalSP4Mp5QIoAjYj4uKRVjXSP+rSAzjSDGfXKkoK8TULTviwa4\ngnB5NL6O9L27g/zRwDQ3L+0CGhjeWq+km7ZbfTbar8ZjkHDO47eb5RnGPDkDmf4oNA5tTLuGwIeB\nN9jSGCzvcD8jgZpCUcHGFxVIOlMakJRHycmljRAROHkpX0oOns6Hy1blCJ7Gs5RooIT4ejRQy2iF\n6JaMBmT0zyP9JdcB+LIT89cA5heZqxD7BF1XfSz69fcDwLWDZ3cdiAZURM5bkaFTNIBu0cC5HEc8\nImBNXqdoIG01wy/L89GIEA1cOaKBlLHd0vbx9XTMO7pG8CSefXNrdKtoya4dt1Emhm8qQ/sKgs+C\nDu7A2P7/8TUn673W9I3b+GteT32/SZEbP5e6rdY3lWfmVynJOd0g8JwxGnSoUG2uBshIN4h8j5Ix\norEciRaGce4/M+/lQNA0xzNB3GTux/YEZMYIa0hcQ/PeCD+K9JfJHJyeIv7OGzFMl3fGf24tmiBj\nDACQBiAxY5NxHh1GmkQDCRMntGHmME5bZWyECGQOrcwRPFFdSWsCErKfovFTrAfQnT6Z7ZciAh7F\n0/vY80TRf+GZtg04LtyOeHdM2wSB6u1T9H2ICDIw7ISIQJApNEXTVf9uwxhu8naxaJ8xsOIW8wiv\nQvsvRz/aoScoXC+ZViroLLeSiDGQ2M38ETaYOLPM8OWUgWqtYeD0LRgJ3AzRXeMMegZPCOi+PgR2\nRLTD4S/qiQho6QkArJwUEdR6xrLzSkhwh7C8B9CmPEejC4xR86gtj9uAo5Gu23aV4R7zXgntPD/y\nEflPvsO1c5j0ka/npB8YvdHc+2FH5fu3T4/KZX0Wf+mzimN6aY4jWYvTmJOV8i2Z+MbRdEmyrEzi\n+fdLj1OllHfLKJ3wZL1/Jq3KEZSTxfTEsBYR7A79lowWERzljovHx4M3/h1CedK/UflLWx1VTHXs\n5zAMqVpgZBB89rVllu+woFyPkY/jDIzOvRJve6VztCmCGHFtw4Rl0D3bur1bv6PGXjPOcxyEV49G\nyxWsnE8Ro16T1wiOs4dxXf5lI8boK7yeAHWyiS7p/TNpVY7gaSEi8EcBR5TMRQSljysnUfqmOmll\n8H0xNU22vp6WYziUo2By6DeofXeYu4LYdzsG7e+YqIDhy1VUMA4zttO2YXPMjY94BSReGzfNELf2\ntfRzPHP5NH5LpvT3We/kac0OoKY5Bs3jENrO++mykWihmyNw6JhBK3MEj6Eu7aCVkdh/Loid9o0P\nL3GHygCM5EokIPFqsvv+zfh6KMb3cL3Tr3PmZErfEVVLevaNB4PNtQ2pMv6lDUdDU0UAolGe0+6V\nseRa+jkeLh3E8Un6JF5LxiPr1TFH95qoJc3hNYQnSw8J/C1G38Nzt1NDz4KE9qUcv1ZempZdtspJ\nF510vJsooCDqww6ZMRKv8uJ5j/Jpvp3Ls4s6AGSK0Ct0f2y7Gl8DR0fAGd6IAfcg/aUjgB7ofy7y\nP5cj8OiZq/+c1GrM7pIT4PhO7ASAlTmCp/EY9mi87P+XS0Mci8tNETmH5mt9klw56FMQeM5AHpKI\n4ncTdL4h1wZ6zzhuD7xB8yUNU+XWbwzz4XqC5IU2zqByfKMxqjaLl7ZreX5vdBDp8/J4F4MjUUJp\nb/199KVTQ4B8P+ekuSkNjxGcszagyZ5jbcDLM5NW5ghKraF5kQC7s6cg70k+vkL6wxjp73aUV0fx\nUwRPrx3oXUPumoH3tkWMezQK8DqCpVNBvdcIJF6NX5Pxykd19RqrF0lGMqrDSx4DeY41AY53ZZHB\nqhzBU3gctLyChPDZ0hI0xz6kw0GkhDwcUjzDOIc/DJu9AwAOZQSqyKAcdip6hyrfDuwXUkfXh7/U\nzSmvRK4P/2fyj+MZn+zy91ltkX7pcz1vr6wmo/W1XC/R1tLeq7+XTE+9vQ4499BrybT2S+1VsG7y\nt7ZFr2fSqhzB0zcHyuRic0C99TFVh5HSvj0f+a6rvmG3GS261nIAqkNLaW/Ab050pupA0qHt5jAR\nvQZ/nTFFY5ZMK4/UJh0+4ngH8AZBOkSlHa6S+rhn4tXp6ffyRPgA+dn0HmdNupeiJa3PHN1RWS+/\nh8/iWeiZrdIRHPezpxGqHypkPuzPeN8U9jr2bSfof4ToM/a59giapwje+2+kz/i/Z5/VFulv+Vzf\nfwTtz40MWnmkNq2PzjWiz9PfyttTtgdpyPUZpc+StXip7khU8IzQXlMdFTxj8G4cfBaPR8cMsg9v\np/RmAK8B8Nmc88sPba8H8OcBfO7A9mM5519iZH8TwO9h/3V8lHN+pTbWM0/vfz2oRuvXE8TPo/qB\nIHzU1xqipwjPQuFzr3tFB3PbIhFCS7uG+DX0akUKljzgL70QLdHQirqXQuu3MQqgtBQUPVVEsIHP\nKCf4zjRYPGeMCN4C4KcAvJW0vzHn/EZDdgDwMOf8u57JPP3UY4cI4FBQbZTjr1D9IRoYofpSTqEg\ntRrR09y2F6UvHQ0svUbA6Zf6LV7vZy0amLNOoI0ptXlz+hq6b4kG5kYLEV0WnTMi8BhILirosU7g\nQfWWbGSNwLvO0LI+ULc9o4w1g0xHkHN+d0rppUyXZ806tG/gmaceHyN7kvPfI3sG1VNkVyMlrY9e\n944OWtA/J8Mh1zlt3LhAfG1Ay5e3rBtYfYAPxXuRfhRRtyJwT5TTSret8JwH9c6hy9pAE81R/7qU\n0p8B8E8B/FDO+fcYngzgH6aUdgB+Ouf8Jk3h8OTjhw9bHdlbKD8SAVAkH4kWtP+pPg/6b40QpDYL\n5S+9PhDp03R6+CM8HMruvUZg6fT095LpoXdNu4Q0uWcg7+qx5C20TfseCePQ8SW+TYCH6+/4N2l1\nBH8LwF/OOeeU0k8AeCOAP8fwfVvO+TMppa/C3iF8OOf8blHrX/uJ4wv2LQ+BP/bw+LJxe+nn9tWO\nQOKdc10cgXStyXB6tbEkPouXiyBoO70Hi79c0/vw9NU6IfD0Ol3cwsvxe+UKj2V0l0D4S54h6HFe\ngNNp0dwTxtE6RJGTyJ4zCS2Hxz79rv0/jncGpUx/DYRj2qeGfrEsFnv7CN/rAXxBWldIKWW8L+vG\nVfqs8dE+en0qZ+CROXWb1E6dldWu6aeGT+uj/VSvp9/LI/HSe7T4S5sHqVuvmlePl06dMupd+mBO\nwTiPHk2eO2EcOYzmOYVsHSrz9P9MQvb+7qlC3ohgFOyklF6Uc/6tw+X3AvjQRCClZwPY5Jy/mFJ6\nDoDvBPDj6ihPY4rmT4H6JV7rWkL8VL927ZGZ2yYh+khU4Gnv2Rftl3gifByvxS/JeOS8enroXgtF\njHdkTablxLAmG0X/3Fx7lY+woogO5Nk++jYADwG8IKX0cQCvB/AdKaVvwt7E/CaAHzzwvhjAm3LO\nrwHw1QD+XkopH8b52ZzzO9TBvnz4v3cUMAf1U0PdYuw9qZ1TRAASoo9EC7S97puD/pdE/ktGCJyM\nZ5yIHolOXTJiDkV/UOZRJ92WLsnwSjIc8ud4OfTP8VGET3ms/k7kSg2dglJKGb+UbeSvIf0logCK\n+r3rAHQ87fqcbXPaaR93754+Tae3n/L04JP4PTKabIueOWOcg3qtG0TTNxEd0eqjXn5PKspTrM4T\nHfz906aGTkNPYR7y90QFHt6e1608XrkWfZ7IwNOuRQBRhK/JeuS9PC28HL9HRpNt0ROhU+G73usC\nwOnWBiKoX+L3RAitawOPAv0zaX2OANDRvjcqmMNbrgHbGHuvPZHCqdqkdmktQUPz51of8PJoC9wW\nr8XPzSMq69Ul0ToC+jFFnQM1pp61AWmMOkrQ9HBbMTU5KZdPeT2ov+YrPJIeq78TrdMRFLp2fKbX\n1KBHDoHRA1LRa894dTQi8SzRJrXT52XpWOIQWY9+L08Lbwt/b/ke1DoHbxmFVvIg2x6WqueBMTrn\nHofGPDz3oejczdFu7nCYZoStPklnuaZrBJbx19YUrINlGo80X0lOauvBK81Pmi+VifZ5+lt4NL4o\nLzB9Fl65Fv6lDo+dkzwOhePRnJi3vAOnQyvvwKHw5OTdMJ8lHkmP1n+Xy1DfOIIe6H9OZLDEdSvP\n3DapfQ2RABcdRfo9Y7TwtfL3kl3jOBYtXTqi0FyLtVQZCU/U1KOo3L2ICOo1Avp/L/QvIfXo/1Q+\ngvxbUP4SEUIm/6x2Tb8WJWh9tK1XJGCh6R6RgBQdeMeM9q+dIouXrT8WIzk+7w/AaMidu675tYgA\nAh9XboLKXiICQteQo4Gl0P9dXheIov5ou9ZnofkeawLS/bXyFdKinyVlvbSWKKCmpdcRCq2xqJz3\n3uei/QWt9bocwVNoXwughicSDVjRgRf9nzsa4MaweFvy/62In+uzZGn7OdYFJJ7e6wK3PSLw0Nyf\nlbQQPSCwPUFhAAAgAElEQVQXnePQtiYj8XsL00nRRf3TKUnpp59p/wqKzi1D9a4hDeUvtWuox7UH\n/feOECLIX0KsLWsA594ltFRE4B3/QnE6tTOsKte76R5FAiccIkDXGBtJajAjfSDXQ9WuyXqvtXMB\nkgzHY+nhdGn6ovL1c9HaNf3ePi5a8fRb+imPxXeh+0Mt34MCAjyF5Er7RuG5NnjoeNH+DrQuR/A0\nxi8895nrg5OXXhdj57326ofBQx2TJtfador2OX2cA5TkNZ5bQbd68iekJY4pz6CoE6GHwDQebvfQ\n9cz+GbQuR/AIMurnDP8pIwJtjaAV7fdu4+Zu8Wrtp4wKzkqthtrjwVppjvOYm9eas/Lby5h7SnVG\n6AROJvIdfwS95HTdz0298zu0LkegRQQQ+ug1zXVHIoSWawndtyJ7Do17EDoXvUi8Le2S/lqGkzsp\nRQaXcmIeOS2/1UJex7G2yKLnr8dQKvssuTEiZovu94w4hIWdh/Vn1yKAuq8DrcsRaBEBvY5ECB5Z\n7dqrnyIBjofT29rGjWHxRlC+hvAXJY/B47yqlzeq3zOGR0eh+7oSLZmbqEWjcDqig64eW3LWL8VQ\n6uw8yldo4a/MuhwBjQhArpeIAlrXAaJRgcSDxraWCEBC8ycDmp6BtIUFygPYMXI0XKEGPxqDa3mB\ntSH6U1OxZp6f77LoEcPvQfwUSlsRS/198NSh9kRAK1sLwdocQYkGI1FBFLWvJQrwInsPSj87otcM\nnIXc637tdJqmQ+KzDO8c1M/p7gXbzrE+0GoKeho1bg5zSpl6U1f1Ku9c3fX6ivVDCNZ4p3MY63IE\n15geIKP/txw263VdI3DOCVAZiYci+TmHyLhxu5NmmKzTYnNPk3nRv+fUmaTXKyPJSosnLdT7lJqX\ntNNWFvU62cSh/OgagxfBF6qfn3UfRbem13Paawc9AirORPtO9XUS63IE3jLU1vXcBWDOsFrje2S8\nqJ+Tldq6kfSl07zMAB1FWxP2FKfx3nT0Vzpaf9VjyT/CQr9DuChtsGy1uVYTdYWYk9zC79A8c3rg\n0KG9O5r8FXqnGdflCMo7FkHrp4gEKMqnuX4reuB4JN1ainkWaYaeQ7SlnZvMWqKAyFpBNGLgxojK\nWvOI9K2RNtjPOVIJLqq/PP/oGkIE6QM+tF/PyYoKnjbG1XRo8p7x47ROR6Ch756RgWX0PeP1jgRm\nU3RFWEKhGvLVkKuVp/KgXi8yjqDzcmw9Skvk3W4j8vdSHSH0em7FTEVfmCvsHYh3Ht4o4oGDz+Kx\nxtLko9GOTetzBK1of2nkP2dNAOj8d+MMmmSwSruE+Lm9pJKX8iB9KYqw+i39HJ/HMLREAtw2sla6\nSxGARHX1tDpCkKq/eak4lPJDJZGIYBOU2x74PVXxLLS/NXissTbQf1C5f3pyfY5grcXlWpH/LIrs\nQ48ie4nfWhPgnIql05pLofIALeQe3fNaxo3I9I4ErPWU207F4Pc0KWWrZ3EI0edX8uxzUoFzdEq6\nSmpH6reiif7fo3U5gnqxWEPs1OCuYV1gFnkRvoSETxkNSFtArf3zlpdcSzRAeee+dFb0c5uJQ/w7\noT1KdQ484sBrlO1JwZVxHkGfc7knTaf0S/Pc/KS//1YZY5loAFibI/Dm9+fuCmrZNdTt4JU3hy+h\nWI63tEn74r38ZVzpRrWHYKFuz0ExL3L36Kp56/xehKKH0Tykvch9d4L0I81A9jIh2hgeZ1wf6KrT\nVC3jcXxlp46ks0QuO/DzpemziHzCcZfRMlHlOh0BBY/adYRXuqZrAVr2w02cAg45S2hZauMMicQL\nhV9qB6PLo0+Ss2S9OiifpS+qVxunV/infanWnDbSdqj0WviuTVFCDPluEXP2tbOw+OuIQaLiJCRd\nZX7SwrUmX2S58fuZ73U6Avo5iui1n6Kk/Zx8E0WQvgflS8g8GhFw7VJqSBu36IPSv9aooNYbkSlE\nkcFSyL0l7PfOpSVV07pbxyI6l5aUR43UIxFDvSV1ru6C4qW5WykyTV6TtcaN0zodwdIRwOzvtYT2\nufoTnnoTgH/PqcQrIXUuvLHQvRQSlT5pVxBgl4Cw/gDRyMD7MlDkEKUlIgTPWKfUydVC7knF3JS5\nlPGi6L/osHL29biF3/rb1+hd0l2bTWnu9TylKECS12Q1uXZapyOgn+dEBNx1mDzIPAtt1KhyKF0y\nvjvwixXcF0sKayR+LSrYVTyUrNyZJuuR9+rheCMGNKLf0lHoLiwA96B6wbZOK0VfRE5PdL3AGrcg\nb0t/4dMWg+szFNI6gTanIq/J3qfFYgvF0yhfy/k3EYegrXy9lP/37jGVeCPynDOwEH5rzt8qEOdB\n/Z7tpJauqE5OJiqn6SnX98UpcCd+W3drFSpmKaqHmjPtb0p3JUn6rxx89QE6Df1LcyprFpx8ORTH\nvXP9Thev0xEAsbz/7Ejdg/gldO6JCjgnIu1R9W5X6hkZtJ4bsBZXpDELRbx2S4Tg5Zdke6ZoJGe2\nRofBHXSSSjy05lq1UtAep1zm4o04KL+G/Gu+pdE/J6/JWuO20XodAb1ePN3DGTWaI+SMJpdH9CL7\nSGU5zzh1u8e5WfoLWZ527kEyQI8youO1zqH3mEuSdT9zXm3pb91qeMr2x5qi0RiV95w7qQu3WfyF\nV4ss6zlI70Sth+qoDTwnr8nWUVJ/8LAuR0Dv0dr9EyJqAKXtRNqiIIf4JRTvaZNuai6vxq+hfG3d\nQFsotmStuVKeos8i7yE0KpOrz63E/X3WRIuWqiW0ET4XyvA7Ug7pA/b9lB+lqSMNacxS5tnipZVV\nJWSvHWKrdXAGvpansnVFVMkx9KH1OYKuyL8F9XO1JSgS8pR4iOT3I2i/VykJLTqwDLb1RzlXYbmI\nXkpd9hALtLZoIkJWOWWgbU2gRsecLksPnZfGT3+eUgM0NIrQUL/EY0UOmnwtSw+Q1XJ3fddQd+Rv\nGXYO9VuFhqQcv6dAEedYuIVmiVfS28KvyVhynq2gmu6ap/7f4o2kJ1oiBk4ewXF7kRVhWdSj1EOL\nE6vXE6TCap5UoKZDK8NA71v621E0z+ktDkuLHuhvGWion5OXDLwmF/n9BJvW5wiayDL+Vq6fQ8Ce\nyMCL5L18PdYNgLa1gLnrAFa/d1ExanhavjSeFNXSc7jNZJmN+m/t+ZtLUYflyLl5SOOV3TdUP9VN\nja+0W6fmocidppQ01A+Mv/Ot0cI8WpcjCJFl/Gm+n8vtWyfNOJ5W1M/xSYY/ivglQyTxW1tKNUdi\noXIPCi9je6OA6BpA6wtSxpmLxL3jSNQz+tB+0L0VUUYcH5e7p6Q5fw35Sjn7+n+Nl9NNv28luqER\nRv03ono4A17Lt6L+Wu4urxGYNNf4e5A/t+PGs3uIk5PQveecgsQr6dX4NRnPeoBmFD1nATyGI2Jc\nvGUo5oyhEZcWXCtpf9uow6mNmdd0eJy4potDvlrJawmISNFDrVva3UTHpzy1Hk4HNfxStEDXTWq9\nNGq6y2sELFlG00L23MEvbvGX7v33LPRSIyqdGZCiDe+agOQoovyajNXn2fPvQePeFFEU2UcOn3Gy\nkTMNraTtxlqSpKggun5Qz99b4tmzTqGVaeBkyzy4XL00HpfeoXx0d5P0e8y1LqqH6qDRQP3eU8Rf\n2x1tXcKzgB+jFTuCaN7f6rdSQ5xMb+Qvnfzl5KX0TPSMgIb45+wO8hhq704cK+qgOhHgr8dokdPm\n0Evf0jTX0aTq/0g6Sft+cMjZIy/tNuL4uYVeiZfjo+8/Z4Dp2oCWkqLyGuKvIwx6v8VO9XMIK3ME\nFvq30H20H+CjAQv5A7yT4FBSZOeQhO6lhWUtGpDQp1YaQntxrRSPB+l7dhoVPjj4qEyPMhGn3B3k\ndZISeV7fubuGgOM8reiFy6drOqXSypq8tHgrRQ7c/nvLMUj6NKQPTNG+VMSPk5X66nvQFpnn0coc\nQSFqyKix7rE2YEUMHI8nrSTpiiJ/7o/cshtIQ9tzFoU9Bt2L9KN7+CMRxJxxvHMpdOrooPdupZY1\ngJq074SG5GuyyjFwRPP8mjOiC66tTqHokSKbokOLErj7Kn0c2r+G/znGaGWOIJL+kdA9vdYcypy8\nv2cNQTL+ka2j0Uih6Nd2BHHUA/H3XCPwzGmufk0HGsa+7RRZAyimwxtxSH9LydDWZO3CoeN4DDkQ\ndwqcuawjG854l7lrqSctQqA6y3P0PDc/rcwRFIrm/63+ltw/x8NFFhzK5PgkAyXl6VvODkhbTlsO\nkxWdGsr1GGqvYS5ph0g+ew7Cj2xhnUPcGKeIHKiR8KZsvCT93S30TsnakWONHYkWPE5B0yelaApp\nDqHMg5PVooBHwpz6pjBX5ghqI9cj/2+lfzgZT+5/KfQfNf5S+xzkrxlXj7GOGP4I4m5B+ksg+3Os\nJbQQB048VAxOq9Pgvh9W3r8mLqfvHc8apzwDzSHUCN46v6A5BCllpKV3JEdS5tR/t1ChlTmCQqeI\nAHrtEIoY9UikwM3RateM9JzqolZO3mukrV1KkTnN1d977LtC0t+Ry3V7ib4LEYNWI3VrfA+qB/T8\nfc2zg+5ctBRNSRlx91pshuZIpOhgGWewMkcQdQDWGoBHp2eRmNMzl0/i1ZBytOCcJmPJeYyqB2l7\ndza0IP65uyaWWDyeQ1aUtBwitIlzjq0Lly174iny9/Jq+uscvURa2gbYf4e056AZb80ZSH3LOIOV\nOYJCNHdNv4RcP7ebh15bC8XeMwJeB9BjrUDK47ei/1OcDfAi85Y1gVakvjTKL89tyTFaHZ9U2nku\nUcfVYk7KPUWKqF3DH6E8gp0y8iy+aga4vNM9nYH2jvV3BuaTTCm9OaX02ZTSB6q216eUPplSeu/h\n36sF2VenlD6SUvpoSulHfFN6hLExusb45aL99HoHPgqgkQN3mIwuFnN6uCiFW2d4hKlR4MbNAm8Z\nizPM9Jlo89PmWWinjAUcn4WWIir3zN0LNw/pHig9gk8vpfI9iIzl0flI+LfmdFKZGzfvnlQ/7yiV\nZ+sl7v2UyLPd2AMytPlZGw+s6Ds6Xl9KOesPKKX07QC+COCtOeeXH9peD+ALOec3KnIbAB8F8CoA\nnwbwHgCvzTl/RODPwOerllNFAa0LwUstGGuoPJoy0mS0OdRzsV4OzwsUTfu0vAC9DbEnOrqLFFnY\n9eproQji9W6l9NybNV9Lhyav3ZMkJ8mU9u9Cznn2H8z8K+Wc351SeinTZQ3+SgC/lnP+GACklH4O\nwPcAYB3BmKKLwZ48vyfF4z0t7E0XSQYwugjcsgZgpXDmLB5b8oWWPjOARhmOToe+1k318+yRftDy\n4BpF0h9lUdczF0vnAD1R4tFx+2jOLxu8LqX0/pTSz6SUnsf0fw2AT1TXnzy0GRR1Ahzq9ywyczw9\nnYBURkJyGFEnYKWBJCfA3Tvt15yAJ4QvRt3jBErKIrqVtDUFUdNSKZK7QjRd2kqnOJjXa9F/rem9\nZak1bvtbAP5yzjmnlH4CwBsB/Ln50/kr1ec/gX1QUZNl4Lk0heVYOB6OT0LJXn1a2oj7EreUk1hD\nFBB56VsM8Fyj0muL6X2hspA699ewtMXUC/npgwA+hJ6/TgY0OoKc8/9bXb4JwC8ybJ8C8LXV9UsO\nbQr9cPXZEwnQa2vHjyd/L60HcLt/5q4H9EwFzVkL6HGC2JpDi74W3T3HvNCe1rTFVqJe6xqWgV3C\nmUXM8L95+FfSU2/rMgOvWxkV+U4pvajq+17sXRSl9wD4+pTSS1NKjwF4LYC3+4aLOgHOkLU4CqlN\n2rnjaYueItbSKVokoNEcJ2DJR3iA8zgBbifZhfzUA32uIRrwGFzrXrX+1oXiqMwZzhGklN4G4CGA\nF6SUPg7g9QC+I6X0TdjD5N8E8IMH3hcDeFPO+TU5511K6XUA3oH903tzzvnD9pSiTgCw00EZviP3\nc9cEOJKcg1QULpoOKnISWXKWE/CkbyKplqhB7oFGL6mgdupR3Kwl8RAxdh79PXYMaXPSHESL8zid\nEwAc20dPRfvto7+FuCPoVU6Cyx1zfNK2wkj+P5re0fLammG19k9bRtmbT19zNABcFoNbqYfRuQ9O\nQHOW2vqKNKYkw7WfaPvoaSnqBDgjZ0UHHE+kTTrly5GE+r28ln6L1uHkz08Jl2cRodvgALyRikfn\nHCdw6nMDy9DKHEGUuBO9lDwGwOsElpK/kExXmB8VXOH+HhDzUq+DZC0mJbIryVtawlMO23Im2jPR\n6gtpP8OpzT/qBPqZ7xU7gpZw3mP0pR958VAP436bFi036OvQygsSfQY9nEF5AS+7h/bk/WF5D7UW\nn4s6n15pIMsBWI5Jk9fGl+YvyUjtfX+Upox0S8jz8t6G8L+HUfOSlRLx9HvSKpF78urkxuhRLTTh\niLDuy5mC4gR7Go/WstTRH67xmqgeDkBD8pa85jx6RQHL/EwlcKscwW2j6EtnhaiSod1CTnlcQY+s\ntoper46az+sMrtCGzMtL0FKyWtJXv1h3IVroXSuoUOSH6TnZiKmJIN5zRwCaY5PktOexfBro9NpP\nTh6k6TF+Em3AGyAuhaLlFiWS5m/JaPQAuiH3GPAH8FVwLAbeY6QLMm8xvvWL1PM3BepogaMyz1NH\nnnN/VL6F5vwQDRBf3IygXY9ujz5LjxUdRHcD9XIAUcdq0y1yBEumVLaYGnIpP861S7wSWpfuRXNS\nmkG3no3HGVgG1Yviy5fUWw10jkMApoh+yQXhW/S6uKmg/B5nBaI6To38PXpa0b/meKQxNYMedRjz\naMXfbMt4ccQZUgmtU6PnNfCc07Dk6Vjaoqlm1Jd0BuVLrOkoRtuzA6denI1ECIU8EQing/tKr+2X\nyE5Nc9I6kq5otBDNb3tRr0evZfwtHdZBsmj+X3Io0jzu9a4hD3HGjRp5znBzDkOKCrjTt5zRLV8G\nqkPKsZc/eiRiAGxnoBleD/L2pHciO3DqFzrymwH1V3Pu1k/pJS9/19u8rbQ28L3XBubsLIrubIls\nIfWYrTl5f0u+ZddQL/Tff8eQNoOVEDXqnIHkDD/l4Ywnp4vj4/RJRlxKEUnGu7zEUmQgIdnyBZEc\njCZbkLe2Y6b+0noihMJnoe4aPUXSOBLqaokaaiovVPQ1WOLXyM71KrYi/Jqi8tEctzcy8PyojEat\nyF9D+L1SP9JaQR9amSOgBpND09TYcjyckefaOGPudRqSEZcQtYTGizHljLqVrtGiA0u27JjxRAiA\nbfyiCJ6+JC3Glfv6Fh1Lovy5hvOUVL4HPebbup6wxIExj05L1xLnBTTjL+mT5iHNv/86wQprDQF9\nfpEM8BeMixSW8/4mgaRX47f65tQVOuePzs+p9bPkds7bXINoyZIDc/BhdF7esTx6PWmTVtRvyUbz\n+JLDiEYK39ql1tBKHQGHKFt+enIuH2dcpZSGhIIlo2kh8VaH4EHjHuMayctHF2PLAvqc9Mqafyx+\nrZTIv1Zq+aGaSBTlRfvWfViH16xFYuuUsBQRaOcK5iJ/yvvNd9kRALwhPEdk0IMX0HPiS0QIpb/H\nbw1Y86DUuo1ziQNdt/2AmIeWWEBsjTii0YR3nB7lJaznZFUYjZ4b0BxR5HyAFiX80bvuCAA+VeEx\nyhyaXaqkdLR0tIZktdSMlbbpUVLai7KLoY9GAXMQfPk7nSoK8PxWA6DPJ5IXPxWVsVpLFRQUHplz\nZMwN+V/js9C+tmXWKruh6ZcQvLY+IMlIDobTw83pTjqCz8OH3k+xHtCDV+PXZCw5C517ULAX4UcQ\ndWvOvTdqv825/xbqXVaidY1gqTpCPUpJW8+od6nppdcGCv+/cVcdQSHrB2ikNk90wPGVNu+agBQh\nSEhZS3to6NpC8la/B4kXHm/11sh3pvC3fs/Ksz/V9zSS0uLuK4ruT/UTjrWtaB2zIOjoVlHveB5H\n4tFn6Wk9IdwiFz08JvFLEcIfueuOAPAhcs7ASr9h7DXynIOReCUDr+XJNaegoWMrh+7JsXvy92Uh\n15smgkOnJNcjGqj/rrf5cFiEqLHoEREUY9OyPdQrVwygh7dHyQjrnlpLSkty2iIyl4qKlpqmEcLL\n7qIj+D3wCJcz9Kf4qUkJbUuoOBolaGNYcpasp987DuWLRgQIynDywPl3CUWczLl/rL1G7XPm0rom\n4B3Ti/Bh8Fk5f8+aQRTta4fxtHYpEqD8nujg6+6qIyjEIX/6InIIWGrzyAK8UYyifg2Zt0YKFpKP\nRAPW37xEBN7vRutunyVKPNDoYB3f73aiKYFeawFzSlNEI4dIxOBZ7/CgfK3fQvqSbBTtazKRiILj\nLXz/2l10BF/C1KBwKNuTMpLQsHdtQEqNSBGFZqi1bZyasZp7SMyL4COlGlqMa6/yzfXzXcf39rzE\nGYweuiKLxUUuEgH0KhlhRROaDs14a6VHtAhJGi+yHhDlfclddQQAb1Tn7ADi2pZE/udC/d78vzcq\ngDEex996fqDQUvv+PaeqlyTLYC1FtWFqGSeC5gufdw3AG1l4CsRph8usw2c9kb4kIz2XSGTAtb3o\nLjqCJ2Gjfw69coZNigi8TkHi1RZRpW2LS6J+D5r3lpeI5P9bUf4Su4C4VOBtJcko9NYdjR6KrPeE\nsIaqOV7PbqE5EYLloFpOBGtrB5HIgOP36njhXXQEzxyuuPRL71IRUoonGimUfDolrUSD5kys0g7W\nbh5v0bWyBuA9NCXdJ0e13jlIexA+X4hfWGzVUxt5L0WL2XkL1lnlJTwHzjQd2jy0BWftEBnnYKTI\nQIpepDlLUUThfd5ddgQ1neLwWPQUsGastbRPa0qoV6kI70GrllPAcw9xnfMQWO901LmL+s4tSBeV\nj9QR8u4o8jxDqyRE7wNkLWUmpHpB3rMC1gG059xFR1Abn55lIyRDq/FyUYG0k6alnEQP5O9B/d5q\nonDoA8ZoP+osav6eOXlt++0aydoGucQ4LUXiIrL1AqqnPITF54k4rChDWzvQZCU5baFYQ/QtSJ/O\nlbYX3sfuoiOQUhWt6wRAn4Jx0bISLRGDNZZX3qMDiO/+uZSPOD31LjXdGqlEq42esnSEpUeLDKwU\nktQnjdfjB+i5Nk5v4dvcVUdQiJuXNwKQkHDkdLF2mExD+VLJaSs6mFtGAgZP0RNZF/DopPqjMpS4\n+d3n9QEJIfbQF1kw9qL9msd7UMyzJqCVfdB0aPKarLReUKICbisnFy1Ih84iRei09Yl01x0BJU+k\nAMzfFSTp1fiXKCxnyXp1ePVE9NXU63cB7hPSn0u9isy1RAiRKOWUkUHvonFaX6+fmfRGBNrf+t45\ngprmRAsaL7cG0LI2oEUA2k6johNKfz1XjSK/NmaNR/lb1wg4maXRfuSEdC9aurwEzV97duNoOrTS\nC5qMRZ4idUtHBZqshO41nZbMXJRfxvDSvXYElLzpHm7njmToJWOubYssRl7qAzN+3d+jhIQ0t8hc\nWvVKcoV6rhec+4CYRh5H0PMQ2dwyFPV8IwfBPPyeA2meg2WeshFSv9bXUlJCkpHGkbaX9qCLIxAo\nUspBigK0baMSmtcOSlm5eU/xOEk35fFECpYuib/18FhkvBbdHJ1jXcGDkntGC/T99yL7Wr4lKij/\ne3cHeU5Ve3ZSWT8UI530lWStaEGSkdYCtMNsHake6vouOoJtXgDIcfcXjQJKX0sUIDkODyr3RgHe\nCMC7U6hO/US+H3Sxf+4f0+P4bgstXWJiTlSwVERQG8w5EcGcInOtW0g1dK/NpbPRt6pSPHUXHcEV\nM5dFSs9EooYyiYhDKX1W+WltLp4SDt4yD9EoYG75iEI9vlvedNjayHo3exw66xERFPKWgijkKQlh\n6Y3wtJSO0GQtuROh+6LSG2TQti/2cQTnPgI5pno2xY5ybbOJPrdiTOhTt3YxXTMypY9O3qPT0+/l\nifBFeSW5K8z33HQe9fO9bVFATUtGBGuMBrx8c3L/lnzr2oDWPpMkhE+H9LZ1ovU6gh2mNyy1zQaF\nLehtB3m72LWisxgzTdaak6Wj5ssOvsjY2nwA4DH0W+uhhn9D+tdI54gEgNjaw5xIwFP4zZpPvc5g\nIW+rX9tueaaKrxKSb0X9WlsnWp8jKO9/fdNcW2mn1C1qAPgnLUUP0oRqshYw80GvVXrC4vHqovPy\n8nNyCMpyRJ8f/XrepohgaSNE9bcUjCtz8e78KfweXk/5CE+kIG3H9MxnYWMv1ZSjQ0d4I2132hF4\nIwFgD2DpHXBtpb0LRcNnOgntBfFGA9nQQ/k86LmM3XJqtciWb2xPtF7/4bS5nSpCiLx9S5SP7qG7\ndxRQyLM3XosCvOMtaOApSVORfkueu31vReo5bTNpXY6Ae4jSLk5gOvuBabPau+4y1L6gWZhEPZnH\njAlFooHyML2HylpQfdkV9RjR1YNotCHRqaMEz7bJQr0KyM0tMxGNAGoei99b2C5SouJEht46Mza3\nnduUJI3prUIRKfQaoPU5AnqTEojm2ncBXq29u4MAfCG1B+lbUUVNjwK8Ud09xptLc+fbQqd+ZeYW\nnmvJIXjH9FqlEyJ4SlpFaWt92tO+RBQgydO2jl/FdTkCDhBvmDYO4Q8B3tL+APxRAk4PqrbFzirN\nSTtRGgA8Xn328M/ZoVOPV7ctRXS+S5MHSc+hHsXlaCLaO+eN8NmiMxp4T1WJSF+Pdi0yoLy922bS\nuhyBNyKIInxpfUCrIi3ZMA2Edl2olsj78tX5+ug3pnjVuXQpJGdTr1LTra/yGY25RNataP0tfT3a\nI5FB77YOtC5HwKH8SJTA8UbbS59UUWGj9BWQqlWbhtK/CLW86L0SkSXdxdF9KC3tyZn3oBUa85o8\nAYcHxVv9Uu6eC4x6RganWAegbZ0zout2BEAM5UtrBNF2Sb+nzxMVeHYwLXKi2ku9DIu1yLjWswBz\naOVGuRd5d67MQfdWf0ue/xxrA5HfsfesAyxgtdflCLicfWLapPw+x6u1a+i+fFmifXUmRrJzmnyE\nB2NQEO8AACAASURBVFi2rtvidE+M5lqJoljvlnxP9QWLR9p8pZ0F0/SXPu4+JJ2SjMXvOSym8Vpt\nnKw2RgdalyN4DFM0ze2WlKoil7QMV6VAqjHH6SkyLX31GFJk4Nmp6V23ja7vck711jmQCwGwjSVH\nktGbw+s532XxePsjfZpOS0ZyktZBr9pYtx4U49osnpm0LkfwAEeUXhPXNijtHMq38v7o3Ofpv8Ix\nVa4Z4xoFaHwPMDXsHiOvOckISYVW5+i8TeTdbj9Hp4XcJTlvjTiLt3dkEEX+mqwnYujRPhf9c/Kt\nUUInMh1BSunNAF4D4LM555eTvh8C8FcBvDDn/DuM7G8C+D3szcOjnPMrzdlwEYB0horLx2soXyoW\n2or+6405nAH0oHVPBHHl5ON0evgluYisV18v3ecij0Gf84K2oHxJ1nN0xcPvTUn0Qv1W/1LIP9pe\nnouH9w5EBG8B8FMA3lo3ppReAuBPAviYIjsAeJhz/t3QbDik/wB8gTkJdZe1SMm4ce2ajNZXdBbi\neDzov+bR+MqzsPhq/lbUX3Y3atW2W8lTM++UtNSxBI/hbNHjiQ6Wjgxao4JW1F/3tawHrCEimIP+\nzxUR5JzfnVJ6KdP11wH8MIC3K+LSkpA+GwmFc5FBXaiOUksEYEUUUt+V0V/LazxevvKsIpECpz8i\np+mo23vtDF1X4tImb0pljj6Pnkhk4F0H8EQFUn6cm9ttiggieX8uStAihx4RQgdqetVSSt8N4BM5\n5w+mpM4oA/iHKaUdgJ/OOb8pNBsO7UuRgdQOyL/LosnUqFuS4/po/xwejs/iRYCfynGyHnkPLZli\nilLrS9QjYpgbHVgG1CPjdVgajEsGzykiAqmvV6SwxFoAlZ8TIXSMYMOOIKX0LAA/hn1a6KZZYP+2\nnPNnUkpfhb1D+HDO+d2i8l9/w/Glf8FD4IUP+VO+JU3jba9TLpIMl/bZVm2aTqm/1smh5CumT+Oj\nc5SQ95XCY6F1+kNAvc6DRSpjnJKWKlWkGdCII2qpPCEdhvLMxzNv6R68B7+s/t5lI6R7lNpbS0RI\nxp9zltIBMdpG5/jpdwGfehd6k+unKg+poV/MOb88pfSNAN4J4MnDFF8C4FMAXplz/pyi4/UAvpBz\nfqPQn/Ef5ng5iJ7lI6xic9Yhr54HxSIHylp/nGeJ4npnPQh3QuqdtpqD8KJz8fJ7+JY8MNbSp6VN\nbkvJCO4H5+h3o/D85Gl/qvLGL+WcPwTgRTcdKf0GgFfQBeGU0rMBbHLOX0wpPQfAdwL4cXM2Enov\n7RIyj0QBaOyzUP4SUYAH0c+JAGqSIgAtMtDmEhljDdQzOuhRXaJHNKBFIN5ogEO4Gs+pogFNZ48S\nEon803iXaPPKdSDP9tG3AXgI4AUppY8DeH3O+S0VS8bhMaWUXgzgTTnn1wD4agB/L6WUD+P8bM75\nHeHZSKUkJGMSjRCsPgs1e1D1ElFAK/Lu8tOeHeZxm0lCi3OoNcKIyHkXGT33Z/FYEU5rpHCuonLe\nqKF1p5B33IU2ULhSQ6eglFLGnz3MhUO1EnKPIHoLhXuiBKnfigQoj8bnGU/j98pYvD3R+9przPUy\n7GuJBDSZSySgt9+GSKBc/6XTpoZOQ/UPXfVaE9AQsIWO564JSNtaW/m84/aWO5W+NdNtjQIi/Ldx\nTUDrO3cRuVbUz7UtXHxuXY7ginzW1goi7XD2UQQc7aefgfGBrEIc0n6g9HFrEVq/1Cb9tVsjgui3\nZw3rA72NeY8IQOK15DmU7o0CJITviRY0hE6vo1HCHLRPP9f3wvFEED+3a2pOxBCNFs69RnBSehzx\nkhFcX+uhMOtwlufwVpRH4+N4Lf65ch49PfSulTbw5dHPeUDMkovwe06qLlk6Ym6fVyZygIxrq78X\n3BmKaNuKyksAa3MEDzA9FJWrPq699HFF5SwZVHJSv3Toq0bwUlkI6aCWFA1w+ix+wF8+QivpEK1f\n5NV7F8h66bSDV159Hh2WkdV4JQenGVKJ79QHxLi+GunTSKelnXMglvEv1y1t0rgeHtrWgdblCOqI\ngJZQsNpB2jUZTq7urw0rPWClyWs8kh5Nl8TvkZHkJFnPNyGa1ulZbmIOeVE+pV6oP6KvJWrwlpSQ\njKCmz+N8Ikje6m9F+t52zrhL43I6lo4CWnlm0rocQYkINAQfjRSsPikiKGRFBF4eb+TA6eTmxsl4\nC8N5EPyciqW3hbwvUhTxa7p7I/8ov2QUJR5pvloE4I0OON1WBED7WtsjkcHSbT14ZtK6HMEV+d+D\n7j2oPxoRcLK0vzbsrWhf4+X+MnPy/prB5hxIyzdjzU6hF3qK1vWJ6Ikif03GGyVYfFyKI9Lfe32g\ndwRQt7dEC1rbKaKCTrQuR/BE9dmL/LVcvob6LUTvze9zPJIejY/j5fg9Mposp6OFbvMC8dyXSDPY\npywep8n0yv17eDwIn5P19EUjgx5rA3NSQh59UTnrugOtyxG0rBGUPm8uP4L6uf4WHo2P47X4PXKW\nrKaDI678922nOUa7RW8L4tfkuNTOnAghkvu3+nv39YoYWqKCpdcFSlskSrjTEQHdNSQhX3qQSZLx\nygH8KWWKtrkDVJRHOhzGIXftIJmE9K1DcN4dPNHDYHd9Z5BGvQ6TtbxtS5SQWPrgmPa8aL7bozMq\nIz2HuSUlligaF+FZsNzEuhxBOVk8gC+kxuXk6z7poJUmRz8Xqreils8c6udKRnDOQTokJhlY6ZDY\n3KJu3LOI0Bp2AUWphOOt1FLewcsTPZAm3Yv3QJp1GIy2ceNph5pqfk2W9s05SCbpknQuXSKCm1s9\nR0lO49lkct03N7tOR1DTGktHeHkifK38vWTPqXtpWvJbPkd3VNaL+L26LR4rEtLkNRTfEhVI7RH0\nH+E9WUSQGZkMpAxsyW/bXo2vN1fDHf0xwCcwReoe9K/J9OjnkD/35LjCb5GyDh7Er/2WQuSvGS1K\nt65vio+i6ZzooTDNKLdEEZHowBsZePgs1B+NCnr3tX7mELb0eTH0nxm5HekbxtfISJt8+L++BjaH\n/9NBxzPoQ+t6vbXZtJaRnov+vT/eEkHM2tpArzFOqeu2U6+3QEO6Pedw24vI3fqIoDLuG9J2xWwx\nJEgeKWOzHW9j3BKe7dV4e+QV1XGQuZuOoN4+6okMon2efinvT8kqBKfJD9AXX6UxvX8tj+Na+i/f\nI2btWFSrSXdkbaFlDUDqj+jqsQ7AXWs5fOu6BenX13PWCKwogOvnEPyo7WCUU67ahklbOrQlguzT\nDbI/IvnNDcI/viibmyiAov9hfJ3G1z1ovY6g5QdmrH7Pj7J4kLI3SvDq6yGzhI77SD12CLVGBr2j\nAs+awhzE33sdQOtr3vGTmfa6zYHgsc/H0zaK4lPK2JC8/gTJJ2C7peifXNN+pjTClVhGoI3W5Qjq\nQm4U5XO1f2gft7PH2ldPn6e0z98qBCfxWecAIuWiKWnppQe4P4fGljwgVvq9EQRHS5SH0Hit/f8e\nHhoNaOcV5paDqPvos2bba2Re81Qv07b6XNo3ecKbNhmpIHq1bZi0bbfHF+wGpW/HqD2lKhogSD4d\nXv5NVec+IR9uNd9c4+a6lul7jH9djuDx6nP0UFhrv6c+j2SwoyUfWktEaH/zpUpL9xxvLvU4PNPr\nUJhXX8uhMm8JiVY+bsE3Kf2avGSwrT7x8FgWDo9lYLubtm+rL2LVnqr2Yrg3TNu+/XrStj3w3hj9\nzYD9r+3qbZsNkcNwY7SPZ8B2ZLr5xgkc/dqOXF8bOvrQuhwBV2IC4Pfg05IRNVklIeaUjLD4NN6I\nQ/GME9XBUd8I87zEIVcPeWRaHYMWTUQiiKWiAQvthwx+VvoMtA4A28y2bxjkLbVzyB3Q0bvZVuf8\nbxB6MeDl/1wtOYxR+xTFFz15wrMVZKTrXrQyR/CoKilBEpwTlE++vXOjAM3ZaHzSeBKvxm/JeWQj\nujRaY1poTmQQcRCnjgbmFJPjeC3jzumZ9GeljyxAOZA6Ur5Z9AQoMtcR+7g9Vzl0D+pnIoQN01Yt\nqmnoe4rYdxM5CcUf2+ZHApvOL+iqHMHmiWcw7IpfvQbobzLvqoThpK+CH5nZ7kENO/d7zxMeYaKR\nKEDit2Q0Oa98q945JM2pdyzLUYujaDXqHh091gbcfOTBs86AoMgtuR4h+oFEMvnQv1e6Gcnmm1QJ\n0ng3yxjBD2x7jdb3unT+MWqvFnRvHM1QGcxKx83ia902zblT1C9HAukQCWRWrkbt45y/FQmkSgev\n806vETx4/GkMu/GUcnW/eyeRxP68Y96QUaRg9APAsAEyE8tzEcQgWIBLJNCXbkMk0JQCYv5ILG8e\nI22Jl/IwuhLhqdFznfumfYCE0A9DC4gcKY+2SNbIfsMge0DhZxD6FG3zKL1uH7ddT9oo+t7n+/Oo\nLRIJxNYELmsEePDYNbhIYFeh/Uz6hl3tGx9NjPgwpEqGiTIA5BFPaZQiBvK2Z0D8s+wEy8A5moms\n0c/Nz0NrNPSt1OQgsv0WaUb9RofgaSnKtsZl+XnezYR3yredoPxB7U/kPqbyuaFvjNjH+92rSIBB\n6+N2MjemfTvKlR934RyNuK5Xbzui76MjoFHFsb1eA+B4xtdld9Bw8/nII8/v2H+HI4LHHn80us6H\nZ/AAwDBs2T4AGIZppDCQ7yyV3+ugMrxlyURXzglZMOaU99huWC1r7ceSn/DPxAzR8eZSjyJa0QM2\n1piGv06SfJoi5j2N0e5RDzM0p3uTJ/6B8iVm7AnPpuYfRwF07pz+m74RKh6nVsZpEUnGas9saqfm\ntXTbbWOdieGLRQYDkyrajeZqrzNk4iDyjXOjawy9aF2OYPsURl+PLbC7gX3XoMZ+N0HcNFog0QAD\ns4YhVQ6hOCLeiF5fc5ZBNri769qY1gtstpHeTcbiToi1G/uBvZeaHhn9p6G0HQSj6iGfYzge6Inr\n2V5pzmQql1ImOXGZl5/XmG/DPh8D5U+QtryN7go80t/LCVHACB2P+6SxrpqjAE63rtcTGdRRxRGN\n2/rlSEGep6XDuu5Bq3IEj1fGp/wZHuARMoFmN33bcj3uvwm8tuV6im4Hgpgzt8BM+K4O1VHzIBvR\nmv/BY8f70WT244/TXg+YSqyDocOj94a4Sq+3gFqO1XMoPKo3KTokWUmG4090IVfiY3J7dOGQokXa\nT7cejhdKrb7M9lHkXhvROq3ia58uiG5GEcZ0vrUOf9tu0rYhY08XfDldx0hhc/M/r+coUw6Padcl\nCqgXmPNhH8Bdjgjw9Oh6NzHgJCKYxO5SP4P0D13DzeMW9FRTmM5nOmZNOWvG2wdzp5GBRP235Yyj\npZmUOGQ6l/xOgUfOMX36/HnZLbsow/NesS+3hsK9PP4IYMPCqj1xeXB+zEgEQPmt9rFzoKh9gx2D\n5Dk+rU1G7xvHWgCNApaJCu7w9tHHD7X06mgAKD7wSMd+oM6c0f7jr1oy6wOol2hKG2906Xw03pof\nCZMFTfqaqfIArpRUPXdfUaKR0ampR+EsDiXzfNMc+2Q+CtLStuxJCI2T8fK2oH9ez3kiAG2s6eGq\ncXsarQPokYGG2qdtPHLnZMt99Y4K5HHG9+yV70GrcgS9I4Lh5hE/EnkA4Jo1zrzJ4KMCmR/Yf/V0\nB+A3xuN7juTxhTGWrPAZpnanwCPl+Bi6HmWtoDPyl+cyB/1P+1vWALS+/pHBudYGPFGBrEdG8MWI\nSzuEfNeXiADjiGB/PbVmNN1TMmvjtj0eGf92Pb9WUHgfEOOrRQZc1DGeo98K0yjHQz0ihlOTF91T\n8h6wsXZbtKB+KdLwRwmZ5bXQPx3X0qGhf9o/1n1Eu1RuQ4xztdFbkIlEAJmM1RIBcG08cufWFHRd\nO1YPzePTdYCjDHc9jK43Cv8djgiembwiAzFmVn8rDyCjfck0DMrjs8wSN6cWPVO9tcO7jbWo27bG\neWIqz7F8SU9dRMwjk4T74Hi5eU22iTbwTPuvxX6KMKWtnXVfIg5i3D4w7XQ76FS/taVUO8hF2+lh\nsxZ9ezl526ikh2419RShG+ukh8/GB+rueERAU0PjzOM09bMBjLz9QFJDEuIfAmj/iPSnqRktQgDG\n9ySRNEeLIhHGWqnNCfhktoZb1fRIkcLctQE+OpjKTqOB8TW9Nx39y3l8GkuPDfc4YmhF/3JUUKPu\nYcJb32M0Kqi3hNSom0fqHjRf6xkmMhtWRrsemHuVI4A7HxEARxQ8Tv1M0zo1eVB/BPFLaF8zJRbK\nz9B/mMyjQ9J7F6hl2dpbfEvTTZGtV5YiUY2fQ3BcAtITRSyF/CWEr40pywws+t8nLL1RwbQkQ42M\n63ZvqQkLue95ekUEQ3XNIXzrmo8AaHTTg1blCJ7AlyeImiJoX2TA7ySaIn45OoiifQvpW5HCnsfH\n16J7zdSKbDzRg6X7iApjsnN2A/lQ/5SHzpWOZ60Z+KMC2sdHBRJCr/XxO4V49M/r9kQP3rbxOO0R\nwViPVTr6OFa9hpBRO4r62tolRKOJHrQqR7A31DHkz6P527km0PpnbYkizk/tx+Q9kYOF8jU9vdYD\nOP6l1gMoz6nWA+hY0jgcyq/bTxERlPTLmMcuPx2JInz5f891re8YRdHooRetyhHQNQKK1XYE+Zc/\ny5SHhuw86r8KlG3IAAZ18VU2T3tZH2q3d7rzNN1Kuz46voxROb+LHB8U4kjX5dmy2cI/LbsAcAWm\nqL5pVDDdDxfp54q0TfvG30JfKQm7nRZdo+3+AnHH9vH9DqO2aMG58dbRPOKxi81Rvdx1Hl0nhl8+\njDZUvLF3wkOrcgRP4Kmbz73WBQAe+UdRv/bYM8ZLZ1H56RzaEf7a1gva3BoQRT3WOFZOdal1AI63\nD+Lfgf61z7MWwKNqzw4hqd2D1m3ZebuBSptWaE7Sa+1Wml5r0YS1ltCHVuYIjhGB55CXt+QDj5b5\nx7hrRP38fP3y03kske7p9dVZxt3EtsTZcxiXFvDL8+hdlvFGBdbhLoA7bMbp4RE6EI0IxusAR0M/\nlrELwo375DLSvoNfRx3SePVYtk7+kJiHRy8D4YsGxg77yoX++eiD19WHVuUIyq4hYLq7xhMh7PmY\nAnMMX1mumbZjMo6mu5C1G8ijg6O1IfyeFHFLkV0Seu5fdwARmanrb9v108KTmO0J40XVsbxUrrnW\n648U6iJtFEHz7UXOXh/IkzZavrmWj5eeHv8OgKSf0z2ODsbztEpM29d7B3AsK4HRdVJke9DKHMFT\n1RWX+59+9SnJufIIL89vy+iyNfkiiHb966B5bsyHeHyRQVQ2slawXEQw7eeKpclzkRD5tG9J1J+Q\nR7uAbB1SNKCPdYwAY2sDclShyejXx9QYH03IPzhT39sRlo7LS0yfxVxalSOoU0PSFtCaeB5M2vbt\nHH6TF2e18gzWcuReXjf0pQq0R5eHWg+hzaGeB1q8JSI8Y2rF5bQSFvyWUHmbnlXcjZ/LVN9026de\nbM4qNNdSZM5bYK69iJzEO95rL43n3y46Li3dq6icpVc77LW/PqL58frJuByFvFV0mOj2vjMeWpUj\neNYoIpizTjDl2/PGSkIfi9ZppPe3If/bgvi91PaF7RERyLuIomsEPD+/dbR/7t/q58ov8+Nrawp2\nTl5qp8bNGtf/s5FtBeV86wK7CU/0h2YSBtA8/5Zc6+sA06jGjkzo3ObTqhxBvUYAAFfMKyyXkaB8\nPEKWzYJs9H0RAO+Q6t9/mYv+13xwbG504EE33oJ0kq5oielIkTlPNNCjrHQE7XNF1Hi5HSsjlZqW\nitNpCH3/uRfq18cpqDkZci0RxVTGLhxX90uRANXP6dpUv7+7yRlpd0cdwbPxpGuHj4yyuSggVjba\nXgPQjXlGS82f/hHArjFw3LBoeC7F9dloZ6m1gTnlpD1rCH0jgX1/PPffFgksmf+n8j0jAe9agqyX\nQ/6zI4Fc3dswjK63O3qdb6a3ne4cnk2rcgSPkTLUhSi6f4zhOfL6IgaNH/AWh/P8yMw8FL/mKGBP\nc4+6245nblG5XqUivJGAZ+3AKicxvhe+3ISkX+8bRJkWxO9ZE0BlBOs9+VypBBn1T8teTNG0N98/\nqHpqGaqX6uCvM9IBvad8QP4FzecDmh8OjmvYG/h0uLXNAJRfLU0ANntkub8uvqo8hn4BwbocweOj\nxeIpUXQvmY+98eUXhyXStnV6zdySZSSi5N2m6inF0IN8MY//EJmkr6XU9H4rpuRExtTv0Ji11TSr\nOrQidpqc50BX3U6fjTS+b+tm3c6ViuAK0k0PUPHbSKeHz6j+lrLUk8NiO9J/Ta8Hcj2e53ZXjZcB\n1Nf0CBMN9OrrzhXmTUeQUnozgNcA+GzO+eWk74cA/FUAL8w5/w4j+2oA/y32VvnNOef/WhvrcbJY\nDEwR8dKRwN7p2gvEnlIQPcpCrysi6FPoqkexuEItUYAX2cd4x5ENv2ago3WM+vXdRbH1gNo4TVE0\nLzNF5nJ0cEyjjMs6zMn9TyMXjo+OIe8EGpehTqhQOoBtzqPrzZCBoVovYZB7fb0tSP1wzaJ3muEb\nqr5W3k7kiQjeAuCnALy1bkwpvQTAnwTwMU4opbQB8DcAvArApwG8J6X0D3LOH5EGehaeYu/Pe5gM\niJWTOOpvLyZ31NHvh2Y+9K7fxjc+fIFz5PWRd8XDQu8fetfv4Bsf/iumXqkMhCannWL2oPspn+cw\n2ZRn2s8j7fe96wv45ofPFvW3RQF8YbrxKeNYNKHr5g6U6VFAaavH+Mfvusa3Pbxio4AJYt/tjsa5\nzM2J2m9kqpx8AqaGmi570K/jHKR/TXgHhXcGmY4g5/zulNJLma6/DuCHAbxdEH0lgF/LOX8MAFJK\nPwfgewCIjoBGBEtFA1OZ424l784eb1npUv5K5pjSr77rc/jmh88ReXufGeARbmTOPK9cqsEzPvCR\nd30O3/LwWSN92p4vWb8/AgB8Of49n57nB/RSzXQsrf9D7/odfPvDDenjkf543PG8ort/6nMUFjqP\n5/2FKCBXSH3IQM5I1cLp//XOHb7zjx//oumA3FP1CDbDAa3n/cQ2Q2XEb+SgG3QOqXPXGj/VX69B\nZ+aak6V6qUPoQE1rBCml7wbwiZzzB1MSDdLXAPhEdf1J7J2DSE8wqSHAV3q60Lkigr2ePuUjrnAt\nPou1UMQNRY/CF9375zCuSGuh/1r+KBP74RlPBCCtrVgIn/J4c/3lO1HL6gidOjNdxpPzL+1SKQoL\noZf2Oq9+bJtGShSp7/UBDx4BT3y5cu6Uj35FphuVYqjdwx+Rp1FAhJf2daKwI0gpPQvAj2GfFrpp\n7jGZ2vgd9xYcST8QNp2Cfyvo+PxCrKSz/9a9em+DIwD8iH9P/m9tSRU8wCM8G0+G9PBbRmX4NOXn\no47pvXJOwObRykCP+8eyV7jGs/BkhdqzIDeVtQ+G0WjiuOOlUMoZm7xH57WabcmrD5QX45w5qrx6\nNU2aWwcEpF5NPz0CNl9ipm7IsberofZyTVG/l5+ifk/EoEUfVE9HSjnbGg+poV/MOb88pfSNAN4J\n4Ens/6YvAfApAK/MOX+ukvlWAG/IOb/6cP2jALK0YJxS6nxrF7rQhS509ynnPBuIeyOCG8idc/4Q\ngBfddKT0GwBekXP+XSLzHgBff3AinwHwWgDfJw3Q42YudKELXehCcTJzFSmltwH4xwBellL6eErp\nBwjLYTkGSCm9OKX0vwFAznkH4HUA3gHgnwH4uZzzh3tO/kIXutCFLjSfXKmhC13oQhe60N2lxU8r\npZTenFL6bErpA0L/V6aU3p5Sen9K6YMppf+46nt1SukjKaWPppR+ZOm5LkmtzyGl9JKU0i+nlP7Z\nof2/OOnEF6A534lD/yal9N6UkrR1+dbQzPfjeSml/yWl9OHD9+OPn2ziC9DMZ/EXU0ofSil9IKX0\nsymlxzgdt4Ucz+L5KaVfSCn9Skrpn6SUvqHqi9vNnPOi/wB8O4BvAvABof+/BPBXDp9fCOC3sV+7\n2AD4fwC8FPsfAHs/gD+89HxX+BxeBOCbDu3PBfCrt/k5zHkWVf9fBPA/AXj7ue/lnM8CwH8P4AcO\nn68AfOW57+cczwLA7wfw6wAeO/T9HQDff+77WfhZ/DcA/qvD5z8E4J2Hz012c/GIIOf8bgB0IXnE\nAuArDp+/AsBv55yvUR1Iyzk/AlAOpN1Kan0OOeffyjm//6DjiwA+jP0ZjVtLM74T5UT7dwH4mUUn\neSJqfRYppa8E8G/nnN9y0HOdc/7/lp3tsjTnewFgC+A5KaUrAM/GvprBrSXHs/gGAL984P1VAH8g\npfRVaLSbayhk8zcAfENK6dMAfgXAXzi0cwfSbrUBNEh6DjeUUvoD2KOE//OkMzs9ac+inGi/L4tb\n0rP41wF8PqX0lkOa7KcPZ3zuMrHPIuf8aQB/DcDHsd/K/i9zzu882yxPQ78C4HsBIKX0SgBfi/1W\n/ia7uQZH8KcAvC/n/PsB/FsA/mZK6blnntM5SH0Oh88/D+AvHCKDu0zss0gp/QfYFz98P6RThHeP\npO/FFYBXAPibOedXYH+u50fPN82TkPS9eD72qPel2KeJnptS+tNnnOcp6CcB/L6U0nsB/GcA3ofp\nGWU3rcER/ACAXwCAnPO/APAbAP4w9p79ayu+cnDtrpL0HHAId38ewP+Yc/4HZ5vh6Uh6Ft8G4LtT\nSr8O4H8G8B0ppbeKWu4GSc/ik9iXefmnB76fx94x3GWSnsW/D+DXc86/k/fb1n8BwJ842yxPQDnn\nL+Sc/2zO+RU55/8IwL+K/TpJk908lSPQ0NvHsP9DIqX01QBehv0N3RxIO+wAeC3kAne3hVqeAwD8\nbQD/POf83y0+w9NR+FnknH8s5/y1Oeevw/778Ms55+8/yWyXpZZn8VkAn0gpvezA9yoA/3zpiZ6A\nWt6RjwP41pTSE2lf/OxV2K+l3XYSn8Vhx9iDw+c/D+AfHTIFbXbzBKvfb8N+4eZp7P9gPwDgpR26\npAAAANFJREFUBwH8J4f+FwP43wF84PDv+yrZV2O/S+bXAPzouVfyz/EcsEfBO+xX/98H4L0AXn3u\n+znXd6LS8e/ibuwamvN+/NHDi/9+7FHw8859P2d8Fq/H3vh/AMD/AODBue9n4WfxrQfb+GHso8Hn\nVbJhu3k5UHahC13oQvec1rBGcKELXehCFzojXRzBhS50oQvdc7o4ggtd6EIXuud0cQQXutCFLnTP\n6eIILnShC13ontPFEVzoQhe60D2niyO40IUudKF7ThdHcKELXehC95z+f4g9QxUVt5HjAAAAAElF\nTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f3d180f22e8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "k1, k2 = (1.8,1.9)\n",
    "A1, A2 = (14,16)\n",
    "kvalues = np.linspace(k1,k2,1000)\n",
    "Avalues = np.linspace(A1,A2,1000)\n",
    "kvalues, Avalues = np.meshgrid(kvalues,Avalues)\n",
    "ERR = np.sum([(f(X[i],kvalues,Avalues) - Y[i])**2 for i in range(0,len(X))],axis=0)\n",
    "ERR = ERR / ERR.max()\n",
    "plt.imshow(ERR,extent=[k1,k2,A1,A2],aspect='auto')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "### 2. A tour of scikit-learn\n",
    "\n",
    "[Machine learning](https://en.wikipedia.org/wiki/Machine_learning) is a *very large* field study and so we won't (and can't) go through everything. Instead, we'll start off with examples notebooks written in-part by [Gael Varoquaux](https://github.com/GaelVaroquaux) for various introductory courses taught at SciPy conferences in the past few years.\n",
    "\n",
    "Before getting into these example notebooks\n",
    ", take a look at the [scikit-learn](http://scikit-learn.org/) documentation and especially the [tutorial](http://scikit-learn.org/stable/tutorial/basic/tutorial.html).\n",
    "\n",
    "To get started, clone the [GitHub repo](https://github.com/GaelVaroquaux/sklearn_ensae_course):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Cloning into 'sklearn_ensae_course'...\n",
      "remote: Counting objects: 903, done.\u001b[K\n",
      "remote: Compressing objects: 100% (311/311), done.\u001b[K\n",
      "remote: Total 903 (delta 588), reused 903 (delta 588), pack-reused 0\u001b[K\n",
      "Receiving objects: 100% (903/903), 16.50 MiB | 4.56 MiB/s, done.\n",
      "Resolving deltas: 100% (588/588), done.\n",
      "Checking connectivity... done.\n"
     ]
    }
   ],
   "source": [
    "! git clone https://github.com/GaelVaroquaux/sklearn_ensae_course.git"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "1. In the Jupyter main page, navigate to the directory `sklearn_ensae_course` that was cloned above.\n",
    "2. Go to the directory `notebooks` and start with `00_introduction_machine_learning_and_data.ipynb`\n",
    "    * You may get an error about the version of the Jupyter notebook -- click OK\n",
    "    * Follow the notebook and execute each cell as you go\n",
    "    * In the second code cell, you'll get an error: `No module named 'urllib2'` and the error message points to an import in the file `sklearn_ensae_course/notebooks/figures/sdss_filters.py`\n",
    "        * Open the file `sklearn_ensae_course/notebooks/figures/sdss_filters.py` and replace `urllib2` with `requests` on lines 9, 30, and 47. (`requests` is the new HTTP package for Python 3.)\n",
    "    * The second cell should work but now we get several warnings: `DeprecationWarning: Passing 1d arrays as data ...`\n",
    "        * Open the file `sklearn_ensae_course/notebooks/figures/sgd_separator.py` and edit line 25:\n",
    "```python\n",
    "p = clf.decision_function(np.array([x1, x2]).reshape(1,2))\n",
    "```\n",
    "        * We need to restart the kernel since we changed something that was already imported into the notebook\n",
    "        * After restarting the kernel (and re-executing the first cell), the notebook should run smoothly ...\n",
    "    * Be sure to change the values of `x_index` and `y_index` in the last cell to plot the different features of the dataset\n",
    "3. Go to the directory `notebooks` and open `03_supervised_classification.ipynb`\n",
    "    * Follow the notebook and execute each cell as you"
   ]
  }
 ],
 "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.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}