!pip install tsfresh
!pip install lightgbm
!pip install --upgrade scipy

!git clone https://github.com/sdsc-bw/Predictive-Maintenance.git
!ls

import pandas as pd
import numpy as np
pd.set_option('display.max_columns', 500)
# Bibliothek speziell für die Visualisierung
import matplotlib.pyplot as plt
#plt.rcParams.keys()
plt.rcParams['ytick.labelsize'] = "large"
plt.rcParams['xtick.labelsize'] = "large"
plt.rcParams['axes.labelsize'] = "large"
plt.rcParams['figure.titlesize'] = "x-large"

from tsfresh.feature_extraction import EfficientFCParameters, MinimalFCParameters
from tsfresh.utilities.dataframe_functions import make_forecasting_frame
from tsfresh import extract_features, select_features

import warnings
warnings.filterwarnings("ignore")

from tqdm import tqdm
import lightgbm as lgb

from sklearn.metrics import mean_squared_error
from math import sqrt

# Die Eingabedateien haben keine Spaltennamen, deshalb werden hier die Spaltennamen definiert.
# Enthalten sind 26 Spalten: Motornummer, Zeitschritt, drei betriebsbereite Sensoreinstellungen und 21 Sensormessungen.

# Motornummer, Zeitschritt
index_columns_names =  ["UnitNumber","Cycle"]
# 3 betriebsbereite Sensoreinstellungen
operational_settings_columns_names = ["Operation Setting "+str(i) for i in range(1,4)]
# 21 Sensormessungen
sensor_measure_columns_names =["Sensor Measure"+str(i) for i in range(1,22)]

input_file_column_names = index_columns_names + operational_settings_columns_names + sensor_measure_columns_names

# Laden der Daten FD001 als Beispiel.
print("Laden der Daten beginnt!")
train_data = pd.read_csv('./Predictive-Maintenance/Turbofan RUL Prediction/data/train_FD001.txt', delim_whitespace=True, header=None)
train_data.columns = input_file_column_names
print("Laden der Daten ist fertig!")
print("Größe der Daten : Rows {}, Columns {}".format(train_data.shape[0], train_data.shape[1]))
train_data.head()

unit = train_data["UnitNumber"].unique()
print("In den Trainingsdaten sind {} Turbofan-Motoren enthalten.".format(len(unit)))

motor_id = 1
Dataframe_id = train_data[train_data["UnitNumber"]==motor_id]
print("Turbofan mit ID {} lief insgesamt {} Zyklen".format(motor_id, Dataframe_id.shape[0]))

motor_id = 2
Dataframe_id = train_data[train_data["UnitNumber"]==motor_id]
print("Turbofan mit ID {} lief insgesamt {} Zyklen".format(motor_id, Dataframe_id.shape[0]))

Zyklen_list = []
for ID in range(len(unit)):
    tmp = train_data[train_data["UnitNumber"]==ID+1]
    Zyklen_list.append(tmp.shape[0])
plt.hist(x=Zyklen_list, bins='auto',alpha=0.7, rwidth=0.85);
plt.xlabel('Lebenszyklus')
plt.ylabel('Frequency');

motor_id = 2
numberOfcolumns = len(sensor_measure_columns_names)
plt.figure(figsize=(15,2.5*len(sensor_measure_columns_names)))
Dataframe_id = train_data[train_data["UnitNumber"]==motor_id]
for i,v in enumerate(range(numberOfcolumns)):
    a = plt.subplot(numberOfcolumns,1,v+1)
    a.plot(Dataframe_id.index.values,Dataframe_id.iloc[:,v+5].values)
    a.title.set_text(sensor_measure_columns_names[v])
    plt.tight_layout()
print("       Zeitliche Verlauf von Turbofan mit id {}. X Axis is the Cycle Number, Y Axis is the Sensor Value ".format(motor_id))
plt.show()

def f1(col):
    # Option 1: Reverse the time evolution, where remaining time of a machine is 1 at the failure.
    return col[::-1]  
train_data['rul'] = train_data[['UnitNumber', 'Cycle']].groupby('UnitNumber').transform(f1)
plt.figure(figsize=(15,5))
motor_id = 1
df_of_id = train_data[train_data["UnitNumber"]==motor_id]
plt.plot(df_of_id.rul, label="RUL")
#plt.plot(df_of_id.Cycle, label="Cycle")
plt.legend()
plt.title("RUL von Turbofan mit ID {}".format(motor_id))
plt.show()

# Hier definieren wir die Fenstergröße mit 5
window_size = 5

engine_id = train_data["UnitNumber"].unique()
feature_x = pd.DataFrame()
RUL = pd.Series()


for index in tqdm(engine_id):
    
    data_frame = train_data[train_data["UnitNumber"]==index]
    add = pd.DataFrame(np.ones((1,data_frame.shape[1])), columns=data_frame.columns)
    data_frame = pd.concat([data_frame,add])
    data_frame.reset_index(inplace=True,drop=True)
    
    df_shift, _ = make_forecasting_frame(data_frame["Sensor Measure1"], kind="x", 
                                             max_timeshift = window_size, 
                                             rolling_direction=1)
    del df_shift["kind"]
    df_shift["Sensor Measure1"] = df_shift["value"]
    del df_shift["value"]
    for column in sensor_measure_columns_names[1:]:
        temp, _ = make_forecasting_frame(data_frame[column], kind="x", 
                                             max_timeshift = window_size, 
                                             rolling_direction=1)
        df_shift[column] = temp["value"]
    # Als Beispiel extrahieren wir nur einige grundlegende statistische Merkmale, wie z.B.: MEAN MAX MIN VARIANZ
    # Sie können die Einstellungen von default_fc_parameters ändern, um die zu extrahierenden Features auszuwählen.
    All_features = extract_features(df_shift, column_id="id", column_sort="time",
                                    default_fc_parameters=MinimalFCParameters(),
                                    impute_function=None, disable_progressbar=True,
                                    show_warnings=False, n_jobs=8)
    
    All_features=All_features.iloc[window_size-1:]
    y = data_frame["rul"].iloc[window_size-1:-1]
    
    feature_x = pd.concat([feature_x, All_features])
    RUL = pd.concat([RUL, y])

RUL.reset_index(drop=True, inplace=True)

def identify_and_remove_missing(Dataframe, Threshold=0.1):
    missing_series = pd.DataFrame( Dataframe.isnull().sum() / Dataframe.shape[0]).rename(columns = {'index': 'feature', 0: 'missing_fraction'})
    record_missing  = missing_series[missing_series["missing_fraction"]>Threshold]
    to_drop = list(record_missing.index)
    Dataframe = Dataframe.drop(columns = to_drop)
    return Dataframe
print("Größe vor der Verarbeitung ", feature_x.shape)
feature_x = identify_and_remove_missing(feature_x)
print("Größe nach der Verarbeitung ", feature_x.shape)

def identify_and_remove_unique_columns(Dataframe):
    unique_counts = Dataframe.nunique()
    record_single_unique = pd.DataFrame(unique_counts[unique_counts == 1]).reset_index().rename(columns = {'variable': 'feature', 0: 'nunique'})
    to_drop = list(record_single_unique['feature'])
    Dataframe = Dataframe.drop(columns = to_drop)
    return Dataframe
print("Größe vor der Verarbeitung ", feature_x.shape)
feature_x = identify_and_remove_unique_columns(feature_x)
print("Größe Nach der Verarbeitung", feature_x.shape)

selected_features = select_features(feature_x, RUL.values.reshape(-1))
print("Größe vor der Auswahl ", feature_x.shape)
feature_x = identify_and_remove_missing(feature_x)
print("Größe nach der Auswahl ", feature_x.shape)

def feature_sigmoid_transform(Featrues_DataFrame):
    """
    transform the operations into range(0,1)
    return transformed feature Frames , the median and IQR for each feature(columns) for new feature transformation
    this is excuted only when there are tsfresh fearues
    """

    # Interquartile range for each features
    Interquartile0_25 = np.array(Featrues_DataFrame.quantile(0.25))
    Interquartile0_75 = np.array(Featrues_DataFrame.quantile(0.75))
    IQR = Interquartile0_75 - Interquartile0_25

    # Median for each features
    median = np.array(Featrues_DataFrame.median())

    feature_columns = Featrues_DataFrame.columns
    # drop the features which has fast 0.0 Interquartile
    position = np.where( IQR == 0 )[0]

    Featrues_DataFrame.drop(Featrues_DataFrame.columns[position], axis=1, inplace = True)

    median = np.delete(median,position)
    IQR = np.delete(IQR,position)

    # sigmoid transformation
    transformed_features = (1+np.exp((-Featrues_DataFrame + median)/(IQR*1.2)))**(-1)

    return median, IQR, transformed_features

median, IQR, transformed_features = feature_sigmoid_transform(selected_features)

# Modellinitialisierung
model = lgb.LGBMRegressor(random_state=12)
# Modelltraining
model.fit(transformed_features, RUL)

prediction = model.predict(transformed_features)
plt.figure(figsize=(12,5))
plt.plot((prediction[:2000]), label="Prediction")
plt.plot((RUL[:2000]), label="Reale RUL")
plt.legend()
plt.show()

rms = sqrt(mean_squared_error(RUL.values.reshape(-1), prediction))
print("Der RMSE auf die Trainingsdaten lautet :", rms)

print("Laden der Testdaten beginnt!")
test_data = pd.read_csv('./Predictive-Maintenance/Turbofan RUL Prediction/data/test_FD001.txt', delim_whitespace=True, header=None)
test_data.columns = input_file_column_names
print("Laden der Daten ist fertig!")
print("Größe der Daten : Rows {}, Columns {}".format(test_data.shape[0], test_data.shape[1]))
test_data.head()

# Reale RUL
RUL_FD001 = pd.read_table("./Predictive-Maintenance/Turbofan RUL Prediction/data/RUL_FD001.txt", header=None, delim_whitespace=True)

engine_id = test_data["UnitNumber"].unique()
window_size = 5
feature_x = pd.DataFrame()
feature_columns = transformed_features.columns

for index in tqdm(engine_id):
    data_frame = test_data[test_data["UnitNumber"]==index]
    # Hier werden nur die Informationen zur letzten Fenstergröße verwendet. 
    data_frame = data_frame.iloc[-window_size:]
    data_frame = data_frame[sensor_measure_columns_names]
    data_frame["id"] = 1
    data_frame["time"] = range(data_frame.shape[0])
    data_frame.reset_index(inplace=True,drop=True)



    All_features = extract_features(data_frame, column_id="id", column_sort="time", #column_value="value", 
                                    default_fc_parameters=MinimalFCParameters(),
                                    impute_function=None, disable_progressbar=True,
                                    show_warnings=False, n_jobs=8)

    All_features = All_features[feature_columns]
    

    feature_x = pd.concat([feature_x, All_features])
    

# Merkmalstransformation
transformed_features = (1+np.exp((-feature_x + median)/(IQR*1.2)))**(-1)


prediction = model.predict(transformed_features)
prediction = np.round(prediction)

rms = sqrt(mean_squared_error(RUL_FD001.values.reshape(-1), prediction))
print("Der RMSE auf die Testdaten lautet :", rms)

feature_names = list(transformed_features.columns)
feature_importance_values = model.feature_importances_
# Notierung der Wichtigkeitkeint der Merkmale
feature_importances = pd.DataFrame({'feature': feature_names, 'importance': feature_importance_values})

# Sortierung der Merkmale nach ihrer Wichtigkeit
feature_importances = feature_importances.sort_values('importance', ascending = False).reset_index(drop = True)

# Normalisierung die Merkmals-Wichtigkeiten, um bis zum Wert Eins zu addieren
feature_importances['normalized_importance'] = feature_importances['importance'] / feature_importances['importance'].sum()
feature_importances['cumulative_importance'] = np.cumsum(feature_importances['normalized_importance'])

plot_n = 15
plt.figure(figsize = (10, 6))
ax = plt.subplot()
ax.barh(list(reversed(list(feature_importances["feature"][:plot_n]))), 
        feature_importances['normalized_importance'][:plot_n], 
        align = 'center', edgecolor = 'k')
ax.set_yticks(list(reversed(list(feature_importances["feature"][:plot_n]))))
ax.set_yticklabels(feature_importances['feature'][:plot_n], size = 12)

# Plot Labeling
plt.xlabel('Normalized Importance', size = 16); plt.title('Feature Importances', size = 18)
plt.show()