#!/usr/bin/env python
# coding: utf-8

# In[ ]:


import dateutil
import math

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import matplotlib.dates as mdates
import plotly.plotly as py
import plotly.graph_objs as go
import plotly.figure_factory as ff
import statsmodels.api as sm

from scipy import stats
from sklearn.metrics import mean_squared_error, mean_absolute_error
from sklearn import linear_model
from sklearn.preprocessing import MinMaxScaler


# # Carga y limpieza de datos

# In[ ]:


df_start = pd.read_excel('../Data/Graficas DustTrack Honeywell 5sept comp.xlsx')
df_start.head()


# In[ ]:


# Eliminar nulos y headers repetidos

df_start = df_start[~df_start['Valor DustTrack'].isna()]
df_start = df_start[df_start.columns.values[:-1]]
sum(df_start['Valor DustTrack'] == 'Valor DustTrack')
df_start = df_start[np.logical_not(df_start['Valor DustTrack'] == 'Valor DustTrack')]

df_start.index = np.arange(len(df_start))


# In[ ]:


# Casting a variables numericas

df_start[['Valor DustTrack', 'Valor Honeywell']] = df_start[['Valor DustTrack', 'Valor Honeywell']].apply(pd.to_numeric)


# In[ ]:


def plot_time_series(df, variables, filename):    
    data = list()
    for variable in variables:
        data.append(go.Scatter(x = df.index,
                                 y = df[variable],
                                 name=variable))
        
    layout = go.Layout(title='DustTrack vs Honeywell',
                       yaxis=dict(title='Pm 2.5 value'),
                       xaxis=dict(title='N° Observation'))
    fig = go.Figure(data=data, layout=layout)
    plot_url = py.iplot(fig, filename=filename)
    return plot_url
    
# Login plotly
py.sign_in('FranciscoJavierOspinaSalazar', 'WLz5dfnwRjfo2cSwlRiL')


# # Visualización inicial de los datos y exploraciones univariadas

# In[ ]:


plot_time_series(df=df_start, variables=['Valor DustTrack', 'Valor Honeywell'], filename='Comparison_sensors')


# ## Boxplot de ambas mediciones

# In[ ]:


def build_boxplot(df, variables, filename):
    if isinstance(variables, str):
        variables = [variables]
    data = list()
    for variable in variables:
        data.append(go.Box(y = df[variable],
                           name=variable))
        
    layout = go.Layout(title='Boxplot',
                       yaxis=dict(title='Pm 2.5 value'))
    fig = go.Figure(data=data, layout=layout)
    plot_url = py.iplot(fig, filename=filename)
    return plot_url
    


# In[ ]:


build_boxplot(df=df_start, variables='Valor DustTrack', filename='Boxplot_DustTrack')


# In[ ]:


build_boxplot(df=df_start, variables='Valor Honeywell', filename='Boxplot_Honeywell')


# In[ ]:


build_boxplot(df=df_start, variables=['Valor DustTrack', 'Valor Honeywell'], filename='Boxplot_Comparison1')


# # Remoción de outliers

# In[ ]:


df = df_start.copy()
df = df[np.logical_not(np.logical_or(df['Valor Honeywell'] > 90, df['Valor DustTrack'] > 492))]


# # Suavizamiento de serie con medias móviles
# 
# La descripción de las médias móviles y su uso para suavizar picos en series de tiempo se encuentra en [http://www.expansion.com/diccionario-economico/media-movil.html](http://www.expansion.com/diccionario-economico/media-movil.html)

# In[ ]:


n_movil_mean = 20
movil_mean = df[['Valor DustTrack', 'Valor Honeywell']].rolling(n_movil_mean).mean()
movil_mean.columns = ['Mm_DustTrack', 'Mm_Honeywell']
movil_mean = movil_mean.loc[n_movil_mean - 1:]


# In[ ]:


## Reescalamiento de los datos con la función Logaritmo


# In[ ]:


movil_mean = movil_mean.assign(Log10_DustTrack=np.log10(movil_mean['Mm_DustTrack']))
movil_mean = movil_mean.assign(Log10_Honeywell=np.log10(movil_mean['Mm_Honeywell']))
movil_mean.index = np.arange(len(movil_mean))
movil_mean.head()


# In[ ]:


plot_time_series(df=movil_mean, variables=['Mm_DustTrack', 'Mm_Honeywell'], filename='serie_movil')


# Con el suavizamiento se suprimieron los picos que superaban 2000 en el eje y

# In[ ]:


plot_time_series(df=movil_mean, variables=['Log10_DustTrack', 'Log10_Honeywell'], filename='log10_serie_movil')


# In[ ]:


build_boxplot(df=movil_mean, variables=['Log10_DustTrack', 'Log10_Honeywell'], filename='Bloxplot_log10')


# ## Partición de la base de datos en traint test

# In[ ]:


size_train = int(len(movil_mean) * .7)
index_train = np.random.choice(np.arange(len(movil_mean)), size=size_train, replace=False)
index_test = np.setdiff1d(np.arange(len(movil_mean)), index_train)
df_train = movil_mean.iloc[index_train]
df_test = movil_mean.iloc[index_test]
df_train = df_train.sort_index()
df_test = df_test.sort_index()


# ## Escalamiento de los datos 

# In[ ]:


scaler_dt = MinMaxScaler()
scaler_dt.fit(df_train['Mm_DustTrack'].values.reshape(-1, 1))
scaled_DustTrack_train = scaler_dt.transform(df_train['Mm_DustTrack'].values.reshape(-1,1))
scaled_DustTrack_test = scaler_dt.transform(df_test['Mm_DustTrack'].values.reshape(-1,1))

scaler_hw = MinMaxScaler()
scaler_hw.fit(df_train['Mm_Honeywell'].values.reshape(-1, 1))
scaled_Honeywell_train = scaler_hw.transform(df_train['Mm_Honeywell'].values.reshape(-1,1))
scaled_Honeywell_test = scaler_hw.transform(df_test['Mm_Honeywell'].values.reshape(-1,1))

df_train = df_train.assign(scaled_DustTrack=scaled_DustTrack_train)
df_train = df_train.assign(scaled_Honeywell=scaled_Honeywell_train)
df_test = df_test.assign(scaled_DustTrack=scaled_DustTrack_test)
df_test = df_test.assign(scaled_Honeywell=scaled_Honeywell_test)
movil_mean = movil_mean.assign()

df_train.head()


# In[ ]:


build_boxplot(df=df_train, variables=['scaled_DustTrack', 'scaled_Honeywell'], filename='Bloxplot_log10')


# # Correlación lineal de Pearson 
# 
# [https://es.wikipedia.org/wiki/Coeficiente_de_correlaci%C3%B3n_de_Pearson](https://es.wikipedia.org/wiki/Coeficiente_de_correlaci%C3%B3n_de_Pearson)
# 
# Debido a que solo son dos variables si la correlación de pearson es alta se pueden hacer modelos que intenten generar la variable **Valor DustTrack** a partir de **Valor Honeywell**.

# In[ ]:


# Coeficinte de correlación sobre las variables transformadas por logaritmos de las
# series suavizadas por las media móviles
stats.pearsonr(df_train['Log10_DustTrack'], df_train['Log10_Honeywell'])[0]


# In[ ]:


stats.pearsonr(df_train['scaled_DustTrack'], df_train['scaled_Honeywell'])[0]


# Tanto en la gráfica de la serie de tiempo como en la correlación, estas son las mejores variables hasta el momento para el modelo. 

# # Construcción de modelos

# In[ ]:


model_min_max_scaler = linear_model.LinearRegression(fit_intercept=True)
model_log_scaler = linear_model.LinearRegression(fit_intercept=True)
model_min_max_scaler.fit(X=df_train['scaled_Honeywell'].values.reshape(-1,1), y=df_train['scaled_DustTrack'].values.reshape(-1,1), )
model_log_scaler.fit(X=df_train['Log10_Honeywell'].values.reshape(-1,1), y=df_train['Log10_DustTrack'].values.reshape(-1,1))


# ## Párametros de los modelos

# In[ ]:


print('Parameters model with log scaler. Coef = {}, Intercept = {}'.format(model_log_scaler.coef_[0][0], model_log_scaler.intercept_[0]))
print('Parameters model with min_max scaler. Coef = {}, Intercept = {}'.format(model_min_max_scaler.coef_[0][0], model_min_max_scaler.intercept_[0]))


# ## Append de predicciones a las series suavizadas

# In[ ]:


pred_min_max_scaler_train = scaler_dt.inverse_transform(model_min_max_scaler.predict(df_train['scaled_Honeywell'].values.reshape(-1,1)).reshape(-1,1))
pred_log_scaler_train = np.power(10, model_log_scaler.predict(df_train['Log10_Honeywell'].values.reshape(-1,1)).reshape(-1,1))

df_train = df_train.assign(pred_min_max_scaler=pred_min_max_scaler_train)
df_train = df_train.assign(pred_log_scaler=pred_log_scaler_train)

pred_min_max_scaler_test= scaler_dt.inverse_transform(model_min_max_scaler.predict(df_test['scaled_Honeywell'].values.reshape(-1,1)).reshape(-1,1))
pred_log_scaler_test = np.power(10, model_log_scaler.predict(df_test['Log10_Honeywell'].values.reshape(-1,1)).reshape(-1,1))

df_test = df_test.assign(pred_min_max_scaler=pred_min_max_scaler_test)
df_test = df_test.assign(pred_log_scaler=pred_log_scaler_test)

df_test.head()


# ## Append de predicciones a las series sin suaizar

# In[ ]:


df_start = df_start.assign(pred_log_scaler=np.power(10, model_log_scaler.predict(np.log10(df_start['Valor Honeywell']).values.reshape(-1,1)).reshape(-1,1)))

scaled_Honeywell = scaler_hw.transform(df_start['Valor Honeywell'].values.reshape(-1,1))
df_start = df_start.assign(scaled_Honeywell=scaled_Honeywell)
pred_min_max_scaler = scaler_dt.inverse_transform(model_min_max_scaler.predict(df_start['scaled_Honeywell'].values.reshape(-1,1)).reshape(-1,1))
df_start = df_start.assign(pred_min_max_scaler=pred_min_max_scaler)
df_start.head()


# # Plots del valor real vs la predicción del modelo (sobre serie suavizada / train)

# In[ ]:


plot_time_series(df=df_train, 
                 variables=['pred_min_max_scaler', 'pred_log_scaler', 'Mm_DustTrack'], 
                 filename='Predictions_01')


# # Plots del valor real vs la predicción del modelo (sobre serie suavizada / test)

# In[ ]:


plot_time_series(df=df_test, 
                 variables=['pred_min_max_scaler', 'pred_log_scaler', 'Mm_DustTrack'], 
                 filename='Predictions_02')


# In[ ]:


plot_time_series(df=df_start, 
                 variables=['pred_min_max_scaler', 'pred_log_scaler', 'Valor DustTrack'], 
                 filename='Predictions_03')


# # Errores de la predicción vs el valor real

# In[ ]:


dc_errors = dict()
for key, value in {'Train Data': df_train, 'Test Data': df_test, 'Original Data': df_start}.items():
    true_val = 'Mm_DustTrack'
    if key == 'Original Data':
        true_val = 'Valor DustTrack'
    dc_errors[key] = [mean_squared_error(value[true_val], value['pred_log_scaler']),
                      mean_squared_error(value[true_val], value['pred_min_max_scaler']),
                      mean_absolute_error(value[true_val], value['pred_log_scaler']),
                      mean_absolute_error(value[true_val], value['pred_min_max_scaler'])]
df_errors = pd.DataFrame(dc_errors)
df_errors.index = ['Log10 transform MSE', 'Min_max transform MSE', 'Log10 transform MAE', 'Min_max transform MAE']
df_errors[['Train Data', 'Test Data', 'Original Data']]


# Nota: Los valores de MAE y MSE corresponden al error absoluto medio y error cuadrático medio

# In[ ]:


def plot_density(df, y_tue, y_pred, group_labels, title, filename):
    hist_data = list()
    for y in y_pred:
        hist_data.append(df[y_tue] - df[y])

    group_labels = group_labels

    # Create distplot with custom bin_size
    fig = ff.create_distplot(hist_data, group_labels, show_hist=False)
    fig['layout'].update(title=title)
    plot_url = py.iplot(fig, filename=filename)
    return plot_url


# In[ ]:


plot_density(df=df_train, 
             y_tue='Mm_DustTrack', 
             y_pred=['pred_log_scaler', 'pred_min_max_scaler'], 
             group_labels=['Log scaler', 'Min_max scaler'], 
             title='Densidad de los residuales train', 
             filename='Residuals_density_01')


# In[ ]:


plot_density(df=df_test, 
             y_tue='Mm_DustTrack', 
             y_pred=['pred_log_scaler', 'pred_min_max_scaler'], 
             group_labels=['Log scaler', 'Min_max scaler'], 
             title='Densidad de los residuales test', 
             filename='Residuals_density_02')


# In[ ]:


plot_density(df=df_start, 
             y_tue='Valor DustTrack', 
             y_pred=['pred_log_scaler', 'pred_min_max_scaler'], 
             group_labels=['Log scaler', 'Min_max scaler'], 
             title='Densidad de los residuales test', 
             filename='Residuals_density_02')


# In[ ]:


df_train = df_train.assign(res_log_scaler=np.abs(df_train['Mm_DustTrack'] - df_train['pred_log_scaler']))
df_train = df_train.assign(res_min_max_scaler=np.abs(df_train['Mm_DustTrack'] - df_train['pred_min_max_scaler']))
df_test = df_test.assign(res_log_scaler=np.abs(df_test['Mm_DustTrack'] - df_test['pred_log_scaler']))
df_test = df_test.assign(res_min_max_scaler=np.abs(df_test['Mm_DustTrack'] - df_test['pred_min_max_scaler']))
df_start = df_start.assign(res_log_scaler=np.abs(df_start['Valor DustTrack'] - df_start['pred_log_scaler']))
df_start = df_start.assign(res_min_max_scaler=np.abs(df_start['Valor DustTrack'] - df_start['pred_min_max_scaler']))


# In[ ]:


step = 0.05
quantile_train = df_train[['res_log_scaler', 'res_min_max_scaler']].quantile(np.arange(0, 1 + step, step))
quantile_train.columns = [['res_log_scaler_train', 'res_min_max_scaler_train']]
quantile_test = df_test[['res_log_scaler', 'res_min_max_scaler']].quantile(np.arange(0, 1 + step, step))
quantile_test.columns = [['res_log_scaler_test', 'res_min_max_scaler_test']]
quantile_start = df_start[['res_log_scaler', 'res_min_max_scaler']].quantile(np.arange(0, 1 + step, step))
quantile_start.columns = [['res_log_scaler_original_data', 'res_min_max_scaler_original_data']]
df_quantile = pd.concat([quantile_train, quantile_test, quantile_start], axis=1)
df_quantile


# # Ecuación final de modelamiento (transformaión Log10)
# 
# La ecuación final que mejor se comporta para el modelamiento de los valores de **Valor DustTrack** a partir de los datos de **Valor Honeywell** es:
# 
# $$y = 10^{1.1675 * log_{10}(x) + 0.4173}$$
# 
# 
# Donde X es el valor de la variable **Valor Honeywell**  y y el valor de la variable **Valor DustTrack** m

# # Conclusiones finales
# 
# - Es posible hacer una reconstrucción no perfecta pero si bastante buena de los datos del sensor **DustTrack** con base en los datos del sensor **Honeywell**.

# In[ ]: