Unsupervised Learning — Clustering

import numpy as np
import matplotlib.pyplot as plt

class KMeans:
    """K-Means clustering from scratch with K-Means++ initialization."""
    
    def __init__(self, k=3, max_iters=100, tol=1e-4, init='kmeans++', random_state=42):
        self.k = k
        self.max_iters = max_iters
        self.tol = tol
        self.init = init
        self.random_state = random_state
        self.centroids = None
        self.labels = None
        self.inertia_ = None
        self.n_iters_ = 0
    
    def _init_centroids(self, X):
        """Initialize centroids using random or K-Means++."""
        rng = np.random.RandomState(self.random_state)
        n = X.shape[0]
        
        if self.init == 'random':
            indices = rng.choice(n, self.k, replace=False)
            return X[indices].copy()
        
        elif self.init == 'kmeans++':
            centroids = [X[rng.randint(n)]]
            
            for _ in range(1, self.k):
                # Compute squared distances to nearest centroid
                dists = np.array([
                    min(np.sum((x - c) ** 2) for c in centroids)
                    for x in X
                ])
                # Probability proportional to D(x)^2
                probs = dists / dists.sum()
                # Weighted random selection
                idx = rng.choice(n, p=probs)
                centroids.append(X[idx])
            
            return np.array(centroids)
    
    def _assign_clusters(self, X):
        """Assign each point to nearest centroid."""
        # Vectorized distance computation
        # X: (n, d), centroids: (k, d)
        # Using broadcasting: (n, 1, d) - (1, k, d) -> (n, k, d)
        diffs = X[:, np.newaxis, :] - self.centroids[np.newaxis, :, :]
        distances = np.sum(diffs ** 2, axis=2)  # (n, k)
        return np.argmin(distances, axis=1)
    
    def _update_centroids(self, X, labels):
        """Update centroids as cluster means."""
        new_centroids = np.zeros_like(self.centroids)
        for j in range(self.k):
            cluster_points = X[labels == j]
            if len(cluster_points) > 0:
                new_centroids[j] = cluster_points.mean(axis=0)
            else:
                # Handle empty cluster: reinitialize randomly
                new_centroids[j] = X[np.random.randint(X.shape[0])]
        return new_centroids
    
    def _compute_inertia(self, X, labels):
        """Compute within-cluster sum of squares (WCSS)."""
        inertia = 0
        for j in range(self.k):
            cluster_points = X[labels == j]
            if len(cluster_points) > 0:
                inertia += np.sum((cluster_points - self.centroids[j]) ** 2)
        return inertia
    
    def fit(self, X):
        """Fit K-Means to data X."""
        self.centroids = self._init_centroids(X)
        
        for i in range(self.max_iters):
            # Assignment step
            self.labels = self._assign_clusters(X)
            
            # Update step
            new_centroids = self._update_centroids(X, self.labels)
            
            # Check convergence
            shift = np.sum((new_centroids - self.centroids) ** 2)
            self.centroids = new_centroids
            self.n_iters_ = i + 1
            
            if shift < self.tol:
                break
        
        self.inertia_ = self._compute_inertia(X, self.labels)
        return self
    
    def predict(self, X):
        """Predict cluster labels for new data."""
        return self._assign_clusters(X)
    
    def fit_predict(self, X):
        """Fit and return labels."""
        self.fit(X)
        return self.labels


# --- Demo ---
np.random.seed(42)
# Generate 3 blobs
X = np.vstack([
    np.random.randn(50, 2) + [0, 0],
    np.random.randn(50, 2) + [5, 5],
    np.random.randn(50, 2) + [10, 0],
])

km = KMeans(k=3, init='kmeans++')
labels = km.fit_predict(X)

print(f"Converged in {km.n_iters_} iterations")
print(f"Inertia (WCSS): {km.inertia_:.2f}")
print(f"Centroids:\n{km.centroids}")

# Elbow method
inertias = []
K_range = range(1, 10)
for k in K_range:
    km_temp = KMeans(k=k)
    km_temp.fit(X)
    inertias.append(km_temp.inertia_)

print("\nElbow Method - Inertias:")
for k, inertia in zip(K_range, inertias):
    bar = '█' * int(inertia / 50)
    print(f"  k={k}: {inertia:8.1f} |{bar}")

import numpy as np
from collections import deque

class DBSCAN:
    """DBSCAN clustering from scratch."""
    
    def __init__(self, eps=0.5, min_pts=5):
        self.eps = eps
        self.min_pts = min_pts
        self.labels = None
        self.core_points_ = None
        self.n_clusters_ = 0
    
    def _get_neighbors(self, X, point_idx):
        """Find all points within eps of X[point_idx]."""
        distances = np.sqrt(np.sum((X - X[point_idx]) ** 2, axis=1))
        return np.where(distances <= self.eps)[0]
    
    def fit(self, X):
        """Fit DBSCAN to data X."""
        n = X.shape[0]
        self.labels = np.full(n, -1)  # -1 = unvisited/noise
        visited = np.zeros(n, dtype=bool)
        cluster_id = 0
        core_mask = np.zeros(n, dtype=bool)
        
        # Precompute distance matrix for efficiency
        # (for small datasets; for large ones, use spatial index)
        dist_matrix = np.sqrt(
            np.sum((X[:, np.newaxis] - X[np.newaxis, :]) ** 2, axis=2)
        )
        
        for i in range(n):
            if visited[i]:
                continue
            visited[i] = True
            
            # Find neighbors
            neighbors = np.where(dist_matrix[i] <= self.eps)[0]
            
            if len(neighbors) < self.min_pts:
                # Mark as noise (may become border later)
                self.labels[i] = -1
                continue
            
            # Core point — start new cluster
            core_mask[i] = True
            self.labels[i] = cluster_id
            
            # BFS expansion
            queue = deque(neighbors.tolist())
            
            while queue:
                q = queue.popleft()
                
                if not visited[q]:
                    visited[q] = True
                    q_neighbors = np.where(dist_matrix[q] <= self.eps)[0]
                    
                    if len(q_neighbors) >= self.min_pts:
                        core_mask[q] = True
                        queue.extend(q_neighbors.tolist())
                
                # Assign to cluster if not yet assigned
                if self.labels[q] == -1:
                    self.labels[q] = cluster_id
            
            cluster_id += 1
        
        self.core_points_ = core_mask
        self.n_clusters_ = cluster_id
        return self
    
    def fit_predict(self, X):
        """Fit and return labels."""
        self.fit(X)
        return self.labels


# --- Demo ---
np.random.seed(42)
# Two crescent moons (non-spherical clusters)
from sklearn.datasets import make_moons
X, y_true = make_moons(n_samples=200, noise=0.1, random_state=42)

db = DBSCAN(eps=0.2, min_pts=5)
labels = db.fit_predict(X)

print(f"Number of clusters found: {db.n_clusters_}")
print(f"Noise points: {np.sum(labels == -1)}")
print(f"Core points: {np.sum(db.core_points_)}")
print(f"Cluster sizes: {[np.sum(labels == i) for i in range(db.n_clusters_)]}")

import numpy as np

def silhouette_score(X, labels):
    """Compute mean silhouette score from scratch."""
    n = X.shape[0]
    unique_labels = np.unique(labels)
    k = len(unique_labels)
    
    if k <= 1 or k >= n:
        return 0.0  # Undefined for 1 cluster or n clusters
    
    # Pairwise distance matrix
    dist_matrix = np.sqrt(
        np.sum((X[:, np.newaxis] - X[np.newaxis, :]) ** 2, axis=2)
    )
    
    silhouettes = np.zeros(n)
    
    for i in range(n):
        own_cluster = labels[i]
        own_mask = labels == own_cluster
        
        # a(i): mean intra-cluster distance
        if own_mask.sum() > 1:
            a_i = dist_matrix[i][own_mask].sum() / (own_mask.sum() - 1)
        else:
            a_i = 0
        
        # b(i): min mean distance to other clusters
        b_i = np.inf
        for label in unique_labels:
            if label == own_cluster or label == -1:
                continue
            other_mask = labels == label
            if other_mask.sum() > 0:
                mean_dist = dist_matrix[i][other_mask].mean()
                b_i = min(b_i, mean_dist)
        
        if b_i == np.inf:
            silhouettes[i] = 0
        else:
            silhouettes[i] = (b_i - a_i) / max(a_i, b_i)
    
    return silhouettes.mean()


# --- Demo ---
np.random.seed(42)
X = np.vstack([
    np.random.randn(30, 2) + [0, 0],
    np.random.randn(30, 2) + [5, 5],
    np.random.randn(30, 2) + [10, 0],
])

from kmeans_scratch import KMeans  # from above

for k in range(2, 7):
    km = KMeans(k=k)
    labels = km.fit_predict(X)
    score = silhouette_score(X, labels)
    bar = '█' * int(score * 30)
    print(f"k={k}: Silhouette = {score:.4f} |{bar}")
# Best k should be 3 (highest silhouette)

import tensorflow as tf
import numpy as np

class TFKMeans:
    """GPU-accelerated K-Means using TensorFlow."""
    
    def __init__(self, k=3, max_iters=100, tol=1e-4):
        self.k = k
        self.max_iters = max_iters
        self.tol = tol
    
    @tf.function
    def _assign_step(self, X, centroids):
        """Vectorized assignment using TF ops."""
        # X: (n, d), centroids: (k, d)
        # Expand: (n, 1, d) - (1, k, d) = (n, k, d)
        diffs = tf.expand_dims(X, 1) - tf.expand_dims(centroids, 0)
        distances = tf.reduce_sum(diffs ** 2, axis=2)  # (n, k)
        return tf.argmin(distances, axis=1)  # (n,)
    
    @tf.function
    def _update_step(self, X, labels, k):
        """Update centroids as cluster means."""
        new_centroids = []
        for j in tf.range(k):
            mask = tf.equal(labels, tf.cast(j, labels.dtype))
            cluster_pts = tf.boolean_mask(X, mask)
            new_centroids.append(tf.reduce_mean(cluster_pts, axis=0))
        return tf.stack(new_centroids)
    
    def fit(self, X_np):
        """Fit K-Means on data."""
        X = tf.constant(X_np, dtype=tf.float32)
        n = X.shape[0]
        
        # K-Means++ init
        indices = [np.random.randint(n)]
        for _ in range(1, self.k):
            centroids_so_far = tf.gather(X, indices)
            diffs = tf.expand_dims(X, 1) - tf.expand_dims(centroids_so_far, 0)
            dists = tf.reduce_min(tf.reduce_sum(diffs**2, axis=2), axis=1)
            probs = dists / tf.reduce_sum(dists)
            idx = np.random.choice(n, p=probs.numpy())
            indices.append(idx)
        
        self.centroids = tf.Variable(tf.gather(X, indices))
        
        for i in range(self.max_iters):
            labels = self._assign_step(X, self.centroids)
            new_centroids = self._update_step(X, labels, self.k)
            
            shift = tf.reduce_sum((new_centroids - self.centroids) ** 2)
            self.centroids.assign(new_centroids)
            
            if shift.numpy() < self.tol:
                print(f"Converged at iteration {i+1}")
                break
        
        self.labels_ = labels.numpy()
        return self


# --- Demo ---
np.random.seed(42)
X = np.vstack([
    np.random.randn(100, 2) + [0, 0],
    np.random.randn(100, 2) + [6, 6],
    np.random.randn(100, 2) + [12, 0],
]).astype(np.float32)

model = TFKMeans(k=3)
model.fit(X)
print(f"TF Centroids:\n{model.centroids.numpy()}")
print(f"Cluster sizes: {[np.sum(model.labels_ == i) for i in range(3)]}")


# --- Deep Clustering with Autoencoder ---
class DeepCluster(tf.keras.Model):
    """Autoencoder-based deep clustering."""
    
    def __init__(self, latent_dim=10, n_clusters=3):
        super().__init__()
        self.encoder = tf.keras.Sequential([
            tf.keras.layers.Dense(128, activation='relu'),
            tf.keras.layers.Dense(64, activation='relu'),
            tf.keras.layers.Dense(latent_dim, activation='relu'),
        ])
        self.decoder = tf.keras.Sequential([
            tf.keras.layers.Dense(64, activation='relu'),
            tf.keras.layers.Dense(128, activation='relu'),
            tf.keras.layers.Dense(784, activation='sigmoid'),  # MNIST
        ])
        self.cluster_centers = tf.Variable(
            tf.random.normal([n_clusters, latent_dim])
        )
    
    def call(self, x):
        z = self.encoder(x)
        x_hat = self.decoder(z)
        return x_hat, z

print("\nDeep Clustering architecture defined for high-dimensional data.")

import numpy as np
import matplotlib.pyplot as plt
from sklearn.cluster import KMeans, DBSCAN, AgglomerativeClustering
from sklearn.mixture import GaussianMixture
from sklearn.preprocessing import StandardScaler
from sklearn.pipeline import Pipeline
from sklearn.metrics import (
    silhouette_score, silhouette_samples,
    davies_bouldin_score, calinski_harabasz_score,
    adjusted_rand_score
)
from sklearn.datasets import make_blobs, make_moons
from scipy.cluster.hierarchy import dendrogram, linkage
import warnings
warnings.filterwarnings('ignore')

# =========================================
# 1. K-MEANS with Elbow & Silhouette
# =========================================
print("=" * 60)
print("1. K-MEANS CLUSTERING")
print("=" * 60)

X, y_true = make_blobs(n_samples=300, centers=4, cluster_std=0.8,
                        random_state=42)
scaler = StandardScaler()
X_scaled = scaler.fit_transform(X)

# Elbow Method
inertias = []
sil_scores = []
K_range = range(2, 10)

for k in K_range:
    km = KMeans(n_clusters=k, init='k-means++', n_init=10, random_state=42)
    km.fit(X_scaled)
    inertias.append(km.inertia_)
    sil_scores.append(silhouette_score(X_scaled, km.labels_))

print("\nElbow Method Results:")
for k, inertia, sil in zip(K_range, inertias, sil_scores):
    bar_i = '█' * int(inertia / 20)
    bar_s = '█' * int(sil * 40)
    print(f"  k={k}: Inertia={inertia:7.1f} |{bar_i}")
    print(f"       Silhouette={sil:.4f} |{bar_s}")

# Best K-Means
best_k = list(K_range)[np.argmax(sil_scores)]
print(f"\nOptimal k by Silhouette: {best_k}")

km_final = KMeans(n_clusters=best_k, init='k-means++', n_init=10, 
                   random_state=42)
labels_km = km_final.fit_predict(X_scaled)

print(f"Final Inertia: {km_final.inertia_:.2f}")
print(f"Silhouette: {silhouette_score(X_scaled, labels_km):.4f}")
print(f"Davies-Bouldin: {davies_bouldin_score(X_scaled, labels_km):.4f}")
print(f"Calinski-Harabasz: {calinski_harabasz_score(X_scaled, labels_km):.1f}")
print(f"ARI (vs true): {adjusted_rand_score(y_true, labels_km):.4f}")


# =========================================
# 2. HIERARCHICAL CLUSTERING
# =========================================
print("\n" + "=" * 60)
print("2. HIERARCHICAL CLUSTERING")
print("=" * 60)

# Test all linkage types
for link in ['ward', 'complete', 'average', 'single']:
    hc = AgglomerativeClustering(
        n_clusters=4, linkage=link
    )
    labels_hc = hc.fit_predict(X_scaled)
    sil = silhouette_score(X_scaled, labels_hc)
    ari = adjusted_rand_score(y_true, labels_hc)
    print(f"  {link:10s} linkage: Silhouette={sil:.4f}, ARI={ari:.4f}")

# Dendrogram with scipy
Z = linkage(X_scaled[:50], method='ward')  # Subset for readability
print(f"\nLinkage matrix shape: {Z.shape}")
print("  Columns: [idx1, idx2, distance, n_points]")
print(f"  First 3 merges:")
for row in Z[:3]:
    print(f"    Merge clusters {int(row[0])} & {int(row[1])}, "
          f"dist={row[2]:.3f}, size={int(row[3])}")


# =========================================
# 3. DBSCAN
# =========================================
print("\n" + "=" * 60)
print("3. DBSCAN CLUSTERING")
print("=" * 60)

X_moons, y_moons = make_moons(n_samples=300, noise=0.08, random_state=42)
X_moons_scaled = StandardScaler().fit_transform(X_moons)

# Parameter search
print("\nParameter sensitivity:")
for eps in [0.2, 0.3, 0.5, 0.7]:
    for min_pts in [3, 5, 10]:
        db = DBSCAN(eps=eps, min_samples=min_pts)
        labels_db = db.fit_predict(X_moons_scaled)
        n_clusters = len(set(labels_db)) - (1 if -1 in labels_db else 0)
        n_noise = (labels_db == -1).sum()
        if n_clusters >= 2:
            sil = silhouette_score(X_moons_scaled[labels_db != -1],
                                    labels_db[labels_db != -1])
        else:
            sil = -1
        print(f"  eps={eps}, min_pts={min_pts}: "
              f"clusters={n_clusters}, noise={n_noise}, sil={sil:.3f}")

# Best DBSCAN
db_best = DBSCAN(eps=0.3, min_samples=5)
labels_db = db_best.fit_predict(X_moons_scaled)
print(f"\nBest DBSCAN: {len(set(labels_db)) - 1} clusters found")
print(f"Noise points: {(labels_db == -1).sum()}")
print(f"ARI vs true: {adjusted_rand_score(y_moons, labels_db):.4f}")


# =========================================
# 4. GAUSSIAN MIXTURE MODEL
# =========================================
print("\n" + "=" * 60)
print("4. GAUSSIAN MIXTURE MODEL (EM)")
print("=" * 60)

# BIC/AIC for model selection
print("\nModel selection via BIC/AIC:")
for k in range(2, 8):
    gmm = GaussianMixture(n_components=k, covariance_type='full',
                           n_init=3, random_state=42)
    gmm.fit(X_scaled)
    labels_gmm = gmm.predict(X_scaled)
    sil = silhouette_score(X_scaled, labels_gmm)
    print(f"  k={k}: BIC={gmm.bic(X_scaled):8.1f}, "
          f"AIC={gmm.aic(X_scaled):8.1f}, Silhouette={sil:.4f}")

# Best GMM
gmm_best = GaussianMixture(n_components=4, covariance_type='full',
                             random_state=42)
gmm_best.fit(X_scaled)
labels_gmm = gmm_best.predict(X_scaled)
proba = gmm_best.predict_proba(X_scaled)

print(f"\nBest GMM (k=4):")
print(f"  Converged: {gmm_best.converged_}")
print(f"  Iterations: {gmm_best.n_iter_}")
print(f"  ARI vs true: {adjusted_rand_score(y_true, labels_gmm):.4f}")
print(f"\nSoft assignments (first 5 points):")
for i in range(5):
    probs_str = ', '.join(f'{p:.3f}' for p in proba[i])
    print(f"  Point {i}: [{probs_str}] → Cluster {labels_gmm[i]}")

print(f"\nMixing weights: {gmm_best.weights_.round(3)}")


# =========================================
# 5. COMPREHENSIVE COMPARISON
# =========================================
print("\n" + "=" * 60)
print("5. ALGORITHM COMPARISON")
print("=" * 60)

algorithms = {
    'K-Means (k=4)': KMeans(n_clusters=4, random_state=42),
    'Hierarchical (Ward)': AgglomerativeClustering(n_clusters=4, linkage='ward'),
    'DBSCAN (0.5, 5)': DBSCAN(eps=0.5, min_samples=5),
    'GMM (k=4)': GaussianMixture(n_components=4, random_state=42),
}

print(f"\n{'Algorithm':<25} {'Clusters':>8} {'Silhouette':>10} "
      f"{'DB-Index':>10} {'ARI':>8}")
print("-" * 65)

for name, algo in algorithms.items():
    if hasattr(algo, 'fit_predict'):
        labels = algo.fit_predict(X_scaled)
    else:
        algo.fit(X_scaled)
        labels = algo.predict(X_scaled)
    
    mask = labels >= 0
    n_clusters = len(set(labels[mask]))
    
    if n_clusters >= 2:
        sil = silhouette_score(X_scaled[mask], labels[mask])
        db = davies_bouldin_score(X_scaled[mask], labels[mask])
    else:
        sil = db = float('nan')
    ari = adjusted_rand_score(y_true, labels)
    
    print(f"{name:<25} {n_clusters:>8d} {sil:>10.4f} {db:>10.4f} {ari:>8.4f}")