G_ECM056_ECM062_ARI082_ECM036_ARI021.py

import pandas as pd
import numpy as np
import seaborn as sns
import matplotlib.pyplot as plt
from sklearn.cluster import KMeans
from sklearn.decomposition import PCA
from sklearn.preprocessing import StandardScaler
from sklearn.model_selection import train_test_split
from sklearn.ensemble import RandomForestClassifier, RandomForestRegressor
from sklearn.metrics import classification_report, confusion_matrix, mean_squared_error, mean_absolute_error, r2_score

train_file_fd1 = 'data/train_FD001.txt'
train_file_fd2 = 'data/train_FD002.txt'
train_file_fd3 = 'data/train_FD003.txt'
train_file_fd4 = 'data/train_FD004.txt'
test_file_fd1 = 'data/test_FD001.txt'
test_file_fd2 = 'data/test_FD002.txt'
test_file_fd3 = 'data/test_FD003.txt'
test_file_fd4 = 'data/test_FD004.txt'
rul_file_fd1 = 'data/RUL_FD001.txt'
rul_file_fd2 = 'data/RUL_FD002.txt'
rul_file_fd3 = 'data/RUL_FD003.txt'
rul_file_fd4 = 'data/RUL_FD004.txt'

columns = ['unit_number', 'time_in_cycles'] + [f'operational_setting_{i}' for i in range(1, 4)] + [f'sensor_measurement_{i}' for i in range(1, 27)]

def read_train_data():
    train_df1 = pd.read_csv(train_file_fd1, sep=' ', header=None, names=columns).dropna(axis=1, how='all')
    train_df2 = pd.read_csv(train_file_fd2, sep=' ', header=None, names=columns).dropna(axis=1, how='all')
    train_df3 = pd.read_csv(train_file_fd3, sep=' ', header=None, names=columns).dropna(axis=1, how='all')
    train_df4 = pd.read_csv(train_file_fd4, sep=' ', header=None, names=columns).dropna(axis=1, how='all')

    final_train_df1 = train_df1.dropna(axis=1, how='all')
    final_train_df2 = train_df2.dropna(axis=1, how='all')
    final_train_df3 = train_df3.dropna(axis=1, how='all')
    final_train_df4 = train_df4.dropna(axis=1, how='all')
   
    final_train_df1.to_excel('Train_df1.xlsx', index=False)
    final_train_df2.to_excel('Train_df2.xlsx', index=False)
    final_train_df3.to_excel('Train_df3.xlsx', index=False)
    final_train_df4.to_excel('Train_df4.xlsx', index=False)

    return [final_train_df1, final_train_df2, final_train_df3,final_train_df4]

def read_test_data():
    test_df1 = pd.read_csv(test_file_fd1, sep=' ', header=None, names=columns).dropna(axis=1, how='all')
    test_df2 = pd.read_csv(test_file_fd2, sep=' ', header=None, names=columns).dropna(axis=1, how='all')
    test_df3 = pd.read_csv(test_file_fd3, sep=' ', header=None, names=columns).dropna(axis=1, how='all')
    test_df4 = pd.read_csv(test_file_fd4, sep=' ', header=None, names=columns).dropna(axis=1, how='all')

    final_test_df1 = test_df1.dropna(axis=1, how='all')
    final_test_df2 = test_df2.dropna(axis=1, how='all')
    final_test_df3 = test_df3.dropna(axis=1, how='all')
    final_test_df4 = test_df4.dropna(axis=1, how='all')

    print(final_test_df1)
    print(final_test_df2)
    print(final_test_df3)
    print(final_test_df4)

    final_test_df1.to_excel('Test_df1.xlsx', index=False)
    final_test_df2.to_excel('Test_df2.xlsx', index=False)
    final_test_df3.to_excel('Test_df3.xlsx', index=False)
    final_test_df4.to_excel('Test_df4.xlsx', index=False)

    return [final_test_df1, final_test_df2, final_test_df3,final_test_df4]


def read_rul_validation_data():
    rul_fd1 = pd.read_csv(rul_file_fd1, sep=' ', header=None).dropna(axis=1, how='all')
    rul_fd2 = pd.read_csv(rul_file_fd2, sep=' ', header=None).dropna(axis=1, how='all')
    rul_fd3 = pd.read_csv(rul_file_fd3, sep=' ', header=None).dropna(axis=1, how='all')
    rul_fd4 = pd.read_csv(rul_file_fd4, sep=' ', header=None).dropna(axis=1, how='all')
    return rul_fd1, rul_fd2, rul_fd3, rul_fd4


def staging_by_kmeans_clustering(dataset_no, train_df, no_of_stages):
    degradation_stage_map = {
        0: "Normal",
        1: "Slightly degraded",
        2: "Moderately degraded",
        3: "Critical",
        4: "Failure"        
    }
    staging_df = train_df.copy()
    sensor_cols = [col for col in staging_df.columns if 'sensor_measurement' in col]
    scaler = StandardScaler()
    scaled_sensor_data = scaler.fit_transform(staging_df[sensor_cols])

    kmeans = KMeans(n_clusters=no_of_stages, random_state=50)
    staging_df['degradation_stage'] = kmeans.fit_predict(scaled_sensor_data)

    staging_df['degradation_stage_name'] = staging_df['degradation_stage'].map(degradation_stage_map)

    print(staging_df)

    pca = PCA(n_components=2)
    pca_result = pca.fit_transform(scaled_sensor_data)

    plt.figure(figsize=(8,5))
    plt.scatter(pca_result[:, 0], pca_result[:, 1], c=staging_df['degradation_stage'], cmap='viridis', s=5)
    cbar = plt.colorbar(label='Degradation Stage')
    cbar.set_ticks([0, 1, 2, 3, 4])
    cbar.set_ticklabels([
        degradation_stage_map[i] for i in range(5)
    ])
    plt.xlabel('PCA Component 1')
    plt.ylabel('PCA Component 2')
    plt.title(f'Data Set No: {dataset_no}: KMeans Clustering')
    plt.grid(True)
    filename=f'data_set_no_{dataset_no}_kmeans.jpg'
    plt.savefig(filename)
    plt.show()

    return scaler, kmeans, staging_df


def train_stage_classifier_and_confusion_matrix(dataset_no, staging_df):
    clf_df = staging_df.copy()
    feature_cols = [col for col in clf_df.columns if 'operational_setting' in col or 'sensor_measurement' in col]

    X = clf_df[feature_cols]
    y = clf_df['degradation_stage']
    X_train, X_val, y_train, y_val = train_test_split(X, y, test_size=0.2, stratify=y, random_state=42)

    clf = RandomForestClassifier(n_estimators=100, class_weight='balanced', random_state=42)
    clf.fit(X_train, y_train)
    y_pred = clf.predict(X_val)
    n_classes=len(np.unique(y))

    print("Classification Report:"); 
    print(classification_report(y_val, y_pred))

    cm = confusion_matrix(y_val, y_pred, labels=list(range(n_classes)))
    plt.figure(figsize=(6,5))
    sns.heatmap(cm, annot=True, fmt='d', cmap='Blues', xticklabels=range(5), yticklabels=range(5))
    ticks = np.arange(n_classes)
    plt.xticks(ticks, ticks)
    plt.yticks(ticks, ticks)
    plt.xlabel('Validation Predicted') 
    plt.ylabel('Validation Actual')
    plt.title(f'Data Set No {dataset_no}: Confusion Matrix')
    filename=f'data_set_no_{dataset_no}_confusion_matrix.jpg'
    plt.savefig(filename)
    plt.show()

    return clf


def train_random_forest_rul(data_set_no, train_df):
    train_df = train_df.copy()
    train_df = train_df.reset_index(drop=True)
    train_df['RUL'] = train_df.groupby('unit_number')['time_in_cycles'].transform('max') - train_df['time_in_cycles']
    features = [col for col in train_df.columns if 'sensor_measurement' in col or 'operational_setting' in col]
    X = train_df[features]
    y = train_df['RUL']
    
    X_train, X_val, y_train, y_val = train_test_split(X, y, test_size=0.2, random_state=42)

    rf_model = RandomForestRegressor(n_estimators=100, random_state=42)
    rf_model.fit(X_train, y_train)
    y_pred = rf_model.predict(X_val)

    sorted_idx = np.argsort(y_val)
    plt.figure(figsize=(10,6))
    plt.plot(np.arange(len(y_val)), y_val.values[sorted_idx], label='Actual RUL', color='green')
    plt.plot(np.arange(len(y_pred)), y_pred[sorted_idx], label='Predicted RUL', color='red', linestyle='--')
    title=f'Random Forest - Actual vs Predicted RUL for data_set_no_{data_set_no} (Validation)'
    plt.title(title)
    plt.xlabel('Sample Index')
    plt.ylabel('RUL')
    plt.legend()
    plt.grid(True)
    filename=f'data_set_no_{data_set_no}_actual_vs_predicted_rul.jpg'
    plt.savefig(filename)
    plt.show()

    return rf_model


def print_model_accuracy(data_set_no, rul_true_values, pred_rul):

    mse = mean_squared_error(rul_true_values, pred_rul)
    rmse = np.sqrt(mse)
    mae = mean_absolute_error(rul_true_values, pred_rul)
    r2 = r2_score(rul_true_values, pred_rul)

    print(f" Model Accuracy Information:")
    print(f" Data Set NO: {data_set_no}: RUL Prediction Metrics")
    print(f"  R2 Score: {r2:.2f}")
    print(f"  RMSE: {rmse:.2f}")
    print(f"  MAE: {mae:.2f}")


def evaluate_rul_on_test(data_set_no, test_df, rul_df, rul_model):
    print(rul_df)
    print(test_df)  
    df = test_df.copy()
    feature_cols = [col for col in df.columns if 'operational_setting' in col or 'sensor_measurement' in col]

    df['unit_number'] = df['unit_number'].astype(int)
    first_cycles = df.groupby('unit_number').first().reset_index()
    X_test = first_cycles[feature_cols]

    pred_rul = rul_model.predict(X_test)
    rul_true_values = rul_df.iloc[:, 0].values

    min_len = min(len(pred_rul), len(rul_true_values))
    pred_rul = pred_rul[:min_len]
    rul_true_values = rul_true_values[:min_len]

    return data_set_no, rul_true_values, pred_rul



def compute_and_plot_risk_score(data_set_no, test_df, scaler, kmeans, clf, rul_model):
    df = test_df.copy()
    sensor_cols = [col for col in df.columns if 'sensor_measurement' in col]
    feature_cols = [col for col in df.columns if 'operational_setting' in col or 'sensor_measurement' in col]

    scaled_data = scaler.transform(df[sensor_cols])
    df['degradation_stage'] = kmeans.predict(scaled_data)

    X = df[feature_cols]
    X_full = df[feature_cols]

    failure_probs = clf.predict_proba(X)[:, 4]

    predicted_time_left = rul_model.predict(X_full)
    predicted_time_left = np.clip(predicted_time_left, 1, None)  
    raw_risk_scores = failure_probs * predicted_time_left

    min_score = raw_risk_scores.min()
    max_score = raw_risk_scores.max()
    normalized_risk_scores = (raw_risk_scores - min_score) / (max_score - min_score + 1e-6)

    plt.figure(figsize=(10,6))
    plt.plot(normalized_risk_scores, marker='x', color='blue', label='Min-Max Normalized Risk Score')

    high_risk_indices = []
    high_risk_scores = []
    for idx, score in enumerate(normalized_risk_scores):
        if score > 0.7:
            high_risk_indices.append(idx)
            high_risk_scores.append(score)
            print(f"High Risk at sample {idx} | Normalized Risk Score = {score:.3f}")

    plt.scatter(high_risk_indices, high_risk_scores, marker='o', color='red', s=50, label='High Risk (> 0.7)')

    plt.title(f'Data Set No: {data_set_no} - Min-Max Normalized Risk Score')
    plt.xlabel('Sample Index')
    plt.ylabel('Normalized Risk Score')
    plt.grid(True)
    plt.legend()
    filename=f'data_set_no_{data_set_no}_risk_score.jpg'
    plt.savefig(filename)
    plt.show()

    urgency_risk_scores = failure_probs / (predicted_time_left + 1e-6)

    plt.figure(figsize=(10,6))
    plt.plot(urgency_risk_scores, marker='x', color='blue', label='Urgency-Based Risk Score')

    print(f"\nRisk Alerts for Data Set No: {data_set_no} (Urgency-Based Inversion):")
    for idx, score in enumerate(urgency_risk_scores):
        if score > 0.03:  
            print(f"High Urgency Risk at sample {idx} | Urgency Risk Score = {score:.4f}")

    plt.title(f'Data Set No: {data_set_no} - Urgency-Based Risk Score')
    plt.xlabel('Sample Index')
    plt.ylabel('Urgency Risk Score')
    plt.grid(True)
    plt.legend()
    filename=f'data_set_no_{data_set_no}_urgency_rick_score.jpg'
    plt.savefig(filename)
    plt.show()


def main():

    print('start of the program')
    data_set_number = ['FD001', 'FD002', 'FD003', 'FD004']
    train_dfs = read_train_data()
    test_dfs = read_test_data()
    rul_dfs = read_rul_validation_data()

    for index, (train_df, test_df, rul_df) in enumerate(zip(train_dfs, test_dfs, rul_dfs), start=0):
        data_set_no = data_set_number[index]
        no_of_stages = 5
        print(f"Phase 1: Clustering & classification for: {data_set_no}")
        scaler, kmeans, staging_df = staging_by_kmeans_clustering(data_set_no, train_df, no_of_stages)

        print(f"Phase 2: Training stage classifier for: {data_set_no}")
        clf = train_stage_classifier_and_confusion_matrix(data_set_no, staging_df)

        print(f"Phase 3: Computing time-to-next failure prediction for:{data_set_no}")
        rf_model = train_random_forest_rul(data_set_no, train_df)

        print(f"Phase 4: Computing risk scores & plotting distributions for:{data_set_no}")
        data_set_no, rul_true_values, pred_rul = evaluate_rul_on_test(data_set_no, test_df, rul_df, rf_model)
        compute_and_plot_risk_score(data_set_no, test_df, scaler, kmeans, clf, rf_model)
        print(f" Printing Model Errors")
        print_model_accuracy(data_set_no, rul_true_values, pred_rul)



    print("FD001 + FD003")
    combine_train_1_3 = pd.concat([train_dfs[0], train_dfs[2]], ignore_index=True)
    combine_test_1_3 = pd.concat([test_dfs[0], test_dfs[2]], ignore_index=True)
    combine_rul_1_3 = pd.concat([rul_dfs[0], rul_dfs[2]], ignore_index=True)

    scaler_1_3, kmeans_1_3, staging_df_1_3 = staging_by_kmeans_clustering("FD001+FD003", combine_train_1_3, 5)
    clf_1_3 = train_stage_classifier_and_confusion_matrix("FD001+FD003", staging_df_1_3)
    rf_model_1_3 = train_random_forest_rul("FD001+FD003", combine_train_1_3)
    data_set_no, rul_true_values, pred_rul = evaluate_rul_on_test("FD001+FD003", combine_test_1_3, combine_rul_1_3, rf_model_1_3)
    compute_and_plot_risk_score("FD001+FD003", combine_test_1_3, scaler_1_3, kmeans_1_3, clf_1_3, rf_model_1_3)
    print_model_accuracy(data_set_no, rul_true_values, pred_rul)

    print("FD002 + FD004")
    combine_train_2_4 = pd.concat([train_dfs[1], train_dfs[3]], ignore_index=True)
    combine_test_2_4 = pd.concat([test_dfs[1], test_dfs[3]], ignore_index=True)
    combine_rul_2_4 = pd.concat([rul_dfs[1], rul_dfs[3]], ignore_index=True)

    scaler_2_4, kmeans_2_4, staging_df_2_4 = staging_by_kmeans_clustering("FD002+FD004", combine_train_2_4, 5)
    clf_2_4 = train_stage_classifier_and_confusion_matrix("FD002+FD004", staging_df_2_4)
    rf_model_2_4 = train_random_forest_rul("FD002+FD004", combine_train_2_4)
    data_set_no, rul_true_values, pred_rul = evaluate_rul_on_test("FD002+FD004", combine_test_2_4, combine_rul_2_4, rf_model_2_4)
    compute_and_plot_risk_score("FD002+FD004", combine_test_2_4, scaler_2_4, kmeans_2_4, clf_2_4, rf_model_2_4)
    print_model_accuracy(data_set_no, rul_true_values, pred_rul)


    print("FD001 + FD002 + FD003 + FD004")
    combine_train_1_2_3_4 = pd.concat([train_dfs[0], train_dfs[1], train_dfs[2], train_dfs[3]], ignore_index=True)
    combine_test_1_2_3_4 = pd.concat([test_dfs[0], test_dfs[1], test_dfs[2], test_dfs[3]], ignore_index=True)
    combine_rul_1_2_3_4 = pd.concat([rul_dfs[0], rul_dfs[1], rul_dfs[2], rul_dfs[3]], ignore_index=True)

    scaler_1_2_3_4, kmeans_1_2_3_4, staging_df_1_2_3_4 = staging_by_kmeans_clustering("FD001+FD002+FD003+FD004", combine_train_1_2_3_4, 5)
    clf_1_2_3_4 = train_stage_classifier_and_confusion_matrix("FD001+FD002+FD003+FD004", staging_df_1_2_3_4)
    rf_model_1_2_3_4 = train_random_forest_rul("FD001+FD002+FD003+FD004", combine_train_1_2_3_4)
    data_set_no, rul_true_values, pred_rul = evaluate_rul_on_test("FD001+FD002+FD003+FD004", combine_test_1_2_3_4, combine_rul_1_2_3_4, rf_model_1_2_3_4)
    compute_and_plot_risk_score("FD001+FD002+FD003+FD004", combine_test_1_2_3_4, scaler_1_2_3_4, kmeans_1_2_3_4, clf_1_2_3_4, rf_model_1_2_3_4)
    print_model_accuracy(data_set_no, rul_true_values, pred_rul)



    print('End of the program')

if __name__ == '__main__':
    main()