%matplotlib inline

import numpy as np
import pandas as pd

import matplotlib as mpl
import matplotlib.dates as mdates
from matplotlib import pyplot as plt
import seaborn as sns

# Set random
np.random.seed(42)

import sys
sys.path.append('../')

from prediction import (datareader, complete_data, cleanup, bikes_probability,
                        time_resampling)

%load_ext watermark

%watermark -d -v -p numpy,pandas,matplotlib -g -m -w

def plot_features_station(result, station, features_to_plot='paa', nb_row=350, draw_type='-'):
    """Plot available bikes and bike stands for a given station"""
    
    data = result[result.station == station].tail(nb_row).copy()
    
    fig, ax = plt.subplots(figsize=(18,5))
    plt.plot(data.index, data.probability, draw_type, label='probability', alpha=0.8)
    plt.plot(data.index, data[features_to_plot], draw_type, label=features_to_plot, alpha=0.6)
    
    
    ax = plt.gca()
    # set major ticks location every day
    ax.xaxis.set_major_locator(mdates.DayLocator())
    # set major ticks format
    ax.xaxis.set_major_formatter(mdates.DateFormatter('\n\n\n%a %d.%m.%Y'))
    # set minor ticks location every one hours
    ax.xaxis.set_minor_locator(mdates.HourLocator(interval=1))
    # set minor ticks format
    ax.xaxis.set_minor_formatter(mdates.DateFormatter('%H:%M'))
    plt.setp(ax.xaxis.get_minorticklabels(), rotation=45)

    plt.legend(loc='best')

def plot_features_station_train_test(train, test, station, features_to_plot='paa', nb_row=350, draw_type='-'):
    """Plot available bikes and bike stands for a given station for a train / test dataset"""
    
    train = train[train.station == station].tail(nb_row).copy()
    test = test[test.station == station].copy()
    
    
    fig, ax = plt.subplots(figsize=(18,5))
    plt.plot(train.index, train.probability, draw_type, label=' Train probability', alpha=0.8)
    plt.plot(train.index, train[features_to_plot], draw_type, label='Train ' + features_to_plot, alpha=0.6)
    
    plt.plot(test.index, test.probability, draw_type, label=' Test probability', alpha=0.8)
    plt.plot(test.index, test[features_to_plot], draw_type, label='Test ' + features_to_plot, alpha=0.6)
    
    
    ax = plt.gca()
    # set major ticks location every day
    ax.xaxis.set_major_locator(mdates.DayLocator())
    # set major ticks format
    ax.xaxis.set_major_formatter(mdates.DateFormatter('\n\n\n%a %d.%m.%Y'))
    # set minor ticks location every one hours
    ax.xaxis.set_minor_locator(mdates.HourLocator(interval=1))
    # set minor ticks format
    ax.xaxis.set_minor_formatter(mdates.DateFormatter('%H:%M'))
    plt.setp(ax.xaxis.get_minorticklabels(), rotation=45)

    plt.legend(loc='best')

"""
Explanation: Experiment and anlyse some features creation
End of explanation
"""


DATAFILE = '../data/lyon.csv'

raw = datareader(DATAFILE)

"""
Explanation: Read and prepare the data
End of explanation
"""


df_clean = cleanup(raw)

df_clean.head()

"""
Explanation: Cleaning dataset
End of explanation
"""


df = (df_clean.pipe(time_resampling)
      .pipe(complete_data)
      .pipe(bikes_probability))

df.head()

df.shape

df.info()

"""
Explanation: Create dataset for learning
End of explanation
"""


# params of learning dataset creation
start = pd.Timestamp("2017-08-01T02:00:00") # Tuesday
predict_date = pd.Timestamp("2017-09-22T09:00:00") # wednesday
# predict the next 30 minutes
freq = '1H'
# number of predictions at 'predict_date'.
# Here, the next 30 minutes and the next hour (30 minutes + 30 minutes).
# If you want to predict the next 3 hours, every 30 minutes, thus set periods=6
periods = 1

"""
Explanation: Create Target learning & analyse mechanics
End of explanation
"""


from prediction import get_summer_holiday

df = get_summer_holiday(df.copy())

df.head(2)

df.tail(2)

"""
Explanation: Here we don't care about DataFrame Sort (timeserie). We only do some matching and ratio
get_summer_holiday
Create a bool (0 classic day / 1 holidays)
End of explanation
"""


from prediction import get_public_holiday

df = get_public_holiday(df.copy(), count_day=5)

df[df.ts >='2017-08-14 23:50:00'].head()

"""
Explanation: get_public_holiday
Count day before and after special holiday (like assomption on 15/08)
End of explanation
"""


from prediction import cluster_station_lyon

"""
Explanation: cluster_station_lyon
End of explanation
"""


df = cluster_station_lyon(df.copy(), path_file='../data/cluster_lyon_armand.csv')

df.head()

"""
Explanation: Cluster of station by ativite (mean on bike by hours of day)
You can find the process of clustering in file ../clustering-Lyon-Armand.ipynb
End of explanation
"""


from prediction import cluster_station_geo_lyon

"""
Explanation: cluster_station_geo_lyon
End of explanation
"""


df = cluster_station_lyon(df.copy(), path_file='../data/station_cluster_geo_armand.csv')

df.head()

"""
Explanation: Cluster of Lyon's station by geography position
You can find the process of clustering in file ../Clustering-Lyon-geo-Armand.ipynb
End of explanation
"""


from prediction import mapping_hours

df['hours_binned'] = df.hour.apply(mapping_hours)

df.head()

"""
Explanation: Binned hour
End of explanation
"""


from prediction import get_statio_ratio_open_by_time

df_temp_1 = get_statio_ratio_open_by_time(df.copy())

df_temp_1.head()

"""
Explanation: get_statio_ratio_open_by_time
Ratio of open station by time
End of explanation
"""


from prediction import get_statio_cluster_geo_ratio_open_by_time

"""
Explanation: get_statio_cluster_geo_ratio_open_by_time
End of explanation
"""


df_temp_2 = get_statio_cluster_geo_ratio_open_by_time(df.copy())

df_temp_2.head()

"""
Explanation: Ratio of open station (on geography cluster) by hours
End of explanation
"""


data = df.sort_values(['station', 'ts']).set_index(["ts", "station"])
observation = 'probability'
label = data[observation].copy()
label.name = "future"
label = (label.reset_index(level=1)
         .shift(-1, freq=freq)
         .reset_index()
         .set_index(["ts", "station"]))
result = data.merge(label, left_index=True, right_index=True)

result.reset_index(level=1, inplace=True)
if start is not None:
    result = result[result.index >= start]

"""
Explanation: Now we sorting our DataFrame to create timeserie features ( /!\ order is important )
End of explanation
"""


data.head(15)

"""
Explanation: Step by step :
First step is to sort dataset by station and time ('ts')
End of explanation
"""


label[6:11]

"""
Explanation: Creating 'label' to shift probability to one hour.
probability at 2017-07-09 01:00:00 become futur at 2017-07-09 00:00:00
End of explanation
"""


result[result.station == 1001][['station', 'bikes', 'stands', 'probability', 'future']].head(15)

"""
Explanation: Merging label and data to create the learning dataset with target shifting
End of explanation
"""


# Original learning dataset :
result.head()

"""
Explanation: Creation of features
End of explanation
"""


from prediction import create_shift_features

df_temp_3 = create_shift_features(result.copy(), features_name='bikes_shift_'+str(freq.replace('H', 'bin')), feature_to_shift='bikes', 
                                features_grp='station', nb_shift=periods)

df_temp_3[['station', 'bikes', 'bikes_shift_1bin']].head(15)

"""
Explanation: create_shift_features
End of explanation
"""


from prediction import create_cumul_trend_features

# Need to use df_temp with 'bikes_shift_1bin' values
df_temp_4 = create_cumul_trend_features(df_temp_3, features_name='bikes_shift_'+str(freq.replace('H', 'bin')))

df_temp_4[df_temp_4.station == 1001][['station', 'bikes', 'bikes_shift_1bin',
                                     'cumsum_trend_sup', 'cumsum_trend_inf', 'cumsum_trend_equal']].head(8)

"""
Explanation: create_cumul_trend_features
End of explanation
"""


from prediction import get_station_recently_closed

"""
Explanation: get_station_recently_closed
End of explanation
"""


df[254350:254361][['station', 'ts', 'bikes', 'is_open', 'probability']]

df_temp_5 = get_station_recently_closed(result, nb_hours=4)

df_temp_5[['station', 'bikes', 'is_open', 'probability', 'was_recently_open']].tail()

"""
Explanation: Sometime station are closed for maintenance, so they can't be use by users. Trying to catch this information to help the learning process
End of explanation
"""


from prediction import filling_bike_on_geo_cluster

"""
Explanation: If station is open since 4 hours, was_recently_open has a value of 24 (4 * 6 (bin is egal to 10 min))
filling_bike_on_geo_cluster
End of explanation
"""


df_temp_6 = filling_bike_on_geo_cluster(df_temp_3.copy(), features_name='bikes_shift_'+str(freq.replace('H', 'bin')))

df_temp_6.tail()

"""
Explanation: Create a features on filling bike by geo station. This give information if some zone (cluster) are empty or full
End of explanation
"""


from prediction import get_paa_transformation

df_temp_7 = get_paa_transformation(result.copy(), features_to_compute='probability', segments=10)

df_temp_7[df_temp_7.station == 1005][['station', 'bikes', 'probability', 'future', 'paa']].tail(22)

plot_features_station(df_temp_7, station=1001, features_to_plot='paa', nb_row=120, draw_type='-o')

plot_features_station(df_temp_7, station=1001, features_to_plot='paa', nb_row=29, draw_type='-o')

df_temp_7[df_temp_7.station == 1001][['station', 'probability', 'future', 'paa']].tail(15)

df_temp_7[df_temp_7.station == 1001][['station', 'probability', 'future', 'paa']][-26:-16]

df_temp_7[df_temp_7.station == 1001].probability[-26:-16].mean()

"""
Explanation: get_paa_transformation
End of explanation
"""


from prediction import get_sax_transformation

df_temp_8 = get_sax_transformation(result.copy(), features_to_compute='probability', segments=10, symbols=8)

df_temp_8[df_temp_8.station == 1001].tail(22)

plot_features_station(df_temp_8, station=1001, features_to_plot='sax', nb_row=35, draw_type='-o')

"""
Explanation: There is data leak in this features (PAA). At 09:40, there is a probability of 0.062 (target is the same one hour later). But PAA is going to mean the next 9 values. So PAA will se the increase in the futur (0.187 / 0.187 / 0.125 - 30 min later) and PAA will be highter. Algorithm will catch it as a win information but can't see it in a production vision.
get_sax_transformation
End of explanation
"""


# Original
result.head()

from prediction import create_rolling_mean_features

df_temp_9 = create_rolling_mean_features(result, 
                                     features_name='mean_6', 
                                     feature_to_mean='probability', 
                                     features_grp='station', 
                                     nb_shift=6)

df_temp_9[df_temp_9.station == 1001].tail(15)

plot_features_station(df_temp_9, station=1001, features_to_plot='mean_6', nb_row=40, draw_type='-o')

"""
Explanation: As PAA, SAX transformation give data leak.
tranform signal with rolling mean
As PAA & SAX can give data leak, we will mean our probability on x bin (rolling mean)
End of explanation
"""


from prediction import create_rolling_std_features

df_temp_10 = create_rolling_std_features(result, 
                                     features_name='std_9', 
                                     feature_to_std='probability', 
                                     features_grp='station', 
                                     nb_shift=9)

plot_features_station(df_temp_10, station=4012, features_to_plot='std_9', nb_row=35, draw_type='-o')

"""
Explanation: Here there is no leak of information in the future. You only take past informations to give context to our algorithm 
Rolling Standard Deviation
Sometime station's bike don't move too much, and sometime it's crazy time. By given this indicator, we want to help our algorithm with context awareness
End of explanation
"""


from prediction import create_rolling_median_features

df_temp_11 = create_rolling_median_features(result, 
                                     features_name='median_6', 
                                     feature_to_median='probability', 
                                     features_grp='station', 
                                     nb_shift=6)

plot_features_station(df_temp_11, station=4012, features_to_plot='median_6', nb_row=40, draw_type='-o')

"""
Explanation: Rolling on median
End of explanation
"""


def create_bool_empty_full_station(df):
    """
    Create a bool features "warning_empty_full"
    If bike <= 2 --> 1
    If Proba >= 0.875 --> 1
    else --> 0
    """
    
    df['warning_empty_full'] = 0
    df.loc[df['bikes'] <= 2, 'warning_empty_full'] = 1
    df.loc[df['probability'] >= 0.875, 'warning_empty_full'] = 1
    
    return df

df_temp_12 = result.copy()

df_temp_12.head()

df_temp_12['warning_empty_full'] = 0
df_temp_12.loc[df_temp_12['bikes'] <= 2, 'warning_empty_full'] = 1
df_temp_12.loc[df_temp_12['probability'] >= 0.875, 'warning_empty_full'] = 1

feature_event_to_plot='warning_empty_full'
features_event_value=1

df_temp_13[df_temp_13[feature_event_to_plot] == features_event_value].head()

def plot_event_station(result, station, feature_event_to_plot='bike', features_event_value=1,
                       nb_row=350, point_type='*'):
    """Plot available bikes and bike stands for a given station"""
    
    data = reswult[result.station == station].tail(nb_row).copy()
    
    fig, ax = plt.subplots(figsize=(18,5))
    plt.plot(data.index, data.probability, '-', label='probability', alpha=0.8)
    plt.plot(data[data[feature_event_to_plot] == features_event_value].index, 
            data[data[feature_event_to_plot] == features_event_value].probability,
            point_type, markerfacecolor='k',
            label=feature_event_to_plot, alpha=0.6)
    
    
    ax = plt.gca()
    # set major ticks location every day
    ax.xaxis.set_major_locator(mdates.DayLocator())
    # set major ticks format
    ax.xaxis.set_major_formatter(mdates.DateFormatter('\n\n\n%a %d.%m.%Y'))
    # set minor ticks location every one hours
    ax.xaxis.set_minor_locator(mdates.HourLocator(interval=1))
    # set minor ticks format
    ax.xaxis.set_minor_formatter(mdates.DateFormatter('%H:%M'))
    plt.setp(ax.xaxis.get_minorticklabels(), rotation=45)

    plt.legend(loc='best')

plot_event_station(df_temp_12, station=10101, feature_event_to_plot='warning_empty_full',
                      features_event_value=1, nb_row=350, point_type='*')

"""
Explanation: Create bool for pretty empty or full station
End of explanation
"""


# to have same value
date = predict_date
print ('date : ' + str(date))
cut = date - pd.Timedelta(freq.replace('T', 'm'))
stop = date + periods * pd.Timedelta(freq.replace('T', 'm'))
print ('cut : ' + str(cut))
print ('stop : ' + str(stop))

train = result[result.index <= cut].copy()
mask = np.logical_and(result.index >= date, result.index <= stop)
test = result[mask].copy()

print('train shape : ' + str(train.shape))
print('test shape : ' + str(test.shape))

train.head()

"""
Explanation: Split learning dataset on train test (avoid data leak feature)
The train test split
We split our dataset to create on date :
    - A trainning Dataset
    - A test Dataset
End of explanation
"""


from prediction import create_mean_by_sta_day_binned_hours

train_temp_1, test_temp_1 = create_mean_by_sta_day_binned_hours(train.copy(), test.copy(), 
                                     features_name='proba_mean_by_sta_day_binned_hour', 
                                     feature_to_mean='probability', 
                                     features_grp=['station', 'day', 'hours_binned'])

plot_features_station_train_test(train_temp_1, test_temp_1, station=1036, features_to_plot='proba_mean_by_sta_day_binned_hour', 
                      nb_row=450, draw_type='-o')

"""
Explanation: Create KPI probability group on binned hour / month / day / is_open ==1
We need to create our binned hours mapping
End of explanation
"""


df_temp7 = result.copy()
df_temp7['ts'] = df_temp7.index
df_temp7 = df_temp7.sort_values(['station', 'ts'])

df_temp7.head()


df_temp7['prob_shit'] = df_temp7.groupby(['station'])['probability'].apply(lambda x: x.shift(1))
df_temp7['prob_diff'] = np.abs(df_temp7['prob_shit'] - df_temp7['probability'])

df_temp7.tail()

df_temp7[df_temp7.prob_diff >= 0.5].tail(6)

plot_features_station(df_temp7, station=1036, features_to_plot='prob_diff', nb_row=300, draw_type='-')

df_temp7[(df_temp7.station == 1036) & (df_temp7.ts >= '2017-09-26 02:30:00')].head(6)

"""
Explanation: Detect anomalie in probability du to some re stock by humain
End of explanation
"""


df_temp7['ano'] = 0
df_temp7.loc[df_temp7['prob_diff'] > 0.5, 'ano'] = 1

df_temp7[(df_temp7.station == 1036) & (df_temp7['ano']==1)][['prob_diff','day', 'hour', 'minute', 'bikes', 'prob_diff']].tail(50)

df_temp7[df_temp7['ano']==1].station.value_counts()

"""
Explanation: In 10 min at 03:00, 13 bikes have been loaded here. It's impossible to predict it as it's bike company who reload this station
End of explanation
"""


lyon_meteo = pd.read_csv('../data/lyon_weather.csv', parse_dates=['date'])
lyon_meteo.rename(columns={'date':'ts'}, inplace=True)

lyon_meteo.head()

"""
Explanation: Weather feature
Exact weather
End of explanation
"""


lyon_forecast = pd.read_csv('../data/lyon_forecast.csv', parse_dates=['forecast_at', 'ts'])
lyon_forecast['delta'] = lyon_forecast['ts'] - lyon_forecast['forecast_at']

lyon_forecast.tail()

lyon_forecast[(lyon_forecast.rain_3h >= 1) & (lyon_forecast.delta == '1H')].tail()

lyon_forecast[(lyon_forecast.ts >= '2017-09-14 14:00:00') & (lyon_forecast.delta == '1H')].head(15)

"""
Explanation: Forcast weather
End of explanation
"""


