Fisher_Vector_feature.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "nnkKq8DV3ThO"
   },
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import datetime\n",
    "import matplotlib.pyplot as plt\n",
    "import scipy\n",
    "import numpy as np\n",
    "import torch\n",
    "import random\n",
    "import scipy\n",
    "from torch import nn, optim\n",
    "from torch.autograd import Variable\n",
    "from IPython.display import Audio\n",
    "from scipy.io import wavfile\n",
    "from sklearn.mixture import GaussianMixture\n",
    "from collections import Counter\n",
    "from collections import defaultdict\n",
    "import math\n",
    "import copy\n",
    "import random\n",
    "import pickle"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "htG5DLW1kFz8"
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "lpUvpUlKEt2E"
   },
   "outputs": [],
   "source": [
    "make_dataset = False\n",
    "noise = 1.5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "53FeW83gQ20S"
   },
   "outputs": [],
   "source": [
    "\n",
    "file = 'Noise = 1.5'\n",
    "if (file[0] == 'N'):\n",
    "    df = pd.read_csv(file + '.csv')\n",
    "    vals = []\n",
    "    for row in df.iterrows():\n",
    "      vals.append(row[1])\n",
    "    vals = torch.Tensor(vals)\n",
    "    vals = vals.squeeze(1)\n",
    "    \n",
    "else:\n",
    "  df = pd.read_csv(file + '.csv')\n",
    "\n",
    "  class Value(object):\n",
    "\n",
    "      \n",
    "      def __init__(self, date, val):\n",
    "          self.date = date\n",
    "          self.value = val\n",
    "          \n",
    "  unordered_time_series = defaultdict(float)\n",
    "\n",
    "\n",
    "  for row in df.values:\n",
    "      if (math.isnan(row[1]) == False ):\n",
    "          \n",
    "          last_date = row[1]\n",
    "        \n",
    "        \n",
    "          unordered_time_series[row[0]] += row[1] \n",
    "  #     else:\n",
    "  #       unordered_time_series[last_date] += row[2]\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "  # Separating the data by date\n",
    "\n",
    "  time_series = []\n",
    "\n",
    "  for key in unordered_time_series.keys():\n",
    "      year = int(key.split('-')[0])\n",
    "      month = int(key.split('-')[1])\n",
    "      day = int(key.split('-')[2].split(' ')[0])\n",
    "      hour = int(key.split(' ')[1].split(':')[0])\n",
    "      date = datetime.date(year, month, day)\n",
    "      time_series.append(Value(date, unordered_time_series[key]))\n",
    "      \n",
    "      \n",
    "      \n",
    "      \n",
    "      \n",
    "  # Sorting the data\n",
    "  time_series.sort(key=lambda r: r.date)\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "  # Creating the data list\n",
    "\n",
    "  vals = [] # list which contains the data\n",
    "  days = []\n",
    "  months = []\n",
    "  weekdays = []\n",
    "  value2index = defaultdict(int)\n",
    "  index = 0\n",
    "  for i in time_series:\n",
    "      months.append(int(str(i.date).split(' ')[0].split('-')[1]))\n",
    "      days.append(int(str(i.date).split(' ')[0].split('-')[1]))\n",
    "      weekdays.append(int(i.date.weekday()))\n",
    "      vals.append(i.value)\n",
    "      \n",
    "      value2index[i.value] = index\n",
    "      index += 1\n",
    "  vals_backup = vals\n",
    "  print (len(vals))\n",
    "  vals = torch.Tensor(vals)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "QoPqgOnPu9ny"
   },
   "outputs": [],
   "source": [
    "# Rescaling the data and separating it in train/test/valid datasets\n",
    "  \n",
    "class RMSELoss(nn.Module):\n",
    "    def __init__(self):\n",
    "        super().__init__()\n",
    "        self.mse = nn.MSELoss()\n",
    "        \n",
    "    def forward(self,yhat,y):\n",
    "        return torch.sqrt(self.mse(yhat,y))\n",
    "\n",
    "\n",
    "      \n",
    "\n",
    "test_proportion = 0.2\n",
    "validation_proportion = 0.2\n",
    "\n",
    "val_break = round(len(vals)*(1-(test_proportion + validation_proportion)))\n",
    "test_break = round(len(vals)*(1-test_proportion))\n",
    "\n",
    "\n",
    "values = vals[0 : val_break]\n",
    "values_valid = vals[val_break : test_break]\n",
    "values_test = vals[test_break:(len(vals))]\n",
    "      \n",
    "      \n",
    "      \n",
    "vals_normed = torch.Tensor(vals)\n",
    "values = torch.Tensor(values)\n",
    "values_test = torch.Tensor(values_test)\n",
    "values_valid = torch.Tensor(values_valid)\n",
    "\n",
    "\n",
    "\n",
    "max_vals = max(vals)\n",
    "\n",
    "vals_normed = (vals_normed/max_vals + 1)/2 - 0.5 \n",
    "values = values/max_vals\n",
    "values_test = values_test/max_vals\n",
    "values_valid = values_valid/max_vals\n",
    "\n",
    "\n",
    "loss_fct = RMSELoss()\n",
    "\n",
    "values_test = values_test.data.tolist()\n",
    "values = values.data.tolist()\n",
    "values_valid = values_valid.data.tolist()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 283
    },
    "colab_type": "code",
    "id": "y5dA8bntQmxi",
    "outputId": "bda13f1e-a5f6-4e71-b1a5-19a131369dec"
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.PathCollection at 0x16678168518>"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJztfX/MZsV13jNAWSn1ktiwSRExXmi2QWtUitkiozaWKiMHowRouq5sWgc3kWiquC1yI5cVpYrY/GOqNrTGauQ4KLhKYye0oaiSgyFNW7Vg7MUBG7oLu1BHJmzZDxLxfW7L59JO//juXc8O8+OcmTNz73vf80ifvve9d+bMuXNnnjnnzI/XWGuhUCgUimXhrKkVUCgUCoU8lNwVCoVigVByVygUigVCyV2hUCgWCCV3hUKhWCCU3BUKhWKBUHJXKBSKBULJXaFQKBYIJXeFQqFYIM6ZquALLrjA7t27d6riFQqFYiXx5JNPvmqt3ZNLNxm57927F0eOHJmqeIVCoVhJGGP+kJJOwzIKhUKxQCi5KxQKxQKh5K5QKBQLhJK7QqFQLBBK7gqFQrFAKLkrFArFAqHkrlAoFAuEkrtCoVAsEEruCgURG1vbU6ugUJCh5K5YFDa2tpuQ8MbWNu588Jlm8hUKaSi5KxaDja1tfPKBp/HJB54WJ+A9u3fh8E2XA8Bpklco5ozJzpZRKGqwsbWNPbt3nXFtz+5duPvgFac/S8sf5R6+6fJq+b0Qew7F8qGWu2Ll4IZIfOzZvStKZlRrOyTfvTYXssw9z6jzsZObnTRSzAkkcjfGXGeMec4Yc8IYc3si3UFjjDXGHJBTUbFK4BBoKUqs59SA4F8LyZ+DxT7qOcb973zwGTx+4tVo+j27d+G2a/fhnkePaxhpDZEld2PM2QA+A+CDAPYD+IgxZn8g3W4Afw/AE9JKKlYDKQINpTt2cpNEOjELnYMYOcd0joVkKLq1gFtndz74DADglmvehb/7m38Qtcw3trZx2YXnTT4oKaYBxXK/GsAJa+2L1trvAvgCgBsD6Q4DuBvAG4L6KVYIIQKNEfNt1+7Dp3732BmTnzGrmjqBmUsjHUPn6Obn42LU0yXrH/mh3fjRP7M7q9ue3bvUcl9DUMj9IgDfdr6/NFw7DWPMlQDeaa3994K6KSZALQn4xB4jv8suPA93H7wCdx+84jT5uGnH/9QBo4Ro3TJKUBseooSI/PL8/3f+xP5g2MXVjVI3Sv7LA4XcTeCaPX3TmLMA/DKAf5AVZMytxpgjxpgjGxsbdC0VXTC6/DUd3c2bIz938jNFRpQBgzMIpORwwSX22JJKN47O0YkSdsm9h9SAo6S/wrDWJv8AXAPgYef7IQCHnO/fD+BVAN8a/t4A8DKAAym5V111lVXMB6c237B/+/NH7NGXX6+WcWrzjeA9rqySeyld/GtcnWoQ08e/J6VT6l3E0ofqhyND0QcAjtgMb1trYXbSxmGMOQfA8wDeD+CPAHwNwM3W2mcj6f8jgF+w1iZ/IPXAgQNWf0O1P46d3MRlF54XvMdd5hdKH7t254PPdJ/Yo+onIbc2X4mutfdD6YG3eko6GTsvGGOetNZmVyRmwzLW2jcBfBzAwwCOAvgta+2zxpi7jDE31Kuq6IGNrW0cO7mJW+77anR1hcTSwjlt/KGueOGgZhI1FzrhlEPRY5RJnYweV+Gk9FKsECjmfYs/Dcv0g+te+2GXGpd7CnedE2ZoVUZJiCkX3gjdy4XIQvd9OZzQyhzCL6nQ1Bz0mwNADMvoDtU1gGs5uyGZ2knF3lYddSNSzXNxLGIqOBOa7rXU5qPQ/ZAcStmlz1WL0ORt7GwgqQnwtQJlBGjxp5Z7G0hOXPZEjUUesk5rJ4Z7po9Zq6HJV/dz7cR1zQR4bbuJla2Wex5Qy339ELMAU+Ccw1J7RklN/HgEJaafs3pLygihZPliaD1/KNbt6uEvVUyVx7G+ubt2qffddKHPqbLH5bGxzW8KBigjQIs/tdzbILb0r9aiO/ry6/bqX3qk2CLOxX6lrbKa5Y7cpZbcuvXTU+rUf6+tlymGLOrU91B+t364yzK56dcJIFruSu4zBbczpeSUEEHJBF+JzJb5SyYyU+vRfVAmNCmyS0m6BalRw15cWa1CXeu4Fl/JfYWR20wyFyKQBpcYU3kodRWzvinkW0rUNYNm6xVAJXFwSR1LwPWclgAl9xVHjMR6WyoU60tCFykre7z+sfueIIdXciGW3MSm+7k2fBHL4z6TtAcz3j/68uvkOqaU07KN9ghNzRVK7jNDKGaaSpOSIxWyyZVPiZtKdLAcqXF0tnbHEs6Re8w6TR3BwLXiqXXjl50asGL3c8ilT8nllFXjbZaEYloQ+9wHCyX3GcFvjB+774m3kA/VcpWM2eZkj9f9z6kwhlS5uXsxWTmCzsnNkSflOscTCF2n3PfTcCdlS+5z5VO8TUoarq41A98qeANK7jODpOVeko8ii1p+aHCq0alETs7LKJXr3pPwRkLXJTwdVwZlJZMkeZbIpwzeEoRa8pzUCfS5QMl9pphLwynt7CN5pjpDznKVWHUTGlwkLbDS/JLhB44MCcvdTdfq2SXySMvkemlTg0ruuolJCJzDmabaQs3Zlu7n83/iLbTBxpcN4C33QoeXcX5EgnosQO1hZZxt++N36vul1jlHRuykT265Y7qSuqM+F/XIA4lNaKUb4ygbuWYPygjQ4m9JljvH0imZBJKIlda6wKl83PCHa2XGrPCYzjWx7BS4YSSOXlxQ2tNcrcoUxoluykqmEu8h5FHWHK/QMrRWA2hYpi9aubyUCT7uqX8pQu2NnC7UeYFYPVCJsnaTl2TcOCenN7mUkGMofWy5Za5catuOLVSI6VLTp6bsO0ruDVH7YrkDQWppXok+PqFKWZylaE2IEpZ7iXfUEj2JXWqFVmlMvoVXnJujmIPhE4OSeyPMwWqS0MO34udiyc8RJd5RD/TwEsb7VA+qhQ4tyvKNplYGTguZVHLXCVUm/Mmm2pMSKeVR9AB4kzzuD1PfffAK3H3wijpFBdB6kqpmgu62a/eRJgy5v340/i/5ZSeJyXmqnNCkObccjg5u3UhjfJ/jyaGpc+RrEDoBtCsoI0CLv1W13F3UnpQoiZj7XJOfmo8SDmkxgcZBjXxuXk6cnzPJGJIhAYrVWlMWN4YtUTcUnfwyW1jvrocs1b6hYRkeSit9DsQ+ojYuWhq7j605d89CoW7Dl9JNWkYL4q2Z8ygdpGLXS8NO3MlRKrh1U2PIxAi4Z+iJAyV3BnrH0TmosSZ6EGJMv9CkbaleoRgpV89VA9Xq58irPU8/RH5zmIuoMWSoRJ+S0RtK7gSECKeGhKR0ca+1dEspaSQmbakyU9dyB1vNdXDmgurplBwg1mKw79GOKJDyyqj9f8o2p+QeQWpkdq+1fHmhDpIirhbEPqXFRSVojpvccvDrFc7x21+p7NQg2SqMUqLT3MDt81OFZJXcA6CQRapTS1k+nB+daIU5drZegxtFD84gnyK0kjPqa5AaJKnWPneQSZW9auAMfmq5z4jcra2zSKRe5hwa/xx0oKC03iVIcvyfMwY4seyUpyiF1CDJIa+Sum/tGbSAtJeSS1/77FRyX7t17r0OkmqhgxQkDzBL/cK9hOySQ6xS66ap2LN712k5r31n+/Shab6c3H4Dd927LyN1wFoNYnU1HvpGlVFS99yD6KRRs1+gpI1wy4m1oyagjAAt/npZ7rUWH8WlnvuEUam8kKVHnaNoEV7g5o/pSckzwl1r/djxjSK9T21+b2I8d5RE6p6UNzIHuCEiSS+rpA2WeikjuL+B28tyXzS5hzp1rhO5S8xSL1yazFI6TRHfcwkp1XlCcxRS+lJDHVRZqfmO3LvmrE5JhUW4zzC+h9Kf12uJ0rpw71HCWpQ+6/fD1Lvm6Jp7xik2Miq5D3AbEKWD+Gm4jYFzP4dWMdlcme7nWOfIkWHLwY66EzdG4iF5lNMKpZ6JM/iN5N6iTmu8IsozUN5V6rlCBkYsHeW6W5cUpH7P1k+X0kMaiyb3UsuNmldiFG5tvbYAdzCheDTc/FQ9x/8pC4264WYKq5hiqbb0CikWc6jsmjbCeScjOEcEU3RJ/XaA/z20UYza3lpiseReWomcuFjPVTFTutkpq1t681IoTU1YixKzrXXJOZBqL25svuasoJj81Hf3esrTrdGD+06kCZNjkFDbVe66NBZL7ta+1RoIfXbBjYv1IuWeoz237FYxRP+YVYouPmJW1VSg6k9tV6EBT2LA5daxrwOnLCn0Kn9KI4uLRZP7CN9FTDVeykqHknIlZLVCKdm0GnRSgyy3rJZeRQu5EnVaa3T0qOPUc/YgUKm2O3X+FNaC3K2lWe4tGtzU1klOXo/QSgl6Tz65mLrjtx5YpI2OUlnU0MdUAy0lf08viQtRcgdwHYDnAJwAcHvg/s8B+CaApwD8FwD7czJ771DNNTjOi6hxWaV08HWRktcDLRu/RGgkl7/mOAJp+PpIhapi4aFame7nuS3tdFHbjlbCcgdwNoAXAFwK4FwAT/vkDeA85/MNAH43J3cOp0Jae6aVS+kYbtraI1R9HbhLvmIdpIdFNKXbGiOdXoSRe0e9V1NIW+6ccGdtOZLpeqGnFxKCJLlfA+Bh5/shAIcS6T8C4Es5ua3JnVvZJZtUJGPWOXJPkUaLsE6s/HFp2tQTwTUrOVpPxMbeUY9BuJWVPRVybW0q3fz6if12QQtIkvtBAJ9zvn8UwL2BdD8/WPjfBrAvJ7cluVOs4Fg+/3vP+GquAVNjmTW6pepOcst4DWo60TipWzPBXjJo+3XaaoCUlDcXi1nSiCoph6IH91fHaiBJ7h8KkPunE+lvBnB/5N6tAI4AOHLxxRc3e3gJco91Sq6cUnBcbS7Z5Cyhmz/7+OwsJUk8dnyjarKwZDCNDcySkCSW0j7UAlJtP1dGrl+kDBuup1Oj95RhmbMAvJ6TO7ewTCyOVvLSJJa8tTy3ptRy58qaM2rIds7PLPUMcyH3niHAXL+QssxrZUiS+zkAXgRwiTOh+m4vzT7n809SCp/LhKqL3Mst2VYvEdZpbQmUyooNiHNC6bOUuuhS+rQEtS1PoWtLD6d04Cs1tkJyuLqEIL0U8noAzw8x9TuGa3cBuGH4/M8BPDsshfx9n/xDf63IvWWDlPAGSsrkHNDUu0PWxJFLwyKctJydrNRQWOl77WmFUnSZE1rHq6mGSEyH2nct+Wxrs4nJxVxcSRe1xE6V0bpzUEG1Tko6UcnzlR7Vy7FqYyTR0grtgdQzSnouvYwTqiGSGtRLyyvJH8PiyT32UuZG7qWosQ5Ly5ME1eOozUPVhZKmxhr3z8zhTsbPrc36hBuzQkOEXOLxzDEEVOOVcsvlYNHkPrdGwgE3rNADrSwnqgWU6lBSenAtd06a0TvwCZ5D7JLLWWvSxvKF3lloEKOEtOaK1IA23pcsoxSLJndr610njttdIj+Wt8e63BJ30G3QLTf5+HWQ2owktdsyZkWXWukpgi81OLhkzQk1tZz05ljurZHqD5w+kJNXC7XcC0C1GijpWhGO1OgfI+KUC02RW/rDFdzwi6tbye+KcnWTepeu7rF7LmK/7NXTOGkVXihFrg64csbPsaXM1F9zmrpeqFhLcreWHu/LvegUIbUCVbfxOahb8LkWHvdZxw4UIzyqm967Y5WWRyWCMd24G9YdzFoYDxzjpSdCJFz761cpj8Q1fsZ+QhmMW3hYLbCW5M6xZlJpQo2mZUzULYPS6DmEWKpzzqr2842/HpSSm9I71VnnCKpurSz3UD6/nDmQe46EubJS36196xEZsb5b2qepfbOlJ7B25F5TmZRGI90IUuXWls/VIwbOQEM5JTM0WeWm9yckS2PYc/zFrR6gDM6cQTskN/Sdm79EDpVYY89cW34oT27AaDV3tXbkbm1ZZ201wlJd9lzHGQmwdIlnq4453suFumKdIUTyNZtMRlB/UpFDFnPyIDi6HH359WC7oYZEQu8l126p+rXyhlu+q9DgGarb1H0JrCW550Al1BblUq7HOlKsk5boId3gaskxR+Yh5Ih7bpa7BFIhhlSe1FyI+z90f2x7FAMkFVKkTJj37pNcUIwPirEjASX3CHo3mpIOGfoupXcJmUqVF7rew+JrhR7eXsjDqZXp3st5QdTQWOw9UjwpqndAeSZpuHU+F2NAyX1G4DaKXqQ7FUn6ll7JJPLU5C5dfzHLMEbqkuW6f6H7PrgrXaiT8qH/oftS+x8oaUqP/W4JJfcGqIkPUi0oqmvH6TC5+1NY8BzL3Q8vzMV672W5h65JboX/2H1P2Js/+ziJxNyBWRqx56L0CW7f5Oz9mLqd+VByFwa3QXA3EoUs91g+jqsrMbFYi9oyQnMOU3W4mgG+Vp6bnjpwU2VySKx1W+GWWdK+WjxDrzap5N4ANY2n1NKpsdxja6x92ac236haGtcqz5iv9ZEIfnk5XUoGeG6a8VroQDIKsUutrpLAFJ5hb/T0JqnkfhYUJGxsbWPP7l3k9G7aja1t3PPocWxsbReVGcp32YXnZfNfduF5OHzT5adl3PngM6dlud9f+842brnvqzh2cpOkkytnvJYDp+78fIdvujz6vNw69fP4n1PPNuoSeiduvfppY4ilGfV4/MSr5PcigdDzl8oZ/x87uZmUWVvWiNL2FUKJTpT33R2UEaDF3ypZ7hKjcon77U82Sk8iSVjuLSwWTniAs2nFzROa3/DzcmLAElvqXYzPxHm2knRS+dz87vLd3M7lmpVktXpK6JSS5d6T1B3rGpbhxjk5k1jccmrhuuat3D7KZBVVzxb65DCSCGfzCCd2H2sLqYFSiiBbvOvUd+lyKHVBvS9ZJzlZXGJPhdakB31r15DcY5ZYKN0Yj4xZbi13LpZa8BKkkdMp9FxScwWcjpRKG5LrrvRwvRAp68tPRz1moQY9BnHO/EHN/dK0uTmjGvSSpZZ7JTgE6K+84JALN01KT26+XvCJ3V+VQ60bCpGELJqUpRz6YQj3vpuGswEn9GMbueerud8a1IG05F1S78cGYu7BW9IT/VO/GwmsFblby2uovVZexHRoLZtaF5RORj3MK5Y/9d1/FzlPytWZQio5V9n9/tjxjUWsnW9h7VIGM/+dlfxASsgo8FETv5fyqqbG2pE7FXO0uqRcY85uz9TOuxjJjv+l6ojqNeUsUa4+Me+t1bunypU647zH4VWhcv32UuopxOB6khRPgSufgjkMDGtF7pJk03tkT5WZs1JT6SnlUix3v/zaOuKQN1UeZ07CDd/0ID3Ouwt5MiVlxr738BpT96lhshRGYq/dsVxat3OYT1kbcpcm5B7ETumAoQZcKl8CKassZ1m730smsCm6cd3y1lZ6SVmUeimxhlPyeiLWzkuIuWbHck1d1OaReA9rQ+7WzsNVoiJlNfqEkDuuNXSd23hqLJjYEs2UHlKWe6mcksGSo5NkvDs3OPZ41y1QY7i4Mmp16IFY/6jBWpF7CL3cUS5yIYGc3jGLsNRSLLUiYpurelnHErHp1h5O6Psc5bfWOVReanMTV9aU+SmypOtzLcg9VZm1MWKJBh+TIdEgSywCKas5l3eKDlP7fqXS+vmmDoXkUOsNlEJiY8+cBnmOp1qLxZN7rjLd/yUdv7bB18gocbMpxD53ovGR83JGcNa0l+hQU2+PHd+oKp+DGpKjypFcRtzC65C+z9HB5xrpnakjFk/u1srHJHOyQ9dT+fy01DPYawcSim6h76VoNWBQfl5wfPZQ3UoNaKX5/U1gLQfWHoM39fdp54IS77amjNDnFnW1FuQeQgmBUUnRfWExKzF3zgRFP4nBiJKecwYLVY7kgEHtHLXPIJnPR8815z28spbELtlXQ99bvYPY4NHqfawtuXNwavONrOvkNwj3tDsuMac8jSnCJRSLmNPhpGOYU3kWPUigRfrWckrkcUKMtTtMpQyW1qjVQck9gzGeS5m1D5FyTSzPH93nGg8v8SqkiL31RBcl35SQDClJL8/kDPrSSzV7tL9StKiTEJTcCaixDker35eVI2tOeKZEJ8q10nJLrHkupMM70rJaluM/u1Q9S3lWJe9mSu+rJ9FzBlG13IXQgjRPbe4cM3vg8JfP2A7tHj2b+31TaYQ6roQbG8pTogtF3liPU1vNIVAGSalzYVpZxtbWHeEssXSxVdmhdiTZlqYKT4YgSu4ArgPwHIATAG4P3P8EgP8G4BsAfg/Au3Iye5C7VIOMWdruxCrFcm9NWlOFUThyU4NPiNxb1hn1/VDbkZTlXnvMLfU5Ynkp746LnDebKy/1rkInYUpa7q2MoFKIkTuAswG8AOBSAOcCeBrAfi/NXwHwfcPnvwPgizm5q2K5+43E/Z5rsBwrWhpTWL8l1q2fx11dIlVnsffDWTZbEhYogZQHkHuOULrcYgFK2e5/yv2UIUQJbboeTy5fDebkTUqS+zUAHna+HwJwKJH+SgD/NSe3Jbm3GrFD7h7H2utJ7L2X3UlYt25HdX9ZqVZPzgAsgdr6LwmdUJ8jls6t+1Jiz/1IiltOSA+Op5vKL2UMSLYHyb4oSe4HAXzO+f5RAPcm0t8L4B9F7t0K4AiAIxdffLHIg4YaRMtYG+dFcRqrNKTXI1OtLYlyQpuXSus8lVd6EPS9OmqeWp04A2tq70BNPbj9rmSAr/WcJL087txHTp50/F+S3D8UIPdPR9L+TQBfAbArJ1dqh6q/TpxK7q2INmUZlayN55Y5fpfcjt/CC8kRbskvQJUSYw1iXh1FbqxNSA1qoevSvzhV4jFw00h4hZx87juVMljmarmTwjIArgVwFMAPUgqWPM/df/kSJMB9Ga4eoaNwUx25RYw1JC92kiO1HCmUWGnSJCeB1EQelYxSA7502/Dbp5+2Vn4qHdVyTdVFKE2uD+faUs4QmuP+E0lyPwfAiwAucSZU3+2luXKYdN1HKdQKx9xDL7FmlKf+cLN7b0yfItCUhVCDMX8qFOP+CHSvRpur4x5lUp63ZOAIeRi1Mn3ZY3uqQajNxUi/lOApaajedG45r/s9R8yhcEisjJrn6w3ppZDXA3h+IPA7hmt3Abhh+PwogFcAPDX8PZST2fKXmGoba0lYwCXzVgeEpZA61Gksyz2hUGpASd2nWkQSVqQvdyTG2G9uUvVMpWvZ8SUsd4r8UgOjVXrKe6ix3GPWfIleJWkksBabmFINs6aiS+S5gwI1TCJluY9IkVjNBFHoWajueM5alrQi3TI5li/VJZf2tKTTcnUonejr5flJl1uan2rcpTw5SSye3FPWX42VUKtTSG7uvPGaLf0Uj6WGlHIyS/KG0nFlU8puISP3TDUeTW+UWrFj3inA7duSXmHuWqvz230sntytDVt/HELpFXfOjeilREkJFZSQkSTZSntVnHJayUnVZc9VHb3y9pRZglibb0m2NUZPLdaC3EeE4nCUSm7hQpWSGTVNqQUd+l4yYEiEJCSXaeZkpvL4nyV0axVyqTFGWhgyc/FAWnirnLJr7pdiLcg9ZqlTRuyeDV7Ki6glsZCMGJG5992121IbO0LXauOpJWVR62MObWUulnuorqYk+Sms5xxaeg2LJ3e/g/rXKduoW7uqHEKkNlAuGdQONi65U/LWQFo21RKv9aw4A0vq/lys4RSo3l4KrX/NqWUdUiMCLSdXF0/u1sbDL5JWZiligw83b206CUt4LmRTGpJoidQAWroyI/WdmrfXIJwb9P3rPX6HtZVRRzGWXONSsmwXa0PuEptTWqHE2uXmoYZXcvkl0Upm7QRWzoKv0Y1znYrSs905A0uNYUDpf6H7Ukcalw72JQZRrMyYJ1NaNgVrQe7WTk/gPmosOD8EkpOdKmOU1XJ7e0yXkmegyi4NB0iEE3oiZwGG4G+8y8nPha0kVv6UeK4pHdx6KW1nKZ1iu9NrjciaAcnH2pD7nFBrXaaIsaSBuY041lhrOl7seuwZajdScawkis41z18LihfB0S23l8L/Tnn/El4fdwCl1MGoN7WvUA2u1IBaM4hJGxFK7hOh9gXWdLYYuGfl5MotHWgoP0YeAuXAM27d9LDaY0Ql7UXESClFghIWeWowiRkWKXmUOqC24VRdp/KW1n9ut69kO1NynwmmdvlL3PtQ/hqrckxfs9U9twlMYvJSEiHrkbLksrSslBWaq3fu4J3Sd3zP7u8JU/SnrHCb43t2y+hB7NauOblLVaYEoUm4v7V5WnoTPeRQOs3Ug2gIte2HWkYuHFNCmDW6+pY7BdQYes5ImVs7iHlpNaCS+1lYGDa2tnHng89gY2tbRM6xk5vFcvfs3oXDN10OAG/JO34ukRvLE5OxZ/cusuwW+Ue9SuWE8vl1UKujL1sCo04UHVNlpu6NbSwke7yWKjeWn1qfbjt2845/VJz/tnTdjH3xnkePB+tjY2u7uC+1hF+/UvxEAmUEaPEnsYkp9r3UbfPzhX6MoRSp+KqE5S4h05ctYQW1jG33toRL9EhZxrm48HhP8ifaJOGGU1rtXHbvxbwBqofslzFV3fay3FeS3HPxzBJZnKN6Ofdz+SQbVm19jJ1DusPOkZhS4BI7hbz9e5wfWp/yZyNzZXLacao+KHklliP6vNGL3CXLWDS5WysbzywlWgqJllh1nPKl7oWsHwliD3kTq0T4HNJyv1NWT1DrgELsvdfrcz3FmNHA0X2Mt0t50ClZpfyS61sSWDy510DSfawl71hjoOjAsWTG75z10BLwO/E4cdbqvBrqu6HKql2XLwWKhV8rm3ufS9ASljv3d4BrvHnK81F/sU0t9w4IWail+anpOfc58lOdLtQ4Wx9olILvEZRYbpQyYlvwa72jFqTK0YG6rJGrE8VIiK2h53yXgN9+cqtQagfn1Hdrw+fk9PBK14bcuRU5Vr7kiy9FTAcpogt9lyyLYz2VbHoq1cn/TH3XMbc/9Z6kJrFzGOuasmZdekI4F0qRJFUO/Hedm/togdCgJx3e9LEW5M4lDT99awuWQlq1L76mAdUOcJzf4ORYmy3IgULsV//SI/ax4xtFOxqnILSaNJT8KS8399w9PBofsf42xXxES295LcjdWnpowk/f+qXH3EbJMv0ypSDaAAAXtElEQVQG1MILyVl1ElY3xcUulc1B7YRdaKCdguRCoOpBscpDeebynCG09qhiZbhGpE6oCoNCXC1ePOWFS8jPlVVbbqmL32qLeEl4hypbCr61O9UcR0wnanqOnCmsYy5a9XOqcSVZvpI7Aa3c/94eQSiNlE698tTKnpJ4UtZtr0nsFoM7R45UKKgFWrYB39urCZNSsXhyb91Y5ybTlZ0j9rnEQFuBQzCtQwYU4qg9lTNWrvu55ORHimzOtZryJCeBJdJT8oRCWJIbAGNYNLm3dP1b5JVC6rlLQjXUNJx03LSU/P6zcZaKctsJ9efRuESdegaJ/D651/QRylxRjNhKPYjWoT+JeqDIVct9AsudatGGXs6cYopU/ajhG2oDliDUkvyxZ4vJrLEmY7/vWUpkKYweRc7KG9NR1nSH8pbqRZHFrWv3OaQNAGq5LcvI5ZPij8WTu7V8NzzUGKluVesXFpKXKjNkqZXoFqqvWqKI5fdjz5z3FELNVvaY7NhmnZpBI1buuJQ05S2Mz9NjUjZXfxI6xAarWFoptCzLH7Raz/0sntxjBMdZez3mCZEkpTGEBgcuYtZ2zgIen7HWGklZpdLWlb91vKaxhwaLUJkUGdTBQxqUOpaoK65OIbheTc07y5Xj3m/1zC2IN/ZsnGemYm3IPfbTYhw5Ncv9aq3GlAWekuV2eqpFlbIqYrI5gyS3LOnBQ1IGt25bgtOmKYYGxzNzr4/EThn8S9oDVafSdDE9xmeh1IH7nVIHLQaoxZO7tXKWVa3FPX4vsaJrVzlw64DbEajpJFfoxMinh+Wa06M3YmTkp7H2zGWXKa+POj+Tup7zLKlGhITFXDNp7g7ilDpw01N+SrBFG1oLcrdWhphr5UnuEC31AKaGlA4x8nHDbXN43hE9dEkRsU9OlF22pZZ76F5JPqoMKiQs95Sc2OCW+8H3Vm1jbcjdWrk4aymxptx36cGnF6R1qB2sXNd5LstgW7ndXH0ooZhS2bm0LSzvFmRP8XwkynHvtdpvsjbkzrGaKZ2gloxdS2SKzi/VwSX15srLpSsh9tLt4ZSlilxwLOdSeaVyuM/K6XeU+7VtjxN2kigvp0usrJryRMkdwHUAngNwAsDtgfvvA/B1AG8COEiRKXkqJGXSqwfR+uGD8VovSFlS7n8pcIi9xURmqZUlrQuXfErkufdieWLXSryUkrpNobXlXlqelOFU20/FyB3A2QBeAHApgHMBPA1gv5dmL4A/D+DzPcndWr5V0AIxyz2XliqzVBfu/RYWewlix+5KlEENabhhoFpQXHMpyz3nNVIGF86gIHVolnSdSCNWbzXySiFJ7tcAeNj5fgjAoUjaX+9N7in0stalt8OXWkm15UuRGaWsXL7Hjm+Il0Hd9BTywLjIlUHJW1KmP6lKtdLda5wz+iU8mxhxcver1Nyn5KP+pF4LPVxIkvtBAJ9zvn8UwL2RtElyB3ArgCMAjlx88cXVD0lBbeOgEF6LjhvrlLWx6xK3myqjNJ2PknkUSrqRiCjWF3Wgo7wnbvuosRCpG8Viu7ClSZWK2ABL8QgoRkuNoRGryxp5NZAk9w8FyP3TkbSzstypiL24U5tvkNay9gLHspa0GFMufCu0kO9vwqmFhIcVy8shrtS1EGmGLG6fyDj6Uu6VIFdWbgCtIeISGZLGSAoalmEiRGDjDxJzdoDWlp1KQyUSv8HXNKhSi6XXYFgy2E3pnVAGSqrlX+LJpeqAQl4xQ4f6492S4JB/Lr2ELpS5h7mR+zkAXgRwiTOh+u5I2pUldx9UYqRaWKn8NSEIP7/EYVqUcnPpe+wcDlmgNfJ8udIEzyGcXoMLN29sR/XH7nvC3vzZx8XamZRn5X5uHW5K5RkHP4kBUHop5PUAnh9WzdwxXLsLwA3D578I4CUA/xPAawCezcmcktwl3KsxnWsJlbrnUiQYIzuK/NSBUNINnZqfQrCtDvfiDO61HgFXd6nTGUs22aTaSI1l6tejtEESIvecRyS1BNaVNyvLvdXfVOSeazglVlvOvS5pJDWdjCvr6Muv26vu+nLQ8pKywktALbPV4V7UZ6cQklR7i505XwLpdy2RP6abhB7u4JPTNbQcNzV4SRmMFKwVude4phLyUnlCYZJc/pxVxV1NQrEspS13CqRcb+lYuvuuUpYt1er13yfXWvbR8qRKKS+yZTmU+otd48xl+BPOsbBKbwNobcidUrESrnANQtZCjsxjaUp+cFnK1Q3Jrclbe6KeT8KUdkBNExpA3UlDKlG4eak6TIVe1mdNHYTq3e9TNW0qla7GcpfE2pC7tflKr93uLU2IMcvATxebfCu1UKWfo2ZyMEfuuXfEKT+mS86S9tP6k4YuYXOeZa7EnlqNJTUoSdRB6H2FBmUpT04ynwTWjtxr1gKn5LaytKS2bcfy9V6LTj0PO5Y/J7/kfi5vzEKPwQ1d+RYjxQvpZRmXIvaOcmRJvc7pp6VW/fg/9U5C75BbTo3nUYu1Indry63ZVFqOBcCRO/6XsII4Hau2rNS9WKhIaoWA1DsIue+UgSfnKVCIu8Uuyh5ItaWY3qnrsTKkV8tQ4uM9LXepd7x25J5DiQUpTYgl1qIrJxQmyIUD3HJLLBQ/f2xACT2r1MoOqv6U+5Q64xITVQ9py53q9UkbJ9Jr9UPpWxoEreqGArXcGyHndkqXlVvx4l9P6Zc6C4RSfq3r6bq6lBDEmFdyrTBH31Q6qvVYqmevkFjKWJG0gkPyl4DaQXwqrC25U60Zah5K+tyvMFFcwxRp1hB0DUIWTsp6lyy7BBL109JilAbXcu9J8HOziCne7fh/ruGxEWtJ7lRrpga+BZ0LPVCtgxxh1ugtFc5YpQ6QwxKegQPJnby5fpa7L43cu3S9X0r+ubeJtSR3a8ssd47sUCPJhR5qG3lNuIDa8HP1Rt26Hcs/R8xVr1KUkptUWbl2ygmdxcqIlZ+Tyfm5xFYDkBTWltxDaBGeSMmVig/7eUqtTQkXmurWhj6vi4Xc+hlzZFdy/EHuHhcUDzWXP7fogGqQUK6VyObg1OYb4ucfKbkP4FoL1Pzc69TyKPrMBX5HjJ2XXjN4SOeRQul7rymvxb6IlnrXGiOlz1raBkvS5uSMc2mSXpOSu4NaS7rkeg1ZSQ0QPTDq4G/N58oo8WSm3GXcapVFLwu7h9wa2aWWdyitZFiKq1cpF6Sg5J5Bi47pyqolqxahHTevNEYrpbQDtXwOTl1xBwtJvaWs6Fbvt1famOVdUzcSoZHQIDFF2FHJ3UFvFzpUZmmeFm54zeRoTGaL2KIkKINliTcg7XFIvAtK+KaVlySVttZyj5VXi9h+k5Zl+lByt+kO26LxcCFpnVM9EZ/cawc6N64YI42pyL6EpEtCStxyJb3GVPmxeyW/CMSpG057LZE5tfFAPVJDMhzkYvHkziFGysjKJdIeFparH0cW1zqUstwputWilVUtZX1R6z73rqjGSKys3GqTkkFM6jdSazykUFvP5ZXQ1/2c+0lB93oLb3bR5C5FjNyGMqaR+i1E7s7WnF6p71NCyjKVGoA5pFIyOOUGutwPuHAG6xTppyzG1uQuTbqh47Fz9SJhWITqMvVj4KEypfviosnd2jJ3rkaOn6fVYVgxAlgKOO+jtKPkOj1FJtXz46LU4uR6Frn0uc1G1HJCaWpJNTXYcT0aiXcXOiM+JbflL2VZuwbkToHEyN1SLtV1nyuoBFRjaUkRROwe1yrm3G+FUm/GlxGrb4nd0FzjK9dGQvpyB7xcmpw8qkfU2ihbS3IPVXqtax5DCwKeO5m7oIYOKB2SSgqt6ocjd9RRyuUvySPp9sf6TM2gkfru3xuXz7plUcJJEt4BZ55qRMsf2aFi7cidc3Z46leDuN8ps+FUK6LGlewxMHDrx71eawVLu/qlMnIrgziyptqnQCW02rqmyB/JvYS0a+shdOx27jO1XLXchUAZ8d20IULOjeSxRkc5kIjq9lMsppzlG0tTi7EjSlmbpTqU5pNylWtJ3ZclmY4qi2olS8jPeRjcdeOU8qnXKPMyJW1Hsr2FsDbkziFP9zqVbKnycjqGrpXoncsjTWTu5xqLtYdnMefyS9GCKFq0lZB8/1rp3AanXM7O1lS5tYNNS0Nrbcjd2nBFpizN0AtvPcMdQiuSLHEjQzJCHaWVK6+Io2WdShsBqeu1YZ7StFOFUVoNnlRyPwsriI2t7TO+79m9Cxtb27jzwWfeci+Ub8/uXTh80+XYs3sXAODYyU3cct9XcezkZjulAxjLbymXUi8xGW4djddC12vLUqTRqp1IyPbfufvfbwslZXHbVKiMVLkbW9tN2u3G1nayr3QBZQRo8VdquXNj2JR81qYnV+aCUkug1bO0dLEVq4PUPJX7v0Z+a8tacv5DY+6VO1RL0nNeYOkOWGocrxRzkrHETVaKOoQMLKldoq0Innqfo0fLfrFocuegdKcn5aVzJ2/mQIYxvUOfKbIUihwkjKAp2lqur0wFJXf7PWJvNVk6heUugdiANJcBSLFsrIIRNGIuerigkrvZSdsfBw4csEeOHGkmf5wkue3afbjswvOalbMEjJM//meFohVS7ezYyc1sn+W00zHtUtq2MeZJa+2BXLqVXC1DwThTvW7EXjLj76+IUShaI9bONra2cc+jx7Or3qirW8a0x05uns7D7SPSq2h6gUTuxpjrjDHPGWNOGGNuD9zfZYz54nD/CWPMXmlFS9CCqKgNqqX8WDpdiqhYdfjLB0NtmbPE0DXyDt90OQCw+ohkn+rdP7Pkbow5G8BnAHwQwH4AHzHG7PeS/SyAP7HW/giAXwbwKWlFpwBl7W4oT+kLpOZ1rRAXk6+rVSgE4BJ7rD+kyD+WNrVPIyZHsk/17p8Uy/1qACestS9aa78L4AsAbvTS3Ajg/uHzAwDeb4wxcmrSIO0++e4c5eXELA9qI8zJH93W267dx96woVCsEnL9oWaDXkxOjedcUnZLUMj9IgDfdr6/NFwLprHWvgngdQDnSyhIhbTL47tz7uhPyevqFLO0U3kpeikUS0fOkHKNnFj/ovS7wzddjte+81YOCfHKqoQ/KeQessD9JTaUNDDG3GqMOWKMObKxsUHRjwyOy8Pdylw62sYGiBJdQnopFOsMd/I1Rrg5Ih7vv/adsEccO4ajNLzSdUDIrZUEcA2Ah53vhwAc8tI8DOCa4fM5AF4FdpZZxv56bWLy0fMkPI4uc1xPq1DMHZQNRtRdqK37oBT3QPDgsK8B2GeMucQYcy6ADwN4yEvzEIBbhs8HAfyHQYnmoMa0Y4eGccqIlVPiprkxfIA3g69QKHZAWcZLCXXG0kn2ydlNqNqdGPrHsWOdHwXwW9baZ40xdxljbhiS/RqA840xJwB8AsBblku2ADWm7aYDeGENf7IlVg7npblycictKhSKdmi18i2GcTNVF1DM+xZ/kj/W4f6PgXrGTO5IgZqDkeZ2RoVCsc6gnB0zHl0i9fsFEqEZLPk8dxfUSc/UpOaI2Eidcv2oVrcvW610hWJahJYt+2fT3/Po8WhkwI8IlJTZEos9WyaG3PkSLc+fWMrZFgrFUuH30dy5NMdObuKeR493Dasu/mwZf8MB55yJVGx+zr96o1Ao2iLkmYeuj6BEBEb0XjCxkuS+sbWNTz7wND75wNOsn8mq/Yk4yrEACoVi+XC5gHKUwRQbn1aS3Pfs3oW7D16Buw9eQVptknoR7vVxWWIoX25FzlS/w6pQKPoiRdSpebveK+JWktyB7x0C5H4PgTtihpY95s5zAXbcs/t/5mo9FkChWDhSRJ271xNrMaEaiqXH4uvuBApw5jkxGjNXKBRTY/ETqhykToELpR3vh64rFArFKmClyb3m2M1U/Ct0f1VOglMoFApghcmdS7apjUmpHwNwv+sRAQqFYlWwsuTOWSGTGggkz1pXKBSKuWBlyR3grZCJ/bTWuApGoVAoloSVJvcYXKs+NDnqpzv/bbs0nq5QKBaFlV0KyVmaSEmrSx0VCsUqYNFLIamTqZwTGJXYFQrFkrCS5E5ZuaJLFxUKxTpjJckdoFnaunRRoVCsK1aW3FNITaIqFArFOmCR5K4bjhQKxbpjkeQO6ASpQqFYbyyW3BUKhWKdoeSuUCgUC4SSu0KhUCwQSu4KhUKxQCi5KxQKxQKh5K5QKBQLhJK7QqFQLBCTnQppjNkA8IeF2S8A8KqgOlJQvXiYq17AfHVTvXhYol7vstbuySWajNxrYIw5QjnysjdULx7mqhcwX91ULx7WWS8NyygUCsUCoeSuUCgUC8Sqkvtnp1YgAtWLh7nqBcxXN9WLh7XVayVj7gqFQqFIY1Utd4VCoVAksHLkboy5zhjznDHmhDHm9s5lv9MY8/vGmKPGmGeNMX9/uP6Lxpg/MsY8Nfxd7+Q5NOj6nDHmxxvq9i1jzDeH8o8M195hjHnEGHN8+P/24boxxvyLQa9vGGPe00inH3Xq5CljzKYx5rYp6ssYc58x5pQx5hnnGrt+jDG3DOmPG2NuaaTXPzHGHBvK/h1jzA8M1/caY/63U2+/4uS5anj/JwbdTQO92O9Nur9G9Pqio9O3jDFPDdd71leMG6ZrY9balfkDcDaAFwBcCuBcAE8D2N+x/AsBvGf4vBvA8wD2A/hFAL8QSL9/0HEXgEsG3c9upNu3AFzgXbsbwO3D59sBfGr4fD2ALwEwAN4L4IlO7+5/AHjXFPUF4H0A3gPgmdL6AfAOAC8O/98+fH57A70+AOCc4fOnHL32uuk8OV8FcM2g85cAfLCBXqz31qK/hvTy7v9TAP94gvqKccNkbWzVLPerAZyw1r5orf0ugC8AuLFX4dbak9barw+ftwAcBXBRIsuNAL5grd221v53ACew8wy9cCOA+4fP9wO4ybn+ebuDrwD4AWPMhY11eT+AF6y1qY1rzerLWvufAfxxoDxO/fw4gEestX9srf0TAI8AuE5aL2vtl621bw5fvwLgh1MyBt3Os9Y+bncY4vPOs4jplUDsvYn315Reg/X91wH8ZkpGo/qKccNkbWzVyP0iAN92vr+ENLk2gzFmL4ArATwxXPr44F7dN7pe6KuvBfBlY8yTxphbh2s/ZK09Cew0PgA/OIFeIz6MMzvd1PUF8Otninr7GexYeCMuMcb8gTHmPxljfmy4dtGgSw+9OO+td339GIBXrLXHnWvd68vjhsna2KqReygu1n25jzHmbQD+DYDbrLWbAP4lgD8L4C8AOIkd1xDoq+9fsta+B8AHAfy8MeZ9ibRd69EYcy6AGwD89nBpDvWVQkyP3vV2B4A3AfzGcOkkgIuttVcC+ASAf22MOa+jXtz31vt9fgRnGhDd6yvADdGkER3EdFs1cn8JwDud7z8M4OWeChhj/hR2Xt5vWGv/LQBYa1+x1v5fa+3/A/Cr+F4ooZu+1tqXh/+nAPzOoMMrY7hl+H+qt14DPgjg69baVwYdJ6+vAdz66abfMJH2EwD+xhA6wBD2eG34/CR24tl/btDLDd000avgvfWsr3MA/BSALzr6dq2vEDdgwja2auT+NQD7jDGXDNbghwE81KvwIab3awCOWmv/mXPdjVf/VQDjTP5DAD5sjNlljLkEwD7sTORI6/WnjTG7x8/YmZB7Zih/nG2/BcC/c/T66WHG/r0AXh9dx0Y4w6Kaur4ccOvnYQAfMMa8fQhJfGC4JgpjzHUA/iGAG6y1/8u5vscYc/bw+VLs1M+Lg25bxpj3Dm30p51nkdSL+9569tdrARyz1p4Ot/Ssrxg3YMo2VjNDPMUfdmaZn8fOKHxH57L/MnZcpG8AeGr4ux7AvwLwzeH6QwAudPLcMej6HCpn5BN6XYqdlQhPA3h2rBcA5wP4PQDHh//vGK4bAJ8Z9PomgAMN6+z7ALwG4Puda93rCzuDy0kA/wc71tHPltQPdmLgJ4a/v9VIrxPYibuObexXhrR/bXi/TwP4OoCfdOQcwA7ZvgDgXgwbFIX1Yr836f4a0mu4/usAfs5L27O+YtwwWRvTHaoKhUKxQKxaWEahUCgUBCi5KxQKxQKh5K5QKBQLhJK7QqFQLBBK7gqFQrFAKLkrFArFAqHkrlAoFAuEkrtCoVAsEP8fNN0FksLkbmcAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.scatter(range(vals.size(0))[:2000], vals_normed.data.tolist()[:2000], s = 0.1)\n",
    "# plt.scatter(range(vals.size(0)), vals_normed.data.tolist(), s = 0.1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 1000
    },
    "colab_type": "code",
    "id": "iseAD9gzzVmH",
    "outputId": "3aeb18ce-4fe8-41f3-971a-b47b829d4259"
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "D:\\python\\lib\\site-packages\\ipykernel_launcher.py:21: ComplexWarning: Casting complex values to real discards the imaginary part\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0\n",
      "1000\n"
     ]
    },
    {
     "ename": "KeyboardInterrupt",
     "evalue": "",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mKeyboardInterrupt\u001b[0m                         Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-10-7de65dcc18df>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m()\u001b[0m\n\u001b[0;32m     27\u001b[0m     \u001b[1;32mfor\u001b[0m \u001b[0mcomponent\u001b[0m \u001b[1;32min\u001b[0m \u001b[0mrange\u001b[0m \u001b[1;33m(\u001b[0m\u001b[0mnumber_gaussian_components\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     28\u001b[0m       \u001b[1;32mfor\u001b[0m \u001b[0mj\u001b[0m \u001b[1;32min\u001b[0m \u001b[0mrange\u001b[0m \u001b[1;33m(\u001b[0m\u001b[0mi\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mi\u001b[0m\u001b[1;33m+\u001b[0m\u001b[0mgmm_batch_size\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m---> 29\u001b[1;33m         \u001b[0mq_i\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mgmm\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mpredict_proba\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mgmm_data\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mj\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mreshape\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;33m-\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mreshape\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;36m4\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m     30\u001b[0m         \u001b[0mcomponent_weight\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mgmm\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mweights_\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mcomponent\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     31\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32mD:\\python\\lib\\site-packages\\sklearn\\mixture\\base.py\u001b[0m in \u001b[0;36mpredict_proba\u001b[1;34m(self, X)\u001b[0m\n\u001b[0;32m    357\u001b[0m         \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m_check_is_fitted\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    358\u001b[0m         \u001b[0mX\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0m_check_X\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mX\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;32mNone\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mmeans_\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mshape\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 359\u001b[1;33m         \u001b[0m_\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mlog_resp\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m_estimate_log_prob_resp\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mX\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m    360\u001b[0m         \u001b[1;32mreturn\u001b[0m \u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mexp\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mlog_resp\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    361\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32mD:\\python\\lib\\site-packages\\sklearn\\mixture\\base.py\u001b[0m in \u001b[0;36m_estimate_log_prob_resp\u001b[1;34m(self, X)\u001b[0m\n\u001b[0;32m    467\u001b[0m             \u001b[0mlogarithm\u001b[0m \u001b[0mof\u001b[0m \u001b[0mthe\u001b[0m \u001b[0mresponsibilities\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    468\u001b[0m         \"\"\"\n\u001b[1;32m--> 469\u001b[1;33m         \u001b[0mweighted_log_prob\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m_estimate_weighted_log_prob\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mX\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m    470\u001b[0m         \u001b[0mlog_prob_norm\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mlogsumexp\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mweighted_log_prob\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0maxis\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    471\u001b[0m         \u001b[1;32mwith\u001b[0m \u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0merrstate\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0munder\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;34m'ignore'\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32mD:\\python\\lib\\site-packages\\sklearn\\mixture\\base.py\u001b[0m in \u001b[0;36m_estimate_weighted_log_prob\u001b[1;34m(self, X)\u001b[0m\n\u001b[0;32m    420\u001b[0m         \u001b[0mweighted_log_prob\u001b[0m \u001b[1;33m:\u001b[0m \u001b[0marray\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mshape\u001b[0m \u001b[1;33m(\u001b[0m\u001b[0mn_samples\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mn_component\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    421\u001b[0m         \"\"\"\n\u001b[1;32m--> 422\u001b[1;33m         \u001b[1;32mreturn\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m_estimate_log_prob\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mX\u001b[0m\u001b[1;33m)\u001b[0m \u001b[1;33m+\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m_estimate_log_weights\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m    423\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    424\u001b[0m     \u001b[1;33m@\u001b[0m\u001b[0mabstractmethod\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32mD:\\python\\lib\\site-packages\\sklearn\\mixture\\gaussian_mixture.py\u001b[0m in \u001b[0;36m_estimate_log_prob\u001b[1;34m(self, X)\u001b[0m\n\u001b[0;32m    673\u001b[0m     \u001b[1;32mdef\u001b[0m \u001b[0m_estimate_log_prob\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mself\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mX\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    674\u001b[0m         return _estimate_log_gaussian_prob(\n\u001b[1;32m--> 675\u001b[1;33m             X, self.means_, self.precisions_cholesky_, self.covariance_type)\n\u001b[0m\u001b[0;32m    676\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    677\u001b[0m     \u001b[1;32mdef\u001b[0m \u001b[0m_estimate_log_weights\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mself\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32mD:\\python\\lib\\site-packages\\sklearn\\mixture\\gaussian_mixture.py\u001b[0m in \u001b[0;36m_estimate_log_gaussian_prob\u001b[1;34m(X, means, precisions_chol, covariance_type)\u001b[0m\n\u001b[0;32m    410\u001b[0m         \u001b[0mlog_prob\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mempty\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mn_samples\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mn_components\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    411\u001b[0m         \u001b[1;32mfor\u001b[0m \u001b[0mk\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;33m(\u001b[0m\u001b[0mmu\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mprec_chol\u001b[0m\u001b[1;33m)\u001b[0m \u001b[1;32min\u001b[0m \u001b[0menumerate\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mzip\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mmeans\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mprecisions_chol\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 412\u001b[1;33m             \u001b[0my\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mdot\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mX\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mprec_chol\u001b[0m\u001b[1;33m)\u001b[0m \u001b[1;33m-\u001b[0m \u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mdot\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mmu\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mprec_chol\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m    413\u001b[0m             \u001b[0mlog_prob\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mk\u001b[0m\u001b[1;33m]\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0msum\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0msquare\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0my\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0maxis\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    414\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;31mKeyboardInterrupt\u001b[0m: "
     ]
    }
   ],
   "source": [
    "\n",
    "gmm_range = (torch.Tensor(range(vals_normed.size(0))) - vals.size(0)/2)/vals.size(0)\n",
    "gmm_data = np.array([[float(vals_normed[i]), float(gmm_range[i])] for i in range (gmm_range.size(0))])\n",
    "\n",
    "gmm_batch_size = 800\n",
    "gmm_interval = 50\n",
    "number_gaussian_components = 4\n",
    "Fisher_vector = []\n",
    "with torch.no_grad():\n",
    "  for i in range (0, len(gmm_data) - gmm_batch_size, gmm_interval):\n",
    "\n",
    "    gmm = GaussianMixture(n_components = number_gaussian_components)\n",
    "    gmm.fit(gmm_data[i : i+gmm_batch_size])\n",
    "    u_k = np.zeros((4, 2))\n",
    "    v_k = np.zeros((4, 2))\n",
    "    cov_halfinvs = []\n",
    "\n",
    "    # Calculating Sigma^{-1/2}\n",
    "    for cov in gmm.covariances_:\n",
    "      eigenval, eigenvec = scipy.linalg.eig(cov)\n",
    "      inveigenval = [1/np.sqrt(float(i)) for i in eigenval]\n",
    "      halfinv = np.dot(np.diag(inveigenval), np.transpose(eigenvec))\n",
    "      cov_halfinvs.append(halfinv)\n",
    "    current_fisher = []\n",
    "\n",
    "    # Calculating u_k and v_k\n",
    "    for component in range (number_gaussian_components):\n",
    "      for j in range (i, i+gmm_batch_size):\n",
    "        q_i = gmm.predict_proba(gmm_data[j].reshape(1,-1)).reshape(4)\n",
    "        component_weight = gmm.weights_[component]\n",
    "\n",
    "        u_k_term2 =  np.dot(cov_halfinvs[component] , \n",
    "                            (gmm_data[j]-gmm.means_[component]).reshape(2))\n",
    "        u_k[component] +=  (q_i[component] / np.sqrt(gmm_batch_size * component_weight)) * u_k_term2\n",
    "\n",
    "\n",
    "        v_k_term2 = np.array([i**2 for i in u_k_term2]) - 1    \n",
    "        v_k[component] +=  q_i[component] * v_k_term2\n",
    "\n",
    "      current_fisher.extend(u_k[component].reshape(2))\n",
    "      current_fisher.extend(v_k[component].reshape(2))\n",
    "    Fisher_vector.append(current_fisher)\n",
    "    if (i%1000 == 0):\n",
    "      print (i)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "CRo6G7R4p167"
   },
   "outputs": [],
   "source": [
    "# The actual network\n",
    "class WaveNet(nn.Module):\n",
    "    def __init__(self, mu=1,n_residue=48, n_skip= 256, dilation_depth=10, \n",
    "                 n_repeat=5, nhpp = 50, nhpp_nbr = 0, additional_features = 0):\n",
    "        # n_residue: residue channels\n",
    "        # n_skip: skip channels\n",
    "        # dilation_depth & n_repeat: dilation layer setup\n",
    "        super(WaveNet, self).__init__()\n",
    "        n_skip = int(np.ceil(n_skip/(dilation_depth*n_repeat))*(dilation_depth*n_repeat))\n",
    "        self.nhpp_nbr = nhpp_nbr\n",
    "        self.dilation_depth = dilation_depth\n",
    "        dilations = self.dilations = [2**i for i in range(dilation_depth)] * n_repeat\n",
    "#         self.one_hot = One_Hot(mu)\n",
    "        self.from_input = nn.Conv1d(in_channels=mu, out_channels=n_residue, kernel_size=1)\n",
    "        self.conv_sigmoid = nn.ModuleList([nn.Conv1d(in_channels=n_residue, out_channels=n_residue, kernel_size=2, dilation=d)\n",
    "                         for d in dilations])\n",
    "        self.conv_tanh = nn.ModuleList([nn.Conv1d(in_channels=n_residue, out_channels=n_residue, kernel_size=2, dilation=d)\n",
    "                         for d in dilations])\n",
    "        self.skip_scale = nn.ModuleList([nn.Conv1d(in_channels=n_residue, out_channels=n_skip, kernel_size=1)\n",
    "                         for d in dilations]) \n",
    "        self.residue_scale = nn.ModuleList([nn.Conv1d(in_channels=n_residue, out_channels=n_residue, kernel_size=1)\n",
    "                         for d in dilations])\n",
    "        self.conv_post_1 = nn.Conv1d(in_channels=n_skip, out_channels=n_skip, kernel_size=1)\n",
    "        \n",
    "        if (additional_features > 0):      \n",
    "          self.linear1 = nn.Linear(in_features= n_skip + nhpp_nbr + additional_features, out_features=mu) # + 3, i.e. month, day, weekday\n",
    "        else:\n",
    "          self.linear1 = nn.Linear(in_features= n_skip + nhpp_nbr, out_features=mu)\n",
    "          \n",
    "          \n",
    "        if (nhpp_nbr + additional_features > 0): \n",
    "          self.linear_nhpp = nn.Linear(in_features = nhpp_nbr + additional_features, \n",
    "                                         out_features = nhpp_nbr + additional_features)\n",
    "          \n",
    "        self.dropout = nn.Dropout (p=0.5)\n",
    "        \n",
    "    def forward(self, input):\n",
    "        # Separating histories\n",
    "        history, nhpp = input\n",
    "        \n",
    "        \n",
    "        \n",
    "        wn_history = self.preprocess(history) # Preparing the history for the Wavenet sequence\n",
    "        ######### Wavenet Sequence ##########\n",
    "        \n",
    "        skip_connections = []\n",
    "        loop_cst = 0\n",
    "        # Convolution sequence\n",
    "        for s, t, skip_scale, residue_scale in zip(self.conv_sigmoid, self.conv_tanh, self.skip_scale, self.residue_scale):\n",
    "            wn_history, skip = self.residue_forward(wn_history, s, t, skip_scale, residue_scale)\n",
    "            \n",
    "            skip_connections.append(skip)\n",
    "        # sum up skip connections\n",
    "        wn_history = sum([s[:,:,-wn_history.size(2):] for s in skip_connections]) # If standard Wavenet\n",
    "        ######## End of Wavenet Sequence #########\n",
    "        # \n",
    "        output = self.postprocess(wn_history, nhpp)\n",
    "        \n",
    "        \n",
    "        return output\n",
    "    \n",
    "    def preprocess(self, input):\n",
    "        output = input.unsqueeze(0).unsqueeze(0).cuda()\n",
    "        output = self.from_input(output)\n",
    "        return output\n",
    "    \n",
    "    \n",
    "    def postprocess(self, wn_history, nhpp):\n",
    "        \n",
    "        \n",
    "            \n",
    "        # Discard first nhpp (+1) points as we don't need them - they don't have nhpp history points\n",
    "        wn_history = wn_history[:,:,-wn_history.size(2)+self.nhpp_nbr:]\n",
    "        \n",
    "        start_point = nhpp.size(1) - wn_history.size(2)\n",
    "        \n",
    "        nhpp = nhpp[:,start_point:].unsqueeze(0)\n",
    "  \n",
    "  \n",
    "        # One layer to process the nhp\n",
    "        if (add_nhpp_FF == True):\n",
    "          if (nhpp.size(1) + additional_features > 0):\n",
    "            nhpp = torch.tanh(self.linear_nhpp(nhpp.transpose(1,2)).transpose(1,2))\n",
    "        \n",
    "        output = torch.cat((wn_history, nhpp), dim = 1)\n",
    "        \n",
    "        \n",
    "        output = output.transpose(1,2)\n",
    "\n",
    "#         Dropout\n",
    "        if (add_dropout == True):\n",
    "          if (training_so_dropout == True):\n",
    "            output = self.dropout(output)\n",
    "      \n",
    "        output = self.linear1(output).squeeze(0).squeeze(1) # Last layer: linear - for regression\n",
    "        \n",
    "        \n",
    "        return output\n",
    "    \n",
    "    def residue_forward(self, input, conv_sigmoid, conv_tanh, skip_scale, residue_scale):\n",
    "        output = input\n",
    "        output_sigmoid, output_tanh = conv_sigmoid(output), conv_tanh(output)\n",
    "        output = torch.sigmoid(output_sigmoid) * torch.tanh(output_tanh)\n",
    "        skip = skip_scale(output)\n",
    "        output = residue_scale(output)\n",
    "        output = output + input[:,:,-output.size(2):]\n",
    "        return output, skip\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "xG1lqPHDGqut"
   },
   "outputs": [],
   "source": [
    "# Function which runs the algorithm\n",
    "\n",
    "\n",
    "# classes = 128\n",
    "# values = periodicise(classes, vals)\n",
    "# vals = ground truth, all dataset\n",
    "def run_algorithm_gpu(nhpp, n_residue, n_skip, dilation_depth, n_repeat,\n",
    "                      max_epoch, batch_size, learn_rate, graph_cst, batches, additional_features = 0,\n",
    "                      dbd_prediction = False, add_nhpp_FF = True, add_dropout = True):\n",
    "    print ('nhpp', nhpp[nhpp_nbr].size(1))\n",
    "    print ('additional features', additional_features)\n",
    "    print ('n_residue', n_residue)\n",
    "    print ('n_skip', n_skip)\n",
    "    print ('dilation_depth', dilation_depth)\n",
    "    print ('n_repeat', n_repeat)\n",
    "    print ('max_epoch', max_epoch)\n",
    "    print ('batch_size', batch_size)\n",
    "    print ('learning rate', learn_rate)\n",
    "    print ('add_nhpp_FF', add_nhpp_FF)\n",
    "    print ('add_dropout', add_dropout)\n",
    "    best_valid_epoch = 0\n",
    "    net = WaveNet(mu=1,n_residue=n_residue,n_skip=n_skip,dilation_depth=dilation_depth,\n",
    "                  n_repeat=n_repeat, nhpp = nhpp, nhpp_nbr = nhpp_nbr,\n",
    "                 additional_features = additional_features).cuda()\n",
    "    optimizer = optim.Adam(net.parameters(),lr=learn_rate)\n",
    "#     optimizer = optim.Adam(net.parameters(),lr=learn_rate)\n",
    "#     opti mizer = optim.SGD(net.parameters(), lr=learn_rate, momentum=0.9)\n",
    "    max_data = len(values)-batch_size\n",
    "    \n",
    "    \n",
    "    loss_fct = RMSELoss()\n",
    "    temp_loss = torch.Tensor(0)\n",
    "\n",
    "    loss_save = [] \n",
    "    valid_loss_save = []\n",
    "    valid_loss = 1000000\n",
    "    last_loss = 1000000\n",
    "    \n",
    "    sequence_lowering_loss = 0\n",
    "    \n",
    "    for epoch in range(max_epoch):\n",
    "        optimizer.zero_grad()\n",
    "        loss = 0\n",
    "        training_so_dropout = True\n",
    "        for j in range(int(batches)):\n",
    "            i = random.randint(0, len(values) - batch_size)\n",
    "            batch = torch.FloatTensor(values[i : i + batch_size]).cuda()\n",
    "\n",
    "            nhpp_sample = nhpp[nhpp_nbr][:,:, i : i + batch_size].cuda()\n",
    "#             print (nhpp_sample.size())\n",
    "            x = (batch[:-1], nhpp_sample[0,:,:-1])\n",
    "            logits = net(x)\n",
    "            sz = logits.size(0)\n",
    "            temp_loss = loss_fct(logits, batch[-sz:])\n",
    "            \n",
    "            loss = loss + temp_loss\n",
    "         \n",
    "        loss = loss/batches\n",
    "        loss.backward()\n",
    "        optimizer.step()\n",
    "        loss_save.append(float(loss))\n",
    "        if (epoch % graph_cst == 0):\n",
    "          training_so_dropout = False\n",
    "      \n",
    "          if (dbd_prediction == False):\n",
    "#               print('-------test set, epoch {}--------'.format(epoch+1))\n",
    "              batch = torch.FloatTensor(values_valid).cuda()\n",
    "              nhpp_sample = nhpp_valid[nhpp_nbr]\n",
    "\n",
    "              x = (batch[:-1], nhpp_sample[0,:,:-1])\n",
    "\n",
    "              logits = net(x)\n",
    "              l = logits \n",
    "              sz = logits.size(0)\n",
    "              v_loss = loss_fct(logits, batch[-sz:])\n",
    "#             t_loss = test_generation(nt = net, nhpp = nhpp[nhpp_nbr], plot=False, prediction_window = 5) # For maenners dataset\n",
    "          else:    \n",
    "            diff, t_loss = test_generation(nt = net, nhpp = nhpp_test[nhpp_nbr], plot=False, prediction_window = 30)\n",
    "            \n",
    "          if (v_loss < valid_loss):\n",
    "            best_valid_epoch = epoch\n",
    "            valid_loss = v_loss\n",
    "            best_net = WaveNet(mu=1,n_residue=n_residue,n_skip=n_skip,dilation_depth=dilation_depth,\n",
    "                                n_repeat=n_repeat, nhpp = nhpp, nhpp_nbr = nhpp_nbr,\n",
    "                                additional_features = additional_features).cuda()\n",
    "            best_net.load_state_dict(copy.deepcopy(net.state_dict()))\n",
    "            \n",
    "            valid_loss_save.append(float(valid_loss))\n",
    "            \n",
    "            \n",
    "#           if(t_loss > last_loss):\n",
    "#             sequence_lowering_loss += 1\n",
    "#           else:\n",
    "#             sequence_lowering_loss = 0\n",
    "# #           print (sequence_lowering_loss, last_loss)\n",
    "#           last_loss = t_loss\n",
    "#           if(sequence_lowering_loss > 20):\n",
    "#               print ('epoch finished', epoch)\n",
    "#               break\n",
    "            \n",
    "            \n",
    "    \n",
    "            \n",
    "#             print (t_loss, test_loss)\n",
    "     \n",
    "    return best_net, loss_save, valid_loss_save, best_valid_epoch"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "sKH1irIjrLda"
   },
   "outputs": [],
   "source": [
    "\n",
    "# Function used to create the near history past points (nhpp for short)\n",
    "def past_points_collector(values, gmm_batch_size):\n",
    "  # Precomputing the allocation of each GMM's to each point and additional features\n",
    "  # gmm_batch_size - starting point of allocation\n",
    "  # every gmm_interval allocate the new batch of GMM's\n",
    "  \n",
    "\n",
    "  \n",
    "  nhpp = torch.zeros([1,16,gmm_batch_size]).cuda()\n",
    "\n",
    "  for i in range (gmm_batch_size + 1, len(values)):\n",
    "\n",
    "    nhpp_temp = []\n",
    "\n",
    "    nhpp_temp.extend(Fisher_vector[int(np.floor(i/gmm_interval))])\n",
    "    nhpp_temp = torch.FloatTensor(nhpp_temp).unsqueeze(1).unsqueeze(0).cuda()\n",
    "    nhpp = torch.cat((nhpp, nhpp_temp), dim = 2)\n",
    "\n",
    "      \n",
    "    \n",
    "    \n",
    "    \n",
    "  return nhpp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 1000
    },
    "colab_type": "code",
    "id": "cjocfKKBf8sJ",
    "outputId": "c3a0ae37-febd-4724-b5ee-81db97fa7e5e"
   },
   "outputs": [],
   "source": [
    "# Applying algorithm to the data, takes a long ass time\n",
    "\n",
    "additional_features = 0\n",
    "prediction_train_size = 100\n",
    "\n",
    "import time\n",
    "\n",
    "\n",
    "n_skips = [500]\n",
    "# n_skips = [300, 700]\n",
    "\n",
    "\n",
    "# learn_rates = [1e-3, 3e-3, 5e-3]\n",
    "learn_rates = [1e-3]\n",
    "\n",
    "# nhpp_nbrs = [0, 14, 28, 42, 56, 70, 84, 98] # tanke & tanke_ts & baerchen\n",
    "nhpp_nbrs = [len(Fisher_vector[0])] # maenners\n",
    "\n",
    "# nhpp_nbrs = [0, 42]\n",
    "\n",
    "\n",
    "# dilation_depths = [6, 7] # tanke & tanke_ts & baerchen\n",
    "dilation_depths = [12] # maenners\n",
    "# dilation_depths = [6]\n",
    "\n",
    "training_so_dropout = True\n",
    "add_nhpp_FF = False\n",
    "add_dropout = False\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "# Constants\n",
    "\n",
    "number_of_experiments = 15 #constant saying how many experiments should be run for one setting of hyperparameters\n",
    "# it's used to calculate the mean standard deviation of the test loss in the case we want to make a comparison with \n",
    "# other models\n",
    "\n",
    "\n",
    "\n",
    "max_epoch = 1000\n",
    "batches = 5\n",
    "graph_cst = math.floor(max_epoch/50)\n",
    "print ('graph cst', graph_cst)\n",
    "n_residue = 25\n",
    "n_repeat = 2\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "calc_length = len(nhpp_nbrs)*len(dilation_depths)*len(n_skips)*len(learn_rates) * number_of_experiments\n",
    "calc_count = 0\n",
    "dbd_prediction = False\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "net = defaultdict(lambda: defaultdict(lambda: defaultdict(lambda: defaultdict(int))))\n",
    "nhpp = defaultdict(int)\n",
    "nhpp_test = defaultdict(int)\n",
    "nhpp_valid = defaultdict(int)\n",
    "losses = defaultdict(lambda: defaultdict(lambda: defaultdict(lambda: defaultdict(lambda: defaultdict(\n",
    "                                lambda: defaultdict(lambda: defaultdict(lambda: defaultdict(int))))))))\n",
    "\n",
    "\n",
    "print (len(values))\n",
    "\n",
    "\n",
    "for nhpp_nbr in nhpp_nbrs:\n",
    "  for learn_rate in learn_rates:\n",
    "    nhpp[nhpp_nbr] = past_points_collector(values = values, gmm_batch_size = nhpp_nbr)\n",
    "    nhpp_test[nhpp_nbr] = past_points_collector(values_test, gmm_batch_size = nhpp_nbr)\n",
    "    nhpp_valid[nhpp_nbr] = past_points_collector(values_valid, gmm_batch_size = nhpp_nbr)\n",
    "\n",
    "    for n in range(len(n_skips)):\n",
    "      \n",
    "      n_skip = n_skips[n] - nhpp_nbr - additional_features\n",
    "\n",
    "        #   for n_skip in n_skips:\n",
    "\n",
    "      for dilation_depth in dilation_depths:\n",
    "            batch_size = 2**(dilation_depth+1) + nhpp_nbr + prediction_train_size\n",
    "            total_train_loss = torch.zeros(max_epoch)\n",
    "            total_valid_loss = []\n",
    "            best_valid_train_loss = []\n",
    "            total_test_loss = []\n",
    "    #     for nhpp_nbr in nhpps:\n",
    "\n",
    "\n",
    "            for repeat in range (number_of_experiments):\n",
    "\n",
    "              \n",
    "              start = time.time()\n",
    "              calc_count += 1\n",
    "\n",
    "\n",
    "              net[nhpp_nbr][n_skip][dilation_depth][learn_rate], loss_save, valid_loss, best_valid_epoch = run_algorithm_gpu(nhpp = nhpp, \n",
    "                                                   n_residue = n_residue, \n",
    "                                                   n_skip = n_skip, \n",
    "                                                   dilation_depth = dilation_depth,  \n",
    "                                                   n_repeat = n_repeat,\n",
    "                                                   batch_size = batch_size, \n",
    "                                                   learn_rate = learn_rate, \n",
    "                                                   max_epoch = max_epoch,\n",
    "                                                   graph_cst = graph_cst,\n",
    "                                                   additional_features = additional_features,\n",
    "                                                   dbd_prediction = dbd_prediction,\n",
    "                                                   batches = batches,\n",
    "                                                   add_nhpp_FF = add_nhpp_FF,\n",
    "                                                   add_dropout = add_dropout)\n",
    "          \n",
    "              print ('best validation epoch', best_valid_epoch)\n",
    "              \n",
    "              \n",
    "              # Calculating test loss\n",
    "              batch = torch.FloatTensor(values_test).cuda()\n",
    "              nhpp_sample = nhpp_test[nhpp_nbr]\n",
    "              x = (batch[:-1], nhpp_sample[0,:,:-1])\n",
    "              logits = net[nhpp_nbr][n_skip][dilation_depth][learn_rate](x)\n",
    "              sz = logits.size(0)\n",
    "              t_loss = loss_fct(logits, batch[-sz:])\n",
    "              \n",
    "              \n",
    "              \n",
    "          \n",
    "              total_train_loss += torch.Tensor(loss_save)\n",
    "              total_valid_loss.append(min(valid_loss))\n",
    "              best_valid_train_loss.append(loss_save[best_valid_epoch])\n",
    "              total_test_loss.append(float(t_loss))\n",
    "              \n",
    "              \n",
    "              \n",
    "              end = time.time()\n",
    "              print ('---------------', 'done', (calc_count/calc_length) * 100, \n",
    "                    '% -----------------', 'time left', (end-start) * (calc_length-calc_count))\n",
    "            \n",
    "            \n",
    "            total_train_loss = total_train_loss / number_of_experiments\n",
    "            total_valid_loss = np.mean(total_valid_loss) \n",
    "            best_valid_train_loss = np.mean(best_valid_train_loss) \n",
    "#             total_test_loss = total_test_loss\n",
    "            \n",
    "            \n",
    "            losses[nhpp_nbr][n_residue][n_skip][dilation_depth][n_repeat][prediction_train_size][learn_rate][max_epoch] = [total_train_loss, \n",
    "                                                                                                                np.mean(total_valid_loss),\n",
    "                                                                                                                np.mean(best_valid_train_loss),\n",
    "                                                                                                                np.mean(total_test_loss),\n",
    "                                                                                                                np.std(total_test_loss)]\n",
    "            \n",
    "              #         test_generation(first_points = 0, last_points = 60, prediction_window = 30)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "tzg822XwpLoz"
   },
   "outputs": [],
   "source": [
    "class l_class:\n",
    "  def __init__(self):\n",
    "\n",
    "        self.nhpp_nbr = 0\n",
    "\n",
    "        self.n_residue  = 0\n",
    "\n",
    "        self.n_skip  = 0\n",
    "\n",
    "        self.n_repeat   = 0\n",
    "        \n",
    "        self.prediction_train_size = 0\n",
    "        \n",
    "        self.dilation_depth = 0\n",
    "        \n",
    "        self.learn_rate = 0\n",
    "        \n",
    "        self.max_epoch = 0\n",
    "        \n",
    "        self.loss_train = []\n",
    "        self.loss_train_zero = []\n",
    "        \n",
    "        self.loss_valid = 0\n",
    "        self.loss_valid_zero = 0\n",
    "        \n",
    "        self.loss_test = 0\n",
    "        self.loss_test_zero = 0\n",
    "        \n",
    "        self.best_valid_train_loss = 0\n",
    "        self.best_valid_train_loss_zero = 0\n",
    "        \n",
    "        self.test_std = 0\n",
    "        \n",
    "        self.nbr_of_experiments = 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 262
    },
    "colab_type": "code",
    "id": "ViWVOynzrkJW",
    "outputId": "df169822-248c-4912-afa8-4e9c13a6dbf5"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "test_los and test_std_zero 0.34683617353439333 0.0012378922087593327\n"
     ]
    },
    {
     "ename": "NameError",
     "evalue": "ignored",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-28-7bad6e0f1e1c>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[1;32m     84\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     85\u001b[0m                     \u001b[0mobject_to_write\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mloss_train\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mlos_train\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 86\u001b[0;31m                     \u001b[0mobject_to_write\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mloss_train_zero\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mlos_train_zero\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     87\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     88\u001b[0m                     \u001b[0mobject_to_write\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mloss_valid\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mlos_valid\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mNameError\u001b[0m: name 'los_train_zero' is not defined"
     ]
    }
   ],
   "source": [
    "# Cell which saves the results for later analysis\n",
    "# Showing the results and deciding which ones are good and which ones are not\n",
    "# write_file_title = \"\"\n",
    "if (make_dataset == True):\n",
    "  write_file_title = \"Fisher_\" + \"noise_\"+str(noise)+\"_results.txt\"\n",
    "\n",
    "else:\n",
    "  write_file_title = \"Autoencoder_Wavenet_\"\n",
    "\n",
    "  if (add_nhpp_FF == True):\n",
    "    write_file_title = write_file_title + \"FF_\"\n",
    "  if (add_dropout == True):\n",
    "    write_file_title = write_file_title + \"dropout_\"\n",
    "  write_file_title = write_file_title + \"add_feat = \" + str(additional_features) + \",\" + file + \"_results.txt\"\n",
    "  \n",
    "test_stds = []\n",
    "\n",
    "f = open(write_file_title,\"wb\")\n",
    "object_list = []\n",
    "\n",
    "loss_fct = RMSELoss()\n",
    "min_los = 10000\n",
    "min_diff = 10000\n",
    "los_list = []\n",
    "for nhpp_nbr in losses.keys():\n",
    "  for n_residue in losses[nhpp_nbr].keys():\n",
    "    for n_skip in losses[nhpp_nbr][n_residue].keys():\n",
    "      for dilation_depth in losses[nhpp_nbr][n_residue][n_skip].keys():\n",
    "        for n_repeat in losses[nhpp_nbr][n_residue][n_skip][dilation_depth].keys():\n",
    "          for prediction_train_size in losses[nhpp_nbr][n_residue][n_skip][dilation_depth][n_repeat].keys():\n",
    "            for learn_rate in losses[nhpp_nbr][n_residue][n_skip][dilation_depth][n_repeat][prediction_train_size].keys():\n",
    "              for max_epoch in losses[nhpp_nbr][n_residue][n_skip][dilation_depth][n_repeat][prediction_train_size][learn_rate].keys():\n",
    "                  if( type(losses[nhpp_nbr][n_residue][n_skip][dilation_depth][n_repeat][prediction_train_size][learn_rate][max_epoch]) != int):\n",
    "                    \n",
    "                    # current results\n",
    "                    los_train = losses[nhpp_nbr][n_residue][n_skip][dilation_depth][n_repeat][prediction_train_size][learn_rate][max_epoch][0]\n",
    "                    los_valid = losses[nhpp_nbr][n_residue][n_skip][dilation_depth][n_repeat][prediction_train_size][learn_rate][max_epoch][1]\n",
    "                    best_valid_train_loss = losses[nhpp_nbr][n_residue][n_skip][dilation_depth][n_repeat][prediction_train_size][learn_rate][max_epoch][2]\n",
    "                    test_los = losses[nhpp_nbr][n_residue][n_skip][dilation_depth][n_repeat][prediction_train_size][learn_rate][max_epoch][3]\n",
    "                    test_std = losses[nhpp_nbr][n_residue][n_skip][dilation_depth][n_repeat][prediction_train_size][learn_rate][max_epoch][4]\n",
    "\n",
    "                    print ('test loss and test std', test_los, test_std)\n",
    "                  \n",
    "\n",
    "                    \n",
    "  \n",
    "                    # Writing to file\n",
    "                    object_to_write = l_class()\n",
    "                    \n",
    "                    object_to_write.nhpp_nbr      = nhpp_nbr\n",
    "\n",
    "                    object_to_write.n_residue       = n_residue\n",
    "\n",
    "                    object_to_write.n_skip      = n_skip\n",
    "\n",
    "                    object_to_write.n_repeat   = n_repeat\n",
    "\n",
    "                    object_to_write.prediction_train_size = n_repeat\n",
    "\n",
    "                    object_to_write.learn_rate = learn_rate\n",
    "                    \n",
    "                    object_to_write.dilation_depth = dilation_depth\n",
    "\n",
    "                    object_to_write.max_epoch = max_epoch\n",
    "                    \n",
    "                    object_to_write.loss_train = los_train\n",
    "                    \n",
    "                    object_to_write.loss_valid = los_valid\n",
    "                    \n",
    "                    object_to_write.loss_test = test_los\n",
    "\n",
    "                    object_to_write.best_valid_train_loss = best_valid_train_loss\n",
    "                    \n",
    "                    object_to_write.test_std = test_std\n",
    "                    \n",
    "                    object_to_write.nbr_of_experiments = number_of_experiments\n",
    "\n",
    "                    object_to_write.total_time_duration = total_duration\n",
    "                    \n",
    "                    object_list.append(object_to_write)\n",
    "                    \n",
    "                    \n",
    "pickle.dump(object_list, f)\n",
    "f.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "1gsBPr7QgGB_"
   },
   "outputs": [],
   "source": [
    "losses"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "wiwtgp5nJYRO"
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "accelerator": "GPU",
  "colab": {
   "collapsed_sections": [],
   "name": "Fisher_Vector_full_covariance_Wavenet_first_datasets(1).ipynb",
   "provenance": []
  },
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}