Whole code.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "/mnt/Data/Università_Magistrale/Sullivan/Hadoop/Hadoop_group_mobile_improved\n"
     ]
    }
   ],
   "source": [
    "cd /mnt/Data/Università_Magistrale/Sullivan/Hadoop/Hadoop_group_mobile_improved/\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from sklearn import preprocessing \n",
    "import pandas as pd\n",
    "import datetime as dt\n",
    "from sklearn.metrics import mean_absolute_error,mean_squared_error,median_absolute_error\n",
    "from sklearn.svm import SVR\n",
    "from sklearn.neighbors import KNeighborsRegressor\n",
    "from sklearn.ensemble import GradientBoostingRegressor\n",
    "from sklearn.decomposition import PCA as sklearnPCA ,KernelPCA,FactorAnalysis,IncrementalPCA,FastICA\n",
    "from sklearn.manifold import Isomap,LocallyLinearEmbedding\n",
    "from scipy.stats.mstats import normaltest \n",
    "from scipy.stats import spearmanr\n",
    "from math import *\n",
    "import matplotlib.pyplot as plt\n",
    "import time\n",
    "from sklearn.grid_search import GridSearchCV\n",
    "from sklearn.cluster import MeanShift, estimate_bandwidth\n",
    "import matplotlib.pyplot as plt\n",
    "import urllib2\n",
    "import json"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#Classes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "class Sensor(object):\n",
    "    ID=None\n",
    "    owner=None\n",
    "    days={} #days is a dictionary containing a dataframe with the safecast data for that specific day\n",
    "    daysList=[] #it contains the days of the measurement, it is a list of the keys of days dictionary\n",
    "    dataset=None\n",
    "    latitude=None\n",
    "    longitude=None\n",
    "    stationary=None\n",
    "    def __init__(self,a,date='Captured Time'):\n",
    "        #given a series of measurement it creates a dataframe for every day\n",
    "        df=pd.DataFrame(a)\n",
    "        df=df.sort('Captured Time')\n",
    "        self.latitude,self.longitude,self.ID=df[['Latitude','Longitude','Sensor']].iloc[0].values\n",
    "        i=lambda x: str(x.year) + '-' + str(x.month) + '-' +str(x.day) #I take just year,month and day \n",
    "        try:\n",
    "             dates= df[date].apply(i) \n",
    "        except AttributeError:\n",
    "            df=df.convert_objects(convert_dates='coerce')\n",
    "            dates= df[date].apply(i)\n",
    "        df['Date']=dates\n",
    "        daysList=dates.unique()\n",
    "        self.stationary=Sensor.isStationary(df)\n",
    "        self.days=dict([(day,df[df['Date']==day]) for day in daysList])\n",
    "        self.daysList=daysList\n",
    "    \n",
    "    def apply(self,f):\n",
    "        '''Apply a generic function on historical data'''\n",
    "        self.days.update((x, f(y)) for x, y in self.days.items())\n",
    "        return self\n",
    "    \n",
    "    def addDay(self,a,date='Captured Time'): \n",
    "        ''' It adds another day to the days dictionary\n",
    "        \n",
    "        '''\n",
    "        df=pd.DataFrame(a)\n",
    "        i=lambda x: str(x.year) + '-' + str(x.month) + '-' +str(x.day) #I take just year,month and day \n",
    "        try:\n",
    "             dates= df[date].apply(i) \n",
    "        except AttributeError:\n",
    "            df=df.convert_objects(convert_dates='coerce')\n",
    "            dates= df[date].apply(i)\n",
    "        df['Day']=dates\n",
    "        daysList=dates.unique()\n",
    "        [self.days.update({day:df[df['Day']==day]}) for day in daysList] \n",
    "        [self.daysList.append(day) for day in daysList]\n",
    "        return self\n",
    "    \n",
    "    def cleanAll(self):\n",
    "        '''It cleans all the measurements applying the static method clean to every day\n",
    "        '''\n",
    "        self.days.update((x, Sensor.clean(y)) for x, y in self.days.items())\n",
    "        return self\n",
    "    \n",
    "    @staticmethod\n",
    "    def clean(df):\n",
    "        '''It cleans a single day\n",
    "        '''\n",
    "        from string import strip\n",
    "        columns=['Captured Time','Latitude','Longitude','Value','Unit','ID','Height','Loader ID','Sensor','Distance']\n",
    "        df=df[columns]\n",
    "        #df=df.dropna(1) #empty rows are deleted \n",
    "        df=df.drop_duplicates('Captured Time') #sometimes there are some duplicates\n",
    "        df.index=xrange(0,len(df))\n",
    "        today=dt.datetime.today().strftime(\"%Y-%m-%d %H:%M:%S\")\n",
    "        df=df.convert_objects(convert_dates='coerce')\n",
    "        df=df[df['Captured Time']<=today] #every row with date field incorrect is deleted\n",
    "        df['Unit']=df['Unit'].apply(strip)\n",
    "        df=df[df.Unit=='cpm'] #all the units that are not expressed in cpm are deleted\n",
    "        #I should add some lines to remove special character like \\n and \\t\n",
    "        return df\n",
    "    \n",
    "    @staticmethod\n",
    "    def convertDate(df,date='Captured Time'):\n",
    "        df[date]=0\n",
    "        try:\n",
    "            f = lambda x: str(int(x.Year)) + '-'+ str(int(x.Month)) + '-' + str(int(x.Day)) + ' ' + str(int(x.Hour)) + ':' + str(int(x.Minute)) + ':' + '00'\n",
    "            df[date]=df.apply(f,1)      \n",
    "        except AttributeError:  \n",
    "            diz={0:'00',0.25:'15',0.5:'30',0.75:'45'}\n",
    "            g = lambda x: str(int(x.Year)) + '-'+ str(int(x.Month)) + '-' + str(int(x.Day)) + ' ' + str(int(x.Hour)) + ':' + diz[x.Hour - int(x.Hour)] + ':' + '00'\n",
    "            df[date]=df.apply(g,1)                                                                                                                     \n",
    "        df=df.drop(['Year','Month','Day','Hour'],axis=1)\n",
    "        fmt=\"%Y-%m-%d %H:%M:%S\"\n",
    "        try:\n",
    "            df[date]=df[date].apply(dt.datetime.strptime(date,fmt))\n",
    "        except ValueError:\n",
    "            pass\n",
    "        return df                                                                                                                \n",
    "    \n",
    "    def createDataset(self):\n",
    "        '''It merge all the dataframe in the days dictionary in a single dataframe\n",
    "        '''\n",
    "        tmp=self.days.values()\n",
    "        df = pd.concat(tmp)\n",
    "        self.dataset=df#.sort('Captured Time')\n",
    "        return self.dataset\n",
    "   \n",
    "    def delDay(self,day):\n",
    "        try:\n",
    "            self.days.pop(day)\n",
    "            self.daysList.remove(day)\n",
    "        except KeyError:\n",
    "            print 'The day ' + str(day) + ' is not present'\n",
    "            return self\n",
    "        return self\n",
    "    \n",
    "    @staticmethod\n",
    "    def distance(a1,b1,a2,b2):\n",
    "        '''Evaluates the distance in m between two points with coordinates expressed in\n",
    "        Latitude and Longitude \n",
    "        '''\n",
    "        a1=a1*np.pi/180\n",
    "        a2=a2*np.pi/180\n",
    "        b1=b1*np.pi/180\n",
    "        b2=b2*np.pi/180\n",
    "        return np.arccos(np.cos(a1-a2)*np.cos(b1)*np.cos(b2)+np.sin(b1)*np.sin(b2))*6378*1000    \n",
    "    \n",
    "    def extractDates(self,date='Captured Time',delta=0.25):\n",
    "        '''It applies the extracDate static method on every day \n",
    "        '''\n",
    "        self.days.update((x, Sensor.extractDate(y,date,delta)) for x, y in self.days.items())\n",
    "        return self\n",
    "    \n",
    "    @staticmethod\n",
    "    def extractDate(df,date='Captured Time',delta=0.25):\n",
    "        '''Add two different fields useful to couple with weather data.\n",
    "        \n",
    "        The field 'DAY': year-month-day and the field 'Hour': hour.minutes\n",
    "        \n",
    "        '''\n",
    "        import datetime as dt\n",
    "        fmt=\"%Y-%m-%d\"\n",
    "        i=lambda x: str(x.year) + '-' + str(x.month) + '-' +str(x.day) #I take just year,month and day \n",
    "        try:\n",
    "             dates= df[date].apply(i) \n",
    "        except AttributeError:\n",
    "            df=df.convert_objects(convert_dates='coerce')\n",
    "            dates= df[date].apply(i)\n",
    "        g = lambda x: dt.datetime.strptime(x,fmt)\n",
    "        dates= dates.apply(g)\n",
    "        h=lambda x : str(x).split(' ')[0]#the conversion adds hour,minutes and seconds \n",
    "        dates= dates.apply(h) #I drop it and return a list of string\n",
    "        df['Year']=df[date].apply(lambda x : x.year)\n",
    "        df['Month']=df[date].apply(lambda x: x.month)\n",
    "        df['Day']=df[date].apply(lambda x: x.day)  \n",
    "        tmp=df[date].apply(lambda x: x.to_datetime())\n",
    "        df['Hour']=tmp.apply(lambda x: x.hour)\n",
    "        tmp=df[date].apply(lambda x: x.minute)\n",
    "        f=lambda x: round(round(x/(60*delta))*delta,3)\n",
    "        \n",
    "        df['Hour']=df['Hour']+tmp.apply(f)\n",
    "        df['Hour']=df['Hour'].replace(24,0.00)\n",
    "        \n",
    "        return df\n",
    "    \n",
    "    def getDays(self):\n",
    "        print self.daysList\n",
    "        \n",
    "    @staticmethod\n",
    "    def isStationary(df):\n",
    "        '''It returns True if the measurement in df belong to a stationary detector\n",
    "        '''\n",
    "        l1=df.Latitude.iloc[0]\n",
    "        l2=df.Longitude.iloc[0]\n",
    "        m1=df.Latitude.iloc[len(df)-1]\n",
    "        m2=df.Longitude.iloc[len(df)-1]\n",
    "        if df.Distance.max()>15: #it checks if the distance between two consevutive measurements is more than\n",
    "            #the maximum value of gps spatial inaccuracy\n",
    "            return False\n",
    "        if Sensor.distance(l1,l2,m1,m2)>100: #it checks if the distance between the first and the last point \n",
    "                                                #is too much\n",
    "            return False\n",
    "        if df.Distance.sum()>2*len(df):\n",
    "            return False\n",
    "        return True\n",
    "\n",
    "    def timeSampling(self,day):\n",
    "        '''It returns the time sampling of the measurement in the day indicated\n",
    "        '''\n",
    "        from numpy import median\n",
    "        df=self.days[day]\n",
    "        df=df.clean()\n",
    "        return median([(df['Captured Time'].loc[n]-df['Captured Time'].loc[m]).total_seconds() for n,m in zip(xrange(1,len(df)),xrange(0,(len(df)-1)))])\n",
    "    def to_csv(self,filename):\n",
    "        with open(filename, 'a') as f:\n",
    "               self.dataset.to_csv(f,index=False,float_format = '%.4f',header=False)\n",
    "\n",
    "class Weather(object):\n",
    "    '''The weather info for every day requested are saved in the dictionary historical {'year-month-day:weather df}\n",
    "    '''\n",
    "    lat=None\n",
    "    lon=None\n",
    "    historical={}\n",
    "    stations=None\n",
    "    state=None\n",
    "    icao=None\n",
    "    dataset=pd.DataFrame()\n",
    "    daysUnavailable=[]\n",
    "    daysList=[]\n",
    "    closestStation=None\n",
    "    key=0\n",
    "    def __init__(self,lat,lon):\n",
    "        '''Given latitude and longitude it find the closest weather station\n",
    "        \n",
    "        it will be used after to find weather  informations'''\n",
    "        self.parser=ParseWeather()\n",
    "        self.city,self.country,self.state=self.parser.getLocation(lat,lon)\n",
    "                                                                                                                                                                                                                                                                                    \n",
    "    def addDay(self,a,date='DateUTC'): \n",
    "        '''Add another day to the historical dictionary'''\n",
    "        df=pd.DataFrame(a)\n",
    "        i=lambda x: str(x.year) + '-' + str(x.month) + '-' +str(x.day) #I take just year,month and day \n",
    "        try:\n",
    "             dates= df[date].apply(i) \n",
    "        except AttributeError:\n",
    "            df=df.convert_objects(convert_dates='coerce')\n",
    "            dates= df[date].apply(i)\n",
    "        df['Day']=dates\n",
    "        daysList=dates.unique()\n",
    "        [self.historical.update({day:df[df['Day']==day]}) for day in daysList] \n",
    "        [self.daysList.append(day) for day in daysList]\n",
    "        return self\n",
    "    \n",
    "    def apply(self,f):\n",
    "        '''Apply a function on historical data'''\n",
    "        self.historical.update((x, f(y)) for x, y in self.historical.items())\n",
    "        return self\n",
    "    \n",
    "    @staticmethod\n",
    "    def clean(df):\n",
    "        '''Clean a specific dataframe containing weather informations'''\n",
    "        info=df.copy()\n",
    "        info=info.convert_objects(convert_numeric=True)\n",
    "        pre={'Light Rain':1,'Heavy Rain':1,'Rain':1,'Light Rain Mist':1,   \\\n",
    "        'Heavy Rain Mist':1,'Rain Mist':1,'Light Rain Showers':1,'Heavy Rain Showers':1,   \\\n",
    "        'Rain Showers':1,'Light Thunderstorms and Rain':1,'Heavy Thunderstorms and Rain':1,   \\\n",
    "        'Thunderstorms and Rain':1,'Light Freezing Drizzle':1,'Heavy Freezing Drizzle':1,      \\\n",
    "         'Freezing Drizzle':1,'Light Freezing Rain':1,'Heavy Freezing Rain':1,'Freezing Rain':1, \\\n",
    "        'Light Snow':1,'Heavy Snow':1,'Snow':1,'Light Snow Grains':1,'Heavy Snow Grains':1, \\\n",
    "        'Snow Grains':1,'LightSnow Showers':1,'Heavy Snow Showers':1,'Snow Showers':1,\n",
    "        'Light Ice Crystals':1,'Heavy Ice Crystals':1,'Ice Crystals':1,'Light Ice Pellets':1,  \\\n",
    "        'Heavy Ice Pellets':1,'Ice Pellets':1,'LightIce Pellet Showers':1,'HeavyIce Pellet Showers':1,   \\\n",
    "        'Ice Pellet Showers':1,'LightHail Showers':1,'Heavy Hail Showers':1, \\\n",
    "        'Hail Showers':1,'Light Small Hail Showers':1,'Heavy Small Hail Showers':1, \\\n",
    "        'Small Hail Showers':1}\n",
    "        f=lambda x: pre.get(str(x) , 0)        \n",
    "        info['Conditions']=info['Conditions'].apply(f)\n",
    "        \n",
    "        #cleaning of NaN and other unexpected values\n",
    "        info.PrecipitationIn=info.PrecipitationIn.fillna(value=0)\n",
    "        info['Wind SpeedMPH']=info['Wind SpeedMPH'].fillna(value=0)\n",
    "        info['Wind Direction']=info['Wind Direction'].replace('Calm',0)\n",
    "        info['Wind SpeedMPH']=info['Wind SpeedMPH'].replace('Calm',0)\n",
    "        #windspeedmph contains strings so it is considered as a generic object type, I convert it in float type\n",
    "        info['Wind SpeedMPH']=info['Wind SpeedMPH'].apply(float)\n",
    "        t=info.TemperatureF.copy()    \n",
    "        h=info.Humidity.copy()\n",
    "        s=info['Sea Level PressureIn'].copy()    \n",
    "        d=info['Dew PointF'].copy()\n",
    "        info['PrecipitationIn']=info['PrecipitationIn']+info['Snow']\n",
    "        info=info.drop('Snow',1)\n",
    "        p=info['PrecipitationIn'].copy()\n",
    "        #sometimes the weather informations show unexpected values (as -9999)\n",
    "        t[t < -100] = np.NaN\n",
    "        h[h<0]=np.NaN\n",
    "        s[s<0]=np.NaN\n",
    "        d[d<0]=np.NaN\n",
    "        p[p<0]=0\n",
    "        info['TemperatureF']=t\n",
    "        info['Humidity']=h\n",
    "        info['Sea Level PressureIn']=s\n",
    "        info['Dew PointF']=d\n",
    "        info['PrecipitationIn']=p\n",
    "        return info\n",
    "    \n",
    "    def conditionsOccurred(self,graph=False):\n",
    "        '''It returns the weather conditions occurred in the dataset'''\n",
    "        conditions=self.dataset.Conditions.value_counts()\n",
    "        print conditions\n",
    "        self.conditions=self.dataset.Conditions.value_counts()\n",
    "        if graph:\n",
    "            conditions.plot(kind='barh')\n",
    "        return self\n",
    "    \n",
    "    def createDataset(self):\n",
    "        '''It merges all the dataframe in the historical dictionary in a single dataframe\n",
    "        '''\n",
    "        tmp=self.historical.values()\n",
    "        df = pd.concat(tmp)\n",
    "        self.dataset=df#.sort('DateUTC')\n",
    "        return self.dataset\n",
    "    \n",
    "    @staticmethod\n",
    "    def extractHour(df,date='DateUTC',delta=0.25): \n",
    "        '''It creates a new field hour\n",
    "        \n",
    "        The field contains the hour in the format Hour.quarters (i.e 13.25 are 13 hours and 15 mins)'''\n",
    "        f=lambda x: round(round(x/(60*delta))*delta,3)\n",
    "        try:\n",
    "             hour=df[date].apply(lambda x: x.hour)\n",
    "        except AttributeError:\n",
    "            df[date]=df[date].convert_objects(convert_dates='coerce')\n",
    "            hour=df[date].apply(lambda x: x.hour)    \n",
    "        minute=df[date].dt.minute.apply(f)\n",
    "        df['Hour']=hour+minute\n",
    "        df['Hour']=df['Hour'].replace(24,0.00)\n",
    "        return df\n",
    "    \n",
    "    def extractHours(self,date='DateUTC',delta=0.25):\n",
    "        '''It applies the extractHour static method on every day \n",
    "        '''\n",
    "        self.historical.update((x, Weather.extractHour(y,date,delta)) for x, y in self.historical.items() )\n",
    "        return self\n",
    "    \n",
    "    def getDays(self):\n",
    "        '''It simply prints the days with weather information available in the instance'''\n",
    "        print self.weather.keys()\n",
    "    \n",
    "    def getHistorical(self, date):\n",
    "        '''Given a specific day it extract the weather information from wunderground.com\n",
    "        '''\n",
    "        s=self.state\n",
    "        c=self.city\n",
    "        key=date[:10]\n",
    "        fmt=\"%Y-%m-%d\"\n",
    "        date=dt.datetime.strptime(key,fmt)\n",
    "        day=date.day\n",
    "        date1=date-dt.timedelta(days=1)\n",
    "        date=str(date)\n",
    "        date1=str(date1)\n",
    "        df1=self.parser.getWeather(date,self.city,self.state)\n",
    "        df2=self.parser.getWeather(date1,self.city,self.state)\n",
    "        df1['Day']=df1['DateUTC'].apply(lambda x: x.day)\n",
    "        df2['Day']=df2['DateUTC'].apply(lambda x: x.day)\n",
    "        df1=df1[df1['Day']==day]\n",
    "        df2=df2[df2['Day']==day]\n",
    "        df=df1.append(df2)\n",
    "        df=df.drop('Day',1)\n",
    "        df=Weather.clean(df)\n",
    "        self.historical[key]=df\n",
    "        self.daysList.append(key)\n",
    "        df=Weather.clean(df)\n",
    "        return df\n",
    "    \n",
    "    def timeSampling(self,date='DateUTC'):\n",
    "        from numpy import median\n",
    "        df=self\n",
    "        df=df.clean()\n",
    "        return median([(df[date].loc[n]-df[date].loc[m]).total_seconds() for n,m in zip(xrange(1,len(df)),xrange(0,(len(df)-1)))])\n",
    "    \n",
    "\n",
    "    \n",
    "class Model(object):\n",
    "    '''This class contains method to prediction the background radiation using a dataframe with background\n",
    "     and weather informations\n",
    "    '''\n",
    "    debug={}\n",
    "    outliers=None\n",
    "    reducedDatasets=None\n",
    "    weather_columns=['Humidity','TemperatureF','Sea Level PressureIn','PrecipitationIn','Dew PointF','Conditions','Wind SpeedMPH']\n",
    "    out_columns=['Value']\n",
    "    #model_columns=['Value','PrecipitationIn','Humidity','Dew PointF','Sea Level PressureIn','TemperatureF']\n",
    "    columns=['Captured Time','Humidity','TemperatureF','Sea Level PressureIn','Conditions','PrecipitationIn','Dew PointF','Value','Wind SpeedMPH']\n",
    "        \n",
    "    def __init__(self,df):\n",
    "        from sklearn import preprocessing \n",
    "        self.ModelInputs={}\n",
    "        self.ModelOutput=None\n",
    "        self.prediction=None\n",
    "        self.metrics={}\n",
    "        self.Threats=[]\n",
    "        self.OutputTest={}\n",
    "        self.CorrelationTable=pd.DataFrame()\n",
    "        self.datasetsAvailable=['Dataset']\n",
    "        self.Sensor=df.Sensor.iloc[0]\n",
    "        self.model_columns=['PrecipitationIn','Humidity','Dew PointF','Sea Level PressureIn','TemperatureF']\n",
    "        '''Define a model object '''\n",
    "        df=df[Model.columns]\n",
    "        df=df.convert_objects(convert_dates='coerce')\n",
    "        df=self.clean(df)\n",
    "        t=df['Captured Time'].iloc[0]\n",
    "        f=lambda x: (x-t).total_seconds()\n",
    "        index=df['Captured Time'].apply(f)\n",
    "        #df=df.drop('Captured Time',1)\n",
    "        self.time=index\n",
    "        df.index=index\n",
    "        self.dataset=df\n",
    "        \n",
    "        \n",
    "    def applyOnInputs(self,method,inp,f=None,window=0,percentage=60):\n",
    "        '''It applies a built-in methods or a custom function f to the input variables\n",
    "        \n",
    "        Methods available:  \n",
    "                            'standardize' , it applies the standardization method of sklearn.preprocessing.scale\n",
    "        '''\n",
    "        if not(self.ModelInputs):\n",
    "            self.getInput()\n",
    "        index=int(percentage*len(self.dataset)/100)\n",
    "        d={'Train':self.ModelInputs[inp][:index,:],'Test':self.ModelInputs[inp][index:,:]}\n",
    "        if method=='standardize':\n",
    "            d.update((x, preprocessing.scale(y)) for x, y in d.items())\n",
    "        else:\n",
    "            d.update((x, f(y)) for x, y in d.items())\n",
    "            \n",
    "        #debug\n",
    "        #dataset=pd.DataFrame(self.ModelInputs['Dataset'])\n",
    "        #dataset['Output']=self.ModelOutput\n",
    "        #self.debug['ApplyOnInputs']=dataset\n",
    "        ###\n",
    "        self.ModelInputs[inp]=np.append(d['Train'],d['Test'],axis=0)\n",
    "        return self\n",
    "    \n",
    "    def applyOnOutput(self,method,f=None,window=0,percentage=60):\n",
    "        '''It applies a built-in methods or a custom function f to the output variable\n",
    "        \n",
    "        Methods available:  'movingaverage', it requires the variable window\n",
    "                            'standardize' , it applies the standardization method of sklearn.preprocessing.scale\n",
    "        '''\n",
    "        \n",
    "        if self.ModelOutput==None:\n",
    "            self.getOutput()\n",
    "        index=int(percentage*len(self.dataset)/100)\n",
    "        self.OutputTest['Original']=self.ModelOutput[index:]\n",
    "        #this function it's used to apply some filtering to the output\n",
    "        #for this reason the data are splitted , in this way every filtering technique won't be anticasual\n",
    "        #i.e. a moving average filtering on the train part will consider also some samples from the test part\n",
    "        #that belong ideally to the \"future\"\n",
    "        d={'Train':self.ModelOutput[:index],'Test':self.ModelOutput[index:]}\n",
    "        if method=='movingaverage':\n",
    "            if not(window):\n",
    "                raise ValueError('A value for the window is required')\n",
    "            d.update((x, Model.moving_average(y,n=window)) for x, y in d.items())\n",
    "            \n",
    "        elif method=='standardize':\n",
    "            self.OutputTest['mean']=np.mean(d['Train'])\n",
    "            self.OutputTest['std']=np.std(d['Train'])\n",
    "            d.update((x, preprocessing.scale(y)) for x, y in d.items())\n",
    "        else:\n",
    "            d.update((x, f(y)) for x, y in d.items())\n",
    "        newOutput=np.append(d['Train'],d['Test'])\n",
    "        \n",
    "        #the moving_average could drop some values at the end of the time series, so if this happens the last\n",
    "        #values is repeated to restore the original dimension\n",
    "        check=len(self.ModelOutput)-len(newOutput)\n",
    "        if check>0:\n",
    "            newOutput=np.append(newOutput,newOutput[-check:])\n",
    "        self.ModelOutput=newOutput\n",
    "        '''        #debug\n",
    "        dataset=pd.DataFrame(self.ModelInputs['Dataset'])\n",
    "        dataset['Output']=self.ModelOutput\n",
    "        self.debug['ApplyOnOutputs']=dataset\n",
    "        ###'''\n",
    "        \n",
    "        return self\n",
    "    \n",
    "    def clean(self,dataset):\n",
    "        dataset.Value=dataset.Value.replace(0,np.nan)\n",
    "        #a weighted interpolation is applied on a windows that correspond to a period of 3 hours\n",
    "        #just for the weather conditions\n",
    "        colnames=['Humidity','TemperatureF','Sea Level PressureIn','Conditions','PrecipitationIn','Dew PointF']\n",
    "        dataset[colnames]=dataset[colnames].replace(np.nan,999) \n",
    "        #the rolling apply function require that there are no nan values, so I use a dummy number\n",
    "        dataset[colnames]=pd.rolling_apply(dataset[colnames],13,Model.weightedInterp)\n",
    "        #at the end a linear interpolation it is used on value field and to fulfill the weather conditions in\n",
    "        #the case that some period had no value to interpolate\n",
    "        dataset=dataset.interpolate(method='linear')\n",
    "        dataset=dataset.dropna() #it drops the NaT captured Time\n",
    "        return dataset\n",
    "    \n",
    "    @staticmethod\n",
    "    def clustering(var1,var2):\n",
    "        '''Given two variables it find the clusters according the Meanshift algorithm\n",
    "        The current function is used by the remove_outliers method\n",
    "        '''\n",
    "\n",
    "        X=[var1,var2]\n",
    "        X=np.array(X)\n",
    "        X=X.T\n",
    "        bandwidth = estimate_bandwidth(X, quantile=0.9, n_samples=500) #estimation of bandwidth parameter needed for the \n",
    "                                                                    #clustering\n",
    "        ms = MeanShift(bandwidth=bandwidth, bin_seeding=True)\n",
    "        ms.fit(X)\n",
    "        labels = ms.labels_\n",
    "        tmp=pd.DataFrame(X)\n",
    "        tmp['Label']=labels\n",
    "        return tmp\n",
    "    \n",
    "    def conditionsOccurred(self,graph=False):\n",
    "        '''It returns the weather conditions occurred in the dataset, if the Condition field has not transformed in \n",
    "        a numerical field yet\n",
    "        '''\n",
    "        conditions=self.dataset.Conditions.value_counts()\n",
    "        print conditions\n",
    "        self.conditions=self.dataset.Conditions.value_counts()\n",
    "        if graph:\n",
    "            conditions.plot(kind='barh')\n",
    "        return self\n",
    "    \n",
    "    @staticmethod\n",
    "    def createDataset(sens,printing=False,filename='Stationary_data_with_weather.csv'):\n",
    "        '''This function instantiates the objects Weather and Sensor, use their method to clean and collect informations \n",
    "\n",
    "        Then merge them in a dataset containing weather and radiation information\n",
    "        '''\n",
    "        w=None\n",
    "        s=Sensor(sens)\n",
    "        s=s.cleanAll()\n",
    "        sensor=s.extractDates() \n",
    "        #value of lat and lon needed to instantiate the weather class\n",
    "        lat,lon=sensor.latitude,sensor.longitude\n",
    "        w= Weather(lat,lon)\n",
    "        \n",
    "        for day in sensor.daysList: \n",
    "            w.getHistorical(day)\n",
    "        #the historical weather has a sampling time of 1 hour, so I resample my sensor data every (15 min default)\n",
    "        #taking the median of the value in that period\n",
    "        wea=w.extractHours()\n",
    "        f= lambda x: x.groupby(x.Hour).median() \n",
    "        wea=wea.apply(f)\n",
    "        wea=wea.apply(lambda x: x.drop('Hour',1))\n",
    "        sensor=sensor.apply(f)\n",
    "        #pieces contains a list of dataframe corresponding to a single day of measurements coupled with the weater\n",
    "        #dataframe with all the measurements coupled\n",
    "        try:\n",
    "            pieces=[sensor.days[date].join(wea.historical[date]) for date in wea.daysList  if not(wea.historical[date].empty) ]\n",
    "        except ValueError:\n",
    "            return pd.DataFrame()\n",
    "        #to make the single days well sampled the holes are filled with a linear interpolation method\n",
    "        #the first and the last are skipped because the first piece probably doesn't start at midnight so it would be filled\n",
    "        #with NaN\n",
    "        #for the last is the same, it probably doesn't finish at midnight\n",
    "        filled=[p.reindex(np.arange(0,24,0.25)).interpolate(method='linear') for num,p in enumerate(pieces) if (num!=0 and num!=len(pieces)-1) ]\n",
    "        try:\n",
    "            filled.insert(0,pieces[0])\n",
    "        except IndexError:\n",
    "            return pd.DataFrame()\n",
    "        filled.append(pieces[-1])\n",
    "        try:\n",
    "            \n",
    "            dataset=pd.concat(filled)\n",
    "        except ValueError:\n",
    "            return pd.DataFrame()\n",
    "        #after the median on every hour all the field that were string become NaN or are dropped\n",
    "        dataset.dropna(1,how='all')\n",
    "        \n",
    "        dataset = dataset[np.isfinite(dataset['Sensor'])] \n",
    "        dataset['Hour']=dataset.index\n",
    "        dataset.drop\n",
    "        #in the line below the field Captured Time is recreated\n",
    "        dataset=Sensor.convertDate(dataset)\n",
    "        if printing:\n",
    "            with open(filename, 'a') as f:\n",
    "                   dataset.to_csv(f,index=False,float_format = '%.4f',header=False)\n",
    "        return dataset\n",
    "    \n",
    "    def destandardize(self):\n",
    "        std=self.OutputTest['std']\n",
    "        mean=self.OutputTest['mean']\n",
    "        array=[(pred*std)+mean for pred in self.prediction]\n",
    "        return array\n",
    "    \n",
    "    def dimensionalityReduction(self,nr=5):\n",
    "        '''It applies all the dimensionality reduction techniques available in this class:\n",
    "        Techniques available:\n",
    "                            'PCA'\n",
    "                            'FactorAnalysis'\n",
    "                            'KPCArbf','KPCApoly'\n",
    "                            'KPCAcosine','KPCAsigmoid'\n",
    "                            'IPCA'\n",
    "                            'FastICADeflation'\n",
    "                            'FastICAParallel'\n",
    "                            'Isomap'\n",
    "                            'LLE'\n",
    "                            'LLEmodified'\n",
    "                            'LLEltsa'\n",
    "        '''\n",
    "        dataset=self.ModelInputs['Dataset']\n",
    "        sklearn_pca = sklearnPCA(n_components=nr)\n",
    "        p_components = sklearn_pca.fit_transform(dataset)\n",
    "        fa=FactorAnalysis(n_components=nr)\n",
    "        factors=fa.fit_transform(dataset)\n",
    "        kpca=KernelPCA(nr,kernel='rbf')\n",
    "        rbf=kpca.fit_transform(dataset)\n",
    "        kpca=KernelPCA(nr,kernel='poly')\n",
    "        poly=kpca.fit_transform(dataset)\n",
    "        kpca=KernelPCA(nr,kernel='cosine')\n",
    "        cosine=kpca.fit_transform(dataset)\n",
    "        kpca=KernelPCA(nr,kernel='sigmoid')\n",
    "        sigmoid=kpca.fit_transform(dataset)\n",
    "        ipca=IncrementalPCA(nr)\n",
    "        i_components=ipca.fit_transform(dataset)\n",
    "        fip=FastICA(nr,algorithm='parallel')\n",
    "        fid=FastICA(nr,algorithm='deflation')\n",
    "        ficaD=fip.fit_transform(dataset)\n",
    "        ficaP=fid.fit_transform(dataset)\n",
    "        '''isomap=Isomap(n_components=nr).fit_transform(dataset)\n",
    "        try:\n",
    "            lle1=LocallyLinearEmbedding(n_components=nr).fit_transform(dataset)\n",
    "        except ValueError:\n",
    "            lle1=LocallyLinearEmbedding(n_components=nr,eigen_solver='dense').fit_transform(dataset)\n",
    "        try:\n",
    "            \n",
    "            lle2=LocallyLinearEmbedding(n_components=nr,method='modified').fit_transform(dataset)\n",
    "        except ValueError:\n",
    "            lle2=LocallyLinearEmbedding(n_components=nr,method='modified',eigen_solver='dense').fit_transform(dataset) \n",
    "        try:\n",
    "            lle3=LocallyLinearEmbedding(n_components=nr,method='ltsa').fit_transform(dataset)\n",
    "        except ValueError:\n",
    "            lle3=LocallyLinearEmbedding(n_components=nr,method='ltsa',eigen_solver='dense').fit_transform(dataset)'''\n",
    "        values=[p_components,factors,rbf,poly,cosine,sigmoid,i_components,ficaD,ficaP]#,isomap,lle1,lle2,lle3]\n",
    "        keys=['PCA','FactorAnalysis','KPCArbf','KPCApoly','KPCAcosine','KPCAsigmoid','IPCA','FastICADeflation','FastICAParallel']#,'Isomap','LLE','LLEmodified','LLEltsa']\n",
    "        self.ModelInputs.update(dict(zip(keys, values)))\n",
    "        [self.datasetsAvailable.append(key) for key in keys ]\n",
    "        \n",
    "        #debug\n",
    "        #dataset=pd.DataFrame(self.ModelInputs['Dataset'])\n",
    "        #dataset['Output']=self.ModelOutput\n",
    "        #self.debug['Dimensionalityreduction']=dataset\n",
    "        ###\n",
    "        return self\n",
    "    \n",
    "    @staticmethod\n",
    "    def extractMetrics(pred,test_y):\n",
    "        '''It extracts three different metrics: mean absolute error,median absolute error,mean square error\n",
    "\n",
    "        '''\n",
    "        try:\n",
    "            meanae=mean_absolute_error(test_y,pred)\n",
    "        except ValueError:\n",
    "            #sometimes the moving average filter on the output reduce the dimensionality of it\n",
    "            #so some value of the predition is dropped\n",
    "            pred=pred[:len(test_y)-len(pred)]\n",
    "            meanae=mean_absolute_error(test_y,pred)\n",
    "        mae=median_absolute_error(test_y,pred)\n",
    "        mse=mean_squared_error(test_y,pred)\n",
    "        return meanae,mae,mse\n",
    "    \n",
    "    def findCorrelations(self,alfa=5,duringRain=False,minimumLength=500):\n",
    "        '''It discovers if the input variables are correlated with the output making use of Spearman correlation technique\n",
    "        \n",
    "        The alfa parameter define the level of significance of the test,it is expressed in percentage\n",
    "        If the p-value evaluated is less than alfa/100 the Null Hypotesis (there is no correlation between the variables) is refused'''\n",
    "        e=self.dataset\n",
    "        if duringRain:\n",
    "            e=e[e['Conditions']==1]\n",
    "        e=e[Model.weather_columns]\n",
    "        e['Value']=self.dataset.Value.copy()\n",
    "        e=e.apply(preprocessing.scale)\n",
    "        if len(e)<minimumLength:\n",
    "            self.CorrelationTable=pd.DataFrame()\n",
    "            return self\n",
    "        pthresh=alfa/100.0\n",
    "        val=e.Value.values\n",
    "        temp=spearmanr(e.TemperatureF.values,val)\n",
    "        hum=spearmanr(e.Humidity.values,val)\n",
    "        sea=spearmanr(e['Sea Level PressureIn'].values,val)\n",
    "        prec=spearmanr(e.PrecipitationIn.values,val)\n",
    "        dew=spearmanr(e['Dew PointF'].values,val)\n",
    "        df=pd.DataFrame({'Temperature':temp,'Sea Level PressureIn':sea,'PrecipitationIn':prec,'Humidity':hum,'Dew PointF':dew},index=['Pearson coefficients','p-values'])\n",
    "        def test(p,threshold):\n",
    "            if p<threshold:\n",
    "                return 'Reject H0'\n",
    "            else:\n",
    "                return 'Accept H0'\n",
    "        df.loc['Results']=[test(p,pthresh) for p in df.loc['p-values']]\n",
    "        self.CorrelationTable=df\n",
    "        return self\n",
    "        \n",
    "    \n",
    "    def GBregression(self,percentage=60,inp='Dataset',n_estimators=100, learning_rate=0.1,max_depth=1, random_state=0, loss='ls'):\n",
    "        '''It applies the ensamble method of gradient boosting trees'''\n",
    "        X=y=prediction=metrics=None\n",
    "        X=self.ModelInputs[inp] #input dataset\n",
    "        samples=int(percentage*len(X)/100) #evaluating the samples number given the percentage\n",
    "        x=X[:samples,:] #training input set\n",
    "        try:\n",
    "            y = self.ModelOutput[:samples] #training output set\n",
    "        except KeyError:\n",
    "            self.getOutput()\n",
    "            y = self.ModelOutput[:samples]\n",
    "        test_x=X[samples:,:] #testing input set\n",
    "        test_y=self.ModelOutput[samples:] # testing output set\n",
    "        gb=GradientBoostingRegressor(n_estimators=100, learning_rate=0.1,max_depth=1, random_state=0, loss='ls')\n",
    "        model=gb.fit(x,y)\n",
    "        prediction=model.predict(test_x)\n",
    "        self.prediction=prediction\n",
    "        self.OutputTest['Standardized']=test_y\n",
    "        metrics=Model.extractMetrics(prediction,test_y)\n",
    "        return prediction,np.median(metrics)\n",
    "    \n",
    "    def getDatasetsAvailable(self):\n",
    "        self.datasetsAvailable=self.ModelInputs.keys()\n",
    "        return self.ModelInputs.keys()\n",
    "    \n",
    "    def getInput(self):\n",
    "        X=self.dataset[self.model_columns].copy()\n",
    "        self.ModelInputs['Dataset']=X.as_matrix()\n",
    "        return self\n",
    "        \n",
    "    def getOutput(self):\n",
    "        Y=self.dataset.copy()\n",
    "        try:\n",
    "            self.ModelOutput=Y[Model.out_columns].as_matrix()\n",
    "        except KeyError:\n",
    "            self.ModelOutput=self.dataset['Output'].as_matrix() #if the preparare dataset has been called \n",
    "            #the output is 'Output' instead of 'Values\n",
    "        return self\n",
    "    \n",
    "    def insertThreat(self,testPercentage=40,wLength=4,meanP=1.1):\n",
    "        '''Method to simulate and insert a threat in the part of the output series that will be used as test\n",
    "        wLenght: the lenght of the window in which the threat will be inserted\n",
    "        testPercentage: indicates the percentage of the test dataset\n",
    "        meanP: is the mean value of the Poisson distribution from which the \"threat\" is extracted\n",
    "        '''\n",
    "        t=None\n",
    "        testPercentage=testPercentage/100.0\n",
    "        t=pd.DataFrame()\n",
    "        t['Value']=self.dataset.Value.copy()#create a copy of the output\n",
    "        startTest=int((1-testPercentage)*len(t)) #define the first index of the output that will be used as test\n",
    "        s=np.random.random_integers(startTest,len(t)) #find a random index in the test part of the output\n",
    "        values=np.random.poisson(t['Value'].mean()*meanP,wLength) #find random values from poisson distribution with E[x]=m\n",
    "        window=np.arange(s,s+4)*(self.dataset.index[1]-self.dataset.index[0]) #define the window\n",
    "        #the window is cleaned, the values are added and the other values are interpolated to maintain the continuity\n",
    "        t['Value'].loc[window]=values\n",
    "        #t.loc[window[1:-1]]=values\n",
    "        self.ThreatsIndex=t.copy()\n",
    "        self.ThreatsIndex['Value']=0\n",
    "        self.ThreatsIndex.loc[window]=1\n",
    "        \n",
    "        d={'Train':t['Value'].iloc[:startTest],'Test':t['Value'].iloc[startTest:]}\n",
    "        d.update((x, preprocessing.scale(y)) for x, y in d.items())      \n",
    "        self.Threats=np.append(d['Train'],d['Test'])#append the window in which there is the threat \n",
    "        self.dataset.Value=t['Value'].values.copy() #the threat is inserted in the dataset\n",
    "        return self\n",
    "\n",
    "    def KNregression(self,percentage,inp='Dataset',neighbors=5,weights='distance',algorithm='auto',leaf=30):\n",
    "        '''It evaluates a prediction using k-nearest neighbors regression approach\n",
    "        \n",
    "        It returns a tuple: (prediction, median of three different metrics) '''\n",
    "        X=y=prediction=metrics=None\n",
    "        X=self.ModelInputs[inp] #input matrix\n",
    "        samples=int(percentage*len(X)/100) #evaulating the number of samples given the percentage\n",
    "        x=X[:samples,0:] #training input set\n",
    "        y = self.ModelOutput[:samples] # training output set\n",
    "        test_x=X[samples:,:] #testing input set\n",
    "        test_y=self.ModelOutput[samples:] #testing output set\n",
    "        knn=KNeighborsRegressor(n_neighbors=neighbors,weights=weights,algorithm=algorithm, leaf_size=leaf)\n",
    "        try:\n",
    "            model=knn.fit(x,y) #evaluating the model\n",
    "        except ValueError:\n",
    "            return np.nan,9999\n",
    "        prediction=model.predict(test_x) #evaluating of the prediction\n",
    "        self.prediction=prediction\n",
    "        self.OutputTest['Standardized']=test_y\n",
    "        metrics=Model.extractMetrics(prediction,test_y)\n",
    "        return prediction,np.median(metrics)\n",
    "    \n",
    "    @staticmethod   \n",
    "    def moving_average(a, n=3) :\n",
    "        ''' Function that implements a moving average filter\n",
    "            [source]:http://stackoverflow.com/questions/14313510/moving-average-function-on-numpy-scipy    \n",
    "        '''\n",
    "        first=np.array([a[0]])\n",
    "        last=np.array([a[-1]])\n",
    "        a=np.concatenate((first,a,last))\n",
    "        ret = np.cumsum(a, dtype=float)\n",
    "        ret[n:] = ret[n:] - ret[:-n]\n",
    "        return ret[n - 1:] / n\n",
    "    \n",
    "    def plotRadiationWeather(self):\n",
    "        '''It plots the Value field with each weather field separately\n",
    "\n",
    "        The function returns a plot object\n",
    "        '''\n",
    "\n",
    "        df=self.dataset\n",
    "        plt.figure()\n",
    "\n",
    "        stand=df.apply(preprocessing.scale,axis=0) #the data are normalized because they have different units\n",
    "        val=stand['Value'].as_matrix()\n",
    "        prec=stand['PrecipitationIn'].as_matrix()\n",
    "        dew=stand['Dew PointF'].as_matrix()\n",
    "        hum=stand['Humidity'].as_matrix()\n",
    "        press=stand['Sea Level PressureIn'].as_matrix()\n",
    "        temp=stand['TemperatureF'].as_matrix()\n",
    "        plt.subplot(3,3,1)\n",
    "        plt.plot(val,prec,'bo')\n",
    "        plt.ylabel('Precipitation')\n",
    "        plt.xlabel('Background Radiation')\n",
    "        plt.subplot(3,2,2)\n",
    "        plt.plot(val,dew,'ro')\n",
    "        plt.ylabel('Dew Point')\n",
    "        plt.xlabel('Background Radiation')\n",
    "        plt.subplot(3,2,3)\n",
    "        plt.plot(val,hum,'yo')\n",
    "        plt.ylabel('Humidity')\n",
    "        plt.xlabel('Background Radiation')\n",
    "        plt.subplot(3,2,4)\n",
    "        plt.plot(val,press,'go')\n",
    "        plt.ylabel('Sea Level Pressure')\n",
    "        plt.xlabel('Background Radiation')\n",
    "        plt.subplot(3,2,5)\n",
    "        plt.plot(val,temp,'mo')\n",
    "        plt.ylabel('Temperature')\n",
    "        plt.xlabel('Background Radiation')\n",
    "        plt.subplot(3,2,6)\n",
    "        plt.plot(val,prec,'bo')\n",
    "        plt.plot(val,dew,'ro')\n",
    "        plt.plot(val,hum,'yo')\n",
    "        plt.plot(val,press,'go')\n",
    "        plt.plot(val,temp,'mo')\n",
    "        #plt.legend(['Precipitation','DewPoint','Humidity','Sea Level Pressure','Temperature'])\n",
    "        plt.xlabel('Background Radiation')\n",
    "        plt.show()\n",
    "        \n",
    "    def plotDataset(self):\n",
    "        self.dataset.plot(subplots=True)\n",
    "        plt.xlabel('Time')\n",
    "        plt.show()\n",
    "        \n",
    "    def  plotPrediction(self):\n",
    "        '''It creates a figure with two graphs: the real and the predicted output\n",
    "                                                the absolute error between them\n",
    "        '''\n",
    "        predicted=self.prediction\n",
    "        real=self.OutputTest['Standardized']#[abs(len(self.OutputTest['Standardized'])-len(self.prediction)):]\n",
    "\n",
    "        rmse=np.sqrt(mean_squared_error(predicted,real))\n",
    "        plt.figure()\n",
    "        plt.subplot(211)\n",
    "        plt.xlabel('Time')\n",
    "        plt.ylabel('Radiation ')\n",
    "        plt.title('Comparison between real and predicted output, RMSE=' + str(rmse))\n",
    "        plt.plot(predicted,'r')\n",
    "        plt.plot(real,'b')\n",
    "        plt.legend(['Predicted output','Real output'])\n",
    "        plt.subplot(212)\n",
    "        plt.xlabel('Time')\n",
    "        plt.ylabel('Absolute error')\n",
    "        plt.plot(abs(real-predicted),'m')\n",
    "        plt.show()\n",
    "        \n",
    "    def prepareDataset(self,n=1,l=1,w=0):\n",
    "        X=self.dataset[Model.weather_columns].copy()\n",
    "        self.model_columns=Model.weather_columns[:] #this fake slicing provide a copy of the list \n",
    "        values=self.dataset.Value.copy()\n",
    "        output=values.shift(-l).copy()\n",
    "        vfield=[]\n",
    "        for m in xrange(0,n+1): #the n parameter sets how much new fields should be created \n",
    "                #if the present value of the output is at the time t there will be created n columns with\n",
    "                #output from 0,1,2,...t-1 , 0,1,2,...t-2, ....... 0,1,2,...t-n\n",
    "            field='Values-' + str(m)\n",
    "            vfield.append(field)\n",
    "            self.model_columns.append(field)\n",
    "            X[field]=values.shift(m) #the shift function creates the new fields \n",
    "        for k in xrange(1,w+1):\n",
    "            a=X[Model.weather_columns].shift(k)\n",
    "            newfields=[col+'-' +str(w) for col in a.columns]\n",
    "            a.columns=newfields\n",
    "            #[self.model_columns.append(f) for f in newfields]\n",
    "            X=pd.concat([X,a], axis=1)\n",
    "        X['Output']=output\n",
    "        X=X.dropna()\n",
    "        ##debug    \n",
    "        #dataset=X.copy()\n",
    "        #dataset['Output']=output.copy()\n",
    "        #self.debug['getInput']=dataset\n",
    "        ##\n",
    "        self.dataset=X.copy()\n",
    "        return self\n",
    "    \n",
    "    def reduceDataset(self,nr=3,method='PCA'):\n",
    "        '''It reduces the dimensionality of a given dataset using different techniques provided by Sklearn library\n",
    "         Methods available:\n",
    "                            'PCA'\n",
    "                            'FactorAnalysis'\n",
    "                            'KPCArbf','KPCApoly'\n",
    "                            'KPCAcosine','KPCAsigmoid'\n",
    "                            'IPCA'\n",
    "                            'FastICADeflation'\n",
    "                            'FastICAParallel'\n",
    "                            'Isomap'\n",
    "                            'LLE'\n",
    "                            'LLEmodified'\n",
    "                            'LLEltsa'\n",
    "        '''\n",
    "        dataset=self.ModelInputs['Dataset']\n",
    "        #dataset=self.dataset[Model.in_columns]\n",
    "        #dataset=self.dataset[['Humidity','TemperatureF','Sea Level PressureIn','PrecipitationIn','Dew PointF','Value']]\n",
    "        #PCA\n",
    "        if method=='PCA':\n",
    "            sklearn_pca = sklearnPCA(n_components=nr)\n",
    "            reduced = sklearn_pca.fit_transform(dataset)\n",
    "        #Factor Analysis\n",
    "        elif method=='FactorAnalysis':\n",
    "            fa=FactorAnalysis(n_components=nr)\n",
    "            reduced=fa.fit_transform(dataset)\n",
    "        #kernel pca with rbf kernel\n",
    "        elif method=='KPCArbf':\n",
    "            kpca=KernelPCA(nr,kernel='rbf')\n",
    "            reduced=kpca.fit_transform(dataset)\n",
    "        #kernel pca with poly kernel\n",
    "        elif method=='KPCApoly':\n",
    "            kpca=KernelPCA(nr,kernel='poly')\n",
    "            reduced=kpca.fit_transform(dataset)\n",
    "        #kernel pca with cosine kernel\n",
    "        elif method=='KPCAcosine':\n",
    "            kpca=KernelPCA(nr,kernel='cosine')\n",
    "            reduced=kpca.fit_transform(dataset)\n",
    "        #kernel pca with sigmoid kernel\n",
    "        elif method=='KPCAsigmoid':\n",
    "            kpca=KernelPCA(nr,kernel='sigmoid')\n",
    "            reduced=kpca.fit_transform(dataset)\n",
    "        #ICA\n",
    "        elif method=='IPCA':\n",
    "            ipca=IncrementalPCA(nr)\n",
    "            reduced=ipca.fit_transform(dataset)\n",
    "        #Fast ICA\n",
    "        elif method=='FastICAParallel':\n",
    "            fip=FastICA(nr,algorithm='parallel')\n",
    "            reduced=fip.fit_transform(dataset)\n",
    "        elif method=='FastICADeflation':\n",
    "            fid=FastICA(nr,algorithm='deflation')\n",
    "            reduced=fid.fit_transform(dataset)\n",
    "        elif method == 'All':\n",
    "            self.dimensionalityReduction(nr=nr)\n",
    "            return self\n",
    "        \n",
    "        self.ModelInputs.update({method:reduced})\n",
    "        self.datasetsAvailable.append(method)\n",
    "        return self\n",
    "    \n",
    "    def remove_outliers(self):\n",
    "        '''It removes the outliers using the MeanShift clustering techniques\n",
    "        '''\n",
    "        dataset=self.dataset[self.model_columns].copy()\n",
    "        dataset['Value']=self.dataset.Value.copy()\n",
    "        stand=dataset.apply(preprocessing.scale,axis=0) #the data are standardized because they have different units\n",
    "        val=stand['Value'].as_matrix()\n",
    "        prec=stand['PrecipitationIn'].as_matrix()\n",
    "        dew=stand['Dew PointF'].as_matrix()\n",
    "        hum=stand['Humidity'].as_matrix()\n",
    "        press=stand['Sea Level PressureIn'].as_matrix()\n",
    "        temp=stand['TemperatureF'].as_matrix()\n",
    "        l=[Model.clustering(val,b) for b in [prec,dew,hum,press,temp] ]\n",
    "        l1=[a.groupby('Label').count().index[0] for a in l ] #it finds the cluster with most of the data\n",
    "        l2=[a[a['Label']!=lab] for a,lab in zip(l,l1)] #the biggest cluster is removed in every dataframe\n",
    "        outliers=pd.concat(l2,join='inner',axis=1).index #the concat with join='inner' option find the intersection between                                              \n",
    "        #the dataframes, the resulting indexes indicate the outliers\n",
    "        #the indexes in outliers are not expressed in seconds\n",
    "        #so I create a fake index\n",
    "        index=list(xrange(0,len(stand)))\n",
    "        #and I remove the indexes that corresponds to the outliers\n",
    "        [index.remove(a) for a in outliers ] \n",
    "        #using iloc I remove them from the original dataset\n",
    "        self.dataset.Value.iloc[outliers]=np.nan\n",
    "        #the dropped value are replaced using a linear interpolation\n",
    "        self.dataset.Value=self.dataset.Value.interpolate(method='linear')\n",
    "        self.dataset=self.dataset.dropna()\n",
    "        index=self.dataset.index-self.dataset.index[0]\n",
    "        self.dataset.index=index\n",
    "        self.outliers=outliers #the outliers are saved\n",
    "        \n",
    "        #DEBUG\n",
    "        self.debug['Removeoutliers']=dataset\n",
    "        ###\n",
    "        \n",
    "        return self\n",
    "\n",
    "    def SVregression(self,percentage,inp='Dataset',kern='rbf',method='standard',c=2048,eps=0,gamma=0.01,tau=3):\n",
    "        '''Given the dataset of the input X and the dataset of the output Y it find a regression model using\n",
    "        Support vector regression algorithm of sklearn library\n",
    "        \n",
    "        It returns a tuple: (prediction, median of three different metrics)\n",
    "        '''       \n",
    "        \n",
    "        \n",
    "        X=y=prediction=metrics=None\n",
    "        X=self.ModelInputs[inp].copy() #input dataset\n",
    "        samples=int(percentage*len(X)/100) #evaluating the samples number given the percentage\n",
    "        x=X[:samples,:] #training input set\n",
    "        try:\n",
    "            y = self.ModelOutput[:samples] #training output set\n",
    "        except KeyError:\n",
    "            self.getOutput()\n",
    "        y = self.ModelOutput[:samples]\n",
    "        test_x=X[samples:,:] #testing input set\n",
    "        test_y=self.ModelOutput[samples:] # testing output set\n",
    "\n",
    "        #Parameters settings based on \"Selection of Meta-Parameters for support vector regression\" \n",
    "        # Vladimir Cherkassky and Yunqian Ma\n",
    "        if method=='standard':\n",
    "            n=len(y)\n",
    "            std=y.std()\n",
    "            c=tau*std\n",
    "            eps=tau*np.sqrt(log(n)/n)\n",
    "        #regression\n",
    "        svr =SVR(kernel=kern,C=c,epsilon=eps,gamma=gamma)\n",
    "        m=None\n",
    "        try:\n",
    "            m=svr.fit(x,y)\n",
    "        except ValueError:\n",
    "            return np.nan,9999\n",
    "            \n",
    "        #debug\n",
    "        #self.debug['SVR']=self.ModelOutput\n",
    "            \n",
    "            \n",
    "        prediction=m.predict(test_x)\n",
    "        self.prediction=prediction\n",
    "        self.OutputTest['Standardized']=test_y\n",
    "        metrics=Model.extractMetrics(prediction,test_y)\n",
    "        return prediction,np.median(metrics)\n",
    "\n",
    "    \n",
    "    @staticmethod\n",
    "    def weightedInterp(array):\n",
    "        l=int(len(array)/2)\n",
    "        if array[l]!=999:\n",
    "            return array[6]\n",
    "\n",
    "        #other weight function could be inserted using scipy.signal module\n",
    "        a=list(np.arange(1,l+1))\n",
    "        l1=[(n*m,m)  for n,m in zip(array[0:6],a) if n!=999]\n",
    "        a.reverse()\n",
    "        l2=[(n*m,m)  for n,m in zip(array[7:13],a) if n!=999]\n",
    "        try:\n",
    "            num=reduce(lambda x,y: x+y, [x[0] for x in l1+l2])\n",
    "        except TypeError:\n",
    "            return np.nan\n",
    "        den= reduce(lambda x,y: x+y, [x[1] for x in l1+l2])\n",
    "        return num/den\n",
    "    \n",
    "\n",
    "class ParseMap(object):\n",
    "    '''Class that implements usefull methods to parse OpenStreetMap xml files'''\n",
    "    way={}\n",
    "    node={}\n",
    "    coord={}\n",
    "    way_limit={}\n",
    "    way_City={}\n",
    "    way_Street={}\n",
    "    way_coor={}\n",
    "    '''\n",
    "    \n",
    "    #notes:\n",
    "    #the use of the tag_filter seems slower than a simple if-then\n",
    "    #not used at the moment\n",
    "    whitelist = set(('name', 'highway'))\n",
    "    \n",
    "    #unused\n",
    "    def tag_filter(tags):\n",
    "        for key in tags.keys():\n",
    "            if key not in whitelist:\n",
    "                del tags[key]\n",
    "        if 'name' in tags and len(tags) == 1:\n",
    "            # tags with only a name have no information\n",
    "            # how to handle this element\n",
    "            del tags['name']\n",
    "    '''\n",
    "    def ways_stationary(self,ways):\n",
    "        for osmid, tags, refs in ways:\n",
    "            if tags.has_key('building'): \n",
    "                self.way[osmid]=refs\n",
    "                if tags.has_key('addr:city'):  #sometimes the ways have also the city name in tags\n",
    "                    self.way_City[osmid]=tags['addr:city']\n",
    "                else:\n",
    "                    self.way_City[osmid]=None\n",
    "                if tags.has_key('name'): \n",
    "                    self.way_Street[osmid]=tags['name']\n",
    "                else:\n",
    "                    self.way_Street[osmid]=None\n",
    "\n",
    "    def ways(self,ways):\n",
    "        for osmid, tags, refs in ways:\n",
    "            if tags.has_key('highway'): #just the streets are needed \n",
    "                self.way[osmid]=refs\n",
    "                if tags.has_key('addr:city'):  #sometimes the ways have also the city name in tags\n",
    "                    self.way_City[osmid]=tags['addr:city']\n",
    "                else:\n",
    "                    self.way_City[osmid]=None\n",
    "                if tags.has_key('name'): \n",
    "                    self.way_Street[osmid]=tags['name']\n",
    "                else:\n",
    "                    self.way_Street[osmid]=None\n",
    "                    \n",
    "    def nodes(self,nodes):\n",
    "        for idnode,tag,coor in nodes:\n",
    "            lat=coor[1] #it's necessary because the coordinates in the nodes \n",
    "            lon=coor[0] #are (lon,lat) while in the coords are (lat,lon)\n",
    "            self.node[idnode]=((lat,lon), tag)\n",
    "            \n",
    "    def coords(self,coords):\n",
    "          for osm_id, lon, lat in coords:\n",
    "            self.coord[osm_id]=(lat,lon)\n",
    "\n",
    "    def fill_way_coords(self): #return a dictionary: {osmid:[list of nodes coordinates]}\n",
    "        for osmid in self.way.keys():\n",
    "            l=[]\n",
    "            for ref in self.way[osmid]:\n",
    "                try:\n",
    "                    val=self.node[ref][0]\n",
    "                except KeyError:\n",
    "                    val=self.coord[ref]\n",
    "                l.append(val)\n",
    "            self.way_coor[osmid]=l\n",
    "                  \n",
    "    def getRange(self):\n",
    "        for osmid in self.way.keys():\n",
    "            a=self.way_coor[osmid]\n",
    "            c=map(list, zip(*a)) #to unzip a list of tuples [(lat1,lon1),(lat2,lon2)] in [ [lat1,lat2),(lon1,lon2)]\n",
    "            lat=c[0]\n",
    "            lon=c[1]\n",
    "            self.way_limit[osmid]=[min(lat),min(lon),max(lat),max(lon)]\n",
    "            \n",
    "            \n",
    "class ParseWeather(object):\n",
    "    '''Class that implement methods to get the weather informations from wunderground.com\n",
    "    '''\n",
    "    key=0\n",
    "    def __init__(self):\n",
    "        if not(ParseWeather.key):\n",
    "            raise Exception('Key is not present, register at http://www.wunderground.com/weather/api/ to get one')\n",
    "    def getLocation(self,lat,lon):\n",
    "        '''Given latitude and longitude it returns the city,country and state corresponding to the coordinates '''\n",
    "        key=ParseWeather.key\n",
    "        url_template='http://api.wunderground.com/api/{key}/geolookup/q/{latitude},{longitude}.json'\n",
    "        url=url_template.format(key=key,latitude=lat,longitude=lon)\n",
    "        g = urllib2.urlopen(url)\n",
    "        json_string = g.read()\n",
    "        location = json.loads(json_string)\n",
    "        g.close()\n",
    "        diz=location['location']['nearby_weather_stations']['airport']['station'][0]\n",
    "        return diz['city'].replace(' ','_'),diz['country'],diz['state']\n",
    "    \n",
    "    def getWeather(self,date,c,s):\n",
    "        '''Given a date a city and a state it returns a DataFrame '''\n",
    "        k=ParseWeather.key\n",
    "        d=date[:10].replace('-','')\n",
    "        url_template='http://api.wunderground.com/api/{key}/history_{date}/q/{state}/{city}.json'\n",
    "        url=url_template.format(key=k,date=d,state=s,city=c)  \n",
    "        f = urllib2.urlopen(url)\n",
    "        json_string = f.read()\n",
    "        weather = json.loads(json_string) #parsing the json\n",
    "        f.close()\n",
    "        forecast=weather['history']['observations']\n",
    "        l=[]\n",
    "        for n in xrange(0,len(forecast)):\n",
    "            #every cycle define a row containing the weather information for a single hour\n",
    "\n",
    "            tmp=pd.DataFrame(forecast[n]) #definition of the dataframe\n",
    "            col=['utcdate','tempi','dewpti','hum','pressurei','visi','wdire','wspdi','precipi','conds','snow','wdird']\n",
    "            year=tmp.ix['year','utcdate'] #info about the day are extracted\n",
    "            month=tmp.ix['mon','utcdate']\n",
    "            day=tmp.ix['mday','utcdate']\n",
    "            hour=tmp.ix['hour','utcdate']\n",
    "            minute=tmp.ix['min','utcdate']\n",
    "            date= year +'-' + month + '-' + day + ' ' + hour + ':' + minute + ':00'\n",
    "            #the name of the columns are changed\n",
    "            newcol=['DateUTC', 'TemperatureF', 'Dew PointF', 'Humidity',\n",
    "                   'Sea Level PressureIn', 'VisibilityMPH', 'Wind Direction',\n",
    "                   'Wind SpeedMPH',  'PrecipitationIn', 'Conditions','Snow',\n",
    "                  'WindDirDegrees']\n",
    "            tmp=tmp[col]\n",
    "            tmp.columns=newcol\n",
    "            tmp=tmp.head(1)\n",
    "            tmp['DateUTC']=date\n",
    "            tmp.index=[hour]\n",
    "            l.append(tmp)\n",
    "            newdate=date[:10]\n",
    "        df=pd.concat(l) #all the weather info are concatenated in a single dataframe\n",
    "        df=df.convert_objects(convert_dates='coerce')\n",
    "        return df\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#Functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def EvaluateResults(filename,plot=False,printing=False):\n",
    "    columns=['Technique','10','15','20','25','30','35','40','45','50','55','60','Sensor','Measurements','Outliers']\n",
    "    results=pd.read_csv(filename,names=columns)\n",
    "    results.index=results.Technique\n",
    "    results=results.drop('Technique',1)\n",
    "    grouped= results.groupby('Sensor')\n",
    "    a={}\n",
    "    b={}\n",
    "    lowestError={}\n",
    "    for group in grouped:\n",
    "        if group[1].Measurements.iloc[0]>500:\n",
    "             #it finds the technique and the best percentage for every sensor\n",
    "            a[group[0]]=group[1].idxmin().value_counts().index[0] #it finds the technique that results better most of the time\n",
    "            b[group[0]]=group[1].iloc[:,:9].min().idxmin() #it find the best percentage for training set\n",
    "            lowestError[group[0]]=min(group[1].iloc[:,:9].min()) #it find the lowest error for every person\n",
    "    metrics=pd.Series(a.values()).value_counts() #it counts the occurences of every metric\n",
    "    percentage=pd.Series(b.values()).value_counts() #it counts the occurrences of every percentage\n",
    "    errors=lowestError.values() \n",
    "    df=pd.DataFrame(lowestError.values()) #it creates a dataframe with the lowest error for every sensor\n",
    "    df.index=lowestError.keys()\n",
    "    if plot:\n",
    "        plt.figure(1)\n",
    "        metrics.plot(kind='barh',colormap='winter')\n",
    "        plt.title('Dimensionality techniques that gave the best metric')\n",
    "        plt.figure(2)\n",
    "        percentage.plot(kind='barh',color='r')\n",
    "        plt.title('Percentage of training set that gave the best metric')\n",
    "        plt.figure(3)\n",
    "        plt.plot(errors,'go')\n",
    "        plt.plot(errors,'r')\n",
    "        plt.title('Error between predicted and real output,median value:' + str(round(np.median(errors),2)))\n",
    "        plt.show()\n",
    "    if printing:\n",
    "        print 'The dimensionality reduction technique that gives best result in terms of metric achived is ' + metrics.index[0] \n",
    "        print 'The best percentage for training set is ' + percentage.index[0] + '%'\n",
    "        print 'The median of the lowest error achieved by every sensor is ' +str(np.median(errors))\n",
    "        print 'These results are achieved analyzing ' + str(len(a)) + ' sensors, with more than 500 measurements, that corresponds to more than 5 consecutive days'\n",
    "    return metrics,percentage,df"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#Dataset"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The dataset is obtained by cleaning the data from safecast with pig and clustering the resulting dataset with hadoop streaming functionality"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "colnames=['Captured Time','Latitude','Longitude',  \\\n",
    "          'Value','Unit','Location','ID','Height', \\\n",
    "          'Surface','Radiation','Upload Time','Loader ID','Sensor','Distance']\n",
    "#US_cleanedResult contains all the measurements clustered, it was obtained using hadoop\n",
    "data = pd.read_csv('US_results.csv',header=None,names=colnames)\n",
    "data.columns=colnames\n",
    "data=data.sort('Captured Time')\n",
    "# I am looking just for the stationary detectors\n",
    "grouped = data.groupby('Sensor')\n",
    "filtered=grouped.filter(lambda x: len(x)>10) #it drops the group that have less then 10 elements\n",
    "regrouped=filtered.groupby('Sensor')\n",
    "stationary=regrouped.filter(Sensor.isStationary) #it returns just the meassurements that belong to stationary detectors\n",
    "stationarySensor=stationary.Sensor.unique()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#Code to get weather"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The stationary detectors are coupled with weather informations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "\n",
    "d={}\n",
    "start_time = time.time()\n",
    "for n in stationarySensor:\n",
    "    print n\n",
    "    if len(data[data['Sensor']==n])<1000:\n",
    "        continue\n",
    "    sens=data[data['Sensor']==n]\n",
    "    d[n]=Model.createDataset(sens,printing=True)\n",
    "elapsed_time = time.time() - start_time\n",
    "print 'Datasets creation:' + str(elapsed_time)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#Complete Dataset"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The dataset obtained by the previous block is used to test the model developed, the function could be used on every dataset that contains both weather and radiation data "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "p=pd.read_csv('Stationary_data_with_weather.csv')\n",
    "colnames=['Latitude', 'Longitude', 'Value', 'ID', 'Height',\n",
    "       'Loader ID','Sensor', 'Distance',\n",
    "      'TemperatureF', 'Dew PointF', 'Humidity', 'Sea Level PressureIn',\n",
    "       'VisibilityMPH', 'Wind SpeedMPH', 'PrecipitationIn', 'Conditions',\n",
    "       'WindDirDegrees', 'Captured Time']\n",
    "p.columns=colnames"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/usr/lib/python2.7/site-packages/IPython/kernel/__main__.py:416: FutureWarning: comparison to `None` will result in an elementwise object comparison in the future.\n"
     ]
    }
   ],
   "source": [
    "sensor=406\n",
    "r=p[p.Sensor==sensor]\n",
    "model=Model(r)\n",
    "model.dataset\n",
    "model=model.remove_outliers()\n",
    "model=model.prepareDataset(n=4,w=0,l=1)\n",
    "model=model.getInput()\n",
    "model=model.getOutput()\n",
    "model=model.applyOnOutput(method='movingaverage',window=4)\n",
    "model=model.applyOnOutput(method='standardize',percentage=60)\n",
    "pred,error=model.SVregression(60,inp='Dataset',kern='rbf',method='custom',c=2048,eps=0.5,gamma=0.001)\n",
    "model.plotPrediction()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-0.97516874, -1.85757412, -1.56558594, ...,  0.69355992,\n",
       "        0.86026461,  0.58941807])"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pred"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-1.29959832, -0.44613713, -0.44613713, ...,  0.97629818,\n",
       "        0.12283699,  0.97629818])"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model.OutputTest['Standardized']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Defining the threshold level for the threats: according to IAEA a threat is considered as a radiation level higher 1.4 times the average background"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "m=[group['Value'].mean() for name,group in p.groupby('Sensor')]\n",
    "std=[group['Value'].std() for name,group in p.groupby('Sensor')]\n",
    "backgroundMean=np.median(m)\n",
    "backgroundStd=np.median(std)\n",
    "s=1.4*backgroundMean/backgroundStd\n",
    "s"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "According to the results, for the detectors in use to avoid false alarms the threshold should be 24sigma higher than average background"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#Testing of SVR , rbf kernel\n",
    "#Testing of SVR , linear kernel\n",
    "#Testing of KN"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "start_time = time.time()          \n",
    "\n",
    "\n",
    "def printResults(dic,n,filename,printing=True):\n",
    "    dataframe=pd.DataFrame.from_dict(dic)\n",
    "    dataframe['Sensor']=n\n",
    "    dataframe['Total samples']=len(model.dataset)\n",
    "    dataframe['Outliers removed']=len(model.outliers)\n",
    "\n",
    "    if printing:\n",
    "        with open(filename, 'a') as f:\n",
    "                       dataframe.to_csv(f,header=False)#,index=False,header=False)\n",
    "    return dataframe\n",
    "li=list(p.Sensor.unique())\n",
    "li.sort()\n",
    "\n",
    "for dd,l in zip(p.groupby('Sensor'),li):\n",
    "    dic1={}\n",
    "    dic2={}\n",
    "    dic3={}\n",
    "    dic4={}\n",
    "    sensor=dd[0]\n",
    "    print sensor\n",
    "    df=dd[1]\n",
    "    if l<2346:\n",
    "        continue\n",
    "    if len(df)<500:\n",
    "        continue\n",
    "    model=None\n",
    "    model=Model(df)\n",
    "    model=model.remove_outliers()\n",
    "    model=model.prepareDataset(n=4,w=0,l=1)\n",
    "    model=model.getInput()\n",
    "    model=model.getOutput()\n",
    "    for percentage in xrange(10,65,5): #different percentage of training set\n",
    "        model=model.applyOnOutput(method='movingaverage',window=4)\n",
    "        model=model.applyOnOutput(method='standardize',percentage=percentage)\n",
    "        model=model.reduceDataset(method='All',nr=10)\n",
    "        keys=model.getDatasetsAvailable()\n",
    "        li=[model.applyOnInputs(inp=e,method='standardize',percentage=percentage) for e in keys]\n",
    "        model=li[-1]\n",
    "        values1=[model.SVregression(percentage,inp=a,kern='rbf')[1] for a in keys]\n",
    "        values4=[model.KNregression(percentage,inp=b,neighbors=5)[1] for b in keys]\n",
    "        values2=[model.SVregression(percentage,inp=c,kern='sigmoid')[1] for c in keys]\n",
    "        values3=[model.SVregression(percentage,inp=d,kern='linear')[1] for d in keys]\n",
    "\n",
    "        dic1[percentage]=dict(zip(keys,values1))\n",
    "        dic2[percentage]=dict(zip(keys,values2))\n",
    "        dic3[percentage]=dict(zip(keys,values3))\n",
    "        dic4[percentage]=dict(zip(keys,values4))\n",
    "        \n",
    "    filenames=['SVRrbf.csv','SVRsigmoid.csv','SVRlinear.csv','KN.csv']\n",
    "    dictionaries=[dic1,dic2,dic3,dic4]\n",
    "    df=[printResults(dic=d,filename=f,n=sensor) for d,f in zip(dictionaries,filenames)]\n",
    "\n",
    "elapsed_time = time.time() - start_time\n",
    "\n",
    "print 'TestModel time:' + str(elapsed_time)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#Evaluation of results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "a=EvaluateResults('SVRrbf.csv',plot=True,printing=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "b=EvaluateResults('SVRlinear.csv',plot=True,printing=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "c=EvaluateResults('SVRsigmoid.csv',plot=True,printing=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "d=EvaluateResults('KN.csv',plot=True,printing=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#Chi-square Test"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "percentage=60\n",
    "normalityReal={}\n",
    "normalityPred={}\n",
    "normalityDiff={}\n",
    "stdReal={}\n",
    "stdPred={}\n",
    "stdDiff={}\n",
    "chiReal={}\n",
    "chiDiff={}\n",
    "chiPred={}\n",
    "for SensorId in p.Sensor.unique():\n",
    "    r=p[p.Sensor==SensorId]\n",
    "    if len(r)<500:\n",
    "        continue\n",
    "    stdReal[SensorId]=np.nan\n",
    "    stdDiff[SensorId]=np.nan\n",
    "    stdPred[SensorId]=np.nan\n",
    "    model=Model(r)\n",
    "    model=model.applyOnInputs(method='standardize')\n",
    "    model=model.remove_outliers()\n",
    "    model=model.reduceDataset(method='FICAParallel')\n",
    "    model=model.applyOnOutput(method='standardize')\n",
    "    model=model.applyOnOutput(method='movingaverage',window=3)\n",
    "    a,b=model.SVregression(percentage,inp='FICAParallel',kern='linear')\n",
    "    real=model.ModelOutput[len(model.ModelOutput)-len(model.prediction):]\n",
    "    predicted=model.prediction\n",
    "    difference=real-predicted\n",
    "    p1=normaltest(difference)\n",
    "    p2=normaltest(real)\n",
    "    p3=normaltest(predicted)\n",
    "    if p2[1]>0.05:\n",
    "        stdReal[SensorId]=np.std(real)\n",
    "    if p1[1]>0.05:\n",
    "        stdDiff[SensorId]=np.std(difference)\n",
    "    if p3[1]>0.05:\n",
    "        stdPred[SensorId]=np.std(predicted)\n",
    "    normalityReal[SensorId]=p2[1]\n",
    "    normalityPred[SensorId]=p3[1]\n",
    "    normalityDiff[SensorId]=p1[1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plots"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "plt.subplot(311)\n",
    "plt.title('P-value for the real output')\n",
    "plt.plot(normalityReal.values(),'bo')\n",
    "plt.plot(normalityReal.values(),'b')\n",
    "plt.plot([0.98]*len(normalityDiff),'g')\n",
    "plt.plot([0.02]*len(normalityDiff),'g')\n",
    "plt.legend(['p-value','','Threshold = 0.02','Threshold= 0.98'])\n",
    "plt.subplot(312)\n",
    "plt.title('P-value for the predicted output')\n",
    "plt.plot(normalityPred.values(),'ro')\n",
    "plt.plot(normalityPred.values(),'r')\n",
    "plt.plot([0.98]*len(normalityDiff),'g')\n",
    "plt.plot([0.02]*len(normalityDiff),'g')\n",
    "plt.legend(['p-value','','Threshold = 0.02','Threshold= 0.98'])\n",
    "plt.subplot(313)\n",
    "plt.title('P-value for the (real-predicted) output')\n",
    "plt.xlabel('Sensors')\n",
    "plt.plot(normalityDiff.values(),'mo')\n",
    "plt.plot(normalityDiff.values(),'m')\n",
    "plt.plot([0.98]*len(normalityDiff),'g')\n",
    "plt.plot([0.02]*len(normalityDiff),'g')\n",
    "plt.legend(['p-value','','Threshold = 0.02','Threshold= 0.98'])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##Parameters Selection"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "#r=p[p.Sensor==1721]\n",
    "params={}\n",
    "bestinput='Dataset'\n",
    "percentage=50\n",
    "start=time.time()\n",
    "li=p.Sensor.unique()\n",
    "li.sort()\n",
    "for n,l in zip(li,xrange(0,1)):\n",
    "    print(n)\n",
    "    model=None\n",
    "    r=p[p.Sensor==n]\n",
    "    if len(r)<500:\n",
    "        continue\n",
    "        \n",
    "    model=Model(r)\n",
    "    model=model.remove_outliers()\n",
    "    model=model.prepareDataset(n=4,w=0,l=1)\n",
    "    model=model.getInput()\n",
    "    model=model.getOutput()\n",
    "    if bestinput!='Dataset':\n",
    "        model=model.reduceDataset(method=bestinput,nr=10)\n",
    "    model=model.applyOnInputs(inp=bestinput,method='standardize',percentage=percentage)\n",
    "    model=model.applyOnOutput(method='movingaverage',window=4)\n",
    "    model=model.applyOnOutput(method='standardize',percentage=percentage)\n",
    "\n",
    "    X=model.ModelInputs[bestinput] #input dataset\n",
    "    samples=int(percentage*len(X)/100) #evaluating the samples number given the percentage\n",
    "    x=X[:samples,0:] #training input set\n",
    "    y = model.ModelOutput[:samples] #training output set\n",
    "    test_x=X[samples:,:] #testing input set\n",
    "    test_y=model.ModelOutput[samples:]\n",
    "    scores = ['precision', 'recall']\n",
    "\n",
    "    svr = GridSearchCV(SVR(),\n",
    "                       param_grid={\"C\": 2**np.arange(1,10),\n",
    "                                   \"gamma\": np.logspace(-2, 2, 10),\n",
    "                                    \"epsilon\" : [0, 0.01, 0.1, 0.5]})\n",
    "    svr.fit(x, y)\n",
    "    params[n]=svr.best_params_\n",
    "\n",
    "print str(time.time()-start) + 's to complete the evaluation of the best parameters'\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Evaluation of the results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "pp=pd.read_csv('Parameters_for_raw_dataset.csv')\n",
    "pp.index=pp.Results\n",
    "pp=pp.drop('Results',1)\n",
    "plt.figure()\n",
    "plt.title('C Parameters')\n",
    "pp.loc['C'].value_counts().plot(kind='barh')\n",
    "plt.figure()\n",
    "plt.title('Epsilon Parameters')\n",
    "pp.loc['epsilon'].value_counts().plot(kind='barh')\n",
    "plt.figure()\n",
    "plt.title('Gamma Parameters')\n",
    "pp.loc['gamma'].value_counts().plot(kind='barh')\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#Correlation study"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The code uses the Pearson correlation technique to every sensor to see if there is some correlation between every input variable and the output variable"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "start=time.time()\n",
    "test={}\n",
    "for n in p.Sensor.unique():\n",
    "    #print(n)\n",
    "    model=None\n",
    "    r=p[p.Sensor==n]\n",
    "    if len(r)<500:\n",
    "        '''The p-values are not entirely\n",
    "        reliable but are probably reasonable for datasets larger than 500 or so.'''\n",
    "        continue\n",
    "    model=Model(r)\n",
    "    model.findCorrelations()\n",
    "    df=model.CorrelationTable\n",
    "    if df.empty:\n",
    "        continue\n",
    "    test[n]=df.loc['Results']\n",
    "elapsed_time = time.time() - start\n",
    "print 'TestModel time:' + str(elapsed_time)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plots"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "ee=pd.DataFrame(test)\n",
    "\n",
    "plt.subplot(321)\n",
    "plt.title('Significance test for correlation between background radiation and weather data ')\n",
    "plt.ylabel('Value-PrecipitationIn')\n",
    "ee.loc['PrecipitationIn'].value_counts().plot(kind='barh')\n",
    "plt.subplot(322)\n",
    "plt.title('HO: there is no correlation between the two variables')\n",
    "plt.ylabel('Value-Humidity')\n",
    "ee.loc['Humidity'].value_counts().plot(kind='barh')\n",
    "plt.subplot(323)\n",
    "plt.ylabel('Value-Dew PointF')\n",
    "ee.loc['Dew PointF'].value_counts().plot(kind='barh')\n",
    "plt.subplot(324)\n",
    "plt.ylabel('Value-Temperature')\n",
    "ee.loc['Temperature'].value_counts().plot(kind='barh')\n",
    "plt.xlabel('')\n",
    "plt.subplot(325)\n",
    "plt.ylabel('Value-Sea Level PressureIn')\n",
    "ee.loc['Sea Level PressureIn'].value_counts().plot(kind='barh')\n",
    "plt.show()"
   ]
  }
 ],
 "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.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}