{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 02. Lab excercies, Supervised learning, KNN\n",
    "\n",
    "----\n",
    "\n",
    "1. Implement naive K nearest neighbour regression as a function only using python and numpy. The signature of the functions should be:\n",
    "\n",
    "\n",
    "    ```\n",
    "    def knn_regression(x2pred, x_train, y_train, k=10):\n",
    "        \"\"\"Return prediction with knn regression.\"\"\"\n",
    "        .\n",
    "        .\n",
    "        .\n",
    "        return y_pred\n",
    "    ```\n",
    "    \n",
    "    \n",
    "2. Apply the KNN regressor on photometric redshift estimation using the provided photoz_mini.csv file. Use a 80-20% train test split. Calculate the mean absolute error of predictions, and plot the true and the predicted values on a scatterplot.\n",
    "\n",
    "3. Apply the KNN regressor on photometric redshift estimation using the provided photoz_mini.csv file. Use 5 fold cross validation. Estimate the mean and satndard deviation of the MAE of the predictions.\n",
    "\n",
    "4. Repeat 3 with the KNN regression class from sklearn. Compare the predictions and the runtime.\n",
    "\n",
    "5. Implement weighted KNN regression and apply it on the same data. Use 5 fold cross validation. Estimate the mean and satndard deviation of the MAE of the predictions. Plot the true and the predicted values from one fold on a scatterplot.\n",
    "\n",
    "---\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Populating the interactive namespace from numpy and matplotlib\n"
     ]
    }
   ],
   "source": [
    "%pylab inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "from sklearn.model_selection import train_test_split"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1, Write KNN regressor"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def knn_regression(x2pred, x_train, y_train, k=10):\n",
    "    \"\"\"Return prediction with knn regression.\"\"\"\n",
    "    dist = [((x2pred-xi)**2).sum() for xi in x_train]\n",
    "    knn = np.argsort(dist)[:k]\n",
    "    return y_train[knn].mean()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2, Apply it Photoz data with 80%-20% split"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = pd.read_csv('../data/photoz_mini.csv')  # load train data\n",
    "x = df[['u','g','r','i','z']].values  # format x as scipy expects it\n",
    "y = df['redshift'].values  # format y as scipy expects it"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "x_train, x_test, y_train, y_test  = train_test_split(x, y, test_size= 0.2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 652 ms, sys: 6.47 ms, total: 658 ms\n",
      "Wall time: 657 ms\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "yp = [knn_regression(xi, x_train, y_train) for xi in x_test]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0, 1)"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAGDCAYAAAA1cVfYAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3X+0nFV97/HPJ4cDnAByoKCrHn7aIpSgEE2VLvqLagG9NWSJCgirtfXKrS1tQZvbePVKENtaU67VXiqmrZerRQkITWOBxtsal9YaL0mTgEFyV4qIOfgjKkEkBzg5fO8fM5NM5syz53lm5pmZM+f9Wusszsw888ye54T9ffbe3723I0IAAGRZ0O8CAAAGG4ECAJBEoAAAJBEoAABJBAoAQBKBAgCQRKAAACQRKIAesH2t7e/YfsL2x20fljj2VbYfsr3X9gbbJ9e9dovtZ23/uO5npPraubb/j+0f2t5t+w7bP9mL74fhRqAASmb7QkkrJL1K0imSXiTp+oxjj5N0l6T/LulYSZskrWk47IMRcWTdz0z1+WMkra5+xsmSnpT0v7r6ZTAvESgw59h+xPZy2/fbfsr239p+ge17bT9p+59tH1N3/B11d/NftL2o7rXDbP+57Udtf9f2zbbHulzk35D0txGxPSIel3SDpLdkHPt6Sdsj4o6IeFrSSkln2z6j1YdExL3V9/0oIvZK+p+SzuvKN8C8RqDAXHWJpF+V9GJJr5N0r6T/Juk4Vf5d/37dsfdKOk3S8yX9u6Rb6177s+o5zpH005ImJL232Qfa/nnbexI/P59R1kWSttU93ibpBbZ/otWxEfGUpP+oPl/zO9Xupc22L8n4TEn6RUnbE68DuRzS7wIAbfrLiPiuJNn+kqTvRcSW6uO/V6WbR5IUER+v/W57paTHbR8t6UeS3ibppRHxw+rrfyLpU5Le1fiBEfGvksbbKOuRkp6oe1z7/ShJP2hy7O6G556oHitJH5H0zupzF0haY/s7EfHl+jfYfqkqAe/iNsoLHIRAgbnqu3W/TzV5fKQkVQd6/1jSGyUdL+m56jHHSTpM0kJJm23X3mtJI10u648lPa/uce33J3McWzv+SUmKiH+ve/4e27eq0l21P1DY/mlVWlF/EBFf6qzoAF1PGH5vVuWu+tWSjlZloFeqBITvqxJUFkXEePXn6Ig4stmJbP9CQ7ZR488vZJRhu6Sz6x6fLem7EdHYmph1rO0jJP2UsruQovpdasefLOmfJd0QEZ/MeA9QCIECw+4oSc+o0sWzUNKf1F6IiOck/bWkD9l+viTZnqhmKc0SEV9qyDZq/Mm6e/+EpLfaPrM6yP4eSbdkHPv3ks6yfYntw1XpPro/Ih6qlu8Nto+0vcD2BZKulLSuVnZJn5d0U0TcnPsKAS0QKDDsPiHpm5ImJT0oaWPD638kaaekjbZ/pMrd+OndLEBE/JOkD0raUC3LNyVdV3vd9nbbV1SP3a3KQP0fS3pc0islXVZ3uj+ofpc9klZJeltEfKH62n9WJfX2uvqWTje/C+Yns3ERACCFFgUAIKm0QFFdpuB7tr+W8bptf8T2zurEqZeVVRYAQPvKbFHcIumixOuvUWUS1GmSrpL00RLLAgBoU2mBIiK+KOmHiUMulvSJqNgoaZwFzABg8PRzjGJC0rfqHu+qPgcAGCD9nJntJs81TcGyfZUq3VM64ogjXn7GGS3XRwMA1Nm8efP3I+L4dt7bz0CxS9KJdY9PkPRYswMjYrUqyydryZIlsWnTpvJLBwBDxPY3231vP7ue1kn69Wr207mSnoiIb/exPACAJkprUdj+tKRflnSc7V2qzEQdlaTq8gL3SHqtKrNi90r6zbLKAgBoX2mBIiIub/F6SPrdsj4fANAdzMwGACQRKAAASQQKAEASgQIAkESgAAAkESgAAEkECgBAEoECAJBEoAAAJBEoAABJBAoAQBKBAgCQRKAAACQRKAAASQQKAEASgQIAkESgAAAkESgAAEkECgBAEoECAJBEoAAAJBEoAABJBAoAQBKBAgCQRKAAACQRKAAASQQKAEASgQIAkESgAAAkESgAAEkECgBAEoECAJBEoAAAJBEoAABJBAoAQBKBAgCQRKAAACQRKAAASQQKAEASgQIAkESgAAAkESgAAEkECgBAEoECAJBEoAAAJBEoAABJBAoAQBKBAgCQRKAAACQRKAAASQQKAEASgQIAkESgAAAkESgAAEkECgBAEoECAJBUaqCwfZHtHbZ32l7R5PWTbG+wvcX2/bZfW2Z5AADFlRYobI9IuknSaySdKely22c2HPYeSbdHxGJJl0n6q7LKAwBoT5ktildI2hkRD0fEs5Juk3RxwzEh6XnV34+W9FiJ5QEAtOGQEs89IelbdY93SXplwzErJX3O9u9JOkLSq0ssDwCgDWW2KNzkuWh4fLmkWyLiBEmvlfRJ27PKZPsq25tsb9q9e3cJRQUAZCkzUOySdGLd4xM0u2vprZJul6SI+IqkwyUd13iiiFgdEUsiYsnxxx9fUnEBAM2UGSjuk3Sa7VNtH6rKYPW6hmMelfQqSbL9M6oECpoMADBASgsUEbFP0tWS1kv6uirZTdttv8/20uph75T0NtvbJH1a0lsiorF7CgDQR2UOZisi7pF0T8Nz7637/UFJ55VZBgBAZ5iZDQBIIlAAAJIIFACAJAIFACCJQAEASCJQAACSCBQAgCQCBQAgiUABAEgiUAAAkggUAIAkAgUAIIlAAQBIIlAAAJIIFACAJAIFACCJQAEASCJQAACSCBQAgCQCBQAgiUABAEgiUAAAkggUAIAkAgUAIIlAAQBIIlAAAJIIFACAJAIFACCJQAEASCJQAACSCBQAgCQCBQAgiUABAEgiUAAAkggUAIAkAgUAIIlAAQBIIlAAAJIIFACAJAIFACCJQAEASCJQAACSCBQAgCQCBQAgiUABAEgiUAAAkggUAIAkAgUAIOmQfhcA3bF2y6RWrd+hx/ZM6YXjY1p+4elatnii38WaN7j+GGYEiiGwdsuk3nXXA5qanpEkTe6Z0rvuekCSqKx6gOuPYUfX0xBYtX7H/kqqZmp6RqvW7+hTieYXrj+GHYFiCDy2Z6rQ8+gurj+GHV1PQ+CF42OabFIpvXB8rA+lmX/6ef0ZG0Ev0KIYAssvPF1joyMHPTc2OqLlF57epxLNL926/mu3TOq8D3xep664W+d94PNau2Wy5fHvuusBTe6ZUujA2Eir9wFFESiGwLLFE/rT179EE+NjsqSJ8TH96etfwp1lj3Tj+rdT6TM2gl6h62lILFs8QWDoo06vf6rSzzovYyPoFQIFMADaqfQ7GRthbANFlNr1ZPsi2zts77S9IuOYN9l+0PZ2258qszzAoMqq3FOVfrtjI4xtoKjSAoXtEUk3SXqNpDMlXW77zIZjTpP0LknnRcQiSdeUVR5gkLVT6bc7NsLYBooqs+vpFZJ2RsTDkmT7NkkXS3qw7pi3SbopIh6XpIj4XonlAQZWrXIv2h3UzthIP8Y26Oqa28oMFBOSvlX3eJekVzYc82JJsv1lSSOSVkbEPzWeyPZVkq6SpJNOOqmUwgL91quEhF7P+2CJk7mvzDEKN3kuGh4fIuk0Sb8s6XJJf2N7fNabIlZHxJKIWHL88cd3vaDAfNLreTftdHUVnVOCcpXZotgl6cS6xydIeqzJMRsjYlrSN2zvUCVw3FdiuYB5rd1urnYV7eqiBTJ4ygwU90k6zfapkiYlXSbpzQ3HrFWlJXGL7eNU6Yp6uMQyAVBv590U7epqZ04JylVa11NE7JN0taT1kr4u6faI2G77fbaXVg9bL+kHth+UtEHS8oj4QVllAtB7Rbu6mEg4eEqdcBcR90i6p+G599b9HpLeUf0BCiObZnDV/22OHhvV4aMLtGfvdMu/E4tcDh5mZmPOoi97cDX+bfZMTWtsdEQfuvScln+b5ReeftB7JRa57LfMrifbT9r+UdZPLws5X5DpUQwTxwZXJ38bFrkcPJktiog4SpJsv0/SdyR9UpWU1yskHdWT0s0j3B0XR1/24Or0b8Mil4Mlz2D2hRHxVxHxZET8KCI+KumSsgs233B3XFw76yMNo0FsifK3GS55AsWM7Stsj9heYPsKSTMt34VCuDsublg3bCpS8Xdzgb9uBpxmfxtLOv8MJszORXkGs98s6cPVn5D0Zc2eD4EOkelRXK8njrWrSGZWni7I+vMtsDUTBy940M6cg/esfUC3bnx0/9IJnXZ9Lls8oU3f/OFB5wxJd26e1JKTjx24vxHSWgaKiHhElcX8UCIyPdoz6H3ZRceesrog33n7tv2P68/XGCRqirRE126ZPKhCr//cTia5bXhod9fPif5oGShsv1jSRyW9ICLOsv1SSUsj4v2ll24emSt3x4NqUOdTFJ1lnFXBz0ToXXc9oMNHF8w6XzNFWqKr1u+YVaE3lqed60t36vDI0/X015KWS/qYJEXE/dUNhggUXTbod8eDapAyxhor1GbdiVJ2ZZl6z9T0TK4gUbQl2moXvbzXt/G7jy8c1eN7p5ueE3NLnsHshRHxfxue21dGYYB2DErGWLOB5WZLKEvSArvpoHG7g70jdttzDrIqbqvSJZrn+jb77j9+ep9GRw6+AnSnzk15WhTft/1Tqi4RbvsNkr5daqmAAgali6NZhRqqVLiNXTu1sYXGu/MND+1Ofsb42Kie2ffcrLGsTiakNRsfs6Qrzj1JyxZP6No1W5u+r/76Nvvu08+FxsdGdcRhhwxcl2AZBrX7sxvyBIrflbRa0hm2JyV9Q5VJd8BAGJSMsazAFKrc6efJUkoFt7HREa1cukhS8bGsVCXWanwsz/XNKvcTU9Paet0FybINg0Hq/ixDnkAREfFq20dIWhART1aXDgcGwqBkjGVVqBPjY/ryil+RJJ264u6m761VtFnnGLEPajUUqXzyVGKp8bE813dQgnW/DPvS6HnGKO6UpIh4KiKerD73mfKKBBQzKGsD5ZkA2GrGctY5bnzT2W1/n07HcPJc32Gd/JjXoHR/liWzRWH7DEmLJB1t+/V1Lz1P0uFlFwwoIuuOuJf9xnlSnFvdnZeRJt2NSqxVRt58T+8e9hZVquvpdEm/Jmlc0uvqnn9S0tvKLBTQDf3oN85boa5ct117piqpo4ePLph1TDfL16tKbD6ndw9K92dZUqvH/oOkf7D9cxHxlR6WCeiKfvQb523BPLPvuf2/P753Wteu2apN3/yh3r/sJV0v07BXYoNg2FtUeQazf9v21yNijyTZPkbSjRHxW+UWDehM1sS1rOc7lbcFk5VG+3cbH9Xd93871y5wRQx7JTYohrlFlSdQvLQWJCQpIh63vbjEMgFdMdIkFbX2fBnytmBSYwO1mcytFgNsVtm3SoEd1koM5cuT9bSg2oqQJNk+Vmyhijkga8G8rOc7lXfQOO/YwNT0jK7/7HZJrZcT7+Zy40CjPIHiRkn/ZvsG2zdI+jdJHyy3WEDnJjIq5Kzn21XbxyEr/DQGhuUXnp65tEejx/dO728ppFJcO0mBzbMPxSBujoTeaRkoIuITquxo911J35P0+oj4ZNkFAzrVi9z++jv5Zpp93rLFE7ri3JNyB4tad1IztefbTYHN0xKhtYLMQGH7edX/HqvKntmfknSrpO9UnwMGWrsT8YrcPTe7k6+3ICMavH/ZS/ShS8/ZX7bxsdHMc9TGHJqpPd/u1qN5WiKDsugi+ic11vApVeZRbNbBa5rV1jh7UYnlArqi6CBu0bkXrTKonnp2Rss/s00r123XE1MHZzM1lu2c6z+3f25Fvdp7Uimu7abAZrU46r/XsM86RmupeRS/Vv0v6zph3igy92LtlsmmK8M2mp6J/QEglc20Z2p61vlqlX2rFNd2U2CzJuO5Wq5liyc6mrA3zCuqziepJTxelnpjRPx794sD9FeRu+fUznAp9YGnsQVTvyz5RMEU13ZSYJdfeLquXbN11vcIaX8Z222tDPuKqvNJquvpxup/D5e0RNI2Vf4Nv1TSVyX9fLlFA/Jp56416z1F7p476XqpvTdr8l39irNlWrZ4Qte02G+i3dbKsK+oOp+kup7OlyTbt0m6KiIeqD4+S9If9qZ4QFo7d62p9xS5e05tW9pKLfCU1f9fJHhO5AiO7Yz1FN0GFoMrzzyKM2pBQpIi4muSzimvSEB+7WTktLrTzZsplZV+e+W5J+1//zELRzW6IHs70HazlVKKprN2O4249vlZhmVF1fkkzwzrr9v+G0l/p0qr+EpJXy+1VEBC/d1y1hhB6q611V18491zLV228e48b5dM1t392i2T2vvs7O3nO53rUbTLp9trQaVShlmMcG7KEyh+U9LbJf1B9fEXJX20tBIBCY3dRllSd61FxiFadW01VrK1lkx9Jdus2ybre4yPjWrl0kU933+im2tBpT6nHxtKoXMtA0VEPG37Zkn3RAQzbNBXrSa4Sa3vWouMQ7S6Oy8yRlLfsmi2d7YkHXHYIR1XpP3eRCe1JSxBYm5qOUZhe6mkrZL+qfr4HNvryi4YWF+nmdTdat7Z10XGIVrdnWcFktpifjWN4wZZCxN2Y6C339uS9vvz0X15up6uk/QKSV+QpIjYavuU8ooEiRz0LKm71SLppHm7WlrdnWdV7LXF/Oq7plq1hOrP24l+7z/R789H9+UJFPsi4gmXtIY/miMHvblu7taWJ4W01eelUmTr/1Z5WgqWdP4Zxxf+Hs30e/+Jfn8+uitPeuzXbL9Z0ojt02z/pSpLjaNErK/TXLsL/TXKm0La6vNSAar+b5WnpRCS7tw8SRcjBo6jxSYuthdKerekC6pPrZf0/oh4uuSyNbVkyZLYtGlTPz66p877wOe70sWC5rp5fbMW86s/V7Msp6x1ovgbowy2N0fEknbem2xR2B6RdH1EvDsifrb6855+BYn5hAHBcnWzxbZy6aKWf6tmLZN25oAA/ZAco4iIGdsv71VhcAADguXqZgpp3r9VY799VquGmcsYNHkGs7dU02HvkPRU7cmIuKu0UkESA4Jl6uaguNT+yq3dLANQljyB4lhJP5BU32kakggUmLMGocU2CGUA8mg5mD1o5stgditsCAOgiE4Gs1u2KGy/SNKHJZ2rSkviK5KuiYhvtPOB6ByT8bIRQOcn/u7lyjOP4lOSbpf0k5JeqMpYxW1lFgppbHbfXNHltecalnRpbtj/7oMgT6BwRHwyIvZVf2rLjaNPmIzXXN4AOhcrXCrDbNw4lS9PoNhge4XtU2yfbPu/Srrb9rG2jy27gJitjM1uhkGeADpXK1wqw2zcOJUvT6C4VNJ/kbRBlYUB3y7ptyRtlsSoch8wGa+5PAF0rla4VIbZuHEqX8tAERGnJn5e1ItC4mDdWu9oLsnTXZQngM7VCpfKMBs3TuXLM48CA2g+TcbLm+WVZ15Cvzf1aReT87IxH6V8zKPAwOvmAn7NFucbGx2ZEy0yUkDRiVLnUWCwzYfKo1vdRbVrNTU9o5HqVqQTTa7ZoF7T+dSKxGDJM+HuXyTdGBH31D23OiKuKrVkaGmYJ97l2V+6SHdR47WaidjfddMYJIb1mgLtypP1dKqkP7J9Xd1zbTVf0F1zNYOnlTz7Sxftn897rYb1mgKdyBMo9kh6laQX2P6s7aNLLhNymqsZPK1k7S89YsuSjlk4qsMOWaBr12zNPWEu77Ua1msKdCLvzOx9EfE7ku6U9K+Snl9usZDHsKZMZlXKz0XoQ5eeo6enn9OeqelCE+byXqthvaZAJ/IEiptrv0TELZLeIulzJZUHBQxr/niqsm63ayjvtRrWawp0Is+Eu481PN4cEb+V5+S2L7K9w/ZO2ysSx73Bdthm7KOAYZ14l6qs2+0aynuthvWaAp0obR5Fdb/t/yfpVyXtknSfpMsj4sGG446SdLekQyVdHRHJSRLMo5gfslJUuzmnot0yAHPRoM6jeIWknRHxsCTZvk3SxZIebDjuBkkflPSHJZYFc0zWnIFezVAmTRY4IM8YRbsmJH2r7vGu6nP72V4s6cSI+MfUiWxfZXuT7U27d+/ufkkxZ/Sqa4g0WeCAMlsUbvLc/n4u2wskfUiVwfGkiFgtabVU6XrqUvkwR/VihjJpssABZbYodkk6se7xCZIeq3t8lKSzJH3B9iOqbLW6jgFtDALSZIEDygwU90k6zfaptg+VdJmkdbUXI+KJiDguIk6JiFMkbZS0tNVgNobfIOxAR5oscEBpXU8Rsc/21ZLWSxqR9PGI2G77fZI2RcS69BkwH5U9iJw3k4mlq4EDWGYcA6XM9Ne5vMQ40KlO0mPL7HoCCitzEJlMJqA9BAoMlDIHkclkAtpDoMBAKXMQmUwmoD0ECgyUMifUkckEtIetUDFwyppQRyYT0B4CBeYV9p0GiiNQYCix8ivQPQQK9FUZFTorvwLdxWA2+qZWoU/umSq0rWkrzJcAuosWBQrrVisgVaF3cufPfAmgu2hRoJButgLKqtCZLwF0F4EChXSrW2ftlkktcLMtSzqv0JkvAXQXXU8opButgFqrZKbJgpTdqNCbzZc4/4zjtWr9Dl27ZitZUEBBBAoU8sLxsaaruxZpBTRrlUjSiN21Wdj18yXIggI6Q9cTCulGt05W62MmIllxt7uhEVlQQGdoUaCQdpbBaMySGl84qsf3Ts86ztVjm51r7ZZJLb9jm6afq3RXTe6Z0vI7th1UpixkQQGdIVCgsLzLYKzdMqnrP7v9oKAwuWdKowuaD2KHlJkau3Ld9v1Bomb6udDKddtblqUb3WXAfEbXE0pRGxdo1nJorPDrNd7l17qb9kzNPo8k7ZmabtkFRRYU0BlaFPNcWWsiZQ1Yt1J/l99s69JmWg1ML1s8oU3f/KE+/dVvaSZCI7YueTmLAwJ50aKYx8paQkNq3f8/Pjba8i4/b7BpNTC9dsuk7tw8uT8ddyZCd26e7Mr3BOYDAsU8VmY2UKr/f2x0RCuXLmq5QVGRwebUsWQ9AZ2h66kEc2WJ6yLZQEW+09otk9r77L6mr42PjWrl0kX735u6LlmD0FnHZiHrCegMLYouK7M7p9vyronU7Dtdu2ar3rP2gVnvzRrEHh8b1V9ceo62XndB7qDZbBB6dIE1OnJw1lSrgWnWfgI6Q6DosrnUzZE3G6jZdwpJt258dFYAzBpXOOKwQwq3qprtn73qjWdr1RvOLrSnNllPQGfoeuqyudTNkXfyXFbZm8176Pb3z5qzUSTosFc20BkCRZcN0uSuPOMKeSbPpcYKGgPAIH3/euyVDbSPrqcuG5Rujm6MldQmu6UGlBsDwKB8fwDdQ4uiywalm6PT3ePyTHazNCsAtPv950qmGDAfEShKMAjdHJ2OFeSZ7BZqPlZQ9PuzDDgw2Oh6GlJFUl+bLd2dJ6DUVnvt1FzKFAPmIwLFkMozVpAax8gz+FzLeurUXMoUA+YjAsWQajYHoXG+QepOvlmgaaYblTkT4oDBRqCYx1J38o2BZsTN95DoRmVOphQw2BjMblO7WTq9yu7JM0Dcas5Dat9pqXuV+aBkigFojkDRhnazdHqZ3ZMnPfb8M47XrRsfVf02QlmVf7cq86xAOQiZYgCaI1C0od05Cp3ObSii1QBxbY+G+iBhKbmhT6eVOWmwwNzEGEUb2s3S6WV2T6sB4qyF/jY8tLvrZakhDRaYmwgUbWg3S6eX2T2tBog7CVpZcy9aIQ0WmJsIFG1oN0unl9k9rdJjs4LTAjtZ8XeyhhRpsMDcxBhFG9od2O11dk9qTGH5hac3XctpJiI5btDJOEuzzyQNFhh8BIo2tTuwOyjZPbUyvPP2bZqJOOi1VMXfSfcRabDA3ESgmMeWLZ7QtWu2Nn0tq+LvdL+JQQmUAPJjjGLItRp4LjpuwCxqYP6hRTGEapPaJvdMydL+uRLN5i0UHTeg+wiYfwgUQ6ZxUls0vN44/tBOxU/3ETC/ECiGTJ4NhxrHH6j4AaQwRjFk8mQfMW8BQBG0KAZUu6vMZmUl1TDwDKAoAsUA6mR12r3P7pv1fG1Ae4KBZwBtIFAMoHZmPzfbL0KSxsdGtXLpIoIDgLYRKArqxcZD7cx+zhrE3jM1rXfevk3XrNlKiwJAWwgUBfRqP4XU7OesQJUKIrUlOtj/AUA7yHoqoFf7KWTNfj7/jOMzV27Nm8nE/g8AiiJQFNCr/RSylgjf8NDuzEB1/hnHyznPz/4PAIqg66mAThfEK6LZJLisBfwm90zN2tY0pVbeXoy3AJj7aFEU0KxLyJLOP+P4nnx+VkAasVvOxq6pzaNotgHRNWu26pzrP5d7xzoA80OpgcL2RbZ32N5pe0WT199h+0Hb99v+F9snl1meTi1bPKFLXj5xUBdPSLpz82RPKtessYvG/SQajbhS4vpd7lJZUgQMAPVK63qyPSLpJkm/KmmXpPtsr4uIB+sO2yJpSUTstf12SR+UdGlZZeqGDQ/tbrnQXhlq3URT0zMasTUTsT/dtbZSbKPxsVFtve6CpudrNU6xZ2qaDCkAksptUbxC0s6IeDginpV0m6SL6w+IiA0Rsbf6cKOkE0osT1f0akC7Xn03kVRJd611IS1bPKHlF56u0QWzh7KfenZfZqsgz7gKGVIApHIDxYSkb9U93lV9LstbJd1bYnm6ouhGP93QKi132eIJHXn47Mbh9ExkVvTLLzxdoyOt86TIkAJQZqBoVgs17Uy3faWkJZJWZbx+le1Ntjft3r27i0Usrh87vOVpxezZO13ovcsWT+iIQ1v3PLLSLIAy02N3STqx7vEJkh5rPMj2qyW9W9IvRcQzzU4UEaslrZakJUuW5M0CLc1hhyzYf4d/zMJRXfe6ctdSypOWm+eYxnTYPVPNg0sNK80CkMptUdwn6TTbp9o+VNJlktbVH2B7saSPSVoaEd8rsSxdURsrqK9gn55+rvTPzdOKaXXM2i2TWv6ZbQelw6bUZ0gBmN9Ka1FExD7bV0taL2lE0scjYrvt90naFBHrVOlqOlLSHa6kcD4aEUvLKlOn2lnVtV6zCW6186YmveXZrrTVMdd/drumZ1o3xsZGRwgQAA7iaJGDP2iWLFkSmzZt6stnn7ri7qaDLJb0jQ/8p+QVlXd3AAAJhklEQVR7my0DPrrAknVQBd6sou7GDOpTVtyd+drE+Bizs4EhZ3tzRCxp570s4VFAJ0t4NGuNTD83O+xMTc/o2tsrS3UsWzzRdMXaa9Zs1cp127u2z8SXV/xKx+cAMLxYwqOAdjOe1m6ZbDkmUC9CWn7HtoMm2TXaMzWta9ds1Skr7tZ5H/h8y1nU42OjhZ4HgBoCRQFZq7qm7urfs/aBzMX8Uqafi/3dTVlq7ZH65cazrFy6aNakvNEF1sqliwqXDcD8whhFidZumdS1a7Zmruo6usBNu5/qLbDU4pCDpHaxY7VYYP5ijGJArVq/I7n096o3nq133XW/phIptkWChJTexa7Z0uUA0ApdTyVqtfzFqvU7cm82VARrNAHoJloUJVm7ZVILqqu8ZikywF0UazQB6BZaFCWopbS22ieiTKzRBKBbCBQlyEpp7ZXREbNGE4CuoeupS9ZumdS7//4BPfVs/wJEzRGHHsKgNYCuIVB0wXvWPqC/2/ho6Z8zkTEzvFGrVWEBoAi6njq0dsukbi05SFjSX1x6jr684lc0kWPsobZHNgB0Ay2KgtZumdTKddv337XbGbsxdVHowJyI5Reernes2arU4ub9HEQHMHxoURSwdsuklt+x7aCunV7UyccsPLAe07LFEzp6YXp9pjytDgDIixZFAavW72i55EYZnpme0Xkf+Pz+pTcez9j2VGJXOgDdR6AooMwJcil7p5/T3upnT+6ZktW8u2vEZtMhAF1H11MBgzJEHJpdlrHREd34prMJEgC6jkBRQDc7ncbHRg8ae1g4ukBjo/n/HCEVWu4cANpF11OfPDE13XT71MalwJ96Zl/TeRET42PsTAegJwgUfZK1FlPjUuDN9tpmwBpALxEo+qDIWky1oMGGQwD6hUBRwPjYaFvLY9gH5lscs3BU171uUaGKng2HAPQTgaKAlUsXafkd2wrNpRgbHWGgGcCcRtZTAcsWT2jVG88+KNsolag0PjZKkAAw59GiKCjPYLMlXXHuSXr/spf0oYQA0F0Eig4x2Axg2BEouoDBZgDDjDEKAEASgQIAkESgAAAkESgAAEkECgBAEoECAJBEoAAAJBEoAABJBAoAQBKBAgCQRKAAACQRKAAASQQKAEASgQIAkESgAAAkESgAAEkECgBAEoECAJBEoAAAJBEoAABJBAoAQBKBAgCQRKAAACQRKAAASQQKAEASgQIAkESgAAAkESgAAEkECgBAUqmBwvZFtnfY3ml7RZPXD7O9pvr6V22fUmZ5AADFlRYobI9IuknSaySdKely22c2HPZWSY9HxE9L+pCkPyurPACA9pTZoniFpJ0R8XBEPCvpNkkXNxxzsaT/Xf39M5JeZdsllgkAUFCZgWJC0rfqHu+qPtf0mIjYJ+kJST9RYpkAAAUdUuK5m7UMoo1jZPsqSVdVHz5j+2sdlm1YHCfp+/0uxIDgWhzAtTiAa3HA6e2+scxAsUvSiXWPT5D0WMYxu2wfIuloST9sPFFErJa0WpJsb4qIJaWUeI7hWhzAtTiAa3EA1+IA25vafW+ZXU/3STrN9qm2D5V0maR1Dcesk/Qb1d/fIOnzETGrRQEA6J/SWhQRsc/21ZLWSxqR9PGI2G77fZI2RcQ6SX8r6ZO2d6rSkrisrPIAANpTZteTIuIeSfc0PPfeut+flvTGgqdd3YWiDQuuxQFciwO4FgdwLQ5o+1qYnh4AQApLeAAAkgY2ULD8xwE5rsU7bD9o+37b/2L75H6UsxdaXYu6495gO2wPbcZLnmth+03VfxvbbX+q12XslRz/j5xke4PtLdX/T17bj3KWzfbHbX8vawqBKz5SvU73235ZrhNHxMD9qDL4/R+SXiTpUEnbJJ3ZcMzvSLq5+vtlktb0u9x9vBbnS1pY/f3t8/laVI87StIXJW2UtKTf5e7jv4vTJG2RdEz18fP7Xe4+XovVkt5e/f1MSY/0u9wlXYtflPQySV/LeP21ku5VZQ7buZK+mue8g9qiYPmPA1pei4jYEBF7qw83qjJnZRjl+XchSTdI+qCkp3tZuB7Lcy3eJummiHhckiLiez0uY6/kuRYh6XnV34/W7DldQyEivqgmc9HqXCzpE1GxUdK47Z9sdd5BDRQs/3FAnmtR762q3DEMo5bXwvZiSSdGxD/2smB9kOffxYslvdj2l21vtH1Rz0rXW3muxUpJV9repUom5u/1pmgDp2h9Iqnk9NgOdG35jyGQ+3vavlLSEkm/VGqJ+id5LWwvUGUV4rf0qkB9lOffxSGqdD/9siqtzC/ZPisi9pRctl7Lcy0ul3RLRNxo++dUmb91VkQ8V37xBkpb9eagtiiKLP+h1PIfQyDPtZDtV0t6t6SlEfFMj8rWa62uxVGSzpL0BduPqNIHu25IB7Tz/j/yDxExHRHfkLRDlcAxbPJci7dKul2SIuIrkg5XZR2o+SZXfdJoUAMFy38c0PJaVLtbPqZKkBjWfmipxbWIiCci4riIOCUiTlFlvGZpRLS9xs0Ay/P/yFpVEh1k+zhVuqIe7mkpeyPPtXhU0qskyfbPqBIodve0lINhnaRfr2Y/nSvpiYj4dqs3DWTXU7D8x345r8UqSUdKuqM6nv9oRCztW6FLkvNazAs5r8V6SRfYflDSjKTlEfGD/pW6HDmvxTsl/bXta1XpannLMN5Y2v60Kl2Nx1XHY66TNCpJEXGzKuMzr5W0U9JeSb+Z67xDeK0AAF00qF1PAIABQaAAACQRKAAASQQKAEASgQIAkESgANpg+xrbC/tdDqAXSI8F2lCd+b0kIr7f5LWRiJjpfamActCiABrY/m3bW6s/37C9oeH135f0Qkkbaq/Z/rHt99n+qqSfs/1IdTa0bC+x/YXq70dU9wy4r7o3QrPVb4GBQqAAGkTEzRFxjqSfVWVtnP/R8PpHVFkf5/yIOL/69BGq7AHwyoj418Tp363KcjM/q8ryGqtsH9H1LwF0EYECyPZhVSr1z+Y4dkbSnTmOu0DSCttbJX1BlTWHTmq7hEAPDORaT0C/2X6LpJMlXZ3zLU83jEvs04EbscPrTy3pkojY0XEhgR6hRQE0sP1ySX8o6crEfgVPqrKseZZHJL28+vsldc+vl/R7td0Yqyv/AgONQAHMdrWkY1UZrN5q+2+aHLNa0r2NA911rpf0YdtfUqVbquYGVVbzvN/216qPgYFGeiwAIIkWBQAgiUABAEgiUAAAkggUAIAkAgUAIIlAAQBIIlAAAJIIFACApP8PbbiQIf+HkVcAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x432 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "figsize(6,6)\n",
    "plot(y_test, yp,'o')\n",
    "xlabel('z true')\n",
    "ylabel('z predicted')\n",
    "title('mae = %.3f' % np.mean(np.abs(y_test-yp)))\n",
    "xlim(0,1)\n",
    "ylim(0,1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3, Apply it in 5 fold cross validation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.model_selection import KFold"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "maes = 0.0487 +/- 0.0031\n",
      "CPU times: user 3.04 s, sys: 13.8 ms, total: 3.06 s\n",
      "Wall time: 3.08 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "kf = KFold(n_splits=5)\n",
    "maes = []\n",
    "for train_index, test_index in kf.split(x):\n",
    "    x_train, x_test = x[train_index], x[test_index]\n",
    "    y_train, y_test = y[train_index], y[test_index]\n",
    "    yp = [knn_regression(xi, x_train, y_train) for xi in x_test]\n",
    "    maes.append(np.mean(np.abs(y_test-yp)))\n",
    "print 'maes = %.4f +/- %.4f' % (np.mean(maes), np.std(maes, ddof=1))  # note ddof!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4, Use sklearn, compare results and runtime"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.neighbors import KNeighborsRegressor"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "maes = 0.0487 +/- 0.0031\n",
      "CPU times: user 13.6 ms, sys: 1.6 ms, total: 15.2 ms\n",
      "Wall time: 13.9 ms\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "knnr = KNeighborsRegressor(n_neighbors=10)\n",
    "kf = KFold(n_splits=5)\n",
    "maes = []\n",
    "for train_index, test_index in kf.split(x):\n",
    "    x_train, x_test = x[train_index], x[test_index]\n",
    "    y_train, y_test = y[train_index], y[test_index]\n",
    "    knnr.fit(x_train, y_train)\n",
    "    yp = knnr.predict(x_test)\n",
    "    maes.append(np.mean(np.abs(y_test-yp)))\n",
    "print 'maes = %.4f +/- %.4f' % (np.mean(maes), np.std(maes, ddof=1))  # note ddof!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 5, Implement weighted"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "def knn_regression_w(x2pred, x_train, y_train, k=10):\n",
    "    \"\"\"Return predictions with weighted knn regression.\"\"\"\n",
    "    dist = np.array([((x2pred-xi)**2).sum() for xi in x_train])\n",
    "    knn = np.argsort(dist)[:k]\n",
    "    w = dist[knn]\n",
    "    return w.mean() * (y_train[knn] * 1./w).mean()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "maes = 0.1387 +/- 0.0122\n",
      "CPU times: user 3.15 s, sys: 24.9 ms, total: 3.17 s\n",
      "Wall time: 3.21 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "kf = KFold(n_splits=5)\n",
    "maes = []\n",
    "for train_index, test_index in kf.split(x):\n",
    "    x_train, x_test = x[train_index], x[test_index]\n",
    "    y_train, y_test = y[train_index], y[test_index]\n",
    "    yp = [knn_regression_w(xi, x_train, y_train) for xi in x_test]\n",
    "    maes.append(np.mean(np.abs(y_test-yp)))\n",
    "print 'maes = %.4f +/- %.4f' % (np.mean(maes), np.std(maes))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0, 1)"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAGDCAYAAAA1cVfYAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3X2UHFd55/Hvo1ELj2zjEbHJwmBjkRg7FgYLT4AcJRssCDYYZMVAjMGbQFi8eQEWh9VmWDj4hSRW8BICOd6AAhwWwovBJorAZrWbyCygRKxHyAJk0DmKDbbGgAVIBqyxPZKe/aO7pVZPVXVVdd3qqu7f5xwdT3dXV9+u8dyn7stzr7k7IiIicRYNugAiIlJtChQiIpJIgUJERBIpUIiISCIFChERSaRAISIiiRQoREQkkQKFSAnM7Goz+4GZPWRmHzGzx8Uct8TMbjGz75qZm9nzE477jpnt7Xp+tZl93cx+amb3mNlVAb6OjBgFCpHAzOwiYBp4AXAm8DTguoS3fBW4EvhBwjHrgAe7PqcB/APwQeAU4HLgr8zsWXnLLgIKFFJDrbvtdWb2DTN72Mw+bGa/aGZfNLOfmdk/mdmyjuM/23E3/2UzW9Hx2uPM7L+b2X1m9kMz+4CZjRdc5N8DPuzuu9x9P/Au4LVRB7r7Y+7+1+7+VeBw1DFmtpxmILmh66UnAI8HPu5NdwLfBs4t5mvIqFKgkLp6OfBbwNOBlwFfBP4bcCrN/6/f3HHsF4GzgCcCXwc+0fHaX7bOcT7wy8Ak8M6oDzSzXzezAwn/fj2mrCuAnR2PdwK/aGa/kOkbH/M3NL/rXOeT7v5D4FPA68xszMx+DXgqzRaKSG6LB10AkZz+plUxYmZfAR509x2tx/9As5sHAHf/SPtnM7sW2G9mpwA/Bd4APNPdf9J6/S+ATwJv6/7A1l3+RI6yngQ81PG4/fPJwI+znMjMfhtY7O7/EDN+8SngQ8D7Wo//0N3vz1ZckeMpUEhd/bDj57mIxycBmNkY8OfAK4HTgCOtY04FHgcsBbabWfu9BowVXNaf0+wSamv//LMsJzGzE4F3Ay+Jef0c4Gbgt4H/Q7MV9QUze8Ddb8taaJE2dT3JsHs1cCnwQpoDvGe2njfgRzSDygp3n2j9O8XdT4o6kZn9hpn9POHfb8SUYRfQOaD8LOCH7p6pNUGz4j8T+IqZ/QD4HPCk1vjLmcAzgN3uvtndj7j7buA24MUZP0fkOAoUMuxOBh6l2cWzFPiL9gvufgT4O+C9ZvZEADObbM1SWsDdv+LuJyX8+0pMGT4GvN7Mzm0Nsr8D+GhcgVsD7Ce0Hi4xsxOs2eT5FnA6zfGU84H/SLMldT5wP7ADOKs1RdbM7JeAl3L8+IhIZgoUMuw+BnwPmAXuBrZ1vf6nwB5gm5n9FPgn4OwiC+Du/4tml9EdrbJ8D7im/bqZ7TKz13S8ZTfNls4ksLn181Pd/ZC7/6D9D/gJcKT1+LC7/xvw+8D7aY6//F/gVuDDRX4fGT2mjYtERCSJWhQiIpIoWKBoLVPwoJl9K+Z1M7P3m9meVuLUs0OVRURE8gvZovgocHHC6y+mOYvjLOAq4G8DlkVERHIKFijc/cs0B9viXAp8rLXUwDZgwsyeFKo8IiKSzyDHKCZpTulr29t6TkREKmSQmdkW8VzkFKzWUslXAZx44okXnHPOOSHLJSIydLZv3/4jdz8tz3sHGSj20kweansK8EDUge6+AdgAMDU15TMzM+FLJyIyRMzse3nfO8iup03A77ZmPz0PeMjdvz/A8oiISIRgLQoz+xTwfODU1i5c1wANAHf/AHA7zcXN9gAHgdeFKouIiOQXLFC4+xU9Xnfgj0N9voiIFEOZ2SIikkiBQkREEilQiIhIIgUKERFJpEAhIiKJFChERCSRAoWIiCRSoBARkUQKFCIikkiBQkREEilQiIhIIgUKERFJNMj9KEQqZeOOWW7cvJsHDszx5Ilx1l10NmtXatNFEQUKEZpB4m2f+yZz84cBmD0wx9s+900ABQsZeep6EgFu3Lz7aJBom5s/zI2bdw+oRCLVoUAhAjxwYC7T8yKjRIFCBHjyxHim50VGiQKFCLDuorMZb4wd99x4Y4x1F509oBKJVIcGs0U4NmCtWU8iCylQiLSsXTmpwCASQV1PIiKSSIFCREQSKVCIiEgiBQoREUmkQCEiIokUKEREJJEChYiIJFKgEBGRRAoUIiKSSJnZMnCjumHQqH5vqR8FChmoUd0waFS/t9STup4ktY07Zlm1fgvLp29j1fotbNwx2/c5R3XDoFH93lJPalFIKqHugEd1w6C47zc75N9b6kktCkkl1B3wqG4YFPf9DAppqYkUSYFCUgl15z+qGwatu+hsLOJ5B3U/SeUoUEgqoe78166c5IbLzmNyYhwDJifGueGy84Z+QHftykk85rVh73aT+tEYhaSy7qKzjxujgOLu/Ed1w6DJifHIMYlh73aT+lGLQlIZ1Tv/kEa1203qRy0KSW1U7/xD0T7dUhcKFJJKXbKI61LONgVfqQMFCumpLlnEdSlnldQtsMpgaIxCeqpLFnFdylkV7cA6e2AO51hgVR6HdFOLYgRlvYusS/Z0XcpZFUmBVa0K6aQWxYjJcxdZl+zpupSzKpKWESlqLa+yhFiHTI5RoBgxebpn6jKNs8hyDqLiKfszkwJonbqh1IUWngLFiMnTPVOXHIqiyjmIimcQnxkVWDvVZXxHY1PhaYxixDw5Jhv4lPEGq9ZviR23qMs0ziLKOYi++0F8ZmceR9yqtXUY39HYVHhqUYyYqLvIxiLj4ccOqeneMoiKZ1CV3dqVk2ydXs1kjcd3NDYVngLFiInqnjnphMXMHz5+ibpRbroPouIZdGVXl3GoKHUue10oUIyg9l3kvesvYev0ag4cnI88blSb7oOoeAZd2dVlHCpKncteFxqjkNhxi1Ftug9iDaYqrPtUl3GoKHUuex2Ye9yq+NU0NTXlMzMzgy7GUOle+gKad7O6KxMZHma23d2n8rxXLQqpxN2siFSXAoUAarqLSDwNZouISCK1KERqQMuByyAFbVGY2cVmttvM9pjZdMTrZ5jZHWa2w8y+YWYvCVkekTrSWkYyaMEChZmNATcBLwbOBa4ws3O7DnsH8Bl3Xwm8CvgfocojUlday0gGLWSL4jnAHne/x90fAz4NXNp1jAOPb/18CvBAwPKI1JLWMpJBCxkoJoH7Ox7vbT3X6VrgSjPbC9wOvClgeURqadDLe4iEDBQW8Vx3dt8VwEfd/SnAS4CPm9mCMpnZVWY2Y2Yz+/btC1BUkeoa9PIeIiFnPe0FTu94/BQWdi29HrgYwN3/1cxOAE4FHuw8yN03ABugmZkdqsAieYWclaSESBm0kIHiTuAsM1sOzNIcrH511zH3AS8APmpmvwKcAKjJILXSvQRKe1YSEFuZZw0sSoiUQQrW9eTuh4A3ApuBb9Oc3bTLzK43szWtw94KvMHMdgKfAl7rdVt8SkZe1llJmu4qdRM04c7db6c5SN353Ds7fr4bWBWyDDLaykhUyzoraRC72Yn0Q5nZQ2zUs3nzdAnlkXWZ9iyBZdR/h1INWutpSA1798bGHbOsWr+F5dO3sWr9lsjvVVaiWtZZSWmnuw7771DqQ4FiSA1zNm/aCrSsRLWsO6ylDSzD/DuUelHX05Aa5mzetH38Ze7cl2VWUtrprsP8O5R6UaAYUsO8vWnaCnTdRWdH7txXhUS1qMDSPR5xyniDA3ML9zMfht+h1Iu6nobUMGfzpu3jz9olNEhR3WkPP3aIxqLjFzgYlt+h1ItaFENqmLN5s7QU6pKoFtWdNn/YWba0wdIli4fudyj1okAxxOpSSWY1jEEwrjtt/8F5li5ZzHsvP7/W30/qTYFCamnYgmDcmBKEy/8QSUtjFCIVEDWm1EnTYmWQ1KIQKUg/WdSd3WlxLQtNi5VBUYtCpABFZFGvXTnJ1unVTGqjIqkYBQqRAhSZRT3MU5ulntT1JFKAIrOoh3FWF2iBwzpToBApQNGZ8MM2q6uslXwlDHU9iRQgT3dRmhVwh4UWOKw3tShECpC1u2jU7rC1wGG9KVDISArRX56lu2jUdrkb5kUqR4G6nmTkVGFDoDR32MPUNaWZXPWmQCEjpwr95b1WwK1CMCtSnVbylYXU9SQjpwr95b1WwB3Grqlhm8k1StSikJGTdj+LkHrdYVchmIm0qUUhI6cqO98l3WFr8FeqRC0KGTl16C/X4K9UiVoUMpLK7C/PMxV3WJfxkHpSoBgiWkunGEVex34S6zT4K1WhQDEkRi3TN5Sir2Oo2Uu6KZAyaYxiSFQhN2AYFH0dQ8xeGrYcC6k+BYohoemUxSj6OoaYiqubAimbAsWQqEJuwDAo+jqGmL2kmwIpm8YohkRVcgPqrujrWMTspe7xiImlDfYfnF9wXHcw0ziGFEWBYkhoOmUxQlzHfmYvRQ2uNxYZjTFj/rAfPa47mGUZlFdAkV7M3XsfVSFTU1M+MzMz6GKIlGLV+i2RGdoT4w1OfNzi2Mo97n2TE+NsnV599HF3QIFm0KlaAqL0z8y2u/tUnveqRSFSYXHjDg/NzXPXNS/K/L7u54dx8UEpngazRSos7+B62vdpYFzSUKAQqYiojYryzprq9b72Z8V1PGu2nHRSoBCpgLgkOiDXAoZJCx92flYUzZaTbhrMFqmAtIPPIT+r/Xma9TScNJgtUmFppp+WOVYQd06DwoOSDAd1PYkElHZdpjIz65XFL1kpUIgElHZdpjI3KtKmSJKVup5EAkrbpVRmZr2y+CUrBQqRgLLsfV3mRkXaFEmyUNeTSEBZu3micilEBk0tChHCLYyXpZtHuxRKVSlQSK2EqNBDV9Bpu3m07pJUlbqepDZCbQFalR3jtO6SVJUChdRGqAo9riKePTBX6jiB8hukqhQopDZC3XEnVcRFtVrSUH6DVJUChdRGqDvuqAq6U7+tlrQzmdoL+U2MN44+d0KjGn+imo012qrxf6FICqHuuDtXWo2Tt9WSZ1zl0UNHjv68/+B8aS2aOKHGhqQ+FChk4LLecWddcjvNZ199810ALFvaiDwub6sl67hKVQbWO1WxTFIuTY+Vgco6NbXIjOKoz24sMhpjxvzhY8vv99NqyTquUuQ4TFFTiTUbS9SikIEa5N1q1GfPH3FOXLK4sFZL1nGVosZhiuwu0mwsUaCQgRrk3WrcZzw0N8/W6dXcu/4Stk6v7qsFk3VcpahxmCIDsGZjibqeZKCyLJpXx8/OuoRHu4IfM+Owe+4d54oMwFptVhQoJFHRS2Z0n+/Cc07j1u2zx939lnW3uu6is48bowj12WnGVbrHSw67Hy1LnutddBDUarOjLWjXk5ldbGa7zWyPmU3HHPM7Zna3me0ys0+GLI9kU/S0yKjz3bp9lpdfMFnoTKa0QsyiyqvosRp1F0mRgrUozGwMuAn4LWAvcKeZbXL3uzuOOQt4G7DK3feb2RNDlUeyK3qRurjz3fGdfX3v1Zy35ZP3TrnollbRYzXqLpIihex6eg6wx93vATCzTwOXAnd3HPMG4CZ33w/g7g8GLI9kVHTlFWrguuzluUN8XojxEnUXSVFCdj1NAvd3PN7beq7T04Gnm9lWM9tmZhdHncjMrjKzGTOb2bdvX6DiSreip0WGmmZZ9hTbEJ+nriKpspCBwiKe867Hi4GzgOcDVwAfMrOJBW9y3+DuU+4+ddpppxVe0FDqvj5O0ZVXqMqw7Cm2IT6vSuMlVVf3v6s6Ctn1tBc4vePxU4AHIo7Z5u7zwL1mtptm4LgzYLlKMQy7lRXdzx2q37ysKbbtcYnuu52iPk9dRb0Nw99VHYUMFHcCZ5nZcmAWeBXw6q5jNtJsSXzUzE6l2RV1T8AylWZYdisruvJKOl/eAeIyprl2V1Ddqt5NFGqr17INy99V3QQLFO5+yMzeCGwGxoCPuPsuM7semHH3Ta3XXmRmdwOHgXXu/uNQZSqT1sfJpp87xbwtlaTKs/u1hx89FBsk8ibFlWWY7sL1dzUYQRPu3P124Pau597Z8bMDf9L6N1QGmXFcR/3eKWZt+SRVnsCC1+IY9D21N7RhugvX39VgaK2nQDSLJZuy7xSv+/yu2MozqmKNU4cKapjuwvV3NRixLQoz+xkLZykd5e6PD1KiIaGEp2zKvFPcuGOW/QfnI1/LUnnWpYIaprtw/V0NRmygcPeTAVpjCj8APk6zpf0a4ORSSldzmsWSXlnrLgGJ+Q7tyjOqYl22tMHSJYtrV0GVeW3LoL+r8qUZo7jI3Z/b8fhvzexrwLsDlUlGUJl3ikmthnblGVWxXvOyFbWsoHQXLv1KEygOm9lrgE/T7Iq6guYMJZFClXWnGNcVMzHeOO7zq1CxFjWtVXfh0o80geLVwPta/xzYysJ8CJHaiOuKuXbNiqOPq7BY4DBNa5V6s+YM1fqYmprymZmZgZZhWJKXRlmI32FUUt54Yyz3Uhyr1m+JbPmMmXHEXf/vSSZmtt3dp/K8t2eLwsyeDvwt8Ivu/gwzeyawxt3/LM8H1l0Rd3kKNIMXoium6HyFuLGUw62bO7UwpCxp8ij+juaeEfMA7v4NmstxjKR+Vw4tejMgqY6i8xXSTF8NuUquSFuaQLHU3f9f13OHQhSmDvqtDMpeElvKU/Qy6lHJZVHqmDgn9ZImUPzIzH6JVvKdmb0C+H7QUlVYv5XBMGXJyvGKzhruXnp8zKJW7q9n4pzUS5pZT38MbADOMbNZ4F6aSXcjqd/kpWHKks1imMdlOr/bxNIGj1u8iIfm5gv5np1jKXGD5XVNnJP6SBMo3N1faGYnAovc/WetpcNHUr/JS3XIki26Uh/maZ7d323/wXnGG2O89/LzC/lu3b+Ll18wyR3f2TeUAVeqq+f0WDP7urs/u+u57e5+QdCSxajC9Nh+VfnuuugpnhA/zXNyYrzyK6/2kvW7ZfndF/m7qPL/c1KOINNjzewcYAVwipld1vHS44ET8nyYNFU5SzbEktTDPC6T5btlbVkV9bsY5hadlCNpMPts4KXABPCyjn/PBt4QvmgyCCEq9aJnAxWliL2Xs3y3rDPeivpdaKad9Cs2ULj7P7r764CXuvvrOv692d3/pcQySolCVOpV3EOgqHyWLN8ta8Vf1O9imFt0Uo4002P/wMwm2g/MbJmZfSRgmWSAQlTq3dM8JyfG+xrzKEJRd9lZvlvWir+o30VVW3RSH2lmPT3T3Q+0H7j7fjNbGbBMMkChlqSu2rhMkXfZab9b1hlvRf0u6jDTTqotTaBYZGbL3H0/gJk9IeX7pKaqVqmHMIh8ljwVfxG/C+1HIf1KU+G/B/gXM7ul9fiVwJ+HK5JIeIO6yx5UEB6F4C/h9AwU7v4xM5sBVtPcCvUyd787eMlEAtJdtkh6sQl3ZvZ4d/9pq6tpAXf/SdCSxRiGhDsRkbKF2o/ikzTzKLbTWhCw/Xmtx0/L84EiIlIvsYHC3V/a+u/IruskIiLJS3g8O+41AHf/evHFERktWoNJ6iCp6+k9rf+eAEwBO2l2Oz0T+Brw62GLJjLctAaT1EXSEh4XuvuFwPeAZ7v7VGvF2JXAnrIKKDKstAaT1EWaPIpz3P2b7Qfu/i0zOz9gmUSCqFo3j9ZgkrpIEyi+bWYfAv6e5mynK4FvBy2VSMGK7OYpKuCM6m6HUj9pFgV8HbAL+M/AW4C7W8+J1EZR3TxFrToL1VxVVyRKmszsR8zsA8Dt7q7O0yFWta6ZIpWxt0PWa6XscKmLnoHCzNYANwJLgOWt8Ynr3X1N6MJJeYZ9Bk5R3TxFjytoDSapgzRdT9cAzwEOALj7XcCZAcskAzCMM3A6d7B7+NFDNMbsuNe1t4NIOmkCxSF3fyh4SWSgyp6BU8Q2pL3O3zmWcGBuHhyWLW30tXmSxhVkFKWZ9fQtM3s1MGZmZwFvBrQV6pApcwZOGd1cUS2k+SPO0iWL2fHOF+U+b9S4woXnnMaNm3dz9c13aZxBhlKaFsWbgBXAozQXCnyI5uwnqZh+7tLLvFMuo5srZAtp7cpJtk6v5t71l7DuorO5dftsIbOgRKoqsUVhZmPAde6+Dnh7OUWSPPq9Sy9zBk4Z3VxltZCKnAUlUlWJgcLdD5vZBWUVRvIrosIqawZOGZV4WTvYKbtaRkGarqcdZrbJzP6DmV3W/he8ZJJJnSqsMrq51q6c5IbLzmNyYhyjOYj9uMWLuPrmuwodPNcsKBkFaQaznwD8mOZWqG0OfC5IiSSXIu7Sy0q4K6ubq91CCjl4Pqi9t0XKlCYzW8t11EC/FVbZCXdlJpqFHEdQdrWMgjSZ2U8D3gc8j2ZL4l+Bt7j7vYHLJhn0W2GlqUzrusRH6G45ZVfLsEvT9fRJ4Cbgt1uPXwV8GnhuqEJJPv1UWL0q0zov8aFVWkX6k2Yw29z94+5+qPWvvdy4DJFeg7J1XuJD2dQi/UkTKO4ws2kzO9PMnmpm/xW4zcyeYGZPCF1AKUevyrROs6q6dc+Ayrt8h8ioStP1dHnrv/+p6/nfp9myeFqhJZKB6DXGMcjumyLGRjq75drny7rkRl3HaET6Ze716kWamprymZmZQRdj5HSPUUCzxRH6zjzqc6GZF3HNy1bk2p0uz/cY1PcXKYqZbXf3qTzvTdP1JDKw7puosRGA/QfnWXfLzsyJc3nHWuo8RiPSrzRdTyLAYKaBJo2BzB92rvv8rkxlyjvWUucxGpF+KVBIJWzcMct1n9/F/oPzAEyMN7h2zYrYsZG29vFx5+weU8g71qIptjLKenY9mdk/m9lLup7bEK5IMmo27phl3S07j6v0D8zNs+6zO7nwnNMWzMZKe87OjYvaeR9R50szVVZTbGWUpRmjWA78qZld0/FcrgERkSg3bt7N/OGFkyrmjzh3fGcfN1x2HhbxPmi2POLOGTWm0D5f1rEWTbGVUZam6+kA8ALg/Wb2eeDKsEWSUZPUz//AgbmjlfG6z+5k/sixgNJYZFy7ZkWmc7bPl6eC11IdMqrSBApz90PAH5nZa4GvAsuClkpqLWu+QdI4xCIzlk/fxpMnxrn8Oadzx3f2pTpvnVbTFam6NIHiA+0f3P2jZvZN4I/DFUnqLM+aUOsuOpt1t+yM7H463MrzmT0wx63bZ1N399RtNV2RKus5RuHuH+x6vN3dfz/Nyc3sYjPbbWZ7zGw64bhXmJmbmcY+ai5PvsHalZPc+IpnsWzpsfGGqDGJLHkL/Y4pKG9C5Jhg02Nb+23fBPwWsBe408w2ufvdXcedDLwZ+Fqoskh5suQbdHftdGZaL5++LdP5o4RcTVdklITMzH4OsMfd73H3x2guTX5pxHHvAt4NPBKwLFKStFuDxk1fbWdaD3qL0UF/vkiVhAwUk8D9HY/3tp47ysxWAqe7+xeSTmRmV5nZjJnN7Nu3r/iSSmHS5hv06toZdN7CoD9fpEpCZmZHdTMfHa00s0XAe4HX9jqRu28ANkBzUcCCyicBpN1pr1fXzqC3GB3054tUSchAsRc4vePxU4AHOh6fDDwD+JKZAfw7YJOZrXF3LQ9bY2nGBtJMXx103sKgP1+kKkJ2Pd0JnGVmy81sCc0tVDe1X3T3h9z9VHc/093PBLYBIx8kNu6YZdX6LSyfvo1V67dkXh21LtS1I1IfwVoU7n7IzN4IbAbGgI+4+y4zux6YcfdNyWeorlCJWKM0d19dOyL1oY2LMgq5gc2q9Vsiu2MmJ8bZOr26r3OLyGjTxkUlCpmIpbn7IlJFChQZhazMNXdfRKpIgSKjkJW5BnhFpIoUKDIKWZlrzwMRqSJthZpRntk6WWZJFTV3X0tki0hRFChyyFKZD2LKa7+fWcUgU8UyiYwKBYrAkmZJharo+vnMuCAz872fpN40qGijlF8iUkUaowgsbjZU3I5uIT+z18ysjTtmeetndkYGmU9suy92pdfQtDeEyGCpRRFAZzfJIrOju7R1stZx7TviIrtW8mwD2r5rjyordKzm2BLXQgnRRaT8EpHBUouiYN37LCRVvO074l57M2SVZ2ZW1F17L90VddHfo035JSKDpUBRsCwVbruiLbprJc8026SusKj14mFhRR2qi0j5JSKDpa6ngmXpDmlXtEnjGMunb8vVhZN1ZpaxsHsJYMyMK557Ordun12wvlV3RR2qi0gLCIoMlgJFweLGB7or4s6KNu49tN4TepbPjZt3RwYJA97zO89i7cpJpp76hJ4VdZ6xkbS0N4TI4KjrqWBx3SSved4ZsV1BUe/pFnKWT9wdv3MsMK1dOcnW6dXcu/4Stk6vjqy01UUkMpzUoihYnm6StSsnmfneT/jU1+6PHfyGcLN84loCkxlbAuoiEhlOChQBZO0m2bhjllu3zyYGCQg3y2fdRWdH7rGRpyWgLiKR4aNAUQFpZkqF7MJRS0BEkihQVEBSl5JBKRV3d7Boj4coWIiIAkUJemUrx40RLFvaYOmSxaVU3FpPSUTiaNZTYGmylaNmCzXGjJ8/cqiwLOeNO2ZZtX4Ly6dvY9X6LQvOo/WURCSOAkVgaSrgqEzqE5csZv6IJ74PegeA9jG9gpXWUxKROOp6CixtBdw9W2j59G0935e2uyjNsuMhk+VEpN7Uoggs74J2ad6XtrsoTbBSspyIxFGgCCxvBZzmfWlbK2mCjvbrFpE46noKLG+OQpr3pe0uSptQp2Q5EYmiQFGCvBVwr/dlCQAA127axYG5eQBOaKgxKSLpKFAUrJ0zMXtgjrHW7naTgRLmsrZWHj105OjP+w/OK09CRFJRoChQ9yyk9tpNIZPX0rZW0sx8EhGJokBRoKQ1mwZdKYfKkwixR7aIVIsCRYF6Vbr9VMr9Vsgh8iS07IfIaNCIZoF6VboTSxu5zpsms7qXC885bcHe1/3mSWjZD5HRoEBRoF471f38kUO51mrqt0Ju73fRuSCIAS+/oHnX32sJkDha9kNkNChQFKgzaS3K/BHPdbfdb4UcFWgc+MLO7/fVUsmbdS4i9aJAUbD23tLd3Txtee6201bIcQsExn3mgbn5vloqWvZDZDQoUARS5N12mgo5aRwj62emDWZa9kNkNJj32Ke5aqampnxmZma2fl9jAAAOiElEQVTQxeipe0YQNMcFHDIn4G3cMXtcVvWypQ0ueeaTuOM7+47Ognr40UNHX+/U/qyoDO4TGovYfzD6PVunV2f7wiJSaWa23d2n8rxX02MD6cyanj0wdzRIQLZppFEB5+ePHOLmO+9n/vCxhL44DxyYi83gBlItASIio00tihKsWr8lsjJPc+ce9960OrdTjcq/iMrPgOyLGIpItalFMSBpk+DSzFqKO1c/U03b26m2u5eiWjLdS4AoiU5EuqlFkVNUl1CnZUsbXPOyFaxdORnbKhgz44g7p4w3ePixQ0e7kqDZBXTDZecd7bpKo7v1kDRuEdeS6af1IyLVpRbFACSt6wTN1VnX3bITiF4OHI4tGhhVmbenqca9t9t4Y+xoYGpLs51q2teURCcyujQ9Nqc0Fef8YT+6EGDnNNIxi8uyWPgZvZL42qKmpeaZoqskOhHppkCRU9qKs92N007Eu3f9JRxJ2d3X/oz2e+OCxeTEeOT4QZ6EOCXRiUg3BYqceq3r1BbVekgbZLor56yVeJ6EOCXRiUg3jVHk1J0nEedwROshzbjDxHhjQeWcZ//tPNuwau9sEemkQNGHzgp15fX/OzLLecyMjTtmj6t4k5LxoNlKuHbNip6fKSJSBgWKglzzshWxM5ui8hA6K/yq7BJXlXKISLUoj6JAG3fMcvVn7iLqklY1D6EdHOJaNhqfEBkOyqOokLi4W8U8hO6kwe6i59nnW60SkeGjWU85xO37kLSPg0PmHeRC65U0CNkCXBFbtopI9ahFkVHSWki9KtWqrZuUJghkSbRL2rK1Ct9XRPJRiyKjpMowTaWaZQe50HqVN2uinZb/EBlOChQZJVWGF55zWuwWqGnOUbaoBL52+fMk2mn5D5HhpK6njJ48MR6ZYGcGf7/tvtTnWT59WyGDvf0MHudJ4EsSt5Oelv8QqTdNj82o1/LiWfUzBTWqLIOe0qpZTyLVpOmxJeq+C19kFrlMR1pz84d562d2HnfuON2V8MHHDlVu8FiZ4yLDJ+gYhZldbGa7zWyPmU1HvP4nZna3mX3DzP7ZzJ4asjwh9BMkOs9x9c138Y6N34w9JmrqadSSIVCdMRARGQ7BAoWZjQE3AS8GzgWuMLNzuw7bAUy5+zOBW4B3hypPUbor7KI48Ilt98XmHKTJeWjT4LGIFClk19NzgD3ufg+AmX0auBS4u32Au9/Rcfw24MqA5SnEdZ/fVdj4RDcH3vqZnVx9811MLG3gDg/NzccOoEdpLDINHotIoUJ2PU0C93c83tt6Ls7rgS8GLE/fNu6Yje3u6TTeyH9ZD7vjNLdSPTA3f7SbKd2eeJD+QBGRdEK2KKKqrMjeGjO7EpgCfjPm9auAqwDOOOOMosqXWdpEubn5I4V/dtpurs7tV0eVZl6JFCtki2IvcHrH46cAD3QfZGYvBN4OrHH3R6NO5O4b3H3K3adOO+20IIVNoy6DxHUpZwhab0qkeCEDxZ3AWWa23MyWAK8CNnUeYGYrgQ/SDBIPBixLIdIMEo83xli2tFFCaeKN8mB20hIrIpJPsEDh7oeANwKbgW8Dn3H3XWZ2vZmtaR12I3AS8Fkzu8vMNsWcrhJ67ZM9Md7ghsvO45qXrUi1n3YIo54JrfWmRIoXNOHO3W8Hbu967p0dP78w5OcXrdcWpg8/eojrPr+LAwfnOWW8wQmNRRw4OF/INNqJ8QYH5pIH0ifVHx87Q2yUW1ki/dKigBmtXTnJ1unVTE6MLwgA80ec/a3AcGBunkfmj/Ca553R90SkyYlxrl2zgsaihWdqjBl/ffn5fHf9JWydXj3SQQKiW32j3soS6ZfWespp+fRtqVoKY30u8dFYZJx0wmIOHJxnYmmDR+YPH51VtWxpg0ue+STu+M4+zfDpoFlPIgtpracBSJsElzdIGHDKeIOHHzt0NHdj/8F5xhtj/PXl57N25WTiJkqjXDFqvSmRYqnrKadeA9v9mJwY5971l2DWzIvo1F5EsH3XrBk+IhKaWhQ5dQ9sF6Xdn56UBX7YPXGpc83wEZEiqUXRh86B7bTGLH5ou3NXuV6tgrn5w7HnSjvDZ+OOWVat38Ly6dtYtX6LktJEJJJaFAVIewe/bGkjtpVgwNbp1ZnOedid8cZYrh3lNL4hImkpUOTQHh+YPTDHmFnqPImkBQW7WwFpBsvbeRN5ZvgkjW8oUIhIJwWKjN6x8Zt8Ytt9R4ND3KymxiJjyeJFPPxY7yXJo5YGj9p/ulO75ZB3ho8ymEUkLY1RZLBxx+xxQSLO5MQ4N77yWUwsXZLqvEsWL1pQ2a9dOckNl53H5MQ4RrPbamK8gXH8WEZeceMYymAWkW5qUWRw4+bdPYNE51jD1Tffleq8ca2OkPkAUS0WZTCLSBQFigzSdMt03pFn2ZnuHRu/WWqGdef0XmUwi0gSBYoMelX83XfkvcYZOv39tvuO/lzWDCRlMItIGhqjyGDdRWfHLvA3ZrZg3CBqnCHtBVeGdTLlgIiURy2KDNaunGTmez9ZMKA93hiLHVzuvmvfuGOWt6Qcu9AMpGjKAREpl1oUGf3Z2vN47+XnH20lZJ2BtHblZOpMbs1AiqY1rkTKpRZFDv327acZu9AMpHjKAREpl1oUA9A9djE5Mc6Vzzsjdytl1CgHRKRcalEMiGYc5accEJFyKVBI7SgHRKRcChRSS2qRiZRHYxQiIpJIgUJERBIpUIiISCIFChERSaRAISIiiRQoREQkkQKFiIgkUqAQEZFEChQiIpJIgUJERBIpUIiISCIFChERSaRAISIiiRQoREQkkQKFiIgkUqAQEZFEChQiIpJIgUJERBIpUIiISCIFChERSaRAISIiiRQoREQkkQKFiIgkUqAQEZFEChQiIpJIgUJERBIpUIiISCIFChERSaRAISIiiRQoREQkkQKFiIgkUqAQEZFEChQiIpJIgUJERBIpUIiISKKggcLMLjaz3Wa2x8ymI15/nJnd3Hr9a2Z2ZsjyiIhIdsEChZmNATcBLwbOBa4ws3O7Dns9sN/dfxl4L/CXocojIiL5hGxRPAfY4+73uPtjwKeBS7uOuRT4n62fbwFeYGYWsEwiIpJRyEAxCdzf8Xhv67nIY9z9EPAQ8AsByyQiIhktDnjuqJaB5zgGM7sKuKr18FEz+1afZRsWpwI/GnQhKkLX4hhdi2N0LY45O+8bQwaKvcDpHY+fAjwQc8xeM1sMnAL8pPtE7r4B2ABgZjPuPhWkxDWja3GMrsUxuhbH6FocY2Yzed8bsuvpTuAsM1tuZkuAVwGbuo7ZBPxe6+dXAFvcfUGLQkREBidYi8LdD5nZG4HNwBjwEXffZWbXAzPuvgn4MPBxM9tDsyXxqlDlERGRfEJ2PeHutwO3dz33zo6fHwFemfG0Gwoo2rDQtThG1+IYXYtjdC2OyX0tTD09IiKSREt4iIhIosoGCi3/cUyKa/EnZna3mX3DzP7ZzJ46iHKWode16DjuFWbmZja0M17SXAsz+53W/xu7zOyTZZexLCn+Rs4wszvMbEfr7+QlgyhnaGb2ETN7MC6FwJre37pO3zCzZ6c6sbtX7h/Nwe9/A54GLAF2Aud2HfNHwAdaP78KuHnQ5R7gtbgQWNr6+Q9H+Vq0jjsZ+DKwDZgadLkH+P/FWcAOYFnr8RMHXe4BXosNwB+2fj4X+O6gyx3oWvx74NnAt2JefwnwRZo5bM8DvpbmvFVtUWj5j2N6Xgt3v8PdD7YebqOZszKM0vx/AfAu4N3AI2UWrmRprsUbgJvcfT+Auz9YchnLkuZaOPD41s+nsDCnayi4+5eJyEXrcCnwMW/aBkyY2ZN6nbeqgULLfxyT5lp0ej3NO4Zh1PNamNlK4HR3/0KZBRuANP9fPB14upltNbNtZnZxaaUrV5prcS1wpZntpTkT803lFK1ystYnQODpsX0obPmPIZD6e5rZlcAU8JtBSzQ4idfCzBbRXIX4tWUVaIDS/H+xmGb30/NptjK/YmbPcPcDgctWtjTX4grgo+7+HjP7NZr5W89w9yPhi1cpuerNqrYosiz/QdLyH0MgzbXAzF4IvB1Y4+6PllS2svW6FicDzwC+ZGbfpdkHu2lIB7TT/o38o7vPu/u9wG6agWPYpLkWrwc+A+Du/wqcQHMdqFGTqj7pVtVAoeU/jul5LVrdLR+kGSSGtR8aelwLd3/I3U919zPd/Uya4zVr3D33GjcVluZvZCPNiQ6Y2ak0u6LuKbWU5UhzLe4DXgBgZr9CM1DsK7WU1bAJ+N3W7KfnAQ+5+/d7vamSXU+u5T+OSnktbgROAj7bGs+/z93XDKzQgaS8FiMh5bXYDLzIzO4GDgPr3P3Hgyt1GCmvxVuBvzOzq2l2tbx2GG8szexTNLsaT22Nx1wDNADc/QM0x2deAuwBDgKvS3XeIbxWIiJSoKp2PYmISEUoUIiISCIFChERSaRAISIiiRQoREQkkQKFSA5m9hYzWzrocoiUQdNjRXJoZX5PufuPIl4bc/fD5ZdKJAy1KES6mNkfmNldrX/3mtkdXa+/GXgycEf7NTP7uZldb2ZfA37NzL7byobGzKbM7Eutn09s7RlwZ2tvhKjVb0UqRYFCpIu7f8Ddzwd+lebaOH/V9fr7aa6Pc6G7X9h6+kSaewA8192/mnD6t9NcbuZXaS6vcaOZnVj4lxApkAKFSLz30azUP5/i2MPArSmOexEwbWZ3AV+iuebQGblLKFKCSq71JDJoZvZa4KnAG1O+5ZGucYlDHLsRO6Hz1MDL3X1334UUKYlaFCJdzOwC4L8AVybsV/Azmsuax/kucEHr55d3PL8ZeFN7N8bWyr8ilaZAIbLQG4En0BysvsvMPhRxzAbgi90D3R2uA95nZl+h2S3V9i6aq3l+w8y+1XosUmmaHisiIonUohARkUQKFCIikkiBQkREEilQiIhIIgUKERFJpEAhIiKJFChERCSRAoWIiCT6/+4Ld5s4wsWSAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x432 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "figsize(6,6)\n",
    "plot(y_test, yp,'o')\n",
    "xlabel('z true')\n",
    "ylabel('z predicted')\n",
    "title('mae = %.3f' % np.mean(np.abs(y_test-yp)))\n",
    "xlim(0,1)\n",
    "ylim(0,1)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.15"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}