#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Fri Apr  3 00:03:51 2020

@author: assa
"""
#%%
from sklearn.cluster import KMeans
from scipy.spatial.distance import cdist
import numpy as np
import matplotlib.pyplot as plt
from sklearn.preprocessing import StandardScaler
import seaborn as sns; sns.set(style="ticks", color_codes=True)
import pandas as pd 
from sklearn.decomposition import PCA
#from __future__ import division, print_function   
import skfuzzy as fuzz
from sklearn.metrics import mean_squared_error,r2_score,mean_absolute_error
from sklearn.linear_model import LinearRegression
from sklearn.preprocessing import PolynomialFeatures
import statsmodels.api as sm
import seaborn as sns
from scipy import stats
#print(__doc__)
df = pd.read_csv('/home/assa/Documents/Assa/diplomka/fc5_105_all.csv')
#%% DATA
# =============================================================================
#sns.boxplot(x=df['SO2'])
emis_SO2=df[['SO2','WSPD10','TEMP2','PBL','PURB']]
poloha=df[['LAT','LON','HT']]
conc_SO2=df[['SO2_conc']]
data_0_so2=pd.DataFrame()
data_0_pm=pd.DataFrame()
data_0_so2_test=pd.DataFrame()

emis_SO2.plot.box()
sns.boxplot(data=df[['SO2','WSPD10','TEMP2','PBL','PURB','SO2_conc']])
plt.show()

sns.boxplot(data=df[['TEMP2','PBL']])
plt.show()

sns.boxplot(data=df[['PBL']])
plt.show()

sns.boxplot(data=df[['TEMP2']])
plt.show()

sns.boxplot(data=df[['SO2','WSPD10','PURB','SO2_conc']])
plt.show()


#%%Spliting and pre procesing CLUSTERING
# =============================================================================

z = np.abs(stats.zscore(conc_SO2))
print(z)
print(np.where(z > 3))

mask_so2=z>3
a=np.average(conc_SO2['SO2_conc'])
y=conc_SO2.values
y[mask_so2]=np.average(conc_SO2['SO2_conc'])
conc_SO2['SO2_conc']=y

standardscaler_SO2=StandardScaler().fit(emis_SO2.values)
s_f_so2 =standardscaler_SO2.transform(emis_SO2.values)


standardscaler_poloha=StandardScaler().fit(poloha.values)
s_f_poloha =standardscaler_poloha.transform(poloha.values)

standardscaler_conc_so2=StandardScaler().fit(conc_SO2.values)
s_f_conc_so2=standardscaler_conc_so2.transform(conc_SO2.values)

#PCA
pca = PCA(n_components=5)
pcaer_so2 = pca.fit(s_f_so2)
p_c_a_so2=pcaer_so2.transform(s_f_so2)


#%%Korelacia funkcia
# =============================================================================
def heatMap(corr,d,method, mirror=False):
   
   # Create Correlation df fmt="d")
   #corr = df.corr()
   # Plot figsize
   fig, ax = plt.subplots(figsize=(30, 30))
   # Generate Color Map
   colormap = sns.diverging_palette(220, 10, as_cmap=True)
   #sns.set(font_scale=3)
   if mirror == True:
      #Generate Heat Map, allow annotations and place floats in map
      sns.heatmap(corr, cmap=colormap, annot=True, fmt='.1f',vmin=-1, vmax=1,annot_kws={"size": 40})
      #Apply xticks
      plt.xticks(np.arange(len(corr.columns))+0.5, corr.columns[:-1],fontsize=30);
      #Apply yticks
      plt.yticks(np.arange(len(corr.columns))+0.5,corr.columns[1:],fontsize=30)
      #show plot
   else:
      # Drop self-correlations
      dropSelf = np.zeros_like(corr)
      dropSelf[np.triu_indices_from(dropSelf)] = True
      # Generate Color Map
      sns.set(font_scale=3)
      colormap = sns.diverging_palette(220, 10, as_cmap=True)
      # Generate Heat Map, allow annotations and place floats in map
      sns.heatmap(corr, cmap=colormap, annot=True, fmt='.1f', mask=dropSelf,vmin=-1, vmax=1,annot_kws={"size": 40})
      # Apply xticks
      plt.xticks(np.arange(len(corr.columns))+0.5, corr.columns[:-1],fontsize=10);
      # Apply yticks
      plt.yticks(np.arange(1,len(corr.columns))+0.5, corr.columns[1:],fontsize=10)
      #plt.colorbar(fontsize=30)
   # show plot
   plt.savefig('/home/assa/Documents/Assa/diplomka/reg_mod_11_jun/2/assa_cor_{0}_{1}.png'.format(d,method))
   plt.show()


#%% measure function
def q_accuracy_q(y_test,y_test_predict):
    rmse_test = np.sqrt(mean_squared_error(y_test, y_test_predict))
    r2_test = r2_score(y_test, y_test_predict)
    rmae_test = np.sqrt(mean_absolute_error(y_test, y_test_predict))
    return (rmse_test,r2_test,rmae_test)


#%% Cluster division 
# =============================================================================

X_train_1=np.concatenate((s_f_so2,s_f_poloha),axis=1)
#X_train_2=np.concatenate((s_f_so2,s_f_poloha,kmeans_so2),axis=1)   
#X_train_3=np.concatenate((s_f_so2,s_f_poloha,u.T),axis=1)  
y_train_so2=s_f_conc_so2

df_fc_culst=pd.DataFrame(np.concatenate((X_train_1,y_train_so2),axis=1),columns=['SO2','WSPD10','TEMP2','PBL','PURB','LAT','LON','HT','SO2_conc'])
plt.figure(figsize=(16,10))
sns.boxplot(data=df_fc_culst)
plt.xticks(fontsize=10);
 # Apply yticks
plt.yticks(fontsize=10)
plt.show()
# df_fc_culst=df_fc_culst.join(u_data_0_so2,rsuffix='fuzz')
# df_fc_culst_corr=df_fc_culst.corr()
# #heatMap(df_corr_km,k,'km')
# heatMap(df_fc_culst_corr,"m=11",'fc')
#   # transforms the existing feat
#%% Rgresion Aplication 
n_c=5
y_predicted_fc  =np.zeros((len(y_train_so2),n_c))

cluster_id=df['fuzzy_cmeans'].values

rmse_fc_for_k=[]

for k in range(n_c):
    mask_train_fc=cluster_id==k  
    y_fc_so2_tst=y_train_so2[mask_train_fc]
    #X_tst_so2_fc=X_test_1[mask_train_fc]
    #cluster_id

    # Linear Regresion
    #===================================
    y_fc=np.zeros((len(y_train_so2),1))
    poly_features = PolynomialFeatures(degree=2)
  
    X_train_k_fc = X_train_1[mask_train_fc]
    df_fc=pd.DataFrame(np.concatenate((X_train_k_fc,y_fc_so2_tst),axis=1),columns=['SO2','WSPD10','TEMP2','PBL','PURB','LAT','LON','HT','SO2_conc'])
    df_corr_fc=df_fc.corr()
    heatMap(df_corr_fc,k,'fc')
     
    plt.figure(figsize=(16,10))
    sns.boxplot(data=df_fc)
    plt.xticks(fontsize=10);
      # Apply yticks
    plt.yticks(fontsize=10)
    plt.savefig('/home/assa/Documents/Assa/diplomka/reg_mod_11_jun/2/boxplot_{0}.png'.format(k))
    plt.show()
  # transforms the existing features to higher degree features.  
    
#    help_fm= y_predicted_km[:,k] 
#    help_fc= y_predicted_fc[:,k] 
#    
#    help_fm[mask_test_km]= y_km[:,0]
#    help_fc[mask_test_fc]= y_fc[:,0]
#    
#    y_predicted_km[:,k] =help_fm
#    y_predicted_fc[:,k] =help_fc
    X_train_poly_fc = poly_features.fit_transform(X_train_k_fc)
    model_fc= sm.OLS(y_fc_so2_tst, X_train_poly_fc).fit()
    predictions = model_fc.predict(X_train_poly_fc)
    print("FC")
    print(model_fc.summary())
        
    kde='/home/assa/Documents/Assa/diplomka/reg_mod_11_jun/2/rano_fc_{0}.csv'.format(k)
    
    with open(kde, 'w') as f:
        f.writelines(model_fc.summary().as_csv())
    plt.figure(figsize=(8,5))
    p=plt.scatter(x=model_fc.fittedvalues,y=model_fc.resid,color='orange')
    xmin=min(model_fc.fittedvalues)
    xmax = max(model_fc.fittedvalues)
    plt.hlines(y=0,xmin=xmin*0.9,xmax=xmax*1.1,color='red',linestyle='--',lw=3)
    plt.xlabel("Fitted values",fontsize=10)
    plt.ylabel("Residuals",fontsize=10)
    plt.xticks(fontsize=10)
    plt.yticks(fontsize=10)
    #plt.title("Fitted vs. residuals plot",fontsize=18)
    #plt.grid(True)
    #plt.xlim(-2,2)
    plt.savefig('/home/assa/Documents/Assa/diplomka/reg_mod_11_jun/2/rezidua_{0}.png'.format(k))
    plt.show()
    
    plt.figure(figsize=(8,5))
    plt.hist(model_fc.resid_pearson,bins=20,color='orange')
    plt.ylabel('Count',fontsize=10)
    plt.xlabel('Normalized residuals',fontsize=10)
    #plt.title("Histogram of normalized residuals",fontsize=18)
    plt.savefig('/home/assa/Documents/Assa/diplomka/reg_mod_11_jun/2/histogram_{0}.png'.format(k))
    plt.show()
    y_fc[mask_train_fc,0]=model_fc.fittedvalues
    y_predicted_fc[:,k] =y_fc[:,0]
    
    print(k)
    rmse_test_fc = np.sqrt(mean_squared_error(y_fc_so2_tst, model_fc.fittedvalues))

    rmse_fc_for_k.append(rmse_test_fc)    
#
#%%
y_train_fc=np.sum(y_predicted_fc,axis=1 )
rmse_fc,r2_fc,rmae_fc=q_accuracy_q(y_train_so2,y_train_fc)
print("lin fuzzy \n",rmse_fc,r2_fc,rmae_fc)
print ('FC \n', rmse_fc_for_k )



#%%
mask_so2=np.abs(s_f_conc_so2)>1
so2_con_out=df[mask_so2]
so2_con_out.to_csv("promblemove_koncentrácie.csv")

mask_so2_1=np.abs(s_f_conc_so2)>1.5
so2_con_out_1=df[mask_so2_1]
so2_con_out_1.to_csv("velmi_promblemove_koncentrácie.csv")

mask_so2_2=np.abs(s_f_conc_so2)>2
so2_con_out_2=df[mask_so2_2]
so2_con_out_2.to_csv("extremne_promblemove_koncentrácie.csv")

mask_so2_3=np.abs(s_f_conc_so2)>5
so2_con_out_3=df[mask_so2_3]
so2_con_out_3.to_csv("velmi_extremne_promblemove_koncentrácie.csv")