# Import necessary libraries for data manipulation, numerical operations, plotting, and machine learning
import pandas as pd # For data manipulation and analysis (DataFrames)
import numpy as np # For numerical operations, especially array handling
import matplotlib.pyplot as plt # For creating static, interactive, and animated visualizations
from sklearn.cluster import KMeans # K-Means clustering algorithm
# Import preprocessing tools and metrics from scikit-learn
from sklearn.preprocessing import StandardScaler, LabelEncoder # StandardScaler is used for scaling (RobustScaler was used in another version)
from sklearn.metrics import silhouette_score, precision_score, recall_score, f1_score, roc_auc_score, roc_curve, confusion_matrix, classification_report # Evaluation metrics
import argparse # For parsing command line arguments
import os # For interacting with the operating system, like manipulating file paths
import seaborn as sns # For drawing attractive statistical graphics (used for confusion matrix heatmap)
# --- Command line arguments setup ---
# This section defines and parses command-line arguments to control script behavior
parser = argparse.ArgumentParser(description='Anomaly detection using K-Means clustering with visualization.')
# Define arguments with their expected data type, default value, and a help message
parser.add_argument('--timesteps', type=int, default=30, help='Number of timesteps for sequences.') # Length of time window for each sequence
parser.add_argument('--n_clusters', type=int, default=5, help='Number of clusters for K-Means (should match the number of failure types + normal).') # Number of clusters
parser.add_argument('--n_init', type=int, default=10, help='Number of initializations for K-Means.') # Number of times K-Means algorithm will be run with different centroid seeds
parser.add_argument('--transition', action='store_true', help='Use transition data for testing.') # Flag: if present, use transition test data files
# Plotting flags: action='store_true' means the variable becomes True if the flag is present
parser.add_argument('--plot_raw', action='store_true', help='Plot raw data.')
parser.add_argument('--plot_clustered', action='store_true', help='Plot clustered data.')
parser.add_argument('--plot_anomalies', action='store_true', help='Plot detected anomalies (based on clusters).')
parser.add_argument('--plot_misclassified', action='store_true', help='Plot misclassified instances (based on clusters).')
# Parse the arguments provided by the user when running the script
options = parser.parse_args()
# Assign the parsed arguments to variables for easier access throughout the script
n_clusters = options.n_clusters
timesteps = options.timesteps
n_init = options.n_init
#####################################################################################################
# --- Data File Configuration ---
# This section defines the list of data files to be used for training and testing
#####################################################################################################
# Specify the number of failure types data is available for (excluding the normal state, which is class 0)
NumberOfFailures = 4 # So far, we have only data for the first 4 types of failures
# Initialize a nested list to store filenames for training and testing data
# datafiles[0] will store training files, datafiles[1] will store testing files
# Each inner list datafiles[train/test][i] corresponds to class i (0 for Normal, 1 to NumberOfFailures for failure types)
datafiles = [[], []] # datafiles[0] for train data filenames, datafiles[1] for test data filenames
# Populate the inner lists. We need NumberOfFailures + 1 inner lists for classes (0 to 4).
for i in range(NumberOfFailures + 1):
datafiles[0].append([]) # Add an empty list for the current class's training files
datafiles[1].append([]) # Add an empty list for the current class's testing files
# Manually assign the base filenames for each class for the training set
# These filenames are expected to be in a 'data' subdirectory relative to the script
datafiles[0][0] = ['2024-08-07_5_', '2024-08-08_5_', '2025-01-25_5_', '2025-01-26_5_'] # Normal training data files (Class 0)
datafiles[0][1] = ['2024-12-11_5_', '2024-12-12_5_', '2024-12-13_5_'] # Failure Type 1 training data files (Class 1)
datafiles[0][2] = ['2024-12-18_5_', '2024-12-21_5_', '2024-12-22_5_', '2024-12-23_5_', '2024-12-24_5_'] # Failure Type 2 training data files (Class 2)
datafiles[0][3] = ['2024-12-28_5_', '2024-12-29_5_', '2024-12-30_5_'] # Failure Type 3 training data files (Class 3)
datafiles[0][4] = ['2025-02-13_5_', '2025-02-14_5_'] # Failure Type 4 training data files (Class 4)
# Assign base filenames for the testing set based on the --transition flag
if options.transition:
# Use test files that are specified as including transition periods
datafiles[1][0] = ['2025-01-27_5_', '2025-01-28_5_'] # Normal test data files (Class 0)
datafiles[1][1] = ['2024-12-14_5_', '2024-12-15_5_', '2024-12-16_5_'] # Failure Type 1 test data files (with TRANSITION) (Class 1)
datafiles[1][2] = ['2024-12-17_5_', '2024-12-19_5_', '2024-12-25_5_', '2024-12-26_5_'] # Failure Type 2 test data files (with TRANSITION) (Class 2)
datafiles[1][3] = ['2024-12-27_5_', '2024-12-31_5_', '2025-01-01_5_'] # Failure Type 3 test data files (with TRANSITION) (Class 3)
datafiles[1][4] = ['2025-02-12_5_', '2025-02-15_5_', '2025-02-16_5_'] # Failure Type 4 test data files (Class 4)
else:
# Use test files that are specified as not including transition periods
datafiles[1][0] = ['2025-01-27_5_', '2025-01-28_5_'] # Normal test data files (Class 0)
datafiles[1][1] = ['2024-12-14_5_', '2024-12-15_5_'] # Failure Type 1 test data files (Class 1)
datafiles[1][2] = ['2024-12-19_5_', '2024-12-25_5_', '2024-12-26_5_'] # Failure Type 2 test data files (Class 2)
datafiles[1][3] = ['2024-12-31_5_', '2025-01-01_5_'] # Failure Type 3 test data files (Class 3)
datafiles[1][4] = ['2025-02-15_5_', '2025-02-16_5_'] # Failure Type 4 test data files (Class 4)
# --- Feature Definition ---
# Define the list of features (column names) to be extracted from the data files
features = ['r1 s1', 'r1 s4', 'r1 s5']
# Store the original number of features. This is needed later for sequence creation.
n_original_features = len(features) # Store the original number of features
# Dictionaries to store display names (using LaTeX-like format) and units for features, mainly used for plot labels
featureNames = {}
featureNames['r1 s1'] = r'T\_\{evap\}' # Evaporator Temperature
featureNames['r1 s4'] = r'T\_\{cond\}' # Condenser Temperature
featureNames['r1 s5'] = r'T\_\{air\}' # Air Temperature
unitNames = {}
unitNames['r1 s1'] = r'($^o$C)' # Degrees Celsius unit
unitNames['r1 s4'] = r'($^o$C)'
unitNames['r1 s5'] = r'($^o$C)'
# Redundant variable, but kept from original code
NumFeatures = len(features)
#####################################################################################################
# --- Data Loading and Preprocessing (Training Data) ---
# This section loads, cleans, and preprocesses the training data
#####################################################################################################
# List to hold processed DataFrames for training data, organized by class
# Each element dataTrain[i] is a DataFrame containing concatenated data for class i
dataTrain = []
# Loop through each list of filenames for each training class (Class 0, 1, 2, 3, 4)
for class_files in datafiles[0]:
class_dfs = [] # Temporary list to hold DataFrames loaded for the current class
# Loop through each base filename in the current class's list
for base_filename in class_files:
# Construct the full absolute file path assuming data is in a 'data' subdirectory
script_dir = os.path.dirname(os.path.abspath(__file__)) # Get the directory where the current script is located
data_dir = os.path.join(script_dir, 'data') # Define the path to the 'data' subdirectory
filepath = os.path.join(data_dir, f'{base_filename}.csv') # Combine directory and filename
try:
# Read the CSV file into a pandas DataFrame
df = pd.read_csv(filepath)
# Convert the 'datetime' column to datetime objects. Try two formats and coerce errors (set invalid parsing to NaT).
df['timestamp'] = pd.to_datetime(df['datetime'], format='%m/%d/%Y %H:%M', errors='coerce')
# Fill any NaT values from the first attempt by trying a second format.
df['timestamp'] = df['timestamp'].fillna(pd.to_datetime(df['datetime'], format='%d-%m-%Y %H:%M:%S', errors='coerce'))
# Convert the specified feature columns to numeric type. Coerce errors (set invalid parsing to NaN).
for col in features:
df[col] = pd.to_numeric(df[col], errors='coerce')
# Set the 'timestamp' column as the DataFrame index
# Resample the data to a 5-minute frequency ('5Min').
# Select only the specified 'features' columns.
# Calculate the mean for each 5-minute interval.
df = df.set_index('timestamp').resample('5Min')[features].mean() # Set index, resample, and calculate mean for features
# Estimate any remaining missing values (NaN) within the feature columns using linear interpolation
df = df[features].interpolate() # Estimate missing values using linear interpolation
# Append the processed DataFrame to the list for the current class
class_dfs.append(df)
except FileNotFoundError:
# If a file is not found, print a warning and skip processing it
print(f"Warning: File {filepath} not found and skipped.")
# If the list of DataFrames for the current class is not empty, concatenate them into a single DataFrame
if class_dfs:
dataTrain.append(pd.concat(class_dfs))
# Concatenate all DataFrames from all training classes into a single large DataFrame for training the scaler and the model
combined_train_data = pd.concat(dataTrain)
#####################################################################################################
# --- Data Loading and Preprocessing (Test Data) ---
# This section loads, cleans, and preprocesses the testing data
# The process is identical to the training data loading, but done separately for test files
#####################################################################################################
# List to hold processed DataFrames for test data, organized by class
# Each element dataTest[i] is a DataFrame containing concatenated data for class i
dataTest = []
# Loop through each list of filenames for each test class (Class 0, 1, 2, 3, 4)
for class_files in datafiles[1]:
class_dfs = [] # Temporary list to hold DataFrames loaded for the current class
# Loop through each base filename in the current class's list
for base_filename in class_files:
# Construct the full absolute file path
script_dir = os.path.dirname(os.path.abspath(__file__))
data_dir = os.path.join(script_dir, 'data')
filepath = os.path.join(data_dir, f'{base_filename}.csv')
try:
# Read the CSV file into a pandas DataFrame
df = pd.read_csv(filepath)
# Convert the 'datetime' column to datetime objects
df['timestamp'] = pd.to_datetime(df['datetime'], format='%m/%d/%Y %H:%M', errors='coerce')
df['timestamp'] = df['timestamp'].fillna(pd.to_datetime(df['datetime'], format='%d-%m-%Y %H:%M:%S', errors='coerce'))
# Convert the specified feature columns to numeric type
for col in features:
df[col] = pd.to_numeric(df[col], errors='coerce')
# Set the 'timestamp' column as the index and resample to 5-minute frequency, calculating the mean
df = df.set_index('timestamp').resample('5Min')[features].mean() # Set index, resample, and calculate mean for features
# Estimate any remaining missing values (NaN) using linear interpolation
df = df[features].interpolate() # Estimate missing values using linear interpolation
# Append the processed DataFrame to the list for the current class
class_dfs.append(df)
except FileNotFoundError:
# If a file is not found, print a warning and skip processing it
print(f"Warning: File {filepath} not found and skipped.")
# If the list of DataFrames for the current class is not empty, concatenate them into a single DataFrame
if class_dfs:
dataTest.append(pd.concat(class_dfs))
#####################################################################################################
# --- Raw Data Plotting (Optional) ---
# This section plots the unprocessed data for visualization if the --plot_raw flag is set
#####################################################################################################
# Check if the plot_raw argument was set to True
if options.plot_raw:
num_features = len(features) # Get the number of features to determine plot layout
# Create a figure and a grid of subplots. Share the x-axis among all subplots.
fig, axes = plt.subplots(num_features, 1, figsize=(15, 5 * num_features), sharex=True)
# Ensure 'axes' is always an array, even if there's only one feature (and thus only one subplot)
if num_features == 1:
axes = [axes]
# Loop through each feature by index (i) and name (feature)
for i, feature in enumerate(features):
# Loop through each test data DataFrame (representing a different class)
for k, df in enumerate(dataTest):
# Plot the data for the current feature and class over time
# Use f-string to include the class number in the label
axes[i].plot(df.index, df[feature], label=f'Class {k}')
# Set the label for the y-axis using the feature's display name and unit
axes[i].set_ylabel(f'{featureNames[feature]} {unitNames[feature]}')
# Set the title for the subplot using the feature's display name
axes[i].set_title(featureNames[feature])
# Add a legend to identify the classes being plotted in each subplot
# Placing it on the last subplot is common to avoid repetition
axes[num_features - 1].legend(loc='upper right') # Place legend on the last subplot
# Adjust layout to prevent plot elements (like labels) from overlapping
plt.tight_layout()
# Display the plot window
plt.show()
# exit(0) # Uncomment this line if you want the script to stop after showing raw data plots
########################################################################################################
# --- Data Scaling ---
# This section scales the data using StandardScaler
########################################################################################################
# Initialize the StandardScaler. This scaler standardizes features by removing the mean and scaling to unit variance.
scaler = StandardScaler() # Using StandardScaler in this version
# Fit the scaler on the training data and then transform the training data.
# The scaler learns the mean and standard deviation from the combined_train_data[features].
scaled_train_data = scaler.fit_transform(combined_train_data[features]) # Fit on training, transform training
# Transform each test data DataFrame using the scaler fitted on the training data.
# The same scaling parameters (mean, standard deviation) from the training data are applied to the test data.
# A list comprehension efficiently applies the transformation to each DataFrame in dataTest.
scaled_test_data_list = []
for df in dataTest: # Loop through each test DataFrame
scaled_test_data_list.append(scaler.transform(df[features])) # Transform each test DataFrame
# Convert the scaled NumPy arrays back into pandas DataFrames.
# This step is optional but can be helpful for inspection and keeping track of timestamps/column names.
scaled_train_df = pd.DataFrame(scaled_train_data, columns=features, index=combined_train_data.index)
# Create a list of scaled test DataFrames, maintaining original indices and column names.
scaled_test_df_list = [pd.DataFrame(data, columns=features, index=df.index) for data, df in zip(scaled_test_data_list, dataTest)]
############################################################################################################
# --- Sequence Creation with Rate of Change Feature Engineering ---
# This section defines a function to create time sequences and adds rate of change features
############################################################################################################
# Function to create time sequences (windows) from data and append rate of change features
# 'data': Input NumPy array (scaled feature values)
# 'timesteps': The length of each sequence (number of time points)
# 'original_features': The number of features in the input data (used for padding) - NOTE: Parameter name is potentially misleading, should be count
def create_sequences_with_rate_of_change(data, timesteps, original_features): # NOTE: original_features parameter likely intended as original_features_count
sequences = [] # List to store the generated sequences
# Iterate through the data to create overlapping sequences.
# The loop runs until the last possible start index for a sequence of length 'timesteps'.
for i in range(len(data) - timesteps + 1):
# Extract a sequence (slice) of 'timesteps' length starting from index 'i'.
sequence = data[i:i + timesteps]
# Calculate the difference between consecutive time points within the sequence.
# np.diff with axis=0 calculates the difference along the rows (time).
# This results in an array with shape (timesteps - 1, number_of_features).
rate_of_change = np.diff(sequence[:timesteps], axis=0)
# Pad the rate of change array to match the original sequence's length ('timesteps').
# np.diff reduces the dimension by 1, so we add a row of zeros at the beginning.
# The padding shape should be (1 row, number of columns equal to original features count).
padding = np.zeros((1, original_features)) # NOTE: This line caused a TypeError in previous debugging if original_features was not the count. It should likely use n_original_features or a correctly passed count.
# Vertically stack the padding row on top of the rate of change array.
# The result has shape (timesteps, number_of_features).
rate_of_change_padded = np.vstack((padding, rate_of_change)) # Stack padding on top of diff
# Horizontally stack the original sequence and the padded rate of change array.
# The resulting combined sequence has shape (timesteps, 2 * number_of_features).
sequences.append(np.hstack((sequence, rate_of_change_padded))) # Concatenate original sequence and padded rate of change
# Convert the list of 3D sequences into a single 3D NumPy array.
return np.array(sequences)
# Create time sequences with rate of change for the scaled training data.
# The output is a 3D array: (number of sequences, timesteps, 2 * n_original_features).
X_train_sequences = create_sequences_with_rate_of_change(scaled_train_df.values, timesteps, n_original_features) # Pass n_original_features (correctly passes count)
# Create time sequences with rate of change for each scaled test data DataFrame.
# This results in a list where each element is a 3D array of sequences for a specific test class.
X_test_sequences_list = [create_sequences_with_rate_of_change(df.values, timesteps, n_original_features) for df in scaled_test_df_list] # Pass n_original_features (correctly passes count)
############################################################################################################
# --- K-Means Clustering Model ---
# This section initializes, trains, and applies the K-Means model
############################################################################################################
# Reshape the training sequences into a 2D array for the K-Means algorithm.
# K-Means expects data in the shape (number of samples, number of features).
# Each sequence (timesteps * total_features) is flattened into a single row for clustering.
n_samples, n_timesteps, n_total_features = X_train_sequences.shape # Get dimensions of the sequence array
X_train_reshaped = X_train_sequences.reshape(n_samples, n_timesteps * n_total_features) # Flatten each sequence
# Initialize the KMeans model.
# n_clusters: The desired number of clusters (set by command line argument).
# random_state=42: Sets the seed for random number generation for initial centroids, ensuring reproducibility.
# n_init=10: Runs the K-Means algorithm 10 times with different centroid initializations and selects the best result based on inertia.
kmeans = KMeans(n_clusters=n_clusters, random_state=42, n_init=n_init) # Initialize KMeans model
# Train (fit) the KMeans model on the reshaped training data.
kmeans.fit(X_train_reshaped)
############################################################################################################################
# --- Predict Clusters for Test Data ---
# This section applies the trained K-Means model to the test data to get cluster assignments
############################################################################################################################
# List to store the predicted cluster labels for each test data class.
cluster_labels_test_list = []
# List to store the reshaped test data arrays (flattened sequences), needed for later evaluation metrics.
X_test_reshaped_list = []
kmeans_models = [] # This variable was declared but not used in the original code (To store kmeans model for each test set)
# Loop through each test data sequence array (one for each test class).
for i, X_test_seq in enumerate(X_test_sequences_list):
# Get the dimensions of the current test sequence array.
n_samples_test, n_timesteps_test, n_total_features_test = X_test_seq.shape
# Reshape the test sequences for prediction, flattening each sequence into a single data point for KMeans.
X_test_reshaped = X_test_seq.reshape(n_samples_test, n_timesteps_test * n_total_features_test)
# Use the trained K-Means model to predict the cluster label for each reshaped test sequence.
labels = kmeans.predict(X_test_reshaped)
# Append the predicted labels for the current test class to the list.
cluster_labels_test_list.append(labels)
# Append the reshaped test data for the current class to the list (needed for Silhouette score calculation later).
X_test_reshaped_list.append(X_test_reshaped) # Append reshaped data
# kmeans_models.append(kmeans) # Append the trained kmeans model (Variable not used)
############################################################################################################################
# --- Plotting Clustered Data (Optional) ---
# This function plots the original data points, colored according to their predicted cluster
############################################################################################################
# Function to plot the original feature data, with points colored based on their assigned cluster ID.
# 'original_data_list': List of original (unscaled) test DataFrames, one per class.
# 'cluster_labels_list': List of predicted cluster label arrays, one for each corresponding test DataFrame.
# 'n_clusters': Total number of clusters used.
# 'features', 'featureNames', 'unitNames': Dictionaries for plotting labels.
def plot_clustered_data(original_data_list, cluster_labels_list, n_clusters, features, featureNames, unitNames):
num_features = len(features) # Number of features to plot (determines number of subplots)
# Create a figure and a set of subplots (one row, num_features columns). Share the x-axis.
fig, axes = plt.subplots(num_features, 1, figsize=(15, 5 * num_features), sharex=True)
# Ensure 'axes' is always an array, even if there's only one feature (and thus only one subplot)
if num_features == 1:
axes = [axes]
# Generate a color map to assign distinct colors to each cluster ID
colors = plt.cm.viridis(np.linspace(0, 1, n_clusters)) # Assign colors to each cluster
# Loop through each original test data DataFrame and its corresponding cluster labels
for k, df in enumerate(original_data_list): # k is the index of the test class/DataFrame
original_indices = df.index # Get the time index from the original DataFrame
# The cluster labels correspond to sequences, which start 'timesteps' points later than the raw data.
# Get the time index corresponding to the end of each sequence for plotting.
time_index = original_indices[timesteps - 1:]
# Loop through each original feature to plot it
for i, feature in enumerate(features): # i is the index of the feature
# Loop through each possible cluster ID
for cluster_id in range(n_clusters):
# Find the indices within the current test data corresponding to the current cluster ID
cluster_indices_kmeans = np.where(cluster_labels_list[k] == cluster_id)[0]
# If there are any data points assigned to this cluster
if len(cluster_indices_kmeans) > 0:
# Plot the original feature data points for this specific cluster ID.
# x-axis: time_index points corresponding to the sequence end
# y-axis: original feature values at those time_index points
# color: color assigned to the cluster
# label: label for the cluster (only show for the first class (k==0) to avoid redundant legends)
# s=10: size of the scatter markers
axes[i].scatter(time_index[cluster_indices_kmeans], df[feature].loc[time_index[cluster_indices_kmeans]], color=colors[cluster_id], label=f'Cluster {cluster_id}' if k == 0 else "", s=10)
# Set the y-axis label and title for the current feature's subplot
axes[i].set_ylabel(f'{featureNames[feature]} {unitNames[feature]}')
axes[i].set_title(featureNames[feature])
# Add a legend to the plot. Place it on the last subplot.
axes[num_features - 1].legend(loc='upper right') # Place legend on the last subplot
# Adjust layout to prevent plot elements (like labels) from overlapping
plt.tight_layout()
# Display the plot window
plt.show()
# Call the plotting function if the --plot_clustered command line flag is provided
if options.plot_clustered:
plot_clustered_data(dataTest, cluster_labels_test_list, n_clusters, features, featureNames, unitNames)
#####################################################################################################
# --- Evaluation and plotting of anomalies and misclassified instances (based on cluster labels) ---
# This section evaluates clustering performance using classification metrics and plots specific instances
#####################################################################################################
# Function to evaluate clustering results and plot anomalies/misclassified instances
# 'kmeans_model': The trained K-Means model.
# 'scaled_test_data_list': List of scaled test data DataFrames.
# 'n_clusters': Number of clusters.
# 'original_test_data_list': List of original (unscaled) test data DataFrames.
# 'true_labels_list': List of arrays containing true class labels for test data.
# 'features', 'featureNames', 'unitNames': Feature information for plotting.
# 'plot_anomalies', 'plot_misclassified': Boolean flags from command line arguments.
def evaluate_and_plot_anomalies(kmeans_model, scaled_test_data_list, n_clusters, original_test_data_list, true_labels_list, features, featureNames, unitNames, plot_anomalies=False, plot_misclassified=False):
# Lists to accumulate true labels, predicted cluster labels, and original sequences for ALL test data
all_y_true_categorical = [] # Stores the true class label (0, 1, etc.) for each sequence across all test data
all_predicted_cluster_labels = [] # Stores the predicted cluster ID for each sequence across all test data
all_original_test_sequences = [] # Stores the original feature values for each sequence (window) across all test data, used for plotting
# Lists to store evaluation metrics calculated per test class DataFrame
inertia_values = [] # Inertia score when the model predicts on each class's data (Note: this likely appends the model's overall inertia repeatedly)
silhouette_scores = [] # Silhouette score calculated for each class's data based on predicted clusters
# Loop through each test data class (scaled data, original data, and true labels)
for i, (scaled_test_df, original_test_df, y_true_categorical) in enumerate(zip(scaled_test_data_list, original_test_data_list, true_labels_list)): # i is the class index
# Create sequences with rate of change for the current scaled test data DataFrame
X_test_sequences = create_sequences_with_rate_of_change(scaled_test_df.values, timesteps, n_original_features) # Create sequences for the current class
# If no sequences could be generated for this class (e.g., data too short), print warning and skip
if X_test_sequences.size == 0:
print(f"Warning: No test sequences generated for class {i}. Skipping evaluation for this class.")
continue # Skip to the next iteration (next class)
# Get the number of sequences generated for the current class
n_samples_test = X_test_sequences.shape[0]
# Reshape these sequences for prediction by the K-Means model (flatten each sequence)
X_test_reshaped = X_test_sequences.reshape(n_samples_test, -1)
# Predict the cluster label for each reshaped sequence using the trained model
cluster_labels_predicted = kmeans_model.predict(X_test_reshaped)
# Append the model's inertia (this seems like a repeated value of the final training inertia)
inertia_values.append(kmeans_model.inertia_) # Check if this is the intended calculation for per-class evaluation
# Calculate the Silhouette score for the current class's data based on its predicted clusters
# Only calculate if there's more than one unique predicted label and at least one sample
if len(np.unique(cluster_labels_predicted)) > 1 and len(cluster_labels_predicted) > 0:
silhouette_scores.append(silhouette_score(X_test_reshaped, cluster_labels_predicted))
else:
silhouette_scores.append(np.nan) # Append NaN if silhouette cannot be calculated
# Get the original time indices corresponding to the *end* of each sequence for this class
original_indices = original_test_df.index[timesteps - 1:]
# Collect true labels, predicted cluster labels, and the actual original sequences for this class
# Loop through the generated sequences (represented by index j) and their corresponding true labels
for j, label in enumerate(y_true_categorical[timesteps - 1:]): # Iterate over true labels aligned with sequence ends
all_y_true_categorical.append(label) # Add the true label
all_predicted_cluster_labels.append(cluster_labels_predicted[j]) # Add the predicted cluster label
# Determine the start and end indices in the original DataFrame for the current sequence (j)
# The sequence at index j in the sequences array corresponds to data starting at index 'start_index' in the original DataFrame
start_index = original_test_df.index.get_loc(original_indices[j]) - (timesteps - 1)
end_index = start_index + timesteps
# Extract and add the original (unscaled) feature values for this specific sequence
all_original_test_sequences.append(original_test_df[features].iloc[start_index:end_index].values) # Append original sequence data
# Convert the accumulated lists across all test classes into NumPy arrays
all_y_true_categorical = np.array(all_y_true_categorical) # Array of true labels for all sequences
all_predicted_cluster_labels = np.array(all_predicted_cluster_labels) # Array of predicted cluster IDs for all sequences
all_original_test_sequences = np.array(all_original_test_sequences) # 3D array of original sequence data for all sequences
# Print overall evaluation metrics based on the accumulated data
print("\nEvaluation Metrics:")
# Print the mean of the recorded inertia values (Note: Check if this average of the model's final inertia is the desired metric here)
print(f"Inertia (final): {np.mean(inertia_values):.4f}")
# Print the mean of the calculated silhouette scores per class (ignoring NaNs)
print(f"Average Silhouette Score (valid cases): {np.nanmean(silhouette_scores):.4f}")
# --- Cluster Analysis: Map Cluster IDs to Dominant True Labels ---
# This maps each cluster ID to the true class label that appears most frequently among the sequences assigned to that cluster.
cluster_dominant_label = {} # Dictionary: {cluster_id: dominant_true_label}
for cluster_id in range(n_clusters): # Loop through each possible cluster ID
# Find the indices of all sequences that were predicted to belong to the current cluster_id
indices_in_cluster = np.where(all_predicted_cluster_labels == cluster_id)[0]
# If the cluster is not empty (contains sequences)
if len(indices_in_cluster) > 0:
# Get the true labels of all sequences that fall into this cluster
labels_in_cluster = all_y_true_categorical[indices_in_cluster]
# If there are actual labels in this subset (should be true if indices_in_cluster > 0)
if len(labels_in_cluster) > 0:
# Count the occurrences of each true label in this cluster and find the index of the most frequent one
dominant_label = np.argmax(np.bincount(labels_in_cluster))
cluster_dominant_label[cluster_id] = dominant_label # Assign the dominant true label to the cluster ID
else:
cluster_dominant_label[cluster_id] = -1 # If for some reason no labels were found, mark as -1
else:
cluster_dominant_label[cluster_id] = -1 # If the cluster is empty, mark as -1
# --- Generate Predicted Labels for Classification Evaluation ---
# Create a new array of predicted labels, where each sequence's predicted cluster ID is mapped to the cluster's dominant true label.
# This allows treating the clustering result as a classification output for evaluation.
# Use .get(cluster_id, -1) to handle cases where a cluster_id might not be in cluster_dominant_label (e.g., if cluster was empty).
predicted_labels_numeric = np.array([cluster_dominant_label.get(cluster_id, -1) for cluster_id in all_predicted_cluster_labels])
# --- Classification Evaluation (using mapped labels) ---
# Evaluate the performance by comparing the true labels with the predicted labels derived from cluster dominance.
# Only include instances where a dominant label was successfully assigned (not -1).
valid_indices = predicted_labels_numeric != -1 # Indices of sequences that have a valid predicted numeric label
# Proceed with evaluation metrics and confusion matrix if there are valid instances and more than one true class is present in these instances.
if np.sum(valid_indices) > 0 and len(np.unique(all_y_true_categorical[valid_indices])) > 1:
print("\nEvaluation Results (Clusters vs True Labels):")
# Print a detailed classification report (Precision, Recall, F1-score, Support for each class, and overall metrics).
print(classification_report(all_y_true_categorical[valid_indices], predicted_labels_numeric[valid_indices]))
# Compute the Confusion Matrix: rows are true labels, columns are predicted dominant labels.
cm = confusion_matrix(all_y_true_categorical[valid_indices], predicted_labels_numeric[valid_indices])
# Create a figure for the confusion matrix plot.
plt.figure(figsize=(8, 6))
# Use seaborn heatmap to visualize the confusion matrix.
# annot=True displays the values in the cells.
# fmt='d' formats the values as integers.
# cmap='Blues' sets the color map.
sns.heatmap(cm, annot=True, fmt='d', cmap='Blues')
plt.xlabel('Predicted Cluster (Dominant True Label)') # Label for the x-axis
plt.ylabel('True Label') # Label for the y-axis
plt.title('Confusion Matrix (Clusters vs True Labels)') # Title of the plot
plt.show() # Display the plot
else:
print("\nCould not perform detailed evaluation (not enough data or classes with assigned dominant labels).")
#################################################################################################
# --- Plotting Anomalies (Optional) ---
# This section plots time series data for sequences identified as anomalies based on clustering
#################################################################################################
# Check if the --plot_anomalies command line flag is provided
if plot_anomalies:
print("\nChecking anomaly data:")
# Define "anomaly clusters" as those whose dominant true label is greater than 0 (i.e., corresponds to any failure type).
anomaly_clusters = [cluster_id for cluster_id, label in cluster_dominant_label.items() if label > 0]
# Find the indices of all sequences that were assigned to any of the "anomaly clusters".
anomaly_indices = np.where(np.isin(all_predicted_cluster_labels, anomaly_clusters))[0]
# If any anomaly sequences were found
if len(anomaly_indices) > 0:
# Determine how many anomaly sequences to plot (up to 5, or fewer if less were found).
num_anomalies_to_plot = min(5, len(anomaly_indices))
colors = ['red', 'green', 'blue'] # Define a simple list of colors to cycle through for features
# Randomly select a few anomaly sequences and plot their original data.
for i in np.random.choice(anomaly_indices, num_anomalies_to_plot, replace=False): # Select random indices without replacement
# Print shape and first few values of the original sequence being plotted for debugging/information
print(f"Shape of all_original_test_sequences[{i}]: {all_original_test_sequences[i].shape}")
print(f"First few values of all_original_test_sequences[{i}]:\n{all_original_test_sequences[i][:5]}")
# Create a new figure for each individual anomaly plot.
plt.figure(figsize=(12, 6))
# Plot each original feature within the selected sequence over the timesteps.
for j, feature in enumerate(features): # j is feature index, feature is feature name
# Plot the feature values (y-axis) against the sequence timestep index (x-axis from 0 to timesteps-1).
# Use colors cycling through the defined list.
plt.plot(np.arange(timesteps), all_original_test_sequences[i][:, j], label=feature, color=colors[j % len(colors)])
# NOTE: true_label and predicted_cluster_for_title variables are not defined in this block before use in the title. This will cause a NameError.
# You need to add lines like:
# true_label = all_y_true_categorical[i]
# predicted_cluster_for_title = all_predicted_cluster_labels[i]
plt.title('Detected Anomalies (based on cluster dominance)') # Title of the plot
plt.xlabel('Time Step') # Label for the x-axis (timestep within the sequence)
plt.ylabel('Value') # Label for the y-axis (feature value)
plt.legend() # Add a legend to identify which line corresponds to which feature
plt.show() # Display the plot window
else:
# If no sequences were assigned to anomaly clusters, print a message.
print("No anomalies detected based on cluster dominance.")
#################################################################################################
# --- Plotting Misclassified Instances (Optional) ---
# This section plots time series data for sequences that were "misclassified" based on cluster dominance
#################################################################################################
# Check if the --plot_misclassified command line flag is provided
if plot_misclassified:
print("\nChecking misclassified data:")
# Find the indices of all sequences where the true label does NOT match the dominant label of their assigned cluster.
misclassified_indices = np.where(all_y_true_categorical != predicted_labels_numeric)[0]
# If any misclassified instances are found
if len(misclassified_indices) > 0:
# Determine how many misclassified sequences to plot (up to 5, or fewer if less were found).
num_misclassified_to_plot = min(5, len(misclassified_indices))
colors = ['red', 'green', 'blue'] # Define a simple list of colors to cycle through for features
# Randomly select a few misclassified sequences and plot their original data.
for i in np.random.choice(misclassified_indices, num_misclassified_to_plot, replace=False): # Select random indices without replacement
# Print shape and first few values of the original sequence being plotted for debugging/information
print(f"Shape of all_original_test_sequences[{i}]: {all_original_test_sequences[i].shape}")
print(f"First few values of all_original_test_sequences[{i}]:\n{all_original_test_sequences[i][:5]}")
# Create a new figure for each individual misclassified plot.
plt.figure(figsize=(12, 6))
# Plot each original feature within the selected sequence over the timesteps.
for j, feature in enumerate(features): # j is feature index, feature is feature name
# Plot the feature values (y-axis) against the sequence timestep index (x-axis from 0 to timesteps-1).
# Use colors cycling through the defined list.
plt.plot(np.arange(timesteps), all_original_test_sequences[i][:, j], label=feature, color=colors[j % len(colors)])
# NOTE: true_label and predicted_label are defined here, which fixes the NameError in this block.
true_label = all_y_true_categorical[i] # Get the true label for the current sequence
predicted_label = predicted_labels_numeric[i] # Get the predicted numeric label (dominant cluster label) for the current sequence
# Set the title for the misclassified plot, indicating true label and the dominant label of the predicted cluster.
plt.title(f'Misclassified Instance (True: {true_label}, Predicted Cluster: {predicted_label})') # Title including true label and dominant predicted label
plt.xlabel('Time Step') # Label for the x-axis
plt.ylabel('Value') # Label for the y-axis
plt.legend() # Add a legend to identify features
plt.show() # Display the plot window
else:
# If no misclassified sequences were found, print a message.
print("No misclassified instances found based on cluster dominance.")
# Return the arrays of true labels and predicted numeric labels (based on cluster dominance)
# These can be used for further analysis or saving results
return all_y_true_categorical, predicted_labels_numeric
#####################################################################################################
# --- Main Execution Flow ---
# This is the main part of the script that calls the functions to run the analysis
#####################################################################################################
# Create a list of true class labels for the test data.
# This list will contain arrays, where each array corresponds to a test class and contains the true label for every data point in that class.
# The labels are the index of the class (0 for Normal, 1 for Failure 1, etc.).
true_labels_list = []
for i, df in enumerate(dataTest): # Loop through each test DataFrame (each class)
# Create a numpy array filled with the class index 'i', with a length equal to the number of rows in the DataFrame.
true_labels_list.append(np.full(len(df), i))
# Call the main evaluation and plotting function.
# Pass the trained model, scaled/original test data, number of clusters, true labels, feature info, and plotting options.
# This function performs the prediction on test data, calculates evaluation metrics, and handles plotting based on flags.
y_true_final, y_pred_final = evaluate_and_plot_anomalies(kmeans, scaled_test_df_list, n_clusters, dataTest, true_labels_list, features, featureNames, unitNames, plot_anomalies=options.plot_anomalies, plot_misclassified=options.plot_misclassified)
#####################################################################################################
# --- Final Evaluation Metrics (on combined test data) ---
# This section calculates and prints overall evaluation metrics after processing all test data
#####################################################################################################
# Calculate and print the final Inertia and Silhouette Score for the combined test data.
# Check if there is any reshaped test data available (i.e., if any test data files were processed).
if X_test_reshaped_list:
# Vertically stack all the reshaped test data arrays from different classes into a single array.
# This array contains all flattened sequences from all test data.
X_test_combined_reshaped = np.vstack(X_test_reshaped_list)
# Concatenate all the predicted cluster labels from different classes into a single array.
all_cluster_labels_test = np.concatenate(cluster_labels_test_list)
# Print a header for the final evaluation metrics.
print("\nK-Means Model Evaluation on Combined Test Data:")
# Print the final Inertia of the trained K-Means model on the training data.
# Note: This inertia value is from the model fitting process, not a specific calculation on the combined test data.
print(f"Inertia: {kmeans.inertia_:.4f}")
# Calculate and print the Silhouette Score for the combined test data based on the predicted cluster labels.
# This metric evaluates how well-separated the clusters are based on the data points within them.
# Only calculate if there is more than one unique predicted cluster label and at least one data point.
if len(np.unique(all_cluster_labels_test)) > 1 and len(all_cluster_labels_test) > 0:
silhouette = silhouette_score(X_test_combined_reshaped, all_cluster_labels_test)
print(f"Silhouette Score: {silhouette:.4f}")
else:
# If Silhouette score cannot be calculated, print a message.
print("Silhouette Score: Not applicable for single cluster.")
else:
# If no test data sequences were available at all, print a message.
print("\nNo test data sequences available to evaluate Inertia and Silhouette Score.")