CSE572 Project 3 - Cluster Validation Project

import pandas as pd
import numpy as np

import matplotlib.pyplot as plt
from matplotlib.pyplot import get_cmap

from scipy.integrate import simpson

from sklearn.impute import KNNImputer
from sklearn.cluster import KMeans, DBSCAN
from sklearn.decomposition import PCA

from sklearn.preprocessing import StandardScaler

insulin_data = pd.read_csv('InsulinData.csv', low_memory=False)
cgm_data = pd.read_csv('CGMData.csv', low_memory=False)

insulin_data

#insulin_data['BWZ Carb Input (grams)']

cgm_data

#insulin_data.info()

cgm_data.info()

def preprocess_date(date_string):
    date_string = date_string.strip()
    if len(date_string.split()) == 1:
        return date_string
    elif len(date_string.split()) == 2:
        return date_string.split()[0]

insulin_data['Date'] = insulin_data['Date'].apply(preprocess_date)

insulin_data['datetime'] = pd.to_datetime(insulin_data['Date'] + ' ' + insulin_data['Time'])

cgm_data['datetime'] = pd.to_datetime(cgm_data['Date'] + ' ' + cgm_data['Time'])

insulin_data[['Date', 'Time', 'datetime']]

cgm_data[['Date', 'Time', 'datetime']]

# This datetime series contains the start of all non NaN and  positive carb input
# data from the column BWZ Carb Input (grams)
meal_times = \
    insulin_data[~insulin_data['BWZ Carb Input (grams)'].isna() & \
    (insulin_data['BWZ Carb Input (grams)'] > 0)]['datetime']
meal_times

meal_times.index

meal_carb_data = []

for i in meal_times.index:
    meal_carb_data.append(insulin_data['BWZ Carb Input (grams)'].iloc[i])

len(meal_carb_data)

meal_carb_data[0:10]

len(meal_carb_data)

max_meal_carbs = np.max(meal_carb_data)
max_meal_carbs

min_meal_carbs = np.min(meal_carb_data)
min_meal_carbs

bin_size = 20 # grams

number_of_bins = int(np.round(((max_meal_carbs - min_meal_carbs) / bin_size), decimals = 0))
number_of_bins

# Create 6 bins in the range of the data
bin_boundaries = np.linspace(min_meal_carbs, max_meal_carbs, 6) 
print("Bin boundaries:", bin_boundaries)

# Initialize meal data list
meal_glucose_data = []
meal_carb_for_analysis_data = []

# Loop through each meal start time
for index, meal_time in meal_times.items():
    # A meal windows is 2 hours long
    meal_window = insulin_data[(insulin_data['datetime'] > meal_time) & 
                               (insulin_data['datetime'] <= (meal_time + pd.Timedelta(hours=2)))]
    #print(meal_window)

    # If that meal window is empty, get corresponding CGM_
    # data that starts at tm-30 and ends at tm+2hrs
    if meal_window.empty:
        # Step 4 condition (a): No meal between tm and tm+2hrs, use this stretch as meal data
        meal_cgm_data = cgm_data[(cgm_data['datetime'] >= (meal_time - pd.Timedelta(minutes=30))) & 
                                 (cgm_data['datetime'] <= (meal_time + pd.Timedelta(hours=2)))]
    else:
        # Step 4 check for condition (b) or (c)
        next_meal_time = meal_window.iloc[0]['datetime']
        if next_meal_time < (meal_time + pd.Timedelta(hours=2)):
            # condition (b): There is a meal at tp between tm and tm+2hrs, use tp data
            meal_time = next_meal_time
            meal_cgm_data = cgm_data[(cgm_data['datetime'] >= (meal_time - pd.Timedelta(minutes=30))) & 
                                     (cgm_data['datetime'] <= (meal_time + pd.Timedelta(hours=2)))]
        # Condition (c): There is a meal exactly at tm+2hrs, use tm+1hr 30min to tm+4hrs
        elif next_meal_time == (meal_time + pd.Timedelta(hours=2)):
            meal_cgm_data = cgm_data[(cgm_data['datetime'] >= (meal_time + pd.Timedelta(hours=1, minutes=30))) & 
                                     (cgm_data['datetime'] <= (meal_time + pd.Timedelta(hours=4)))]

    # Check if the extracted data has the correct number of points (30 for meal data)
    if len(meal_cgm_data) == 30:
        #print(f"Meal Time Data {meal_time} has {len(meal_cgm_data)} number of points.")
        meal_glucose_data.append(meal_cgm_data['Sensor Glucose (mg/dL)'].values)
        meal_carb_for_analysis_data.append(insulin_data['BWZ Carb Input (grams)'].iloc[index])

len(meal_glucose_data)

len(meal_carb_for_analysis_data)

meal_bin_indices = np.digitize(meal_carb_for_analysis_data, bin_boundaries, right=False)
meal_bin_indices = meal_bin_indices - 1
print("Bin indices:", meal_bin_indices)

len(meal_bin_indices)

min(meal_bin_indices)

max(meal_bin_indices)

meal_data_matrix = np.vstack(meal_glucose_data)
meal_data_matrix

meal_data_matrix.shape

# 10%
dump_thresold = 0.1

knn_imputer = KNNImputer(n_neighbors=7)

cleaned_meal_glucose_data = []
cleaned_meal_bin_indices = []
cleaned_meal_carb_for_analysis_data = []

for index, row in enumerate(meal_data_matrix):
    
    proportion_missing_data = np.isnan(row).mean()

    if proportion_missing_data > dump_thresold:

        # Ignore this row
        #print(f"Ignoring row with {proportion_missing_data*100:.2f}% missing data")
        continue

    else:
        # The data is good, use it
        cleaned_meal_glucose_data.append(row)
        cleaned_meal_carb_for_analysis_data.append(meal_carb_for_analysis_data[index])
        cleaned_meal_bin_indices.append(meal_bin_indices[index])

len(cleaned_meal_glucose_data)

len(cleaned_meal_bin_indices)

cleaned_meal_glucose_data_matrix = \
    np.array(cleaned_meal_glucose_data)

cleaned_meal_glucose_data_matrix.shape

cleaned_imputed_meal_glucose_data_matrix = \
    knn_imputer.fit_transform(cleaned_meal_glucose_data_matrix)

#cleaned_imputed_meal_glucose_data_matrix[0]

def compute_first_derivative(glucose_data, delta_t=5):
    return np.diff(glucose_data) / delta_t

def compute_second_derivative(glucose_values, delta_t=5):
    second_derivative = np.zeros(len(glucose_values) - 2)
    for i in range(1, len(glucose_values) - 1):
        second_derivative[i - 1] = \
            (glucose_values[i + 1] - 2 * glucose_values[i] + glucose_values[i - 1]) / (delta_t ** 2)
    return second_derivative

def compute_AUC(glucose_values, delta_t=5):
    return simpson(glucose_values, dx=delta_t)

def feature_extractor(data_matrix):

    features = []

    for row in data_matrix:
        mean_glucose = np.mean(row)
        std_glucose  = np.std(row)
        max_glucose  = np.max(row)
        min_glucose  = np.min(row)
        
        first_derivative = compute_first_derivative(row)
        mean_first_derivate = np.mean(first_derivative)

        if len(row) > 2:
            second_derivative = compute_second_derivative(row)
            mean_second_derivative = np.mean(second_derivative)
        else:
            mean_second_derivative = 0

        AUC = compute_AUC(row)

        time_to_peak = np.argmax(row)

        features_vector = [
            mean_glucose,
            std_glucose,
            max_glucose,
            min_glucose,
            mean_first_derivate,
            mean_second_derivative,
            AUC,
            time_to_peak
        ]

        features.append(features_vector)
    
    return np.array(features)

meal_glucose_features_matrix = \
    feature_extractor(cleaned_imputed_meal_glucose_data_matrix)
#len(meal_gluxose_features_matrix[0])

len(cleaned_imputed_meal_glucose_data_matrix)

len(meal_glucose_features_matrix)

len(cleaned_meal_carb_for_analysis_data)

len(cleaned_meal_bin_indices)

cleaned_meal_carb_for_analysis_data[0:5]

cleaned_meal_bin_indices[0:5]

scaler = StandardScaler()
scaled_meal_glucose_features_matrix = \
    scaler.fit_transform(meal_glucose_features_matrix)

meal_glucose_features_matrix[0:5]

scaled_meal_glucose_features_matrix[0:5]

meal_bin_indices_column_vector = \
    np.array(cleaned_meal_bin_indices).reshape(-1, 1)

meal_bin_indices_column_vector[0:5]

#combined_feature_bin_data = \
#    np.hstack((meal_gluxose_features_matrix, meal_bin_indices_column_vector))
combined_feature_bin_data = meal_glucose_features_matrix
#combined_feature_bin_data = scaled_meal_glucose_features_matrix

combined_feature_bin_data[0:5]

number_of_bins

cleaned_meal_bin_indices[0:5]

max(cleaned_meal_bin_indices)

min(cleaned_meal_bin_indices)

# Set the number of clusters for K-means
k = int(number_of_bins)
y_bins = np.array(cleaned_meal_bin_indices)

# Initialize and fit KMeans
kmeans = KMeans(n_clusters=k, random_state=123)
y_kmeans = kmeans.fit_predict(scaled_meal_glucose_features_matrix)

y_kmeans[0:5]

y_kmeans.shape

# Set the number of clusters for K-means
k = int(number_of_bins)
y_bins = np.array(cleaned_meal_bin_indices)

# Initialize and fit KMeans
kmeans = KMeans(n_clusters=k, random_state=123)
y_kmeans = kmeans.fit_predict(scaled_meal_glucose_features_matrix)

# Get cluster labels and centroids
kmeans_cluster_labels = kmeans.labels_
centroids = kmeans.cluster_centers_

pca = PCA(n_components=2)
X_pca_kmeans = pca.fit_transform(scaled_meal_glucose_features_matrix)

# Plotting the clusters found by k-means
plt.figure(figsize=(10, 5))

# Get the viridis color map and extract six colors from it
viridis = get_cmap('viridis', k)
colors = [viridis(i) for i in range(k)]

# Subplot 1: Scatter plot by k-means clusters
plt.subplot(1, 2, 1)
for i in range(k):
    cluster_points = X_pca_kmeans[y_kmeans == i]
    plt.scatter(cluster_points[:, 0], cluster_points[:, 1], 
                color=colors[i], label=f'Cluster {i + 1}', alpha=0.6, marker='o')

# Plotting the centroids after PCA transformation
centroids_pca = pca.transform(centroids)
plt.scatter(centroids_pca[:, 0], centroids_pca[:, 1], 
            color='red', edgecolor='black', marker='X', s=100, label='Centroids')

plt.title("K-Means Clusters")
plt.xlabel('PCA Component 1')
plt.ylabel('PCA Component 2')
plt.legend()

plt.show()

# Set the number of clusters for K-means
k = int(number_of_bins)
y_bins = np.array(cleaned_meal_bin_indices)

# Initialize and fit KMeans
kmeans = KMeans(n_clusters=k, random_state=123)
y_kmeans = kmeans.fit_predict(scaled_meal_glucose_features_matrix)

# Get cluster labels and centroids
kmeans_cluster_labels = kmeans.labels_
centroids = kmeans.cluster_centers_

pca = PCA(n_components=2)
X_pca = pca.fit_transform(scaled_meal_glucose_features_matrix)

# Plotting the clusters found by k-means
plt.figure(figsize=(10, 5))

# Subplot 1: Scatter plot by original bins
plt.subplot(1, 2, 1)
plt.scatter(
    np.arange(0, y_bins.shape[0]),
    y_bins,
    c=y_bins,
    cmap='rainbow',
    marker='o')
plt.title("Original Bins")
plt.xlabel("Sample Index")
plt.ylabel("Bin Label")

# Subplot 2: Scatter plot by k-means clusters
plt.subplot(1, 2, 2)
plt.scatter(
    np.arange(0, y_kmeans.shape[0]),
    y_kmeans,
    c=y_kmeans,
    cmap='rainbow',
    marker='o')
plt.title("K-Means Clusters")
plt.xlabel("Sample Index")
plt.ylabel("Bin Label")

plt.show()

# Compute SSE
SSE_kmeans = np.float64(kmeans.inertia_)
SSE_kmeans

# Set the number of clusters for K-means
k = int(number_of_bins)
y_bins = np.array(cleaned_meal_bin_indices)

# Initialize and fit KMeans
kmeans = KMeans(n_clusters=k, random_state=123)
y_kmeans = kmeans.fit_predict(scaled_meal_glucose_features_matrix)

# Get cluster labels and centroids
kmeans_cluster_labels = kmeans.labels_
centroids = kmeans.cluster_centers_

pca = PCA(n_components=2)
X_pca = pca.fit_transform(scaled_meal_glucose_features_matrix)

# Plotting the clusters found by k-means
plt.figure(figsize=(10, 5))

# Get the viridis color map and extract six colors from it
viridis = get_cmap('viridis', k)
colors = [viridis(i) for i in range(k)]

# Subplot 1: Scatter plot by k-means clusters
plt.subplot(1, 2, 1)
for i in range(k):
    cluster_points = X_pca[y_kmeans == i]
    plt.scatter(cluster_points[:, 0], cluster_points[:, 1], 
                color=colors[i], label=f'Cluster {i + 1}', alpha=0.6, marker='o')

# Plotting the centroids after PCA transformation
centroids_pca = pca.transform(centroids)
plt.scatter(centroids_pca[:, 0], centroids_pca[:, 1], 
            color='red', edgecolor='black', marker='X', s=100, label='Centroids')

plt.title("K-Means Clusters")
plt.xlabel('PCA Component 1')
plt.ylabel('PCA Component 2')
plt.legend()

# Get axis limits to keep the same scale for the second plot
x_lim = plt.xlim()
y_lim = plt.ylim()

# Subplot 2: Scatter plot by DBSCAN clusters
plt.subplot(1, 2, 2)

eps = 280  # Maximum distance between two points to be considered neighbors
min_samples = 15  # Minimum number of points in a neighborhood to form a dense region

dbscan = DBSCAN(eps=eps, min_samples=min_samples)
y_dbscan = dbscan.fit_predict(meal_glucose_features_matrix)

# Create a color palette for the clusters found by DBSCAN, including noise (-1) as black
unique_labels = set(y_dbscan)
dbscan_colors = [plt.cm.viridis(float(i) / (max(unique_labels) if max(unique_labels) > 0 else 1))
                 for i in unique_labels]

# Color each point by DBSCAN cluster, use black for noise points
for label in unique_labels:
    if label == -1:
        # Noise points are labeled as -1
        plt.scatter(X_pca[y_dbscan == label, 0], X_pca[y_dbscan == label, 1], 
                    color='black', s=10, label='Noise', alpha=0.6, marker='x')
    else:
        plt.scatter(X_pca[y_dbscan == label, 0], X_pca[y_dbscan == label, 1], 
                    color=dbscan_colors[label], label=f'Cluster {label + 1}', alpha=0.6, marker='o')

plt.title("DBSCAN Clusters")
plt.xlabel('PCA Component 1')
plt.ylabel('PCA Component 2')
plt.xlim(x_lim)  # Apply same x-axis limits as K-means plot
plt.ylim(y_lim)  # Apply same y-axis limits as K-means plot
plt.legend()

plt.show()

y_bins = np.array(cleaned_meal_bin_indices)

# DBSCAN Parameters
eps = 280  # Maximum distance between two points to be considered neighbors
min_samples = 15  # Minimum number of points in a neighborhood to form a dense region

# Fit DBSCAN
dbscan = DBSCAN(eps=eps, min_samples=min_samples)
y_dbscan = dbscan.fit_predict(meal_glucose_features_matrix)

# Adjust labels to start from 0 (handle noise points separately)
# Noise points are marked as -1, so we'll make their label `max_label + 1`
adjusted_y_dbscan = y_dbscan.copy()
max_label = y_dbscan.max()

# Replace all noise labels (-1) with a new label which is `max_label + 1`
adjusted_y_dbscan[adjusted_y_dbscan == -1] = max_label + 1

# Plotting the clusters found by DBSCAN and original bins
plt.figure(figsize=(10, 5))

# Subplot 1: Scatter plot by original bins
plt.subplot(1, 2, 1)
plt.scatter(
    np.arange(0, y_bins.shape[0]),
    y_bins,
    c=y_bins,
    cmap='rainbow',
    marker='o')
plt.title("Original Bins")
plt.xlabel("Sample Index")
plt.ylabel("Bin Label")

# Subplot 2: Scatter plot by DBSCAN clusters
plt.subplot(1, 2, 2)
plt.scatter(
    np.arange(0, adjusted_y_dbscan.shape[0]),
    adjusted_y_dbscan,
    c=adjusted_y_dbscan,
    cmap='rainbow',
    marker='o')
plt.title("DBSCAN Clusters")
plt.xlabel("Sample Index")
plt.ylabel("Cluster Label")

plt.tight_layout()
plt.show()

# Get unique clusters excluding noise (-1)
unique_clusters = set(y_dbscan) - {-1}

# Initialize SSE
SSE_dbscan = 0

# Loop over each cluster to compute SSE
for cluster in unique_clusters:
    # Get points in the current cluster
    cluster_points = meal_glucose_features_matrix[y_dbscan == cluster]
    
    # Compute the pseudo-centroid (mean of points in the cluster)
    pseudo_centroid = np.mean(cluster_points, axis=0)
    
    # Compute SSE for this cluster
    sse_cluster = np.sum((cluster_points - pseudo_centroid) ** 2)
    
    # Add the cluster's SSE to the total SSE
    SSE_dbscan += sse_cluster

print(f"Approximate SSE for DBSCAN (excluding noise): {SSE_dbscan}")

df_kmeans_cluster = \
    pd.DataFrame(
        {'Cluster': kmeans_cluster_labels,
         'Bin': np.array(cleaned_meal_bin_indices)}
    )
df_kmeans_cluster

df_kmeans_cluster_bin_matrix = \
    pd.crosstab(
        df_kmeans_cluster['Cluster'], 
        df_kmeans_cluster['Bin'])
df_kmeans_cluster_bin_matrix

#y_dbscan

df_dbscan_cluster = \
    pd.DataFrame(
        {'Cluster': y_dbscan,
         'Bin': np.array(cleaned_meal_bin_indices)}
    )

df_dbscan_cluster = \
    df_dbscan_cluster[df_dbscan_cluster['Cluster'] != -1]
df_dbscan_cluster

df_dbscan_cluster_bin_matrix = \
    pd.crosstab(
        df_dbscan_cluster['Cluster'], 
        df_dbscan_cluster['Bin'])
df_dbscan_cluster_bin_matrix

def entropy(entries):
    entries = entries[entries > 0]
    return -np.sum(entries * np.log2(entries))

k_means_entropy_values = \
    df_kmeans_cluster_bin_matrix.apply(lambda row: entropy(row / row.sum()), axis=1)
k_means_entropy_values

dbscan_entropy_values = \
    df_dbscan_cluster_bin_matrix.apply(lambda row: entropy(row / row.sum()), axis=1)
dbscan_entropy_values

k_means_cluster_sizes = df_kmeans_cluster_bin_matrix.sum(axis=1)
k_means_cluster_sizes

dbscan_cluster_sizes = df_dbscan_cluster_bin_matrix.sum(axis=1)
dbscan_cluster_sizes

total_data_points = len(meal_glucose_features_matrix)
total_data_points

k_means_total_entropy = \
    np.sum(k_means_cluster_sizes * k_means_entropy_values) / total_data_points
k_means_total_entropy

dbscan_total_entropy = \
    np.sum(dbscan_cluster_sizes * dbscan_entropy_values) / total_data_points
dbscan_total_entropy

k_means_purity_values = \
    df_kmeans_cluster_bin_matrix.apply(lambda row: row.max() / row.sum(), axis=1)
k_means_purity_values

dbscan_purity_values = \
    df_dbscan_cluster_bin_matrix.apply(lambda row: row.max() / row.sum(), axis=1)
dbscan_purity_values

k_means_total_purity = \
    np.sum(k_means_cluster_sizes * k_means_purity_values) / total_data_points
k_means_total_purity

dbscan_total_purity = \
    np.sum(dbscan_cluster_sizes * dbscan_purity_values) / total_data_points
dbscan_total_purity

results = [
    SSE_kmeans,
    SSE_dbscan,
    k_means_total_entropy,
    dbscan_total_entropy,
    k_means_total_purity,
    dbscan_total_purity
]
results

results_matrix_1x6 = np.array(results).reshape(1, 6)
results_matrix_1x6

np.savetxt('Result.csv', results_matrix_1x6, delimiter=',', fmt='%10.5f')