JKJanosko commited on
Commit
2829e3a
1 Parent(s): dfffde9

Add files via upload

Browse files
Assignment2/index.ipynb ADDED
@@ -0,0 +1,333 @@
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
1
+ {
2
+ "cells": [
3
+ {
4
+ "cell_type": "markdown",
5
+ "metadata": {
6
+ "id": "Mz7n1DKlNfDZ"
7
+ },
8
+ "source": [
9
+ "# Bike Rides and the Poisson Model\n",
10
+ "\n",
11
+ "To help the urban planners, you are called to model the daily bike rides in NYC using [this dataset](https://gist.github.com/sachinsdate/c17931a3f000492c1c42cf78bf4ce9fe/archive/7a5131d3f02575668b3c7e8c146b6a285acd2cd7.zip). The dataset contains date, day of the week, high and low temp, precipitation and bike ride couunts as columns. \n",
12
+ "\n"
13
+ ]
14
+ },
15
+ {
16
+ "cell_type": "markdown",
17
+ "metadata": {
18
+ "id": "fB-BnQuhNfDc"
19
+ },
20
+ "source": [
21
+ "## Maximum Likelihood I \n",
22
+ " \n",
23
+ "The obvious choice in distributions is the [Poisson distribution](https://en.wikipedia.org/wiki/Poisson_distribution) which depends only on one parameter, λ, which is the average number of occurrences per interval. We want to estimate this parameter using Maximum Likelihood Estimation.\n",
24
+ "\n",
25
+ "Implement a Gradient Descent algorithm from scratch that will estimate the Poisson distribution according to the Maximum Likelihood criterion. Plot the estimated mean vs iterations to showcase convergence towards the true mean. \n",
26
+ "\n",
27
+ "References: \n",
28
+ "\n",
29
+ "1. [This blog post](https://towardsdatascience.com/the-poisson-process-everything-you-need-to-know-322aa0ab9e9a). \n",
30
+ "\n",
31
+ "2. [This blog post](https://towardsdatascience.com/understanding-maximum-likelihood-estimation-fa495a03017a) and note the negative log likelihood function. \n"
32
+ ]
33
+ },
34
+ {
35
+ "cell_type": "code",
36
+ "execution_count": 285,
37
+ "metadata": {
38
+ "colab": {
39
+ "base_uri": "https://localhost:8080/",
40
+ "height": 766
41
+ },
42
+ "id": "LAz1foK9NfDd",
43
+ "outputId": "1034c2ad-6856-46d3-bfc5-7915646979b4"
44
+ },
45
+ "outputs": [
46
+ {
47
+ "output_type": "stream",
48
+ "name": "stdout",
49
+ "text": [
50
+ "lambda* from gradient descent : 2680.0420560747625\n",
51
+ "Log Likelyhood for Lambda* : 3952001.7180423485\n"
52
+ ]
53
+ },
54
+ {
55
+ "output_type": "display_data",
56
+ "data": {
57
+ "text/plain": [
58
+ "<Figure size 1440x720 with 4 Axes>"
59
+ ],
60
+ "image/png": "iVBORw0KGgoAAAANSUhEUgAABZcAAALICAYAAAADj84wAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAAsTAAALEwEAmpwYAACM7klEQVR4nOzdfXxkZXnw8d9FiBIVDMhqSWBdVIxFsLt1pbZUq/gS6uuKLYpW0VK3tlq1+qRl1QpifaDG90erhUrBqiiVNaJAIyq+C7qYlQUhCohKFmErhBeNdAnX88c5WSYhyc4smTlJ5vf9fOazM/c5Z851ZjP3mbnmPtcdmYkkSZIkSZIkSY3YreoAJEmSJEmSJElLj8llSZIkSZIkSVLDTC5LkiRJkiRJkhpmclmSJEmSJEmS1DCTy5IkSZIkSZKkhplcliRJkiRJkiQ1zOSyFpWIODEiPnEftr8jIh5R3j8jIv554aKT1E4i4qUR8aUmPXfL+qeIeEVEfKsV+1qqIuKCiDi26jh2VaP/x/f1XCs1U0RkRDyqznVXlevv3uy4FqOIeFJEjNY8vi4inr4Lz7Nju4h4c0T8e3m/Za9vRHwtIv6q2fvZSQw7+saIWFl+r+ioMiZJS8di6Mca1czvOw3E8NGI+KcqY9B9Z3JZ05Qd4i0Rcf86129Z0iIinhIRd5cf9O6IiLGIeHvtOpn5oMy8tgWxfC0ifhsRt0fEbRFxaUQcX+/rVj5H3V+e7otW7UdaiiLijyPiOxFxa0TcHBHfjognAGTmJzPzmYsgxqZ9UK1JHEz1qzdGxBcj4hkNPEdLzgPN3E9m/mlmnrnQz1uet65f6OdttfLv5MSq41C1djVxudy14nPWXK99Zn4zM/sWcl+Z+X8zc9ElRyLioIj4dERsKz97/yQi/l9E7N+M/WXmz8vvFZP39bl2dh5fiHNxqzmIR4uZ56v6zfy+0+xz2myf5zPz1Zn5jmbtU61hclk7RMQq4ElAAs+rNpo5bS0/6D0I+GPguIhY18wdzjNi4bWZuSewH/Am4MXA+RERzYxH0sKIiL2ALwL/D9gH6AXeDtxZZVwV6S771d8DLgQ+FxGvqDYkVS0inhgRbwF2Lx8/uXwsSS1TJjouAbYCazJzL+Bw4BqK7wOzbbMUR7N7LpbUNM3uF5dov6sFYnJZtV4OXAycAUy7PDgiDoiIjeVogV9FxIci4neBjwJ/WP7KPl6uO+3X+Zm/TkXEByLiFzUjfp+0K8Fm5k+B7wAH1zz3rL+0RcSeEXFRRHwwCo+JiAvLkYqjEXF0zbpnRMRHIuL8iPg18NSdxPHrzPwaRUL+D4Fnl89zWER8NyLGI+KG8jW7X7nsG+XmPyxfuxdFxN7lKIVtUYwe/2LtaIzydbw2itHSP42Il9Ys+8uIuLLcbjgiHj7Xfhp4iaXl7tEAmXlWZk5m5kRmfikzL4NZ+66MiL8tR0vdHhHviIhHRjHy+baIOLvmPX6vX+Xn6Z/mfO9HxDspfvT7UPke/lDZPl8f9pCIOLeM6XvAI+t9QTLzl5n5AeBE4F8iYrfyOY+PiGvK4/5RRLygbJ/rPPDsiBgpY/hF1Ix6jYg9IuIT5blkPCK+HxEPK5c9OCI+VvaZYxHxzxHRMdd+Znktp41UiemXOM+33x3nran/u4h4d/n/8dOI+NOa5zwwIr5RvhZfjogPxy6UmNjJazQ1iu2V5bJbIuLVEfGEiLisjP9D937K+FAUo/CvioinzYj562XMFwL7ztjwvyLil+W234iIxwJk5sXA5cBHKH5A/VPgA40eq5a3+fqwcvnXyvfyd8r37xfKfuqT5d//96MY4FDrWVF85vmfiBis6Ys6yvfm/0TEtZSfuWr29cooPg/dXm7/13PEfP/yfXRITduKiJiIiIdGxL7lcYxH0c9+cyqGBl6XB0fEx8vX5WcR8dYZx/Ge8jh+GhGvjV0oPxHzXB0REb9bPvcx5ePnRMTm8pi+ExGPm2O7Hf1mjZdGxM/LeN9Ss+79I+L9EbG1vL0/aq7ii4hXRcTV5Wt4bkT01Cx7RtlX3Vr2Z/MNzjgR+HZmvjEzrwfIzJsy8/2Z+ena1yIi/jEifgn8Rx1/m3P2jTGjJEjMcX4ql8153og5zuPzmedc3BMR55TH89OIeF1NvIdFxKbyPXVjRLy3ZtnUVVrjUZxTXlHz//fu8v/2xiguTe+a8Xq+KSJuKo/7leWy9cBLgX8oj+kLOzsmaTGoo09o6Hw1Xz8WxXeEr0bxufN/yufonie2Pyqf/9by3z8q218UEZtmrPv3EXFueb+e9/GOfnGW/e74zhJz5A5invNHFJ+9/zEiLgN+HRG7R+PfG86ImishYv5zR0bxmfgnZTwfjigG90XEo6Lo028tX/PPzPV6a+GZXFatlwOfLG/9cc+X7g6K0X0/A1ZRjO77dGZeCbwa+G45mri7zv18H1hNMVLwU8B/RcQejQYbEQdRjFq4eCfrPQT4CsWH0tcBD6AYDfAp4KEUX5j/NSIOrtnsJcA7gT2Bui7DzsyfA5soPkACTAJ/T/FB9Q+BpwF/W6775HKd3ytfu89QvB//A3g4sBKYAKYSSQ8EPgj8aTla+o+AzeWy5wNvBo4CVgDfBM6aZz+SCj8GJiPizIj404jYu45t+oHHA08E/gE4FfgL4ADgEOCYXYhjzvd+Zr6F4j392vI9/NqyP5ivD/sw8FuKqyr+srw1amP53FOXW19D0bc9mGJ09yciYr95zgO/pjindFMkf/4m7rnK5NjyeQ4AHlJuP1EuOwO4C3gUsAZ4JvBX9+F8U2u+/c70B8AoRf/9LuBjUx9cKV7375XPcSLwsl2IBeZ/jWrjOAh4EfB+4C3A04HHAkdHxJ/MWPeaMuYTgI0RsU9NzJeWy97BjB+QgQvK/TwU+AHF54ApWXN/csZjCebpw2q8mOK90kvxg9d3y232Aa6k+Jut9QJgLfD7wPO5px97FfAciv5hLfBnM7a7qVy+F/BK4H0R8fszA87MOyn6udo++2jg65l5E8UVaddTfK56GMXnrEb/9v8fRZ/zCOBPKN7vr6w5jj+l+Dz8+8C6Bp97XuUxDwN/l5lnRcQa4HTgryn6rn8Dzo36y7n9McX54GnA28oEARR90hMpjuP3gMOAt5YxHAGcTPG67kfxPWIqEbwvxev/Vop+6RqKz/RzeTpwTh1x/g7F39TDgfXs/G9zZ31jrTOY5fxUs3zW88Zs5/E6jmPKjnNxFAnmLwA/pHgfPQ14Q0T0l+t+APhAOar7kcDZAFEMOLmA4u9xBcX/1eZym1MofmhfXR5XL/C2mv3/DsXfcC9wHPDhiNg7M0+lOE+8qzym5zZwTFKVFux8VUc/FhR9YA/wuxSfP0+cLajy89p5FN/3HwK8FzivzGN8gaIPOKhmk5dQ9F9Q3/u4tl+c02y5gzrPH8dQfJbtzsy7aPx7Q+1rMee5o8ZzgCcAjyvXm+oH3wF8Cdgb2J+i31OLmFwWUPyiTdHhnJ2Zl1J0CC8pFx9G0SkOlKN0f5uZu1z3MjM/kZm/ysy7MvM9wP25J4GxMz3lL1S3USSGLmH+5G8P8HXgvzLzrWXbc4DrMvM/yhhGKD6w/nnNdp/PzG9n5t2Z+dsGDm8rRedNZl6amReX+7iOoiP+k7k2LF+TczLzN5l5O0Vyu3b9u4FDIqIrM2/IzCvK9lcDJ2fmlWVn/n+B1eWHSUlzyMzbKL4wJ3AasK38dfxh82z2rsy8rXz/XQ58KTOvzcxbKb68rdmFOHb23p9pzj6s/DHwhcDbyv76cmBXaglvLf+d6s/+KzO3ln3iZ4CfUJwb5jqmr2XmlnL9yyh+8Jo6pu0UH04fVY4YvzQzbytf92cBbyhjvwl4H8WH/IUw637nWPdnmXlaFrU2z6T4cPuwiFhJ8WH2bZn5v+W58NxdCWYnr9GUd5Tn3C9RJKPPKkfrjVEkK2r/3m4C3p+Z28v/o1Hg2TUx/1Nm3pmZ36D4olIby+mZeXuZcDsR+L0oRuk9keKD+99QfLC/EHj9rhyvqhURp0cx+vDyOtc/uhxtdEVEfGq+devsw/4jM6+p6Suvycwvl59b/ot7953/kpk3lz/cv597ksBHU/yd/yIzb6b4Aloby3nlfjIzv07xJXOuK+Q+xfT+pfbL+naK9/3Dy/fUNzOz7uRy2Re/GNhQvreuA97DPT9GHU2RCLw+M2+hSA4slCdR9Esvz8wvlm3rgX/LzEvK/u9MihJQT6zzOd+exdU9P6RIbv5e2f5S4KSyX9pGkUR4Wc2y0zPzB2XfsoFitNoqir7+isz8bGZup/g//uU8+9+3dnkUI73Hy5Fvp9WsdzdwQtnXTcz3t1lP31izv3rOT7OeN+Y5pnrUnoufAKzIzJPK88+1FJ9dpmLYDjwqIvbNzDuyuPIEir/rL2dxldb28jXZXP5guh74+/K9djvFd4jaY9pO8f+7PTPPB+6g/u9s0qKzwOerefuxzLw6My8s+5dtFAnjuT7fPxv4SWb+Z/nZ/izgKuC5mfkb4POU58EyyfwYigRvPe/jaf1i469aXeePD5bn5Yny2Bv63jDDfOeOKadk5nj5GeEiisQ6FH3Ww4GevI85KzXO5LKmHEuRJPmf8vGnuOfX+wMoPjDdtRA7ioj/E8Uli7eWl0I8mBmX6M5ja2Z2Z/GrfDfFr43zJU6eDXRRXH4x5eHAH5QfSsfLGF5K8avelF/UGc9MvcDNABHx6CgutfllmQz/v8xznBHxgIj4tygunbwN+AbQHREdmflripFrrwZuiIjzIuIxNcfzgZpjuZnil9LeXTwGqW1k8aPMKzJzf4qRxz0UHw7ncmPN/YlZHj+o0Rjme+/Pscl8fdgKivq4tX3YzxqNiXv6j6n+7OVxz+Vw4xSv1Xz92R9EUYpoW0TcStF3Ta3/nxQj6j4dxWXU74qIzvK4Oin6uKn9/BvFqK2FMNd+Z1P75eA35d0HUfx93FzTBrt4vtjJazSlkb+3sRnJr5+V8fYAt5TnkdplU3F0RMQpUVy+eBtwXblo3yx+IP1nitF6ZObXM/P/NnywWgzOAI6sZ8Xyi+sG4PDMfCzwhp2sX08f1mjfObMPm7oktmeWZbWx/GlEXBzFpbTjFF/+5+qrLgIeUL4XV1F8Of1cuWwQuBr4UhTlNY6f4znmsi9Ff1Yb38+4p2+deRy7+rlzNq8GvpNFybYpDwfeNOO8cQD3vK47U5v4/Q33/H/1cO9j7JltWWbeAfyK4jWYdvxl3zXfa/ArimTt1PofymLE2/spXucp27JmUMhO/jbn7RtnqOf8NNd5476oPRc/nHsG2UzF8GbuSWAfRzF68aooLql/Ttl+AMWgoZlWUFzNeWnN8/132T7lVzO+/9X+30tLzgKfr+btxyLiYVFMQjpW7usTzH0+mtmXwvRzxqe450fWlwBDZT9Tz/t4Wr+4C+o5f0zrvxv93jDDfOeOKXOdk/6BIg/yvfLH8V25elO7yOSyiKImz9HAn5SJ0F9SlHP4vYj4PYrOYmXMXgdutlEcv6bo5KbsSNpGUV/5H8r97V1+MLyV+euszar8NfFTwHyXYp1G0cGeH8Wl5FAcz9fLJPXU7UGZ+Te1T99oPBFxAMXl8t8smz5C8YvjQWUy/M3Mf5xvohgN8Afl+lOXpQRAZg5n5jMoPlxfVR7b1PH89Yzj6crM7zR6DFI7y8yrKBIwh+xk1XpM6wcj4nfmWXfe9z737o/m68O2USQCD6hZf+UuxP8CipGwo1FcBXEa8FrgIWW/ffk88UHRN58LHJCZD6b4gW+qL9uemW/PzIMpSvw8h+Jy8V9QjITYt+a49iqTW3PtZ6Y5zz/z7LcRNwD7RETtPg6Ya+WdmPM12kW95QiWKSspRr3dAOxdcw6cWjblJRRlB55O8WPvqrJ9x3Nl5nWZeeJ9iE0Vy2JU5s21bVHUg/zvKOa/+GbNj9avAj6cxYhashilOZ+d9WG7YmYfNjWC84ZZlhU7Ky7RPQd4N/Cwsq86f644shhhejbFF/ZjgC+Wo77IYrTxmzLzERRzarwxauqY1+F/uGcEVW2sYzXHsX/Nsl3tR2bzaorP7e+rafsF8M4Z540HlKPj7out3PsYt862rOyDHkLxGkz7fyz7rvleg69QlH/bmZnnifn+NnfWN9ba2fmp0bjqteNcXMbw0xn/h3tm5rMAMvMnmXkMRcL7X4DPlsf2C2afe+F/KBJlj615vgdnMaFgM49JqtJCnq921o/9X4r3yaHlvv5inv3M7Eth+jnjQmBFRKymOF9NXWVTz/v4vr5X6zl/7NjHLn5vqDXfuWNeWdSrf1Vm9lCU8fjXmGW+GzWHyWVBUedtkmJivNXl7XcpkqQvp6gteQNwSkQ8MIpJkabqCd0I7B/lJFalzcBR5S+Dj6L4JX3KnhSJj23A7hHxNoq6eA2LiAdRXPJxxU5WfS3Fh7IvlIn0LwKPjoiXRURneXtC3FM/rtE4HhBF3cvPU7xW55eL9gRuA+4ov7D9zYxNb6Sow0fN+hPAeBR1l3bUHyx/+Xx+2bneSXFZ2t3l4o8CG6KcgCmKS5lrS3zM3I8kIIpJ8d4U90yedwDFB7Z567jX6YfAYyNidRQ15U+cZ9053/ulme/hOfuwMlmyETix7JsOZv4aktOUfc1ryxg2ZObdwAMpPghuK9d5JdMT8LOdB/akGOH724g4jHvKLBERT42IQ8tRIrdRJGDuzswbKC5hf09E7BURu5XJrz+ZZz8zbQZeXL4m0+qxzrXfel8bgMz8GUVt/RMj4n4R8YfM/wPn1L73mHEL5nmNdtFDgdeVx/7nFOfx82tifnsZ8x/PiHlPivPKrygS845Mbh+nUtTkfTzwf4B/LdsfTdHHfDuKUcC1I547Z/wt787O+7BdMRDFxEsHUJRimZoz4myKv/P9o6iTXzui+H4Upda2AXdFMaHaM3eyn09RXBn2Uu75sj41edGjyvfqrRSfk+frL+5X+7rUxPrOKCaVfjjwRoqRa1PLXh8RvVFM8PSPO4kTZn/tZ3M7xSj1J0fEVLmN04BXRzFKO8rP88+OiD3r2O98zgLeGsVkiPtS1Pn8RM2yV5bnwftT9C2XZFEi5DyKc+RR5XG8julXEM50IvCkiHhvRPTCjnqnO/vsPuffZh19IzXr7uz8tDMNfRaf41z8PeD2KCbO6oriqpNDIuIJ5TZ/ERErynXHy6e6m6I28tOjKHWzexSTk60u1zuNoi75Q8vn6I17ajgv6DFJFWj2+Wpn/dieFN/Zby37rYF5nut8ivPuS8r36YsocjNfhGKABEVJjkGKMjkXlu339X08m5nv7UbPH7vyvaHWfOeOeUXEn8c9EzTeUsbR0Gd97TqTy4Ii8fAfmfnz8teeX2bmLymK27+U4lem51IUiP85xQQnLyq3/SpFcveXETFVUuN9wP9SdBxnMn1ioGGKkcQ/prjc4bc0dilgTxT11e4ot9+njHFO5SUq68u4P0+RUHgmRWJ6K8VlFf9C8YWkER+KiNspjvP9FKNljiw7eSi+qL2E4kP+adzzxWjKicCZUVwucnT5HF0Uv0BeTPE6TdmN4kvJVoqRR39CmazOzM+V8X86iktuLqeYJGau/Ugq3E4xAc8lEfFrivfd5RSjGu6TzPwxcBLwZYo6Y/PV/Ho/c7/3oZik58+imNX6g+XIuvn6sNdSXB72S4qR2PeaGXoW4+VrsIXiMvI/z8zTy2P5EUWt0O9S9HeHAt+u2Xa288DfAieVfeTbKCcWKv0O8FmKBO+VFHXx/7Nc9nKKBNGPKD4UfpZ7LoWebT8z/RPFCK1bKGp/1taKnW+/jXgpxSStvwL+maJvv3Oe9XspvsjU3h7J/K/RrriEYlK+/6GoIfhnmfmrctlLKP7Wb6b4IvXxmu0+TnE+HaN43RfixxUtclH8QP9HFJMqb6a4xH/qvbY7xd/SUyh+cDst7pnh/nym/y2fyM77sF3xeYqJ1jZTfIH/WNl+GsVnyR9STD65cWqDsm98HcV76RaKv/t5a6Jn5iUUVzz0UNTWnHIQRf99B0Xf96+ZedE8T3UF01+XVwJ/Vz73tRTngE9RTIo0dRxfAi4DRihe17sokthzme21n+u4xoFnAH8aEe/IzE0UI9I/RPHaXA28Yp591eufKRK0l1GcP35QtpGZX6bok8+hGKTySMo6oFmU4ftzilrTv6J4vb/NHMpz6h9QjPb+YdlvfpviHPhP88T3fub/25yvb5xpvvPTzkw7j8+z3nzn4kmKK25WAz8tj+nfKa44geIHhSvK70gfAF6cRd3pn5fP9abyODdzT83sf6T4W7i4/A7xZeqvqfwx4ODy+8VQndtIrdTU81Ud/djbKSZsvZXiPLZx5nPUPNevKN7fbyqf6x+A5+Q9JUuhOIc8nWIuqdpyNfflfTybE6nJHTR6/tjF7w2128957qjDEyi+191Bcf5/fRb16dUCkfXPTSFJkqQaEfEZ4KrMXIjRmlLTRFFX+IuZeUhE7AWMZua9kmMR8VGKUUL/UT7+CnB8Zn6/pQG3kShGWX80M52MWZIkLTmOXJYkSapTFCVIHlleFn0kRb3ioYrDkhqSmbcBP42yjFZ5qevUaMYhilHLU6UHHk0x+lYLpCxr8Kzy8udeilGzn9vZdpIkSYuRyWVJkqT6/Q7wNYrL5T8I/E1mjlQakbQTEXEWxSWqfRFxfUQcR1Hi5biI+CHFJarPL1cfBn4VET8CLgIGakqsaGEExSXTt1CUxbiSojyOJEnSkmNZDEmSJEmSJElSwxy5LEmSJEmSJElq2O5VB1CPfffdN1etWlV1GJK0U5deeun/ZOaKquNoJftoSUuFfbQkLV720ZK0eM3XRy+J5PKqVavYtGlT1WFI0k5FxM+qjqHV7KMlLRX20ZK0eNlHS9LiNV8fbVkMSZIkSZIkSVLDTC5LkiRJkiRJkhpmclmSJEmSJEmS1DCTy5IkSZIkSZKkhplcliRJkiRJbSciDoiIiyLiRxFxRUS8vmz/TERsLm/XRcTmsn1VREzULPtozXM9PiK2RMTVEfHBiIiKDkuSWmr3qgOQJEmSJEmqwF3AmzLzBxGxJ3BpRFyYmS+aWiEi3gPcWrPNNZm5epbn+gjwKuAS4HzgSOCCpkUuSYuEI5clSZIkSVLbycwbMvMH5f3bgSuB3qnl5ejjo4Gz5nueiNgP2CszL87MBD4OrGtW3JK0mJhcliRJkiRJbS0iVgFrKEYeT3kScGNm/qSm7cCIGImIr0fEk8q2XuD6mnWupyZJPWM/6yNiU0Rs2rZt28IdgCRVxOSyJEmSJElqWxHxIOAc4A2ZeVvNomOYPmr5BmBlZq4B3gh8KiL2amRfmXlqZq7NzLUrVqy4r6FLUuWsuSxJkiRJktpSRHRSJJY/mZkba9p3B44CHj/Vlpl3AneW9y+NiGuARwNjwP41T7t/2SZJy54jlyVJkiRJUtspayp/DLgyM987Y/HTgasy8/qa9VdEREd5/xHAQcC1mXkDcFtEPLF8zpcDn2/JQUhSxUwuS5IkSZKkdnQ48DLgiIjYXN6eVS57MfeeyO/JwGURsRn4LPDqzLy5XPa3wL8DVwPXABc0O3hJWgwsiyFJkiRJktpOZn4LiDmWvWKWtnMoSmjMtv4m4JCFjE+SlgJHLkuSJEmSJEmSGmZyWZIkSZIkSZLUMJPLkiRJkiRJkqSGmVyWJEmSJEmSJDXM5LIkSZIkSZIkqWEmlyVJkiRJkiRJDTO5LEmSJEmSJElqmMllSZIkSZIkSVLDTC5LkiRJkiRJkhpmclmSJEmSJEmS1DCTy5IkSZIkSZKkhplcliRJkiRJkiQ1bPeqA5AkSZIkSYvP0MgYGzZexsT2u3dp++tOefYCR6QqvPS07/Lta26uOgxJC+jwR+7DJ1/1hwvyXCaXJakNRcQewDeA+1OcCz6bmSdExIHAp4GHAJcCL8vM/60uUkmSJDXD0MgYJ557BeMT25u2j1XHn2eCeQkZGhnjjZ/ZzK79lCBpKfn2NTfz0tO+uyAJZpPLktSe7gSOyMw7IqIT+FZEXAC8EXhfZn46Ij4KHAd8pMpAJUmStOuGRsYYHB5lbHyCALLqgLToODJZak8L9b43uSxJbSgzE7ijfNhZ3hI4AnhJ2X4mcCImlyVpwcx15ciMdR4OnA6sAG4G/iIzry+XTQJbylV/npnPa1XskpaO2oRyLRPLqvXWoS184uKfVx2GpCXO5LK0iK06/ryGt/GyM9UrIjooSl88CvgwcA0wnpl3latcD/TOse16YD3AypUrmx+s2s6u9H9gH6glYdYrRzLz4pp13g18PDPPjIgjgJOBl5XLJjJzdWtDlrTYOTpZjXrGe7/GT276ddVhSFoGTC5LUpvKzElgdUR0A58DHtPAtqcCpwKsXbvW7y+SVKd5rhypdTBFmSKAi4ChlgQnacmYL5nsBzPtzEtP+66JZUkc/sh9FuR5TC5LUpvLzPGIuAj4Q6A7InYvRy/vD4xVG50kLT8zrxzJzEtmrPJD4CjgA8ALgD0j4iGZ+Stgj4jYBNwFnJKZQ7M8v1eXSMvQXAnlxZxM9oqixeetQ1usryyJwx+5z4JM5gcmlyWpLUXECmB7mVjuAp4B/AvFCLk/Az4NHAt8vrooJWl5mnnlSEQckpmX16zyf4APRcQrKOozjwGT5bKHZ+ZYRDwC+GpEbMnMa2Y8v1eXSMvEYkko7xbwkj9YyT+vO7TFe9ZCGhoZu881lg966AO58I1PWZiAJC0LJpclqT3tB5xZjp7bDTg7M78YET8CPh0R/wyMAB+rMkhJWs5qrhw5Eri8pn0rxchlIuJBwAszc7xcNlb+e21EfA1YQ1EzX9Iy0aqE8m4Bdyf0dncx0N/HujWzTrWhZeSNZ29uaP2/eKI/KEjaOZPLktSGMvMyioTEzPZrgcNaH5EktYd5rhypXWdf4ObMvBvYAJxetu8N/CYz7yzXORx4V0sPQFJTNDOhbBJZUJTDuLuOP6jdAt579Gr/TiTVzeSyJEmS1DpzXTlyErApM88FngKcHBFJURbjNeW2vwv8W0TcXW57Smb+qOVHIGlBDY2MsWHjFia2F9Vv7mtC2WSyZlNPOQxLXkjaFSaXJUmSpBaZ58qRt9Xc/yzw2VnW+Q7g9cnSMlE7WnlXTI1w7ohgMtNksuY0NLLzObpNLEvaVSaXJUmSJElqsqlk8tbxCR7c1cmv//cutk82Nk55KqFsIlmN+MdzLpt3+cP2vJ+JZUm7zOSyJEmSJElNNLP0xfjE9rq3NaGs+2JoZIw777p73nUuecszWhSNpOXI5LIkSZIkSQusdqTybmXpinqZUNZCecvntsy7/C+euLJFkUharkwuS5IkSZK0gGaOVG4ksWxCWQvp1/87Oe/yf15nKX9J943JZUmSJEmSFtDg8OiOxHK9ujo7OPmoQ00qa8HsbCI/Ry1LWggmlyVJkiRJuo9qy2DUM065c7fgQXvszvhvttPjaGU1wdu/cMW8yx21LGkhmFyWJEmSJOk+mFkGYy4dEdydaTJZLXHLb+aeOLKrc7cWRiJpOTO5LEmSJElSA2pHKfd0d/Gb/71rp4lly15oMTn5qMdVHYKkZcLksiRJkiRJdZo5SnlsfGLe9QMcqayW21m9Zf8WJS0Uk8uSJEmSJNWpkcn6eru7+PbxRzQ5IunedlZvWZIWikV2JEmSJEmq09adjFSe0tXZwUB/X5OjkWY3X73l7q7OFkYiablz5LIkSZIkSXOYWV+5+wGdsybuurs6eeD9d9+xnmUwtFid+LzHVh2CpGXE5LIkSZIkSbOYrb5y525BZ0ewfTJ3rNfV2cGJz3usyWQtCtZbltRKlsWQJEmSJGkWs9VX3n538sD77U5vdxdBUVf55KMONWGnRWO+esuWxJC00By5LEmSJEnSLOaqr3zrxHY2n/DMFkcj1We+esuWxJC00EwuS5IkSZLa3szaygP9ffR0dzE2S4K5p7urggil+84R9pIWmmUxJEmSJEltbaq28tj4BElRW3nDxi089TEr6OrsmLZuV2cHA/191QQqSdIiY3JZkiRJktTWZqutPLF9kouu2sbJRx1qfWUtC1F1AJKWJctiSJIkSZLa2ly1lbeOT7BuTa/JZC0ZQyNjcy7LFsYhqX04clmSJEmS1NbmqqFsbWUtNW//whVzLuv171lSEzhyWZIkSZLUNmabuG+gv48NG7dMK41hbWUtRbf8Zvucy/x7ltQMjlyWJEmSJLWFuSbuA6ytrGXPv2dJzeDIZUmSJElSW5hr4r7B4VG+ffwRJt+05AWz11Z2Mj9JzdK0kcsRsUdEfC8ifhgRV0TE28v2MyLipxGxubytblYMkiRJkiRNmW/iPmk5mGvSPifzk9QszRy5fCdwRGbeERGdwLci4oJy2UBmfraJ+5YkSZIkaZqe7i7GZkkkO3GflovdAu6eJZPcEY5dltQcTRu5nIU7yoed5c0fyyRJkiRJlRjo76Ors2NamxP3abkYGhmbNbEMMJmmYyQ1R1Mn9IuIjojYDNwEXJiZl5SL3hkRl0XE+yLi/nNsuz4iNkXEpm3btjUzTEmSJEnSMjE0Msbhp3yVA48/j8NP+SpDI2M7lq1b0+vEfVq2BodH51zW6+h8SU3S1An9MnMSWB0R3cDnIuIQYAPwS+B+wKnAPwInzbLtqeVy1q5d609skiRJkqR5DY2MsWHjlh2T9o2NT7Bh4xaAHQnkdWt6TSZrWZqvdrij8yU1S1NHLk/JzHHgIuDIzLyhLJlxJ/AfwGGtiEGSJEmStLwNDo/uSCxPmdg+Oe+ITmm56H5A56ztD+jczR9UJDVN05LLEbGiHLFMRHQBzwCuioj9yrYA1gGXNysGSZIkSVL7mGvk5nwjOqXlYq6yyvfbvWP2BZK0AJpZFmM/4MyI6KBIYp+dmV+MiK9GxAoggM3Aq5sYgyRJkiSpTfR0dzE2SyK5x3qzagPjE9sbapekhdC05HJmXgasmaX9iGbtU5IkSZLUvgb6+6bVXAbo6uyw3qzaQkcEk7MMX+6IqCAaSe2iqRP6SZIkSZLUKlN1ZQeHR9k6PkFPdxcD/X3Wm1VbmC2xPF+7JC0Ek8uSJEmSpCVlaGRszgTyujW9JpPVdoZGxghgtjRyr2VhJDWRyWVJkiRJ0pIxNDI2rfTF2PgEGzZuATCprLY1ODw6a2I5wLIwkppqt6oDkCRJkiSpXoPDo9NqKgNMbJ9kcHi0ooik6s02kSUUI5n90UVSM5lcliRJkiQtGVvnSKLN1S61g7km7XMyP0nNZnJZkiRJkrRk9MxRP3audqkdOJmfpKqYXJYkSZIkLRkD/X10dXZMa+vq7LCurNra3g/obKhdkhaKE/pJkiRJkpaMqfqxg8OjbB2foKe7i4H+PuvKqq3NNUDZgcuSms3ksiRJkiRpSVm3ptdkslTj1ontDbVL0kKxLIYkSZIkSWpLEXFARFwUET+KiCsi4vVl+4kRMRYRm8vbs2q22RARV0fEaET017QfWbZdHRHHt/I4rEUuqSqOXJYkSZIkLRpDI2OWvFAr3QW8KTN/EBF7ApdGxIXlsvdl5rtrV46Ig4EXA48FeoAvR8Sjy8UfBp4BXA98PyLOzcwfteIgnvqYFXzi4p/P2i5JzWRyWZIkSZK0KAyNjLFh4xYmtk8CMDY+wYaNWwBMMKspMvMG4Iby/u0RcSUw3x/b84FPZ+adwE8j4mrgsHLZ1Zl5LUBEfLpctyXJ5Yuu2tZQuyQtFMtiSJIkSZIWhcHh0R2J5SkT2ycZHB6tKCK1k4hYBawBLimbXhsRl0XE6RGxd9nWC/yiZrPry7a52mfuY31EbIqITdu2LVzid2x8Ytb2rXO0S9JCMbksSZIkSVoU5kqEmSBTs0XEg4BzgDdk5m3AR4BHAqspRja/ZyH2k5mnZubazFy7YsXClKwYGhkj5lhmzWVJzWZyWZIkSZK0KDgpmaoQEZ0UieVPZuZGgMy8MTMnM/Nu4DTuKX0xBhxQs/n+Zdtc7U03ODxKztIewEB/XytCkNTGTC5LkiRJkhaFgf4+ujo7prV1dXaYIFPTREQAHwOuzMz31rTvV7PaC4DLy/vnAi+OiPtHxIHAQcD3gO8DB0XEgRFxP4pJ/85txTHMNbI/sVa5pOZzQj9JkiRJ0qIwlQgbHB5l6/gEPd1dDPT3mSBTMx0OvAzYEhGby7Y3A8dExGqKHO11wF8DZOYVEXE2xUR9dwGvycxJgIh4LTAMdACnZ+YVrTiAnu6uWWsu9zriX1ILmFyWJEmSJC0a69b0mkxWy2Tmt2DWksXnz7PNO4F3ztJ+/nzbNctTH7OCT1z881nbJanZLIshSZIkSZK0RF101baG2iVpIZlcliRJkiRJWqLmqrk8V7skLSSTy5IkSZIkSUtU9wM6G2qXpIVkclmSJEmSJGmJymysXZIWkhP6SZIkSZKabmhkjMHhUbaOT9DT3cVAf58T90kL4NaJ7Q21S9JCcuSyJEmSJKmphkbG2LBxC2PjEyQwNj7Bho1bGBoZqzo0acnr6e5qqF2SFpLJZUmSJElSUw0OjzKxfXJa28T2SQaHRyuKSFo+Bvr76OyIaW2dHcFAf19FEUlqJyaXJUmSJElNtXV8oqF2SQ2aWV/ZesuSWsTksiRJkiSpqbxsX2qeweFRtt89PZu8/e70ygBJLWFyWZIkSZLUVAP9fXR1dkxr6+rs8LJ9aQF4ZYCkKu1edQCSJEmSpOVt3ZpeoBhhuXV8gp7uLgb6+3a0S9p13Q/o5JbfbL9Xu1cGSGoFk8uSJEmSpKZbt6bXZLK0wIZGxrjjt3fdq90J/SS1imUxJEmSJEmSlqDZ6i0DPPB+u/tjjqSWMLksSZIkSZK0BM1VV/nWiXuXyZCkZjC5LEmSJEmStATNVVfZesuSWsXksiRJkiRJ0hI00N9HZ0dMa7PesqRWMrksSZIkSZK0VM0suXzvEsyS1DQmlyVJkiRJkpag2Sb02353Mjg8WlFEktqNyWVJkiRJkqQlaK4J/eZql6SFZnJZkiRJkiRpCXJCP0lVM7ksSZIkSWrY0MgYh5/yVQ48/jwOP+WrDI2MVR2S1Hac0E9S1XavOgBJkiRJ0tIyNDLGho1bmNg+CcDY+AQbNm4BYN2a3ipDk9qPE/pJqpAjlyVJkqQWiYg9IuJ7EfHDiLgiIt4+yzoPj4ivRMRlEfG1iNi/ZtmxEfGT8nZsa6OX7jE4PLojsTxlYvukk4hJLeaEfpKqZnJZkiRJap07gSMy8/eA1cCREfHEGeu8G/h4Zj4OOAk4GSAi9gFOAP4AOAw4ISL2blXgUi0nEZMWB9+LkqpmclmSJElqkSzcUT7sLG8zL2A+GPhqef8i4Pnl/X7gwsy8OTNvAS4EjmxyyNKsnERMWhx8L0qqmsllSZIkqYUioiMiNgM3USSLL5mxyg+Bo8r7LwD2jIiHAL3AL2rWu75sm/n86yNiU0Rs2rZt24LHL0ExiVhXZ8e0tq7ODicRk1rsqY9ZQcxo870oqZWc0E9qgVXHn1d1CJIkaZHIzElgdUR0A5+LiEMy8/KaVf4P8KGIeAXwDWAMmLzXE839/KcCpwKsXbvWaZ3UFFOT9g0Oj7J1fIKe7i4G+vuczE9qoaGRMc65dGza5S8BvPDxvb4XJbWMyWVJkiSpApk5HhEXUZS2uLymfSvlyOWIeBDwwnLdMeApNU+xP/C1lgUszbBujQksqUqzTayZwEVXedWKpNaxLIYkSZLUIhGxohyxTER0Ac8Arpqxzr4RMfU5fQNwenl/GHhmROxdTuT3zLJNktSGnMxP0mJgclmSJElqnf2AiyLiMuD7FDWXvxgRJ0XE88p1ngKMRsSPgYcB7wTIzJuBd5TbfR84qWyTJLUhJ/OTtBhYFkOSJElqkcy8DFgzS/vbau5/FvjsHNufzj0jmSVJbWygv48NG7dMK43hZH6SWs2Ry5IkSZIkSUvMujW9vPDxvXREANAR4WR+klrO5LIkSZIkSdISMzQyxjmXjjGZCcBkJudcOsbQyFjFkUlqJyaXJUmSJEmSlpjB4dFpJTEAJrZPMjg8WlFEktqRyWVJkiRJkqQlZuv4REPtktQMJpclSZIkSZKWmJ7urobaJakZTC5LkiRJkiQtMQP9fXR1dkxr6+rsYKC/r6KIJLWj3asOQJIkSZIkSY1Zt6YXKGovbx2foKe7i4H+vh3tktQKjlyWJEmSJEmSJDXMkcuSJEmSJElLzNDIGBs2bmFi+yQAY+MTbNi4BcDRy5JaxpHLkiRJkiRJS8zg8OiOxPKUie2TDA6PVhSRpHZkclmSJEmSJGmJ2To+0VC7JDWDZTEkSZIkqY0NjYw5IZi0BPV0dzE2SyK5p7urgmgktSuTy2pbq44/r+Ftrjvl2U2IRGq9iDgA+DjwMCCBUzPzAxFxIvAqYFu56psz8/xqopQkSc1mzVZp6Rro75v2/gXo6uxgoL+vwqgktRvLYkhSe7oLeFNmHgw8EXhNRBxcLntfZq4ubyaWJUlaxqzZKi1d69b08sLH99IRAUBHBC98fK8/DElqKZPLktSGMvOGzPxBef924ErAT6GSJLUZa7ZKS9fQyBjnXDrGZCYAk5mcc+kYQyNjFUcmqZ1YFkOS2lxErALWAJcAhwOvjYiXA5soRjffMss264H1ACtXrmxdsJpmV8r7wK6X+LGckCQtP9ZslZau+a48cPSypFZx5LIktbGIeBBwDvCGzLwN+AjwSGA1cAPwntm2y8xTM3NtZq5dsWJFq8KVJEkLbKC/j67Ojmlt1myVlgavPJC0GJhclqQ2FRGdFInlT2bmRoDMvDEzJzPzbuA04LAqY5QkSc21bk0vJx91KL3dXQTQ293FyUcd6qhHaQmY6woDrzyQ1EqWxZCkNhQRAXwMuDIz31vTvl9m3lA+fAFweRXxSZKk1lm3xgnApKVooL+PDRu3TCuN4ZUHklrN5LIktafDgZcBWyJic9n2ZuCYiFgNJHAd8NdVBCdJkiRpflM/Cg0Oj7J1fIKe7i4G+vv8sUhSS5lclqQ2lJnfAmKWRee3OhZJkiRJu8YrDyRVzeSyJEmSJEnSEjM0MuaoZUmVM7ksSZIkSZK0hAyNjE2rtzw2PsGGjVsATDBLaqndqg5AkiRJkiRJ9RscHp02kR/AxPZJBodHK4pIUrsyuSxJkiRJkrSEbB2faKhdkprF5LIkSZIkSdIS0tPd1VC7JDVL05LLEbFHRHwvIn4YEVdExNvL9gMj4pKIuDoiPhMR92tWDJIkSZIkScvNQH8fXZ0d09q6OjsY6O+rKCJJ7aqZI5fvBI7IzN8DVgNHRsQTgX8B3peZjwJuAY5rYgySJEmSJEnLyro1vZx81KH0dncRQG93FycfdaiT+Ulqud2b9cSZmcAd5cPO8pbAEcBLyvYzgROBjzQrDkmSJEmSpOVm3Zpek8mSKte05DJARHQAlwKPAj4MXAOMZ+Zd5SrXA7P2hBGxHlgPsHLlymaGKUmSJEmStKQMjYwxODzK1vEJerq7GOjvM9ksqeWaOqFfZk5m5mpgf+Aw4DENbHtqZq7NzLUrVqxoVoiSJEmSJElLytDIGBs2bmFsfIIExsYn2LBxC0MjY1WHJqnNNDW5PCUzx4GLgD8EuiNiasT0/oA9nyRJkiRJUp0Gh0eZ2D45rW1i+ySDw6MVRSSpXTUtuRwRKyKiu7zfBTwDuJIiyfxn5WrHAp9vVgySJEmSJEnLzdbxiYbaJalZmjlyeT/gooi4DPg+cGFmfhH4R+CNEXE18BDgY02MQZIkSZIk6V4i4oCIuCgifhQRV0TE68v2wYi4KiIui4jP1QycWxURExGxubx9tOa5Hh8RWyLi6oj4YEREM2Pv6e5qqF2SmqVpyeXMvCwz12Tm4zLzkMw8qWy/NjMPy8xHZeafZ+adzYpBkiRJkiRpDncBb8rMg4EnAq+JiIOBC4FDMvNxwI+BDTXbXJOZq8vbq2vaPwK8CjiovB3ZzMAH+vvo6uyY1tbV2cFAf18zdytJ97L7zleRJElaGlYdf17D21x3yrObEIkkSVrsMvMG4Iby/u0RcSXQm5lfqlntYu4p7TmriNgP2CszLy4ffxxYB1zQjLgB1q3pBYray1vHJ+jp7mKgv29HuyS1isllSZIkSZLU1iJiFbAGuGTGor8EPlPz+MCIGAFuA96amd8EeoHra9a5vmybbT/rgfUAK1euvE8xr1vTazJZUuWaWXNZkiRJkiRpUYuIBwHnAG/IzNtq2t9CUTrjk2XTDcDKzFwDvBH4VETs1ci+MvPUzFybmWtXrFixMAcgSRVy5LIkSZIkSWpLEdFJkVj+ZGZurGl/BfAc4GmZmQDlnFF3lvcvjYhrgEcDY8D+NU+7f9nWVEMjY5bFkFQ5Ry5LkiRJkqS2ExEBfAy4MjPfW9N+JPAPwPMy8zc17SsioqO8/wiKifuuLWs33xYRTyyf8+XA55sZ+9DIGBs2bmFsfIIExsYn2LBxC0MjTc9pS9I0JpclSZIkSVI7Ohx4GXBERGwub88CPgTsCVxYtn20XP/JwGURsRn4LPDqzLy5XPa3wL8DVwPX0MTJ/KCYyG9i++S0tontkwwOjzZzt5J0L5bFkCRJkiRJbSczvwXELIvOn2P9cyhKaMy2bBNwyMJFN7+t4xMNtUtSs5hcliRJkqQlzLqrUvvp6e5ibJZEck93VwXRSGpnlsWQJEmSpCXKuqtSexro76Ors2NaW1dnBwP9fRVFJKldmVyWJEmSpCXKuqtSe1q3ppeTjzqU3u4uAujt7uLkow71qgVJLWdZDEmSJElaoqy7KrWvdWt6TSZLqpwjlyVJkiRpiZqrvqp1VyVJUiuYXJYkSZKkJcq6q1L7GhoZ4/BTvsqBx5/H4ad81VrrkiphWQxJkiRJWqKmLokfHB5l6/gEPd1dDPT3eam8tMxNTeY5VXN9ajJPwPe/pJYyuSxJkiRJS5h1V6X2M99knvYHklrJshiSJEmSJElLiJN5SlosTC5LkiRJkiQtIU7mKWmxMLksSZIkSZK0hDiZp6TFwprLkiRJkiRJS4iTeUpaLEwuS5IkSZIkLTFO5ilpMbAshiRJkiRJkiSpYSaXJUmSJEmSJEkNM7ksSZIkSZIkSWqYNZclSZIkSZKWkKGRMSfzk7QomFyWJEmSJElaIoZGxtiwcQsT2ycBGBufYMPGLQAmmCW1nGUxJEmSJEmSlojB4dEdieUpE9snGRwerSgiSe3M5LIkSZIkSdISsXV8oqF2SWomy2JIDVh1/HlVhyBJkiRJamM93V2MzZJI7unuqiAaSe3OkcuSJEmSJElLxEB/H12dHdPaujo7GOjvqygiSe3MkcuSJEmSJElLxNSkfYPDo2wdn6Cnu4uB/j4n85NUCZPLkiRJkiRJS8i6Nb0mkyUtCpbFkCRJkiRJkiQ1zOSyJEmSJEmSJKlhJpclSZIkSZIkSQ2z5rIkSZIkSdISMjQy5oR+khYFk8uSJEmSJElLxNDIGBs2bmFi+yQAY+MTbNi4BcAEs6SWsyyGJEmS1CIRsUdEfC8ifhgRV0TE22dZZ2VEXBQRIxFxWUQ8q2xfFRETEbG5vH209UcgSara4PDojsTylIntkwwOj1YUkaR25shlSZIkqXXuBI7IzDsiohP4VkRckJkX16zzVuDszPxIRBwMnA+sKpddk5mrWxqxJGlR2To+0VC7JDWTI5clSZKkFsnCHeXDzvKWM1cD9irvPxjY2qLwJElLQE93V0PtktRMJpclSZKkFoqIjojYDNwEXJiZl8xY5UTgLyLieopRy39Xs+zAslzG1yPiSXM8//qI2BQRm7Zt29aEI5AkVWmgv4+uzo5pbV2dHQz091UUkaR2ZnJZkiRJaqHMnCxLW+wPHBYRh8xY5RjgjMzcH3gW8J8RsRtwA7AyM9cAbwQ+FRF7zdiWzDw1M9dm5toVK1Y09VgkSa23bk0vJx91KL3dXQTQ293FyUcd6mR+kiphzWVJkiSpApk5HhEXAUcCl9csOq5sIzO/GxF7APtm5k0UNZvJzEsj4hrg0cCm1kYuSaraujW9JpMlLQqOXJYkSZJaJCJWRER3eb8LeAZw1YzVfg48rVznd4E9gG3lth1l+yOAg4BrWxS6JEmSdC+OXJYkSZJaZz/gzDJJvBtwdmZ+MSJOAjZl5rnAm4DTIuLvKSb3e0VmZkQ8GTgpIrYDdwOvzsybKzoOSZIkyeSyJEmS1CqZeRmwZpb2t9Xc/xFw+CzrnAOc09QAJUmSpAZYFkOSJEmSJEmS1DCTy5IkSZIkSZKkhlkWQ5IkSZIkaQkZGhljcHiUreMT9HR3MdDfx7o1vVWHJakNmVyWJEmSJElaIoZGxtiwcQsT2ycBGBufYMPGLQAmmCW1nGUxJEmSJEmSlojB4dEdieUpE9snGRwerSgiSe3M5LIkSZIkSdISsXV8oqF2SWomk8uSJEmSJElLRE93V0PtktRMJpclSZIkSZKWiIH+Pro6O6a1dXV2MNDfV1FEktqZE/pJkiRJkiQtEVOT9g0Oj7J1fIKe7i4G+vuczE9SJUwuS5IkSZIkLSHr1vSaTJa0KJhcliRJkqQWGRoZc7ShJElaNkwuS5IkSVILDI2MsWHjFia2TwIwNj7Bho1bAEwwS5KkJckJ/SRJkiSpBQaHR3cklqdMbJ9kcHi0oogkSZLuG5PLkiRJktQCW8cnGmqXJEla7EwuS5IkSVIL9HR3NdQuSZK02JlcliRJkqQWGOjvo6uzY1pbV2cHA/19FUUkaSkaGhnj8FO+yoHHn8fhp3yVoZGxqkOS1Mac0E+SJEmSWmBq0r7B4VG2jk/Q093FQH+fk/lJqpsTg0pabEwuS5IkSVKLrFvTawJI0i6bb2JQ+xZJVbAshiRJkiRJ0hLgxKCSFhuTy5IkSZIkSUuAE4NKWmxMLkuSJEmSJC0BTgwqabGx5rIkSZIkSdIS4MSgkhYbk8uSJEmSJElLhBODSlpMLIshSZIkSZIkSWqYyWVJkiRJkiRJUsNMLkuSJEmSJEmSGmZyWZIkSZIkSZLUMJPLkiRJkiSpLUXEARFxUUT8KCKuiIjXl+37RMSFEfGT8t+9y/aIiA9GxNURcVlE/H7Ncx1brv+TiDi2qmOSpFYyuSxJkiRJktrVXcCbMvNg4InAayLiYOB44CuZeRDwlfIxwJ8CB5W39cBHoEhGAycAfwAcBpwwlZCWpOWsacnleX79OzEixiJic3l7VrNikCRJkiRJmktm3pCZPyjv3w5cCfQCzwfOLFc7E1hX3n8+8PEsXAx0R8R+QD9wYWbenJm3ABcCR7buSCSpGrs38bmnfv37QUTsCVwaEReWy96Xme9u4r4lSZIkSZLqFhGrgDXAJcDDMvOGctEvgYeV93uBX9Rsdn3ZNlf7zH2spxjxzMqVKxcwekmqRtOSy2UnfEN5//aImPr1T5IkSZIkadGIiAcB5wBvyMzbImLHsszMiMiF2E9mngqcCrB27dpdes6hkTEGh0fZOj5BT3cXA/19rFtjukVSNVpSc3nGr38Ary0L359uDSJJkiRJklSViOikSCx/MjM3ls03luUuKP+9qWwfAw6o2Xz/sm2u9gU1NDLGho1bGBufIIGx8Qk2bNzC0MiC70qS6tL05PLMX/8oit0/ElhNMbL5PXNstz4iNkXEpm3btjU7TEmSJEmS1GaiGKL8MeDKzHxvzaJzgWPL+8cCn69pf3kUngjcWl65PQw8MyL2LgfRPbNsW1CDw6NMbJ+c1jaxfZLB4dGF3pUk1aWZNZdn/fUvM2+sWX4a8MXZtl2IS0UkSZIkSZLmcTjwMmBLRGwu294MnAKcHRHHAT8Dji6XnQ88C7ga+A3wSoDMvDki3gF8v1zvpMy8eaGD3To+0VC7JDVb05LLc/36FxH71RTFfwFwebNikCTNLiIOAD5OMTFJAqdm5gciYh/gM8Aq4Drg6HK2a0mSJGnZycxvATHH4qfNsn4Cr5njuU4HTl+46O6tp7uLsVkSyT3dXc3crSTNqZllMaZ+/TsiIjaXt2cB74qILRFxGfBU4O+bGIMkaXZ3AW/KzIOBJwKviYiDgeOBr2TmQcBXyseSJEmSFoGB/j66OjumtXV1djDQ31dRRJLaXdNGLs/z69/5zdqnJKk+5RUkN5T3b4+IK4Fe4PnAU8rVzgS+BvxjBSFKkiRJmmHdml6gqL28dXyCnu4uBvr7drRLUqs1teayJGnxi4hVwBrgEuBhNaWLfklRNmO2bdYD6wFWrlzZgiglSZIkQZFgNpksabFoZlkMSdIiFxEPoph49Q2ZeVvtsrKe3KwTqmbmqZm5NjPXrlixogWRSpIkSZKkxcbksiS1qYjopEgsfzIzN5bNN0bEfuXy/YCbqopPkiRJkiQtbiaXJakNRUQAHwOuzMz31iw6Fzi2vH8s8PlWxyZJkiRJkpYGay5LUns6HHgZsCUiNpdtbwZOAc6OiOOAnwFHVxOeJEmSJEla7EwuS1IbysxvATHH4qe1MhZJkiRJkrQ0WRZDkiRJkiRJktQwk8uSJEmSJEmSpIaZXJYkSZIkSZIkNczksiRJkiRJkiSpYU7oJ0mSJEmStEQMjYwxODzK1vEJerq7GOjvY92a3qrDktSmTC5LkiRJkiQtAUMjY2zYuIWJ7ZMAjI1PsGHjFgATzJIqYVkMSZIkSZKkJWBweHRHYnnKxPZJBodHK4pIUrszuSxJkiRJkrQEbB2faKhdkprN5LIkSZIkSdIS0NPd1VC7JDWbyWVJkiRJkqQlYKC/j67OjmltXZ0dDPT3VRSRpHbnhH6SJEmSJElLwNSkfYPDo2wdn6Cnu4uB/j4n85NUGZPLkiRJkiRJS8S6Nb0mkyUtGpbFkCRJkiRJkiQ1zOSyJEmSJEmSJKlhJpclSZIkSZIkSQ0zuSxJkiRJkiRJapgT+kmSJEkNioiHAocDPcAEcDmwKTPvrjQwSZIkqYVMLkuSJEl1ioinAscD+wAjwE3AHsA64JER8VngPZl5W2VBSpIkSS1iclmSJEmq37OAV2Xmz2cuiIjdgecAzwDOaXVgkiRJUquZXJYkSZLqlJkD8yy7CxhqXTSSJElStZzQT5IkSWpQRLw+IvaKwsci4gcR8cyq45IkSZJayeSyJEmS1Li/LOsqPxPYG3gZcEq1IUmSJEmtZXJZkiRJalyU/z4L+M/MvKKmTZIkSWoL1lyWJEmSGndpRHwJOBDYEBF7AndXHJMkqQ0MjYwxODzK1vEJerq7GOjvY92a3qrDktSmTC5LkiRJjTsOWA1cm5m/iYiHAK+sNiRJ0nI3NDLGho1bmNg+CcDY+AQbNm4BMMEsqRKWxZAkSZIalJl3AzcCB0fEk4HHAt2VBiVJWvYGh0d3JJanTGyfZHB4tKKIJLU7Ry5LkiRJDYqIfwFeBPwImPqWn8A3KgtKkrTsbR2faKhdkprN5LIkSZLUuHVAX2beWXUgkqT20dPdxdgsieSe7q4KopEky2JIkiRJu+JaoLPRjSJij4j4XkT8MCKuiIi3z7LOyoi4KCJGIuKyiHhWzbINEXF1RIxGRP99PAZJ0hIz0N9HV2fHtLauzg4G+vsqikhSu3PksiRJktS43wCbI+IrwI7Ry5n5up1sdydwRGbeERGdwLci4oLMvLhmnbcCZ2fmRyLiYOB8YFV5/8UU9Z17gC9HxKMzc3LmTiRJy9PUpH2Dw6NsHZ+gp7uLgf4+J/OTVBmTy5IkSVLjzi1vDcnMBO4oH3aWt5y5GrBXef/BwNby/vOBT5elOH4aEVcDhwHfbTQOSdLStW5Nr8lkSYuGyWVJkiSpQZl5ZkTcD3h02TSamdvr2TYiOoBLgUcBH87MS2asciLwpYj4O+CBwNPL9l6gdoTz9WXbzOdfD6wHWLlyZV3HI0mSJO0Kay5LkiRJDYqIpwA/AT4M/Cvw44h4cj3bZuZkZq4G9gcOi4hDZqxyDHBGZu4PPAv4z4io+3N7Zp6amWszc+2KFSvq3UySJElqmCOXJUmSpMa9B3hmZo4CRMSjgbOAx9f7BJk5HhEXAUcCl9csOq5sIzO/GxF7APsCY8ABNevtX7ZJkiRJlXDksiRJktS4zqnEMkBm/piifvK8ImJFRHSX97uAZwBXzVjt58DTynV+F9gD2EZR4/nFEXH/iDgQOAj43n0/FEmSJGnXOHJZkiRJatymiPh34BPl45cCm+rYbj/gzLLu8m7A2Zn5xYg4CdiUmecCbwJOi4i/p5jc7xXlRIBXRMTZwI+Au4DXZObkwh6WJEmSVD+Ty5IktZlVx59XdQjScvA3wGuA15WPv0lRe3lemXkZsGaW9rfV3P8RcPgc278TeOcuxCtJkiQtOJPLkiRJUoMy807gveVNkiRJaksmlyVJkqQ6RcTZmXl0RGyhKFkxTWY+roKwJEmSpEqYXJYkSZLq9/ry3+dUGoUkSZK0CJhcliRJkuqUmTeU//6s6lgkSZKkqplcliRJkuoUEbczSzkMIIDMzL1aHJIkSZJUGZPLkiRJUp0yc8+qY5AkSZIWi50mlyPigcBEZt4dEY8GHgNckJnbmx6dJGle9tGS1FoRsc98yzPz5lbFIkmSJFWtnpHL3wCeFBF7A18Cvg+8CHhpMwOTJNXFPlqSWutSirIYMcuyBB7R2nAkSZKk6tSTXI7M/E1EHAf8a2a+KyI2NzkuSVJ97KMlqYUy88CqY5Akta+hkTEGh0fZOj5BT3cXA/19rFvTW3VYktpYXcnliPhDilFwx5VtHc0LSZLUAPtoSWqhiHhMZl4VEb8/2/LM/EGrY5IktYehkTE2bNzCxPZJAMbGJ9iwcQuACWZJlaknufx6YAPwucy8IiIeAVzU3LAkSXWyj5ak1nojsB54zyzLEjiiteFIktrF4PDojsTylIntkwwOj5pcllSZepLLt2Tm86YeZOa1wOuaF5IkqQH20ZLUQpm5vvz3qVXHIklqL1vHJxpql6RWqCe5/K8RcX/gDOCTmXlrc0OSJDXAPlqSKhARHcCzgVXUfKbOzPdWFZMkaXnr6e5ibJZEck93VwXRSFJht52tkJlPAv4COAC4NCI+FRHPbHpkkqSdso+WpMp8AXgF8BBgz5qbJElNMdDfR1fn9OlVujo7GOjvqygiSapv5DKZ+eOIeCuwCfggsCYiAnhzZm5sZoCSpPnZR0tSJfbPzMdVHYQkqX1M1VUeHB5l6/gEPd1dDPT3WW9ZUqV2mlyOiMcBr6S47O9C4LmZ+YOI6AG+C5i4kKSK2EdLUmUuiIhnZuaXqg5E983QyJiJGklLxro1vfZRkhaVekYu/z/g3ylGwO0o7pOZW8uRcpKk6thHS1I1LgY+FxG7AduBADIz96o2LDViaGSMDRu3MLF9EoCx8Qk2bNwCYPJGkiSpDjtNLmfmn8yz7D8XNhxJUiPsoyWpMu8F/hDYkplZdTDaNYPDozsSy1Mmtk8yODxqclmSJKkO9ZTFOAg4GTgY2GOqPTMf0cS4JEl1sI+WpMr8ArjcxPLStnV8oqF2SZIkTVdPWYz/AE4A3gc8laK2527NDEqSVDf7aOk+WnX8ebu03XWnPHuBI9EScy3wtYi4ALhzqjEz31tdSGpUT3cXY7Mkknu6uyqIRlIVIuJ04DnATZl5SNn2GaCvXKUbGM/M1RGxCrgSGC2XXZyZry63eTxwBtAFnA+83h8gJbWDehIQXZn5FSAy82eZeSLFxFGSpOrZR0tSNX4KfAW4H7BnzU1LyEB/H12dHdPaujo7GOjvm2MLScvQGcCRtQ2Z+aLMXJ2Zq4FzmD5J9jVTy6YSy6WPAK8CDipv055TkparekYu31lOVPKTiHgtMAY8qLlhSZLqZB8tSRXIzLdXHYPuu6m6yoPDo2wdn6Cnu4uB/j7rLUttJDO/UY5IvpeICOBo4Ij5niMi9gP2ysyLy8cfB9YBFyxosJK0CNWTXH498ADgdcA7KDrVY5sZlCSpbvbRkiTdB+vW9JpMljSXJwE3ZuZPatoOjIgR4DbgrZn5TaAXuL5mnevLtnuJiPXAeoCVK1c2JWhJaqWdJpcz8/vl3TsoanlKkhYJ+2hJkiSpaY4Bzqp5fAOwMjN/VdZYHoqIxzbyhJl5KnAqwNq1a63JLGnJm7fmckQcGxE/iIhfl7dNEfHyVgUnSZqbfbQkSZLUHBGxO3AU8Jmptsy8MzN/Vd6/FLgGeDRFabr9azbfv2yTpGVvzpHLEXEs8AbgjcAPgAB+HxiMiMzM/2xJhJKke7GPlqRqRcQKiombVlHzmToz/7KqmCRJC+rpwFWZuaPcRdn335yZkxHxCIqJ+67NzJsj4raIeCJwCfBy4P9VErUktdh8ZTH+BnhBZl5X0/bViHgh8GnAxIUkVcc+WpKq9Xngm8CXgcmKY5Ek7aKIOAt4CrBvRFwPnJCZHwNezPSSGABPBk6KiO3A3cCrM/PmctnfAmcAXRQT+TmZn6S2MF9yea8ZSQsAMvO6iNireSFJkupgHy1J1XpAZv5j1UFIku6bzDxmjvZXzNJ2DnDOHOtvAg5Z0OAkaQmYr+byxC4ukyQ1n320JFXrixHxrKqDkCRJkqo038jl342Iy2ZpD+ARTYpHklQf+2hJqtbrgTdHxJ3Ador+NzPTq0ckSZLUNuZNLrcsCklSo+yjJalCmbln1TFIkiRJVZszuZyZP2tlIJKk+tlHS1I1IuIxmXlVRPz+bMsz8wetjkmSJEmqynwjl++TiDgA+DjwMCCBUzPzAxGxD/AZYBVwHXB0Zt7SrDgkSZKkBfRGYD3wnlmWJXBEa8ORJNWKiL2Bg4A9ptoy8xvVRSRJy1vTksvAXcCbMvMHEbEncGlEXAi8AvhKZp4SEccDxwPOtC1JkqRFLzPXl/8+tepYJEnTRcRfUdTE3x/YDDwR+C7L6Ie/oZExBodH2To+QU93FwP9faxb01t1WJLa2G5zLYiIr5T//suuPHFm3jB1WWBm3g5cCfQCzwfOLFc7E1i3K88vSe3svvbRkqRdExF/vJPle0XEIa2KR5I0zeuBJwA/K38EXAOMVxrRAhoaGWPDxi2MjU+QwNj4BBs2bmFoZKzq0CS1sflGLu8XEX8EPC8iPk0xA/YOjdSTi4hVFJ36JcDDMvOGctEvKcpmzLbNeopLDlm5cmW9u5KkdrFgfbQkqSEvjIh3Af8NXApso7j0+lHAU4GHA2+qLjxJamu/zczfRgQRcf+yRn5f1UEtlMHhUSa2T05rm9g+yeDwqKOXJVVmvuTy24B/oric5L0zltVdTy4iHgScA7whM2+LuCf/kZkZETnbdpl5KnAqwNq1a2ddR5La2IL00ZKkxmTm35dziLwQ+HNgP2CC4iq9f8vMb1UZnyS1uesjohsYAi6MiFuAZTMR9tbxiYbaJakV5kwuZ+Zngc9GxD9l5jt25ckjopMisfzJzNxYNt8YEftl5g0RsR9w0648tyS1s4XooyVJuyYzbwZOK2+SpEUiM19Q3j0xIi4CHkxxpcmy0NPdxdgsieSe7q4KopGkwpw1l6dk5jsi4nkR8e7y9px6njiKIcofA67MzNpRdecCx5b3jwU+32jQkqTCrvbRkiRJ0nIREfvMvAFbgG8BD6o4vAUz0N9HV2fHtLauzg4G+pdN5Q9JS9B8ZTEAiIiTgcOAT5ZNr4+IP8rMN+9k08OBlwFbImJz2fZm4BTg7Ig4juLylKN3JXBJ0n3qoyVJkqTl4lKK0nABrARuKe93Az8HDqwssgU0VVd5cHiUreMT9HR3MdDfZ71lSZXaaXIZeDawOjPvBoiIM4ERikTxnMp6czHH4qc1EqQkaU671EdLrbLq+POqDkGSJC1zmXkgQEScBnwuM88vH/8psK7C0BbcujW9JpMlLSr1JJeh+LXv5vL+g5sTiiRpF3VjHy1JLRERR823vGaeEUlS6z0xM1819SAzL4iId1UZkCQtd/Ukl08GRspi+AE8GTi+qVFJkuplHy1JrfXc8t+HAn8EfLV8/FTgO4DJZUmqztaIeCvwifLxS4GtFcYjScveTpPLmXlWRHwNeELZ9I+Z+cumRiVJqot9tCS1Vma+EiAivgQcnJk3lI/3A86oMDRJEhwDnAB8rnz8jbJNktQkdZXFKD80n9vkWCRJu8A+WpIqccBUYrl0I8UkUpKkimTmzcDrq45DktpJvTWXJUmSJN3jKxExDJxVPn4R8OUK45GkthcRaykmtl5FTb4jMx9XVUyStNyZXJYkSZIalJmvjYgXUNS6Bzg1Mz833zaSpKb7JDAAbAHurjgWSWoL8yaXI6IDuCIzH9OieCRJdbKPlqTKfQe4C0jgexXHIkmCbZlpuThJaqHd5luYmZPAaERYP06SFhn7aEmqTkQcTZFQ/jPgaOCSiPizaqOSpLZ3QkT8e0QcExFHTd2qDkqSlrN6ymLsDVwREd8Dfj3VmJnPa1pUkqR62UdLUjXeAjwhM28CiIgVFDWXP1tpVJLU3l4JPAbo5J6yGAlsrCwiSVrm6kku/1PTo5Ak7Sr7aEmqxm5TieXSr9jJVYGSpKZ7Qmb2VR2EJLWTnSaXM/PrEfFw4KDM/HJEPADoaH5okqSdsY+WpMr8d0QMA2eVj18EnF9hPJIk+E5EHJyZP6o6EElqFztNLkfEq4D1wD7AI4Fe4KPA05obmiRpZ+yjJakamTlQ1vH847Lp1Mz8XJUxSZJ4IrA5In4K3AkEkJn5uGrDkqTlq56yGK8BDgMuAcjMn0TEQ5salSSpXvbRklSdbwPbKep5fq/iWCRJcGTVAUhSu6mnLtydmfm/Uw8iYneKD9CSpOrZR0tSBSLiaIqE8p8BRwOXRMSfVRuVJLW3zPxZZv4MmKD4TDx1kyQ1ST0jl78eEW8GuiLiGcDfAl9obliSpDrZR0tSNd5CMXHUTQARsQL4MvDZSqOSpDYWEc8D3gP0ADcBDweuBB5bZVyStJzVM3L5eGAbsAX4a4qJSt7azKAkSXWzj5akauw2lVgu/Yr6PltLkprnHRR1l3+cmQdSzENycbUhSdLyttORy5l5d0ScSVHPM4HRzPSyEklaBHa1j46I04HnADdl5iFl24nAqyiS1QBvzszzmxK4JC19/x0Rw8BZ5eMXUfzAJ0mqzvbM/FVE7BYRu2XmRRHx/qqDkqTlbKfJ5Yh4NvBR4BqKmVYPjIi/zswLmh2cJGl+96GPPgP4EPDxGe3vy8x3L3igkrTMZOZARLwQOLxsOjUzP1dlTJIkxiPiQcA3gE9GxE3AryuOSZKWtXpqLr8HeGpmXg0QEY8EzgNMLktS9Xapj87Mb0TEquaHJ0nLV2aeA5xTdRySpB2eD/wW+HvgpcCDgZMqjUiSlrl66sLdPpW0KF0L3N6keCRJjVnoPvq1EXFZRJweEXvfx9gkadmKiKMi4icRcWtE3BYRt0fEbVXHJUntLDN/nZmTmXlXZp6ZmR/MzF9VHZckLWdzjlyOiKPKu5si4nzgbIp6nn8OfL8FsUmS5tCkPvojFJOgZPnve4C/nGP/64H1ACtXrtzF3UnSkvYu4LmZeWXVgUhSu4uI2yk+w95rEZCZuVeLQ5KktjFfWYzn1ty/EfiT8v42oKtpEUmS6rHgfXRm3jh1PyJOA744z7qnAqcCrF271kleJbWjG00sS9LikJl7Vh2DJLWrOZPLmfnKVgYiSapfM/roiNgvM28oH74AuHyh9yFJS92MK0c+AwwBd04tz8yNVcQlSWoPQyNjDA6PsnV8gp7uLgb6+1i3prfqsCS1sZ1O6BcRBwJ/B6yqXT8zn9e8sCRJ9djVPjoizgKeAuwbEdcDJwBPiYjVFJcUXgf8dTNilqQlrvbKkd8Az6x5nIDJZUlSUwyNjLFh4xYmtk8CMDY+wYaNWwBMMEuqzE6TyxSjMT4GfAG4u6nRSJIaNcQu9NGZecwszR9boJgkadmaunIkIh7iJFGSpFYaHB7dkVieMrF9ksHhUZPLkipTT3L5t5n5waZHIknaFfbRklSNiyNiM/AfwAWZaf15SVJTbR2faKhdklqhnuTyByLiBOBLTK8n94OmRSVJqpd9tCRV49HA04G/BD4YEWcDZ2Tmj6sNS5K0XPV0dzE2SyK5p3uX5vOWpAVRT3L5UOBlwBHcc8l1lo8lSdWyj5akCpQjlS8ELoyIpwKfAP42In4IHJ+Z3600QEnSsjPQ3zet5jJAV2cHA/19FUYlqd3Vk1z+c+ARmfm/zQ5GktQw+2hJqkBEPAT4C4of+G6kmFz1XGA18F/AgZUFJ0lalqbqKg8Oj7J1fIKe7i4G+vustyypUvUkly8HuoGbmhuKJGkX2EdLUjW+C/wnsC4zr69p3xQRH60oJknSMrduTa/JZEmLSj3J5W7gqoj4PtPreT6vWUFJkurWjX20JFWhb65J/DLzX+baKCL2AL4B3J/is/hnM/OEGeu8D3hq+fABwEMzs7tcNglsKZf93P5ekiRJVaonuXzCzleRJFXEPlqSWigivkBR256IuNfyOpK9dwJHZOYdEdEJfCsiLsjMi2ue4+9r9vd3wJqa7Scyc/WuH4EkSZK0cHaaXM7Mr7ciEElS4+yjJanl3n1fNi5HO99RPuwsb7OOgC4dgz8kSpIkaZHaaXI5Im7nng+896P4APzrzNyrmYFJknbOPlqSWqv2R72I6AJWZuZoI88RER3ApcCjgA9n5iVzrPdwiokBv1rTvEdEbALuAk7JzKFZtlsPrAdYuXJlI6FJkiRJDdltZytk5p6ZuVeZqOgCXgj8a9MjkyTtlH20JFUjIp4LbAb+u3y8OiLOrWfbzJwsS1vsDxwWEYfMseqLKWoyT9a0PTwz1wIvAd4fEY+c5flPzcy1mbl2xYoVdR+TJEmS1KidJpdrZWEI6G9OOJKkXWUfLUktdSJwGDAOkJmbKUYZ1y0zx4GLgCPnWOXFwFkzthkr/70W+BrT6zFLkiRJLVVPWYyjah7uBqwFftu0iCRJdbOPlqTKbM/MW2dM6jdf7WQAImJFue14WVbjGcC/zLLeY4C9ge/WtO0N/CYz74yIfYHDgXfdt8OQJEmSdt1Ok8vAc2vu3wVcBzy/KdFIkhplHy1J1bgiIl4CdETEQcDrgO/Usd1+wJll3eXdgLMz84sRcRKwKTOnSmu8GPh0OQHglN8F/i0i7i63PSUzf7RQByRJkiQ1aqfJ5cx8ZSsCkSQ1zj5akirzd8BbgDuBTwFfAk7a2UaZeRmzlLLIzLfNeHziLOt8Bzh018KVJEmSFt6cyeWIeNtcyyhKe76jCfFIkupgHy1JlTsmM99CkWAGICJOAY6vLiRJkiSpteYbufzrWdoeCBwHPAQwcSFJ1bGPliq26vjzGt7mulOe3YRIVJEXRsRvM/OTABHxIaCr4pgkSZKklpozuZyZ75m6HxF7Aq8HXgl8GnjPXNtJkprPPlqSKvdC4Nyy/vGRwHhmHldxTJIkSVJLzVtzOSL2Ad4IvBQ4E/j9zLylFYFJkuZnHy1JrVf2vVP+ChgCvg28PSL2ycybKwlMkiRJqsB8NZcHgaOAU4FDM/OOlkUlSZqXfbQkVeZSIIGo+ffZ5S2BR1QXmiRJktRa841cfhPF7NdvBd4SEVPtQTFZ1F5Njk2SNDf7aEmqQGYeWHUMkiRJ0mIxX83l3VoZiCSpfvbRklSNiDgiM78aEUfNtjwzN7Y6JkmSJKkq89ZcliRJkjTNnwBfBZ47y7IETC5LkiSpbZhcliRJkuqUmSeU/75y5rKIeGHrI5IkSZKq42XVkiRJ0sJ4X9UBSJIkSa1kclmSJElaGLHzVSRJkqTlw+SyJEmStDCy6gAkSZKkVrLmsiRJklSniNjC7EnkAB7W4nAkSZKkSplcliRJkur3nKoDkCQtnIg4naJvvykzDynbTgReBWwrV3tzZp5fLtsAHAdMAq/LzOGy/UjgA0AH8O+ZeUorj0OSqmJyWZIkSapTZv6s6hgkSQvqDOBDwMdntL8vM99d2xARBwMvBh4L9ABfjohHl4s/DDwDuB74fkScm5k/ambgkrQYmFyWJEmSJEltKTO/ERGr6lz9+cCnM/NO4KcRcTVwWLns6sy8FiAiPl2ua3JZ0rLnhH6SJEmSJEnTvTYiLouI0yNi77KtF/hFzTrXl21ztd9LRKyPiE0RsWnbtm2zrSJJS4rJZUmSJEmSpHt8BHgksBq4AXjPQj1xZp6amWszc+2KFSsW6mklqTKWxZAkSZIaFBFbgJzRfCuwCfjnzPxV66OSJC2EzLxx6n5EnAZ8sXw4BhxQs+r+ZRvztEvSsmZyWZIkSWrcBcAk8Kny8YuBBwC/pJgc6rnVhCVJuq8iYr/MvKF8+ALg8vL+ucCnIuK9FBP6HQR8DwjgoIg4kCKp/GLgJc2IbWhkjMHhUbaOT9DT3cVAfx/r1sxagUOSWsLksiRJktS4p2fm79c83hIRP8jM34+Iv6gsKklSQyLiLOApwL4RcT1wAvCUiFhNcYXKdcBfA2TmFRFxNsVEfXcBr8nMyfJ5XgsMAx3A6Zl5xULHOjQyxoaNW5jYPgnA2PgEGzZuATDBLKkyJpclSZKkxnVExGGZ+T2AiHgCRUIBioSDJGkJyMxjZmn+2DzrvxN45yzt5wPnL2Bo9zI4PLojsTxlYvskg8OjJpclVcbksiRJktS4vwJOj4gHUVwOfRtwXEQ8EDi50sgkScvS1vGJhtolqRVMLkuSJEkNyszvA4dGxIPLx7fWLD67mqgkSctZT3cXY7Mkknu6uyqIRpIKu1UdgCRJkrTURMSDywmdvgJ8JSLeM5VoliSpGQb6++jq7JjW1tXZwUB/X0URSZLJZUmSJGlXnA7cDhxd3m4D/qPSiCRJy9q6Nb2cfNSh9HZ3EUBvdxcnH3Wo9ZYlVcqyGJIkSVLjHpmZL6x5/PaI2FxVMJKk9rBuTa/JZEmLiiOXJUmSpMZNRMQfTz2IiMMBZ1SSJElSW3HksiRJktS4VwMfr6mzfAtwbIXxSJIkSS1nclmSJElqUGb+EPi9iNirfHxbRLwBuKzSwCRJkqQWsiyGJEmStIsy87bMvK18+MZKg5EkSZJazOSyJEmStDCi6gAkSZKkVjK5LEmSJC2MrDoASZIkqZWsuSxJkiTVKSJuZ/YkcgBdLQ5HkiRJqpTJZUmSJKlOmbln1TFIkiRJi4VlMSRJkiRJkiRJDWtacjkiTo+ImyLi8pq2EyNiLCI2l7dnNWv/kiRJkiRJkqTmaebI5TOAI2dpf19mri5v5zdx/5IkSZIkSZKkJmlacjkzvwHc3KznlyRJkiRJkiRVp4qay6+NiMvKshl7z7VSRKyPiE0RsWnbtm2tjE+SJEmSJEmStBOtTi5/BHgksBq4AXjPXCtm5qmZuTYz165YsaJF4UmSJEmSJEmS6tHS5HJm3piZk5l5N3AacFgr9y9JkiRJkiRJWhgtTS5HxH41D18AXN7K/UuSJEmSJEmSFsbuzXriiDgLeAqwb0RcD5wAPCUiVgMJXAf8dbP2L0mSJEmSJElqnqYllzPzmFmaP9as/UmSJEmSJEmSWqfVE/pJkiRJkiRJkpYBk8uSJEmSJEmSpIaZXJYkSZIkSZIkNczksiRJkiRJkiSpYSaXJUmSJEmSJEkNM7ksSZIkSZIkSWqYyWVJkiRJkiRJUsNMLkuSJEmSJEmSGmZyWZIkSZIkSZLUMJPLkiRJkiRJkqSGmVyWJEmSJEmSJDVs96oDkGqtOv68hre57pRnNyESSZIkSZIkSfNx5LIkSZIkSZIkqWEmlyVJkiRJkiRJDTO5LEmSJEmSJElqmMllSZIkSZIkSVLDTC5LkiRJkiRJkhpmclmSJEmSJEmS1DCTy5IkSZIkSZKkhu1edQCSJEmSJEma39DIGIPDo2wdn6Cnu4uB/j7WremtOixJbc7ksiRJkiRJ0iI2NDLGho1bmNg+CcDY+AQbNm4BMMEsqVKWxZAkSZIkSVrEBodHdySWp0xsn2RweLSiiCSpYHJZkiRJkiRpEds6PtFQuyS1isllSZIkSZKkRaynu6uhdklqFZPLkiRJkiRJi9hAfx9dnR3T2ro6Oxjo76soIkkqOKGfJEmSJEnSIjY1ad/g8Chbxyfo6e5ioL/PyfwkVc7ksiRJkiRJ0iK3bk2vyWRJi45lMSRJkqQWiYg9IuJ7EfHDiLgiIt4+yzrvi4jN5e3HETFes+zYiPhJeTu2pcFLkiRJMzhyWZIkSWqdO4EjMvOOiOgEvhURF2TmxVMrZObfT92PiL8D1pT39wFOANYCCVwaEedm5i0tPQJJkiSp5MhlSZIkqUWycEf5sLO85TybHAOcVd7vBy7MzJvLhPKFwJFNC1aSJEnaCZPLkiRJUgtFREdEbAZuokgWXzLHeg8HDgS+Wjb1Ar+oWeX6sm3mdusjYlNEbNq2bduCxi5JkiTVMrksSZIktVBmTmbmamB/4LCIOGSOVV8MfDYzJxt8/lMzc21mrl2xYsV9jFaSJEmam8llSZIkqQKZOQ5cxNylLV7MPSUxAMaAA2oe71+2SZIkSZUwuSxJkiS1SESsiIju8n4X8AzgqlnWewywN/DdmuZh4JkRsXdE7A08s2yTJEmSKrF71QFIkiRJbWQ/4MyI6KAY6HF2Zn4xIk4CNmXmueV6LwY+nZk7JvvLzJsj4h3A98umkzLz5lYGL0mSJNUyuSxJkiS1SGZeBqyZpf1tMx6fOMf2pwOnNyU4SZIkqUGWxZCkNhQRp0fETRFxeU3bPhFxYUT8pPx37ypjlCRJkiRJi5vJZUlqT2dw7wmkjge+kpkHAV8pH0uSJEmSJM3K5LIktaHM/AYws07n84Ezy/tnAutaGZMkSZIkSVparLksSZrysMy8obz/S+Bhc60YEeuB9QArV65sQWhLx6rjz6s6BEmSJEmSWsKRy5Kke8nMBHKe5adm5trMXLtixYoWRiZJkiRJkhYLk8uSpCk3RsR+AOW/N1UcjyRJktRUc0x0PRgRV0XEZRHxuYjoLttXRcRERGwubx+t2ebxEbElIq6OiA9GRFRwOJLUciaXJUlTzgWOLe8fC3y+wlgkSZKkVjiDe090fSFwSGY+DvgxsKFm2TWZubq8vbqm/SPAq4CDytvM55SkZcnksiS1oYg4C/gu0BcR10fEccApwDMi4ifA08vHkiRJ0rI120TXmfmlzLyrfHgxsP98z1Fe9bdXZl5clpf7OE6OLalNOKGfJLWhzDxmjkVPa2kgkiRJ0uL2l8Bnah4fGBEjwG3AWzPzm0AvcH3NOteXbffixNiSlhtHLkuSJEmSJM0QEW8B7gI+WTbdAKzMzDXAG4FPRcRejTynE2NLWm4cuSxJkiRJklQjIl4BPAd4Wlnqgsy8E7izvH9pRFwDPBoYY3rpjP3LNkla9hy5LEmSJEmSVIqII4F/AJ6Xmb+paV8RER3l/UdQTNx3bWbeANwWEU+MiABejpNjS2oTjlyWJEmSJEltqZzo+inAvhFxPXACsAG4P3BhkSvm4sx8NfBk4KSI2A7cDbw6M6cmA/xb4AygC7igvEnSsmdyWZIkSZIktaU5Jrr+2BzrngOcM8eyTcAhCxiaJC0JlsWQJEmSJEmSJDXM5LIkSZIkSZIkqWEmlyVJkiRJkiRJDTO5LEmSJEmSJElqmMllSZIkSZIkSVLDTC5LkiRJkiRJkhpmclmSJEmSJEmS1DCTy5IkSZIkSZKkhplcliRJkiRJkiQ1zOSyJEmSJEmSJKlhJpclSZIkSZIkSQ3bveoApPtq1fHnVR2CJEmSJEmS1HYcuSxJkiRJkiRJapjJZUmSJEmSJElSw0wuS5IkSZIkSZIaZnJZkiRJkiRJktQwk8uSJEmSJEmSpIaZXJYkSZIkSZIkNczksiRJkiRJkiSpYSaXJUmSJEmSJEkNM7ksSZIkSZIkSWqYyWVJkiRJkiRJUsNMLkuSJEmSJEmSGmZyWZIkSZIkSZLUMJPLkiRJkiRJkqSGNS25HBGnR8RNEXF5Tds+EXFhRPyk/HfvZu1fkiRJkiRJktQ8zRy5fAZw5Iy244GvZOZBwFfKx5IkSZIkSZKkJaZpyeXM/AZw84zm5wNnlvfPBNY1a/+SJEmSJEmSpObZvcX7e1hm3lDe/yXwsLlWjIj1wHqAlStXtiA0LaRVx59XdQhta1df++tOefYCRyJJkiRJkqTlrLIJ/TIzgZxn+amZuTYz165YsaKFkUmSJEmSJEmSdqbVI5dvjIj9MvOGiNgPuKnF+5ckSZIkSVpyhkbGGBweZev4BD3dXQz097FuTW/VYUlqc60euXwucGx5/1jg8y3evyRJkiRJ0pIyNDLGho1bGBufIIGx8Qk2bNzC0MhY1aFJanNNSy5HxFnAd4G+iLg+Io4DTgGeERE/AZ5ePpYkSZIkSdIcBodHmdg+Oa1tYvskg8OjFUUkSYWmlcXIzGPmWPS0Zu1TkiRJkiRpudk6PtFQuyS1SmUT+kmSJEmSJGnnerq7GmqXpFYxuSxJkiRJkrSIDfT30dXZMa2tq7ODgf6+iiKSpELTymJIkiRJkiTpvlu3phcoai9vHZ+gp7uLgf6+He2SVBWTy5IkSZIkSYvcujW9JpMlLTqWxZAkSZIkSZIkNczksiRJkiRJkiSpYSaXJUmSJEmSJEkNM7ksSZIkSZIkSWqYyWVJkiRJkiRJUsNMLkuSJEmSJEmSGmZyWZIkSZIkSZLUMJPLkiRJkiRJkqSGmVyWJEmSJEmSJDXM5LIkSZIkSZIkqWEmlyVJkiRJkiRJDTO5LEmSJEmSJElqmMllSZIkSZIkSVLDTC5LkiRJkiRJkhpmclmSJEmSJEmS1DCTy5IkSZIkSZKkhplcliRJkiRJkiQ1zOSyJEmSJEmSJKlhJpclSZIkSZIkSQ0zuSxJkiS1SETsERHfi4gfRsQVEfH2OdY7OiJ+VK7zqZr2yYjYXN7ObV3kkiRJ0r3tXnUAkiRJUhu5EzgiM++IiE7gWxFxQWZePLVCRBwEbAAOz8xbIuKhNdtPZObq1oYsSZIkzc7ksiRJktQimZnAHeXDzvKWM1Z7FfDhzLyl3Oam1kUoSZIk1c+yGJIkSVILRURHRGwGbgIuzMxLZqzyaODREfHtiLg4Io6sWbZHRGwq29e1KGRJkiRpVo5cliRJklooMyeB1RHRDXwuIg7JzMtrVtmd/9/evUdbUpZ3Hv/+pBERFVBaBkFscACHXETtEAxoUCIXMaLGKIwzYm4kjiZRI06TZCmaZBZKTEwyjhGViBlF0Qj2SBRQURMT5CY3ERSwjbQIKDeNRLk880e9BzaHc053nT5n3/r7WWuvU/ut2lXPW5e36jy76t2wB3AgsAvwxSQ/U1W3AU+oqvVJdgc+l+Tyqrp2cP5JjgGOAdh1112XvT6SJEnafJlc3oysWnNm78+sO+HwZYhEkiRJVXVbknOBQ4HB5PL1wJer6i7gm0m+TpdsvqCq1rfPXpfk88BTgGtnzfck4CSA1atXz+5yQ5IkSVoydoshSZIkDUmSle2OZZJsDTwHuGrWZGfQ3bVMkh3ousm4Lsn2SbYaKN8fuHIogUuSJElz8M5lSZIkaXh2Ak5JsgXdjR6nVdUnk7wFuLCq1gJnAQcnuRK4Bzi2qr6f5BeAdye5t332hKoyuSxJkqSRMbksSZIkDUlVXUbXlcXs8jcODBfwuvYanOZfgJ9Z7hglaXOS5GTgecBNVfXTrezRwEeAVcA64CVVdWuSAH8FPBf4EfCKqrq4feZo4I/bbP+0qk4ZZj0kaVTsFkOSJEmSJG2u3k/X9/2gNcBnq2oP4LPtPcBhdH3g70H3w6nvgvuS0W8Cfh7YF3hTku2XPXJJGgMmlyVJkiRJ0mapqr4I3DKr+Ahg5s7jU4AXDJR/oDrnAdsl2Qk4BDinqm6pqluBc3hwwlqSppLJZUmSJEmSpPvtWFU3tOHvAju24Z2Bbw9Md30rm6/8QZIck+TCJBfefPPNSxu1JI2AyWVJkiRJkqQ5tH7wawnnd1JVra6q1StXrlyq2UrSyJhcliRJkiRJut+NrbsL2t+bWvl64PED0+3SyuYrl6SpZ3JZkiRJkiTpfmuBo9vw0cAnBspfns5+wO2t+4yzgIOTbN9+yO/gViZJU2/FqAOQJEmSJEkahSSnAgcCOyS5HngTcAJwWpLfAL4FvKRN/o/Ac4FrgB8BvwZQVbck+RPggjbdW6pq9o8EStJUMrksSZIkSZI2S1V11DyjDppj2gJeNc98TgZOXsLQJGki2C2GJEmSJEmSJKk3k8uSJEmSJEmSpN5MLkuSJEmSJEmSejO5LEmSJEmSJEnqzeSyJEmSJEmSJKk3k8uSJEmSJEmSpN5MLkuSJEmSJEmSejO5LEmSJEmSJEnqzeSyJEmSJEmSJKk3k8uSJEmSJEmSpN5MLkuSJEmSJEmSejO5LEmSJEmSJEnqzeSyJEmSJEmSJKk3k8uSJEmSJEmSpN5MLkuSJEmSJEmSelsx6gAkSeMlyTrgB8A9wN1VtXq0EUmSJEmSpHFkclmSNJdnVdX3Rh2EJEmSJEkaX3aLIUmSJEmSJEnqzTuXJUmzFXB2kgLeXVUnzZ4gyTHAMQC77rrrkMOTJteqNWcOdXnrTjh8qMuTJEmStHnxzmVJ0mwHVNVTgcOAVyV55uwJquqkqlpdVatXrlw5/AglSZIkSdLImVyWJD1AVa1vf28CTgf2HW1EkiRJkiRpHNktxgQa9iO1kjYfSbYBHlJVP2jDBwNvGXFYkiRJkiRpDJlcliQN2hE4PQl054gPVdWnRxuSJEmSJEkaRyaXJUn3qarrgCePOg5JkiRJkjT+7HNZkiRJkiRJktSbyWVJkiRJkiRJUm8mlyVJkiRJkiRJvZlcliRJkiRJkiT1ZnJZkiRJkiRJktSbyWVJkiRJkiRJUm8mlyVJkiRJkiRJvZlcliRJkiRJkiT1ZnJZkiRJkiRJktSbyWVJkiRJkiRJUm8mlyVJkiRJkiRJva0YxUKTrAN+ANwD3F1Vq0cRhyRJkiRJkiRpcUaSXG6eVVXfG+HyJUmSJEmSJEmLZLcYkiRJkiRJkqTeRnXncgFnJyng3VV10uwJkhwDHAOw6667Di2wVWvOXNTn1p1w+NCWJS2HYe77kiRJkiRJmnyjSi4fUFXrkzwWOCfJVVX1xcEJWsL5JIDVq1fXKIKUJG2+/AJQkiRJkqSFjaRbjKpa3/7eBJwO7DuKOCRJkiRJkiRJizP05HKSbZI8cmYYOBi4YthxSJIkSZIkSZIWbxTdYuwInJ5kZvkfqqpPjyAOSZIkSZIkSdIiDT25XFXXAU8e9nIlSZIkSZIkSUtnJH0uS5IkSZIkSZImm8llSZIkSZIkSVJvJpclSZIkSZIkSb2ZXJYkSZIkSZIk9WZyWZIkSZIkSZLUm8llSZIkSZIkSVJvJpclSZIkSZIkSb2ZXJYkSZIkSZIk9WZyWZIkSZIkSZLUm8llSZIkSZIkSVJvJpclSZIkSZIkSb2tGHUA02LVmjNHHYIkSZIkSZIkDY13LkuSJEmSJEmSejO5LEmSJEmSJEnqzeSyJEmSJEmSJKk3k8uSJEmSJEmSpN78QT9JkiTdZ7E/UrzuhMOXOBJJkiRJ4847lyVJkiRJkiRJvZlcliRJkiRJkiT1ZnJZkiRJkiRJktSbyWVJkiRpSJI8LMn5SS5N8tUkb55nupckubJN86GB8qOTfKO9jh5e5JK0eUmyV5JLBl53JHlNkuOTrB8of+7AZ45Lck2Sq5McMsr4JWlY/EE/SZIkaXh+DDy7qn6YZEvgn5N8qqrOm5kgyR7AccD+VXVrkse28kcDbwJWAwVclGRtVd06/GpI0nSrqquBfQCSbAGsB04Hfg34y6r688Hpk+wNHAn8FPA44DNJ9qyqe4YZtyQNm3cuS5IkSUNSnR+2t1u2V82a7LeAd84kjavqplZ+CHBOVd3Sxp0DHDqEsCVpc3cQcG1VfWuBaY4APlxVP66qbwLXAPsOJTpJGiGTy5IkSdIQJdkiySXATXTJ4i/PmmRPYM8kX0pyXpKZBPLOwLcHpru+lc2e/zFJLkxy4c0337wMNZCkzc6RwKkD71+d5LIkJyfZvpXZRkvaLJlcliRJkoaoqu6pqn2AXYB9k/z0rElWAHsABwJHAe9Jsl2P+Z9UVauravXKlSuXJmhJ2kwleSjwfOCjrehdwBPpusy4AXh7n/nZRkuaNiaXJUmSpBGoqtuAc3lw1xbXA2ur6q72aPXX6ZLN64HHD0y3SyuTJC2fw4CLq+pGgKq6sX1JeC/wHu7v+sI2WtJmyeSyJEmSNCRJVs7chZxka+A5wFWzJjuD7q5lkuxA103GdcBZwMFJtm+PYR/cyiRJy+coBrrESLLTwLgXAle04bXAkUm2SrIb3ZeC5w8tSkkakRWjDkDjbdWaM0cdgiRJ0jTZCTglyRZ0N3qcVlWfTPIW4MKqWsv9SeQrgXuAY6vq+wBJ/gS4oM3rLVV1y/CrIEmbhyTb0H0J+NsDxW9Lsg/dj7GumxlXVV9NchpwJXA38KqqumeoAUvSCJhcliRJkoakqi4DnjJH+RsHhgt4XXvNnu5k4OTljFGS1KmqfwceM6vsvy8w/Z8Bf7bccUnSOLFbDEmSJEmSJElSbyaXJUmSJEmSJEm9mVyWJEmSJEmSJPVmclmSJEmSJEmS1JvJZUmSJEmSJElSbyaXJUmSJEmSJEm9mVyWJEmSJEmSJPW2YtQBSJI2T6vWnDm0Za074fChLUtSP4tpCzymJUmSpPHgncuSJEmSJEmSpN5MLkuSJEmSJEmSejO5LEmSJEmSJEnqzeSyJEmSJEmSJKk3k8uSJEmSJEmSpN5MLkuSJEmSJEmSejO5LEmSJEmSJEnqbcWoA1hOq9acOeoQJEmSJEmSJGkqTXVyWZIkSdLGO+Mr6znxrKv5zm138rjttubYQ/biBU/ZedRhSZKwjZY0nkwuS5IkSeKMr6znuI9fzp133QPA+tvu5LiPXw5g8kKSRsw2WtK4ss9lSZIkSZx41tX3JS1m3HnXPZx41tUjikiSNMM2WtK4MrksSZIkie/cdmevcknS8NhGSxpXJpclSZIk8bjttu5VLkkaHttoSePK5LIkSZIkjj1kL7becosHlG295RYce8heI4pIkjTDNlrSuPIH/SRJkiTd94NQJ551Nd+57U4et93WHHvIXv5QlCSNAdtoSePK5LIkSZIkoEtemKiQpPFkGy1pHNkthiRJkiRJkiSpN5PLkiRJkiRJkqTeTC5LkiRJkiRJknozuSxJkiRJkiRJ6s0f9JMkTb1Va84cdQjSSLjvS5IkSVpO3rksSZIkSZIkSerN5LIkSZIkSZIkqTeTy5IkSZIkSZKk3uxzWdImWUx/nutOOHxoy1qsxcYoSZIkSZK0ufDOZUmSJEmSJElSbyaXJUmSJEmSJEm9mVyWJEmSJEmSJPVmclmSJEmSJEmS1JvJZUmSJEmSJElSbyaXJUmSJEmSJEm9mVyWJEmSJEmSJPVmclmSJEmSJEmS1JvJZUmSJEmSJElSbyaXJUmSJEmSJEm9mVyWJEmSJEmSJPVmclmSJEmSJEmS1JvJZUmSJEmSJElSbyNJLic5NMnVSa5JsmYUMUiS5mYbLUmSJEmSNsbQk8tJtgDeCRwG7A0clWTvYcchSXow22hJkiRJkrSxRnHn8r7ANVV1XVX9BPgwcMQI4pAkPZhttCRJkiRJ2igrRrDMnYFvD7y/Hvj52RMlOQY4pr39YZKrZ02yA/C9ZYlwfEx7Hae9fmAd55S3LlMkS2ggxr71e8KSBzNcS9VGD9skHWuTFCtMVrzGunw2GO8w2/YNLGuhWCe9je7toosu+l6Sby3zYsZlfzaOBxqHOMYhBjCO2cY1DtvojTcu23CpTWu9YHrrNq31gumt22LrNW8bPYrk8kapqpOAk+Ybn+TCqlo9xJCGbtrrOO31A+s4Daa9fou1oTZ62CZpO01SrDBZ8Rrr8pmkeCcp1mGoqpXLvYxxWefGMX5xjEMMxmEc42yxbfS0rrtprRdMb92mtV4wvXVbjnqNoluM9cDjB97v0sokSaNnGy1JkiRJkjbKKJLLFwB7JNktyUOBI4G1I4hDkvRgttGSJEmSJGmjDL1bjKq6O8mrgbOALYCTq+qri5jV2DyOvYymvY7TXj+wjtNg2uv3AEvYRg/bJG2nSYoVJiteY10+kxTvJMU6LcZlnRvHA41DHOMQAxjHbMYx+aZ13U1rvWB66zat9YLprduS1ytVtdTzlCRJkiRJkiRNuVF0iyFJkiRJkiRJmnAmlyVJkiRJkiRJvU1kcjnJoUmuTnJNkjWjjmdjJTk5yU1Jrhgoe3SSc5J8o/3dvpUnyV+3Ol6W5KkDnzm6Tf+NJEePoi7zSfL4JOcmuTLJV5P8fiufinomeViS85Nc2ur35la+W5Ivt3p8pP0QGkm2au+vaeNXDczruFZ+dZJDRlSleSXZIslXknyyvZ+qOiZZl+TyJJckubCVTcV+Og0WaEuOT7K+bbdLkjy3la9KcudA+d8OzOtpbVtf07ZjhhFrG/e7Sa5q5W8bKJ/z2MgQzm994x3HddvanJl41iW5ZOAzI1m3fWMd5XrdQLz7JDlvpm1Msm8rH1k7uIhYD0xy+8C6fePAvCbyGnKUkuw1sC4vSXJHktdknva4fWZJzv8Zk2vneeI4MV17eVmS05Ns18qX7dieJ47e22FTj4N54hhqW7dAuzDU/WOBOIa6fywQx1D3jwXiGMtz4STa1ON31Jby2B1HWYL/ocdRku2SfKy1a19L8vRp2GZJXtv2wyuSnJou3zOR2yyjvmaqqol60f3A1LXA7sBDgUuBvUcd10bG/kzgqcAVA2VvA9a04TXAW9vwc4FPAQH2A77cyh8NXNf+bt+Gtx913QbqsxPw1Db8SODrwN7TUs8W5yPa8JbAl1vcpwFHtvK/BV7Zhv8H8Ldt+EjgI21477bvbgXs1vbpLUZdv1l1fR3wIeCT7f1U1RFYB+wwq2wq9tNpeC3QlhwPvH6O6Vcx0LbOGnd+225p2/GwIcX6LOAzwFZt3GPb3zmPDYZ0fltEvGO3bmdN83bgjaNet4uIdWTrdQP7wdkzy6Nr+z4/MDySdnARsR5IO3fNms/EXkOOy6utw+8CT2D+9njJzv+MybXzPHEcDKxow28diGPZju154ui1HZbiOJgrjlnjl72tW6BdGOr+sUAcQ90/FohjqPvHfHEMe/+Y1tdSHL+jfi3VsTuuLzbxf+hxfQGnAL/Zhh8KbDfp2wzYGfgmsPXAtnrFpG4zRnzNNIl3Lu8LXFNV11XVT4APA0eMOKaNUlVfBG6ZVXwE3YFK+/uCgfIPVOc8YLskOwGHAOdU1S1VdStwDnDosge/karqhqq6uA3/APga3UE7FfVscf6wvd2yvQp4NvCxVj67fjP1/hhwUPvW/Qjgw1X146r6JnAN3b49FpLsAhwOvLe9D1NWx3lMxX46DRZoS3pp2+lRVXVedWfMD3D/dl3uWF8JnFBVP27jbmofme/YGMr5bRHxzmnE63YmhgAvAU5tRSNbt4uIdU7DWK8biLeAR7XJtgW+04ZH1g4uItb5TOw15Bg5CLi2qr61wDRLdv4fl2vnueKoqrOr6u729jxgl4XmsRTH9jzrYz7L1h4uFMew2rpx+Z9jvjiGvX8s4rppWfaPSTsXTqCJP48t4bE7dpbof+ixk2RbusTl+wCq6idVdRtTsM2AFcDWSVYADwduYEK32aivmSYxubwz8O2B99eziITDGNmxqm5ow98FdmzD89VzYurfHhN4Ct3dvVNTz/aoyyXATXQH27XAbQMXkIOx3lePNv524DGMcf2adwBvAO5t7x/D9NWxgLOTXJTkmFY2NfvpNJnVlgC8uj2+c/LMoz3Nbu0xtC8keUYr25luu8xY1m00K9Y9gWe0x6a+kOTnBmIai/1pI+OF8Vu3M54B3FhV3xiIaeTrdiNjhTFYr3PE+xrgxCTfBv4cOG4grnFbt/PFCvD0dF1YfSrJT7Uy2+xNdyQPTArN1R4v93oex3P1r9PdATRj2Md2n+2w3Otj6G3duPzPMU/bD0PePzbyumlU62Nsz4UTYqrOY5t47I6jd7Dp/0OPo92Am4G/a8fqe5Nsw4Rvs6paT3f9+G90SeXbgYuYjm02Y2jnxElMLk+t9q1sjTqOpZDkEcA/AK+pqjsGx016Pavqnqrah+4OhH2BJ402oqWV5HnATVV10ahjWWYHVNVTgcOAVyV55uDISd9Pp8Ucbcm7gCcC+9BdBLy9TXoDsGtVPYX2OFqSRz14jkONdQXdI0X7AccCp43Tt9s94h3HdTvjKDZw99Ow9Yh15OsV5oz3lcBrq+rxwGtpd6mMgx6xXgw8oaqeDPwNcMYIwp066focfD7w0VY0X3s8NONwrk7yR8DdwAdb0bCP7ZFvh1mG2taNy/8c88Ux7P2jx3XTspq0c6GGb1yO3aUy5f9Dr6DrbuFd7Vj9d7ouFu4zodtse7o7eHcDHgdswxQ/hbzc22gSk8vrgccPvN+llU2qG2ceEWh/Zx5Dnq+eY1//JFvSnSg+WFUfb8VTV8/2KMi5wNPpHiNY0UYNxnpfPdr4bYHvM9712x94fpJ1dI9aPRv4K6arjjPfVM48+n863RcFU7efTrK52pKqurF9wXMv8B7ao9btsc7vt+GL6J4o2JNueww+iros22iedu964OPtcaPz6e5i2IEx2J/6xDum63amvXkR8JGByUe6bvvEOur1ukC8RwMzwx/l/u4MxnHdzhlrVd1RrQurqvpHYMskCx172jiHARdX1Y0wf3vM8q/nsTlXJ3kF8DzgZe2ftqEf24vYDsu5Poba1o3L/xwLtP2vYIj7R5/rJkazPsbyXDhhpuI8tkTH7rhZqv+hx9H1wPVVNfMUwsfoks2Tvs1+CfhmVd1cVXfRXVPuz3RssxlDOydOYnL5AmCPdL/g+FC6x/PWjjimTbGW7p8j2t9PDJS/PJ39gNvb7exnAQcn2b5903JwKxsL7S639wFfq6q/GBg1FfVMsjL3/9rz1sBz6PqJOhd4cZtsdv1m6v1i4HPt4nItcGS6XxzdDdiD7kcrRq6qjquqXapqFd3x9bmqehlTVMck2yR55Mww3f51BVOyn06D+dqSPLC/rhfSbbeZY3OLNrw73f52XdtOdyTZr83z5dy/XZc1Vrq7JZ/VptmT7scvvsf8x8ZQzm994x3TdQvdBeFVVTX42OzI1m3fWEe5XjcQ73eAX2zDzwZmHl0eWTvYN9Yk/6l9hiT70l3vfp/pu4YctgfccThfe8zyn//H4lyd5FC6x5+fX1U/Gigf6rG9iO2wnMfB0Nq6cfmfY4HrlaHuH32vm1im/WPSzoUTaOLPY0t47I6VJfwfeuxU1XeBbyfZqxUdBFzJhG8zuu4w9kvy8LZfztRr4rfZgOGdE2sMftWw74vulw2/TvfN5h+NOp4ecZ9K9+jPXXTf/vwGXR8tn6X7h+gzwKPbtAHe2ep4ObB6YD6/TvejC9cAvzbqes2q4wF0t9pfBlzSXs+dlnoCPwt8pdXvCu7/pePd6S7IrqG7e2qrVv6w9v6aNn73gXn9Uav31Yzprx8DB3L/L91OTR1bXS5tr6/OtCPTsp9Ow2uBtuTv2za4jO6kuFOb/lfatryE7pH4Xx6Y1+p2vF4L/G8gQ4r1ocD/bcu+GHj2wGfmPDYYwvmtb7zjuG7buPcDvzPHZ0aybvvGOsr1uoH94AC6/uYupesD8Wlt+pG1g4uI9dVt3V5K90NavzDMY2waX3SPin4f2HagbM72uI1bkvM/Y3LtPE8c19D1SzizT878gvyyHdvzxNF7O2zqcTBXHK38/QyprWNM/udYII6h7h8LxDHU/WO+OIa9f0zza1O2zzi8lvLYHdcXm/g/9Di+6LrWubBttzOA7adhmwFvBq5qbc3fA1tN6jZjxNdMaR+WJEmSJEmSJGmjTWK3GJIkSZIkSZKkETO5LEmSJEmSJEnqzeSyJEmSJEmSJKk3k8uSJEmSJEmSpN5MLkuSJEmSJEmSejO5rCWXpJK8feD965Mcv0Tzfn+SFy/FvDawnF9N8rUk584qX5XkziSXJLk0yb8k2auNW53kr9vw8Ulev9xxSlJfttG20ZI2T0l+uAzzXJdkh1EsW5KmRZJ/aX9XJfmvSzzvP5xrWdJSMrms5fBj4EUbc6E5TElW9Jj8N4DfqqpnzTHu2qrap6qeDJwC/CFAVV1YVb+3BHFusanzkKQF2EZvWpy20ZIkSVoyVfULbXAV0Cu5vBHX0A9ILg8sS1oyJpe1HO4GTgJeO3vE7LvaZu5iSHJgki8k+USS65KckORlSc5PcnmSJw7M5peSXJjk60me1z6/RZITk1yQ5LIkvz0w339Ksha4co54jmrzvyLJW1vZG4EDgPclOXEDdX0UcOvAsj45xzJ+K8mnkmyd5L+1Ol2S5N0zSYokP0zy9iSXAk9v9b+y1eXPNxCDJPVhG/3AZdhGS9psJfnlJF9O8pUkn0myYys/PskprY3+VpIXJXlba5M/nWTLgdm8oZWfn+Q/t8/vluRfW/mfDizvEUk+m+TiNu6IIVdZksZO7n+64wTgGe1a9LV9rqGTnJHkoiRfTXJMKzsB2LrN74ODy0rnxHadfXmSlw7M+/NJPpbkqiQfTJLhrhFNmj53CUl9vBO4LMnbenzmycB/AW4BrgPeW1X7Jvl94HeB17TpVgH7Ak8Ezm0XsS8Hbq+qn0uyFfClJGe36Z8K/HRVfXNwYUkeB7wVeBpd8uHsJC+oqrckeTbw+qq6cI44n5jkEuCRwMOBn5+vQkleDTwHeAGwO/BSYP+quivJ/wFeBnwA2Ab4clX9QZLHAO8DnlRVlWS7Da86SerFNhrbaEkC/hnYr7Vnvwm8AfiDNu6JwLOAvYF/BX6lqt6Q5HTgcOCMNt3tVfUzSV4OvAN4HvBXwLuq6gNJXjWwvP8AXlhVd6R7gua8JGurqpa3mpI0EdbQXePO3KBxDBt/Df3rVXVLkq2BC5L8Q1WtSfLqqtpnjmW9CNiH7hp/h/aZL7ZxTwF+CvgO8CVgf7rzhTQnk8taFu2C8QPA7wF3buTHLqiqGwCSXAvMNJqX013Yzjitqu4FvpHkOuBJwMHAz+b+O+62BfYAfgKcPztp0fwc8Pmqurkt84PAM7n/Qnk+1840zu3bvZOAQ+eY7uXAt4EXtETFQXRJkgvaF39bAze1ae8B/qEN30534f2+dpfdg+60k6RNYRsN2EZLEsAuwEeS7AQ8FBhsjz/V2sfLgS2AT7fyy+m+SJxx6sDfv2zD+wO/0ob/nu7LQoAA/yvJM4F7gZ2BHYHvLlWFJGmK9LmG/r0kL2zDj2/TfX+BeR8AnFpV9wA3JvkC3fX3HW3e1wO0mzZWYXJZCzC5rOX0DuBi4O8Gyu6mdceS5CF0F7EzfjwwfO/A+3t54L46+86GortQ/d2qOmtwRJIDgX9fTPAbaS0PrN+gy+m+CdyF7kI9wClVddwc0/5Ha9SpqruT7AscBLwYeDXw7CWOW5LegW30PthGS9q8/Q3wF1W1trXJxw+M+zFAVd2b5K6Bu4sXavfnG57xMmAl8LSWuF4HPGxTKiBJU2yjrqHb+18Cnl5VP0ryeTatbR287r8Hc4faAPtc1rKpqluA0+h+eGnGOro7wwCeD2xJf7+a5CHp+vjcHbgaOAt45Uz/b0n2TLLNBuZzPvCLSXZI16/mUcAXesZyAHDtPOO+Avw2sLY93v1Z4MVJHttifHSSJ8z+UJJHANtW1T/S9Yn65J4xSdIG2UbbRksS3V1w69vw0Yucx0sH/v5rG/4ScGQbftms5d3UEsvPAh7UzkrSZuwHdF27zdjYa+htgVtbYvlJwH4D4+7KA/vJn/FPwEvT9eu8ku4JwfOXpBba7Pjtg5bb2+nu6prxHuAT6X4U6dMs7o61f6Nr9B4F/E5V/UeS99I9qnFx62z+Zro+NOdVVTckWQOcS/eN4JlV9YmNWP5Mf56hexzlNxdYxj8neT1wJl2/nn9M12/oQ4C7gFcB35r1sUfSraOHtWW8biNikqTFsI22jZa0+Xh4kusH3v8F3Z3KH01yK/A5YLdFzHf7JJfR3el2VCv7feBDSf4nMNh2fxD4f62rjQuBqxaxPEmaVpcB97Rr8ffT9V+/ig1fQ38a+J0kX6O7seO8gXEn0f3WysVVNfhl3+nA04FL6Z40eUNVfbclp6Ve4m8nSJIkSZIkSZL6slsMSZIkSZIkSVJvJpclSZIkSZIkSb2ZXJYkSZIkSZIk9WZyWZIkSZIkSZLUm8llSZIkSZIkSVJvJpclSZIkSZIkSb2ZXJYkSZIkSZIk9fb/AS6v/bujpzG5AAAAAElFTkSuQmCC\n"
61
+ },
62
+ "metadata": {
63
+ "needs_background": "light"
64
+ }
65
+ }
66
+ ],
67
+ "source": [
68
+ "\n",
69
+ "import pandas as pd\n",
70
+ "import numpy as np\n",
71
+ "import matplotlib.pyplot as plt\n",
72
+ "import scipy\n",
73
+ "from scipy.stats import poisson\n",
74
+ "\n",
75
+ "\n",
76
+ "#estimationg for mu using maximum log likelihood derivation, from reference 2\n",
77
+ "def getMu(x):\n",
78
+ " sigX=0\n",
79
+ " n=0\n",
80
+ " for i in x:\n",
81
+ " sigX+=i\n",
82
+ " n+=1\n",
83
+ " return sigX/n\n",
84
+ "\n",
85
+ "#factorial function\n",
86
+ "def factorial(x):\n",
87
+ " sum=0\n",
88
+ " for i in range(x):\n",
89
+ " sum+=i\n",
90
+ " return sum\n",
91
+ "\n",
92
+ "#calculates log likelihood\n",
93
+ "def LLK(X, mu):\n",
94
+ " llk=0\n",
95
+ " for x in X:\n",
96
+ " llk+=np.log(mu)*x-mu-np.log((x))\n",
97
+ " return llk\n",
98
+ "\n",
99
+ "#calculates the partial derivitave of the log likelihood function with respect to lambda (LOSS FUNCTION)\n",
100
+ "def LLK_Gradient(X,mu):\n",
101
+ " sum=0;\n",
102
+ " n=0\n",
103
+ " for x in X:\n",
104
+ " sum+=x\n",
105
+ " n+=1\n",
106
+ " return (1/mu)*sum-n\n",
107
+ "\n",
108
+ "#iterative method for finding lambda*\n",
109
+ "def iterative(X):\n",
110
+ " llks=[]\n",
111
+ " mus=[]\n",
112
+ " for i in range(X.min(),X.max(),4):\n",
113
+ " llks.append(LLK(X,i))\n",
114
+ " mus.append(i)\n",
115
+ " return mus,llks\n",
116
+ "\n",
117
+ "#Gradient descent for finding lambda*\n",
118
+ "def GradDesc(X,eta,iter):\n",
119
+ " iteration=[]\n",
120
+ " mu_set=[]\n",
121
+ " mu=np.random.randint(X.min(),X.max())\n",
122
+ " for i in range(iter):\n",
123
+ " iteration.append(i)\n",
124
+ " mu_set.append(mu)\n",
125
+ " #print(LLK_Gradient(X,mu))\n",
126
+ " mu=mu+eta*LLK_Gradient(X,mu)\n",
127
+ " return mu, mu_set,iteration\n",
128
+ "\n",
129
+ "\n",
130
+ "\n",
131
+ "\n",
132
+ "data=pd.read_csv(\"./nyc_bb_bicyclist_counts.csv\")\n",
133
+ "bb_count=data.loc[:,\"BB_COUNT\"].values\n",
134
+ "\n",
135
+ "\n",
136
+ "\n",
137
+ "\n",
138
+ "\n",
139
+ "\n",
140
+ "#print(iterative(bb_count))\n",
141
+ "lam, lam_set,iteration=GradDesc(bb_count,0.7,1000)\n",
142
+ "lamLLK_set=[]\n",
143
+ "for mu in lam_set:\n",
144
+ " lamLLK_set.append(LLK(bb_count,mu))\n",
145
+ "\n",
146
+ "simulated=[]\n",
147
+ "for i in range(len(bb_count)):\n",
148
+ " simulated.append(np.random.poisson(lam))\n",
149
+ "\n",
150
+ "mu_set, muLLK_set=iterative(bb_count)\n",
151
+ "\n",
152
+ "print(\"lambda* from gradient descent : \",lam)\n",
153
+ "print(\"Log Likelyhood for Lambda* : \",(LLK(bb_count,getMu(bb_count))))\n",
154
+ "\n",
155
+ "fig, axis= plt.subplots(1, 4, figsize=(20, 10))\n",
156
+ "\n",
157
+ "axis[0].hist(bb_count,20)\n",
158
+ "axis[0].set_title(\"Actual Biker Dataset\")\n",
159
+ "axis[0].set_xlabel(\"Number of Bikers\")\n",
160
+ "axis[0].set_ylabel(\"Number of Days\")\n",
161
+ "\n",
162
+ "axis[1].hist(simulated,20)\n",
163
+ "axis[1].set_title(\"Simulated Dataset using Lambda*\")\n",
164
+ "axis[1].set_xlabel(\"Number of Bikers\")\n",
165
+ "axis[1].set_ylabel(\"Number of Days\")\n",
166
+ "\n",
167
+ "axis[2].scatter(lam_set,lamLLK_set)\n",
168
+ "axis[2].set_title(\"Lambda vs Log Likelihood Gradient Descent\")\n",
169
+ "axis[2].set_xlabel(\"Lambda\")\n",
170
+ "axis[2].set_ylabel(\"Log Likelyhood (in millions)\")\n",
171
+ "\n",
172
+ "axis[3].scatter(iteration,lam_set)\n",
173
+ "axis[3].set_title(\"Lamda over iterations\")\n",
174
+ "axis[3].set_xlabel(\"iteration\")\n",
175
+ "axis[3].set_ylabel(\"lamda\")\n",
176
+ "\n",
177
+ "fig.tight_layout()\n",
178
+ "\n",
179
+ "\n",
180
+ " \n",
181
+ "plt.show()\n",
182
+ "\n"
183
+ ]
184
+ },
185
+ {
186
+ "cell_type": "markdown",
187
+ "source": [
188
+ "The first figure on the left is a histogram of the provided dataset. The second fingure, second from left, shows a simulated dataset using 𝝺* calculated from the gradient descent, and numpy.random.poisson(). The second figure has less outliers, and forms a more uniform distrobution, but maintains a similar shape to the acutal dataset. \n",
189
+ "\n",
190
+ "The 3rd figure, second from the right, is a scatter plot showing the different log-likelihoods for all lambdas genorated during the gradient descent. The first lambda used is the lowest point on the graph, and all subsequent points move towards the true mean, 𝝺*, and gradually increase in log likelihood. The final point with the highest log likelihood is close to 𝝺*, with approximatly the maximum log likelihood.\n",
191
+ " The rightmost figure is another scatterplot which shows the progression of lamda over the iterations of the gradient descent. The values for lamda converger quickly to the approximate value of 𝝺*."
192
+ ],
193
+ "metadata": {
194
+ "id": "7WZ8JQgwAUnS"
195
+ }
196
+ },
197
+ {
198
+ "cell_type": "markdown",
199
+ "metadata": {
200
+ "id": "fbhG8Vs9NfDe"
201
+ },
202
+ "source": [
203
+ "## Maximum Likelihood II\n",
204
+ "\n",
205
+ "A colleague of yours suggest that the parameter $\\lambda$ must be itself dependent on the weather and other factors since people bike when its not raining. Assume that you model $\\lambda$ as \n",
206
+ "\n",
207
+ "$$\\lambda_i = \\exp(\\mathbf w^T \\mathbf x_i)$$\n",
208
+ "\n",
209
+ "where $\\mathbf x_i$ is one of the example features and $\\mathbf w$ is a set of parameters. \n",
210
+ "\n",
211
+ "Train the model with SGD with this assumption and compare the MSE of the predictions with the `Maximum Likelihood I` approach. \n",
212
+ "\n",
213
+ "You may want to use [this partial derivative of the log likelihood function](http://home.cc.umanitoba.ca/~godwinrt/7010/poissonregression.pdf)"
214
+ ]
215
+ },
216
+ {
217
+ "cell_type": "code",
218
+ "execution_count": 422,
219
+ "metadata": {
220
+ "id": "FoWbZHpfNfDe",
221
+ "colab": {
222
+ "base_uri": "https://localhost:8080/"
223
+ },
224
+ "outputId": "4c5b26b9-1484-48ed-f4b9-52f077427b64"
225
+ },
226
+ "outputs": [
227
+ {
228
+ "output_type": "stream",
229
+ "name": "stdout",
230
+ "text": [
231
+ "Mean squared error for method 1: 692017.6584942794\n",
232
+ "Mean squared error for method 2: 3000035.1120429407\n"
233
+ ]
234
+ }
235
+ ],
236
+ "source": [
237
+ "from numpy.core.multiarray import dot\n",
238
+ "data=pd.read_csv(\"./nyc_bb_bicyclist_counts.csv\")\n",
239
+ "y=data.loc[:,\"BB_COUNT\"].values\n",
240
+ "precip=data.loc[:,\"PRECIP\"].values\n",
241
+ "high=data.loc[:,\"HIGH_T\"].values\n",
242
+ "low=data.loc[:,\"LOW_T\"].values\n",
243
+ "m=len(y)\n",
244
+ "\n",
245
+ "n_epochs=500\n",
246
+ "np.random.seed(42)\n",
247
+ "\n",
248
+ "def makeX(precip,high,low,size):\n",
249
+ " X=[]\n",
250
+ " for i in range(size):\n",
251
+ " X.append([precip[i],high[i]/100,low[i]/100])\n",
252
+ " return np.array(X)\n",
253
+ "\n",
254
+ "\n",
255
+ "X=makeX(precip,high,low,m)\n",
256
+ "\n",
257
+ "w=[1,1,1]\n",
258
+ "\n",
259
+ "#def getGrads(X,y,w,m):\n",
260
+ "# gradients=[0]\n",
261
+ "# for i in range(m):\n",
262
+ "# gradients=gradients+(y[i]-np.exp(w @ X[i].T))*X[i]\n",
263
+ "# return gradients\n",
264
+ "\n",
265
+ "for epoch in range(n_epochs):\n",
266
+ " for iteration in range(m):\n",
267
+ " random_index=np.random.randint(m)\n",
268
+ " xi=X[random_index:random_index+1]\n",
269
+ " yi=y[random_index:random_index+1]\n",
270
+ " #print(yi,xi,w, \"yi-xi-w\")\n",
271
+ " #print(np.dot(w,xi.T),\"wTx\")\n",
272
+ " #print(np.exp(np.dot(w,xi.T)), \"yhat\")\n",
273
+ " gradients=(yi-np.exp(np.dot(w,xi.T)))*xi\n",
274
+ " #gradients=getGrads(X,y,w,m)\n",
275
+ " #print(gradients,\"gradient\")\n",
276
+ " eta=8*10**-5\n",
277
+ " w=w+eta*gradients\n",
278
+ " #print(w,\"new w\")\n",
279
+ "\n",
280
+ "#print(w)\n",
281
+ "\n",
282
+ "lamdas=[]\n",
283
+ "for i in range(m):\n",
284
+ " lamdas.append(np.exp(np.dot(w,X[i])))\n",
285
+ " #print(lamdas[i][0],y[i])\n",
286
+ "\n",
287
+ "\n",
288
+ "\n",
289
+ "MSELLK1 = 0\n",
290
+ "MSELLK2 = 0\n",
291
+ "for i in range(m):\n",
292
+ " MSELLK1+=(lam_set[i]-y[i])**2\n",
293
+ " MSELLK2+=(lamdas[i][0]-y[i])**2\n",
294
+ "MSELLK1=MSELLK1/m\n",
295
+ "MSELLK2=MSELLK2/m\n",
296
+ "\n",
297
+ "print(\"Mean squared error for method 1: \",MSELLK1)\n",
298
+ "print(\"Mean squared error for method 2: \",MSELLK2)"
299
+ ]
300
+ },
301
+ {
302
+ "cell_type": "markdown",
303
+ "source": [
304
+ "The mean squared error is much lower in the first method than the second one. Furthermore, the mean squared error is noticable very high in lambda set using the second method. This could be indicitave of a flaw in the GDS algorithm which optimizes our parameter set, w. It could also indicated that the weather has very little effect on the number of cyclists crossing the bridge. I beleive that the former is the case, and specifically that the gradient function is not working correctly. "
305
+ ],
306
+ "metadata": {
307
+ "id": "fthnVmfDh_2R"
308
+ }
309
+ }
310
+ ],
311
+ "metadata": {
312
+ "kernelspec": {
313
+ "display_name": "Python 3",
314
+ "language": "python",
315
+ "name": "python3"
316
+ },
317
+ "language_info": {
318
+ "name": "python",
319
+ "version": "3.10.9"
320
+ },
321
+ "orig_nbformat": 4,
322
+ "vscode": {
323
+ "interpreter": {
324
+ "hash": "7d6993cb2f9ce9a59d5d7380609d9cb5192a9dedd2735a011418ad9e827eb538"
325
+ }
326
+ },
327
+ "colab": {
328
+ "provenance": []
329
+ }
330
+ },
331
+ "nbformat": 4,
332
+ "nbformat_minor": 0
333
+ }
Assignment2/nyc_bb_bicyclist_counts.csv ADDED
@@ -0,0 +1,215 @@
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
1
+ Date,HIGH_T,LOW_T,PRECIP,BB_COUNT
2
+ 1-Apr-17,46.00,37.00,0.00,606
3
+ 2-Apr-17,62.10,41.00,0.00,2021
4
+ 3-Apr-17,63.00,50.00,0.03,2470
5
+ 4-Apr-17,51.10,46.00,1.18,723
6
+ 5-Apr-17,63.00,46.00,0.00,2807
7
+ 6-Apr-17,48.90,41.00,0.73,461
8
+ 7-Apr-17,48.00,43.00,0.01,1222
9
+ 8-Apr-17,55.90,39.90,0.00,1674
10
+ 9-Apr-17,66.00,45.00,0.00,2375
11
+ 10-Apr-17,73.90,55.00,0.00,3324
12
+ 11-Apr-17,80.10,62.10,0.00,3887
13
+ 12-Apr-17,73.90,57.90,0.02,2565
14
+ 13-Apr-17,64.00,48.90,0.00,3353
15
+ 14-Apr-17,64.90,48.90,0.00,2942
16
+ 15-Apr-17,64.90,52.00,0.00,2253
17
+ 16-Apr-17,84.90,62.10,0.01,2877
18
+ 17-Apr-17,73.90,64.00,0.01,3152
19
+ 18-Apr-17,66.00,50.00,0.00,3415
20
+ 19-Apr-17,52.00,45.00,0.01,1965
21
+ 20-Apr-17,64.90,50.00,0.17,1567
22
+ 21-Apr-17,53.10,48.00,0.29,1426
23
+ 22-Apr-17,55.90,52.00,0.11,1318
24
+ 23-Apr-17,64.90,46.90,0.00,2520
25
+ 24-Apr-17,60.10,50.00,0.01,2544
26
+ 25-Apr-17,54.00,50.00,0.91,611
27
+ 26-Apr-17,59.00,54.00,0.34,1247
28
+ 27-Apr-17,68.00,59.00,0.00,2959
29
+ 28-Apr-17,82.90,57.90,0.00,3679
30
+ 29-Apr-17,84.00,64.00,0.06,3315
31
+ 30-Apr-17,64.00,54.00,0.00,2225
32
+ 1-May-17,72.00,50.00,0.00,3084
33
+ 2-May-17,73.90,66.90,0.00,3423
34
+ 3-May-17,64.90,57.90,0.00,3342
35
+ 4-May-17,63.00,50.00,0.00,3019
36
+ 5-May-17,59.00,52.00,3.02,513
37
+ 6-May-17,64.90,57.00,0.18,1892
38
+ 7-May-17,54.00,48.90,0.01,3539
39
+ 8-May-17,57.00,45.00,0.00,2886
40
+ 9-May-17,61.00,48.00,0.00,2718
41
+ 10-May-17,70.00,51.10,0.00,2810
42
+ 11-May-17,61.00,51.80,0.00,2657
43
+ 12-May-17,62.10,51.10,0.00,2640
44
+ 13-May-17,51.10,45.00,1.31,151
45
+ 14-May-17,64.90,46.00,0.02,1452
46
+ 15-May-17,66.90,55.90,0.00,2685
47
+ 16-May-17,78.10,57.90,0.00,3666
48
+ 17-May-17,90.00,66.00,0.00,3535
49
+ 18-May-17,91.90,75.00,0.00,3190
50
+ 19-May-17,90.00,75.90,0.00,2952
51
+ 20-May-17,64.00,55.90,0.01,2161
52
+ 21-May-17,66.90,55.00,0.00,2612
53
+ 22-May-17,61.00,54.00,0.59,768
54
+ 23-May-17,68.00,57.90,0.00,3174
55
+ 24-May-17,66.90,57.00,0.04,2969
56
+ 25-May-17,57.90,55.90,0.58,488
57
+ 26-May-17,73.00,55.90,0.10,2590
58
+ 27-May-17,71.10,61.00,0.00,2609
59
+ 28-May-17,71.10,59.00,0.00,2640
60
+ 29-May-17,57.90,55.90,0.13,836
61
+ 30-May-17,59.00,55.90,0.06,2301
62
+ 31-May-17,75.00,57.90,0.03,2689
63
+ 1-Jun-17,78.10,62.10,0.00,3468
64
+ 2-Jun-17,73.90,60.10,0.01,3271
65
+ 3-Jun-17,72.00,55.00,0.01,2589
66
+ 4-Jun-17,68.00,60.10,0.09,1805
67
+ 5-Jun-17,66.90,60.10,0.02,2171
68
+ 6-Jun-17,55.90,53.10,0.06,1193
69
+ 7-Jun-17,66.90,54.00,0.00,3211
70
+ 8-Jun-17,68.00,59.00,0.00,3253
71
+ 9-Jun-17,80.10,59.00,0.00,3401
72
+ 10-Jun-17,84.00,68.00,0.00,3066
73
+ 11-Jun-17,90.00,73.00,0.00,2465
74
+ 12-Jun-17,91.90,77.00,0.00,2854
75
+ 13-Jun-17,93.90,78.10,0.01,2882
76
+ 14-Jun-17,84.00,69.10,0.29,2596
77
+ 15-Jun-17,75.00,66.00,0.00,3510
78
+ 16-Jun-17,68.00,66.00,0.00,2054
79
+ 17-Jun-17,73.00,66.90,1.39,1399
80
+ 18-Jun-17,84.00,72.00,0.01,2199
81
+ 19-Jun-17,87.10,70.00,1.35,1648
82
+ 20-Jun-17,82.00,72.00,0.03,3407
83
+ 21-Jun-17,82.00,72.00,0.00,3304
84
+ 22-Jun-17,82.00,70.00,0.00,3368
85
+ 23-Jun-17,82.90,75.90,0.04,2283
86
+ 24-Jun-17,82.90,71.10,1.29,2307
87
+ 25-Jun-17,82.00,69.10,0.00,2625
88
+ 26-Jun-17,78.10,66.00,0.00,3386
89
+ 27-Jun-17,75.90,61.00,0.18,3182
90
+ 28-Jun-17,78.10,62.10,0.00,3766
91
+ 29-Jun-17,81.00,68.00,0.00,3356
92
+ 30-Jun-17,88.00,73.90,0.01,2687
93
+ 1-Jul-17,84.90,72.00,0.23,1848
94
+ 2-Jul-17,87.10,73.00,0.00,2467
95
+ 3-Jul-17,87.10,71.10,0.45,2714
96
+ 4-Jul-17,82.90,70.00,0.00,2296
97
+ 5-Jul-17,84.90,71.10,0.00,3170
98
+ 6-Jul-17,75.00,71.10,0.01,3065
99
+ 7-Jul-17,79.00,68.00,1.78,1513
100
+ 8-Jul-17,82.90,70.00,0.00,2718
101
+ 9-Jul-17,81.00,69.10,0.00,3048
102
+ 10-Jul-17,82.90,71.10,0.00,3506
103
+ 11-Jul-17,84.00,75.00,0.00,2929
104
+ 12-Jul-17,87.10,77.00,0.00,2860
105
+ 13-Jul-17,89.10,77.00,0.00,2563
106
+ 14-Jul-17,69.10,64.90,0.35,907
107
+ 15-Jul-17,82.90,68.00,0.00,2853
108
+ 16-Jul-17,84.90,70.00,0.00,2917
109
+ 17-Jul-17,84.90,73.90,0.00,3264
110
+ 18-Jul-17,87.10,75.90,0.00,3507
111
+ 19-Jul-17,91.00,77.00,0.00,3114
112
+ 20-Jul-17,93.00,78.10,0.01,2840
113
+ 21-Jul-17,91.00,77.00,0.00,2751
114
+ 22-Jul-17,91.00,78.10,0.57,2301
115
+ 23-Jul-17,78.10,73.00,0.06,2321
116
+ 24-Jul-17,69.10,63.00,0.74,1576
117
+ 25-Jul-17,71.10,64.00,0.00,3191
118
+ 26-Jul-17,75.90,66.00,0.00,3821
119
+ 27-Jul-17,77.00,66.90,0.01,3287
120
+ 28-Jul-17,84.90,73.00,0.00,3123
121
+ 29-Jul-17,75.90,68.00,0.00,2074
122
+ 30-Jul-17,81.00,64.90,0.00,3331
123
+ 31-Jul-17,88.00,66.90,0.00,3560
124
+ 1-Aug-17,91.00,72.00,0.00,3492
125
+ 2-Aug-17,86.00,69.10,0.09,2637
126
+ 3-Aug-17,86.00,70.00,0.00,3346
127
+ 4-Aug-17,82.90,70.00,0.15,2400
128
+ 5-Aug-17,77.00,70.00,0.30,3409
129
+ 6-Aug-17,75.90,64.00,0.00,3130
130
+ 7-Aug-17,71.10,64.90,0.76,804
131
+ 8-Aug-17,77.00,66.00,0.00,3598
132
+ 9-Aug-17,82.90,66.00,0.00,3893
133
+ 10-Aug-17,82.90,69.10,0.00,3423
134
+ 11-Aug-17,81.00,70.00,0.01,3148
135
+ 12-Aug-17,75.90,64.90,0.11,4146
136
+ 13-Aug-17,82.00,71.10,0.00,3274
137
+ 14-Aug-17,80.10,70.00,0.00,3291
138
+ 15-Aug-17,73.00,69.10,0.45,2149
139
+ 16-Aug-17,84.90,70.00,0.00,3685
140
+ 17-Aug-17,82.00,71.10,0.00,3637
141
+ 18-Aug-17,81.00,73.00,0.88,1064
142
+ 19-Aug-17,84.90,73.00,0.00,4693
143
+ 20-Aug-17,81.00,70.00,0.00,2822
144
+ 21-Aug-17,84.90,73.00,0.00,3088
145
+ 22-Aug-17,88.00,75.00,0.30,2983
146
+ 23-Aug-17,80.10,71.10,0.01,2994
147
+ 24-Aug-17,79.00,66.00,0.00,3688
148
+ 25-Aug-17,78.10,64.00,0.00,3144
149
+ 26-Aug-17,77.00,62.10,0.00,2710
150
+ 27-Aug-17,77.00,63.00,0.00,2676
151
+ 28-Aug-17,75.00,63.00,0.00,3332
152
+ 29-Aug-17,68.00,62.10,0.10,1472
153
+ 30-Aug-17,75.90,61.00,0.01,3468
154
+ 31-Aug-17,81.00,64.00,0.00,3279
155
+ 1-Sep-17,70.00,55.00,0.00,2945
156
+ 2-Sep-17,66.90,54.00,0.53,1876
157
+ 3-Sep-17,69.10,60.10,0.74,1004
158
+ 4-Sep-17,79.00,62.10,0.00,2866
159
+ 5-Sep-17,84.00,70.00,0.01,3244
160
+ 6-Sep-17,70.00,62.10,0.42,1232
161
+ 7-Sep-17,71.10,59.00,0.01,3249
162
+ 8-Sep-17,70.00,59.00,0.00,3234
163
+ 9-Sep-17,69.10,55.00,0.00,2609
164
+ 10-Sep-17,72.00,57.00,0.00,4960
165
+ 11-Sep-17,75.90,55.00,0.00,3657
166
+ 12-Sep-17,78.10,61.00,0.00,3497
167
+ 13-Sep-17,82.00,64.90,0.06,2994
168
+ 14-Sep-17,81.00,70.00,0.02,3013
169
+ 15-Sep-17,81.00,66.90,0.00,3344
170
+ 16-Sep-17,82.00,70.00,0.00,2560
171
+ 17-Sep-17,80.10,70.00,0.00,2676
172
+ 18-Sep-17,73.00,69.10,0.00,2673
173
+ 19-Sep-17,78.10,69.10,0.22,2012
174
+ 20-Sep-17,78.10,71.10,0.00,3296
175
+ 21-Sep-17,80.10,71.10,0.00,3317
176
+ 22-Sep-17,82.00,66.00,0.00,3297
177
+ 23-Sep-17,86.00,68.00,0.00,2810
178
+ 24-Sep-17,90.00,69.10,0.00,2543
179
+ 25-Sep-17,87.10,72.00,0.00,3276
180
+ 26-Sep-17,82.00,69.10,0.00,3157
181
+ 27-Sep-17,84.90,71.10,0.00,3216
182
+ 28-Sep-17,78.10,66.00,0.00,3421
183
+ 29-Sep-17,66.90,55.00,0.00,2988
184
+ 30-Sep-17,64.00,55.90,0.00,1903
185
+ 1-Oct-17,66.90,50.00,0.00,2297
186
+ 2-Oct-17,72.00,52.00,0.00,3387
187
+ 3-Oct-17,70.00,57.00,0.00,3386
188
+ 4-Oct-17,75.00,55.90,0.00,3412
189
+ 5-Oct-17,82.00,64.90,0.00,3312
190
+ 6-Oct-17,81.00,69.10,0.00,2982
191
+ 7-Oct-17,80.10,66.00,0.00,2750
192
+ 8-Oct-17,77.00,72.00,0.22,1235
193
+ 9-Oct-17,75.90,72.00,0.26,898
194
+ 10-Oct-17,80.10,66.00,0.00,3922
195
+ 11-Oct-17,75.00,64.90,0.06,2721
196
+ 12-Oct-17,63.00,55.90,0.07,2411
197
+ 13-Oct-17,64.90,52.00,0.00,2839
198
+ 14-Oct-17,71.10,62.10,0.08,2021
199
+ 15-Oct-17,72.00,66.00,0.01,2169
200
+ 16-Oct-17,60.10,52.00,0.01,2751
201
+ 17-Oct-17,57.90,43.00,0.00,2869
202
+ 18-Oct-17,71.10,50.00,0.00,3264
203
+ 19-Oct-17,70.00,55.90,0.00,3265
204
+ 20-Oct-17,73.00,57.90,0.00,3169
205
+ 21-Oct-17,78.10,57.00,0.00,2538
206
+ 22-Oct-17,75.90,57.00,0.00,2744
207
+ 23-Oct-17,73.90,64.00,0.00,3189
208
+ 24-Oct-17,73.00,66.90,0.20,954
209
+ 25-Oct-17,64.90,57.90,0.00,3367
210
+ 26-Oct-17,57.00,53.10,0.00,2565
211
+ 27-Oct-17,62.10,48.00,0.00,3150
212
+ 28-Oct-17,68.00,55.90,0.00,2245
213
+ 29-Oct-17,64.90,61.00,3.03,183
214
+ 30-Oct-17,55.00,46.00,0.25,1428
215
+ 31-Oct-17,54.00,44.00,0.00,2727