1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203

import numpy as np
from sklearn.cluster import KMeans
from scipy.spatial.distance import cdist
 
class BisectingKMeans:
    """
    Bisecting K-means clustering algorithm.
    
    Parameters:
    -----------
    n_clusters : int
        Target number of clusters (K)
    n_bisections : int
        Number of bisection trials per split (best one is kept)
    selection : str
        Cluster selection strategy: 'largest', 'highest_sse', 'hybrid'
    max_iter_bisect : int
        Max iterations for each K=2 split
    random_state : int
        Random seed for reproducibility
    """
    
    def __init__(self, n_clusters=8, n_bisections=5, selection='largest',
                 max_iter_bisect=100, random_state=None):
        self.n_clusters = n_clusters
        self.n_bisections = n_bisections
        self.selection = selection
        self.max_iter_bisect = max_iter_bisect
        self.random_state = random_state
    
    def _compute_sse(self, X, center):
        """Compute sum of squared errors for a cluster."""
        return np.sum((X - center) ** 2)
    
    def _bisect_cluster(self, X, n_trials=5):
        """
        Split a cluster into two using K-means with K=2.
        Run multiple trials and keep the best split.
        """
        best_inertia = np.inf
        best_labels = None
        best_centers = None
        
        for trial in range(n_trials):
            km = KMeans(
                n_clusters=2,
                n_init=1,
                max_iter=self.max_iter_bisect,
                random_state=self.random_state + trial if self.random_state else None
            )
            km.fit(X)
            
            if km.inertia_ < best_inertia:
                best_inertia = km.inertia_
                best_labels = km.labels_
                best_centers = km.cluster_centers_
        
        return best_labels, best_centers, best_inertia
    
    def _select_cluster_to_split(self, clusters):
        """
        Select which cluster to split next.
        
        Strategies:
        - 'largest': Split the cluster with most points
        - 'highest_sse': Split the cluster with highest SSE
        - 'hybrid': Balance size and SSE
        """
        if self.selection == 'largest':
            # Select largest cluster
            sizes = [c['size'] for c in clusters]
            return np.argmax(sizes)
        
        elif self.selection == 'highest_sse':
            # Select cluster with highest SSE
            sses = [c['sse'] for c in clusters]
            return np.argmax(sses)
        
        elif self.selection == 'hybrid':
            # Score = size * sse (balance both factors)
            scores = [c['size'] * c['sse'] for c in clusters]
            return np.argmax(scores)
        
        else:
            raise ValueError(f"Unknown selection strategy: {self.selection}")
    
    def fit(self, X):
        """Fit bisecting K-means to data."""
        if self.random_state is not None:
            np.random.seed(self.random_state)
        
        n_samples, n_features = X.shape
        
        # Initialize: all points in one cluster
        all_indices = np.arange(n_samples)
        initial_center = X.mean(axis=0)
        initial_sse = self._compute_sse(X, initial_center)
        
        clusters = [{
            'indices': all_indices,
            'center': initial_center,
            'sse': initial_sse,
            'size': n_samples
        }]
        
        # Split history for building hierarchy
        self.split_history_ = []
        
        # Repeatedly bisect until K clusters
        while len(clusters) < self.n_clusters:
            # Select cluster to split
            split_idx = self._select_cluster_to_split(clusters)
            cluster_to_split = clusters[split_idx]
            
            # Get data for this cluster
            cluster_data = X[cluster_to_split['indices']]
            
            # Bisect
            labels, centers, inertia = self._bisect_cluster(
                cluster_data, 
                n_trials=self.n_bisections
            )
            
            # Create two new clusters
            mask_0 = labels == 0
            mask_1 = labels == 1
            
            indices_0 = cluster_to_split['indices'][mask_0]
            indices_1 = cluster_to_split['indices'][mask_1]
            
            new_cluster_0 = {
                'indices': indices_0,
                'center': centers[0],
                'sse': self._compute_sse(X[indices_0], centers[0]),
                'size': len(indices_0)
            }
            
            new_cluster_1 = {
                'indices': indices_1,
                'center': centers[1],
                'sse': self._compute_sse(X[indices_1], centers[1]),
                'size': len(indices_1)
            }
            
            # Record split for hierarchy
            self.split_history_.append({
                'parent_idx': split_idx,
                'parent_sse': cluster_to_split['sse'],
                'child_sses': [new_cluster_0['sse'], new_cluster_1['sse']],
                'n_clusters_after': len(clusters) + 1
            })
            
            # Replace parent with children
            clusters.pop(split_idx)
            clusters.append(new_cluster_0)
            clusters.append(new_cluster_1)
        
        # Store final results
        self.cluster_centers_ = np.array([c['center'] for c in clusters])
        self.labels_ = np.zeros(n_samples, dtype=int)
        for k, c in enumerate(clusters):
            self.labels_[c['indices']] = k
        
        # Compute inertia
        self.inertia_ = sum(c['sse'] for c in clusters)
        
        return self
    
    def predict(self, X):
        """Predict cluster labels for new data."""
        distances = cdist(X, self.cluster_centers_)
        return np.argmin(distances, axis=1)
    
    def fit_predict(self, X):
        """Fit and return labels."""
        self.fit(X)
        return self.labels_
 
 
# Example usage
np.random.seed(42)
 
# Create well-separated clusters
X = np.vstack([
    np.random.randn(200, 2) * 0.5 + [0, 0],
    np.random.randn(200, 2) * 0.5 + [5, 0],
    np.random.randn(200, 2) * 0.5 + [2.5, 4],
    np.random.randn(200, 2) * 0.5 + [7.5, 4],
    np.random.randn(200, 2) * 0.5 + [2.5, -4],
])
 
# Compare bisecting K-means with standard K-means
bisecting = BisectingKMeans(n_clusters=5, selection='largest', random_state=42)
bisecting.fit(X)
 
from sklearn.cluster import KMeans
kmeans = KMeans(n_clusters=5, random_state=42, n_init=10)
kmeans.fit(X)
 
print("Bisecting K-means vs Standard K-means:")
print(f"  Bisecting inertia: {bisecting.inertia_:.2f}")
print(f"  Standard inertia:  {kmeans.inertia_:.2f}")
print(f"  Improvement: {(1 - bisecting.inertia_ / kmeans.inertia_) * 100:.2f}%")

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51

import numpy as np
from sklearn.metrics import silhouette_score
 
def compare_selection_strategies(X, K, n_trials=5):
    """
    Compare different cluster selection strategies.
    """
    strategies = ['largest', 'highest_sse', 'hybrid']
    results = {}
    
    for strategy in strategies:
        inertias = []
        silhouettes = []
        
        for trial in range(n_trials):
            bkm = BisectingKMeans(
                n_clusters=K,
                selection=strategy,
                n_bisections=3,
                random_state=42 + trial
            )
            bkm.fit(X)
            
            inertias.append(bkm.inertia_)
            silhouettes.append(silhouette_score(X, bkm.labels_))
        
        results[strategy] = {
            'mean_inertia': np.mean(inertias),
            'std_inertia': np.std(inertias),
            'mean_silhouette': np.mean(silhouettes),
            'std_silhouette': np.std(silhouettes)
        }
        
        print(f"{strategy:15}: inertia={results[strategy]['mean_inertia']:.2f} "
              f"(±{results[strategy]['std_inertia']:.2f}), "
              f"silhouette={results[strategy]['mean_silhouette']:.4f}")
    
    return results
 
# Create test data with varying cluster sizes
np.random.seed(42)
X = np.vstack([
    np.random.randn(500, 2) * 1.0 + [0, 0],      # Large, spread
    np.random.randn(100, 2) * 0.3 + [5, 5],      # Small, tight
    np.random.randn(300, 2) * 0.5 + [10, 0],     # Medium
    np.random.randn(50, 2) * 2.0 + [-5, 5],      # Small, very spread
])
 
print("Comparing selection strategies on imbalanced data:")
print("=" * 60)
results = compare_selection_strategies(X, K=4, n_trials=5)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78

import numpy as np
from sklearn.cluster import KMeans
from sklearn.metrics import adjusted_rand_score, normalized_mutual_info_score
import time
 
def comprehensive_comparison(X, y_true, K, n_trials=10):
    """
    Comprehensive comparison of bisecting vs standard K-means.
    """
    results = {
        'bisecting': {'inertia': [], 'ari': [], 'nmi': [], 'time': []},
        'standard': {'inertia': [], 'ari': [], 'nmi': [], 'time': []}
    }
    
    for trial in range(n_trials):
        # Bisecting K-means
        start = time.time()
        bkm = BisectingKMeans(n_clusters=K, selection='highest_sse',
                               n_bisections=5, random_state=trial)
        bkm.fit(X)
        results['bisecting']['time'].append(time.time() - start)
        results['bisecting']['inertia'].append(bkm.inertia_)
        results['bisecting']['ari'].append(adjusted_rand_score(y_true, bkm.labels_))
        results['bisecting']['nmi'].append(normalized_mutual_info_score(y_true, bkm.labels_))
        
        # Standard K-means
        start = time.time()
        km = KMeans(n_clusters=K, n_init=10, random_state=trial)
        km.fit(X)
        results['standard']['time'].append(time.time() - start)
        results['standard']['inertia'].append(km.inertia_)
        results['standard']['ari'].append(adjusted_rand_score(y_true, km.labels_))
        results['standard']['nmi'].append(normalized_mutual_info_score(y_true, km.labels_))
    
    print("Comparison: Bisecting vs Standard K-means")
    print("=" * 60)
    
    for method in ['bisecting', 'standard']:
        print(f"\n{method.upper()}:")
        print(f"  Inertia:  {np.mean(results[method]['inertia']):.2f} "
              f"± {np.std(results[method]['inertia']):.2f}")
        print(f"  ARI:      {np.mean(results[method]['ari']):.4f} "
              f"± {np.std(results[method]['ari']):.4f}")
        print(f"  NMI:      {np.mean(results[method]['nmi']):.4f} "
              f"± {np.std(results[method]['nmi']):.4f}")
        print(f"  Time:     {np.mean(results[method]['time']):.4f}s "
              f"± {np.std(results[method]['time']):.4f}s")
    
    # Improvement summary
    inertia_improvement = (1 - np.mean(results['bisecting']['inertia']) / 
                           np.mean(results['standard']['inertia'])) * 100
    print(f"\nBisecting inertia improvement: {inertia_improvement:.2f}%")
    
    return results
 
 
# Generate hierarchical data (where bisecting should excel)
np.random.seed(42)
 
# Create 2-level hierarchical structure
# Top level: 2 super-clusters
# Bottom level: 4 clusters each
X_list = []
y_list = []
cluster_id = 0
 
for super_center in [[0, 0], [20, 20]]:
    for offset in [[-2, -2], [2, -2], [-2, 2], [2, 2]]:
        center = np.array(super_center) + np.array(offset) * 2
        points = np.random.randn(100, 2) * 0.5 + center
        X_list.append(points)
        y_list.extend([cluster_id] * 100)
        cluster_id += 1
 
X_hierarchical = np.vstack(X_list)
y_true = np.array(y_list)
 
results = comprehensive_comparison(X_hierarchical, y_true, K=8, n_trials=10)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148

import numpy as np
from collections import defaultdict
 
class BisectingKMeansWithHierarchy(BisectingKMeans):
    """
    Extended bisecting K-means that explicitly tracks hierarchy.
    """
    
    def fit(self, X):
        """Fit and build explicit tree structure."""
        super().fit(X)
        
        # Build tree from split history
        self._build_tree()
        
        return self
    
    def _build_tree(self):
        """Construct explicit tree structure from splits."""
        # Tree nodes: {node_id: {'parent': id, 'children': [ids], 'cluster_idx': k}}
        self.tree_ = {}
        
        # Start with root
        self.tree_[0] = {
            'parent': None,
            'children': [],
            'level': 0,
            'indices': None  # Will be filled
        }
        
        # Track current leaves
        current_leaves = [0]
        node_counter = 1
        
        # Replay splits to build tree
        for split in self.split_history_:
            split_node = current_leaves[split['parent_idx']]
            
            # Create two children
            child_1 = node_counter
            child_2 = node_counter + 1
            node_counter += 2
            
            self.tree_[split_node]['children'] = [child_1, child_2]
            self.tree_[child_1] = {
                'parent': split_node,
                'children': [],
                'level': self.tree_[split_node]['level'] + 1
            }
            self.tree_[child_2] = {
                'parent': split_node,
                'children': [],
                'level': self.tree_[split_node]['level'] + 1
            }
            
            # Update leaves
            current_leaves.pop(split['parent_idx'])
            current_leaves.append(child_1)
            current_leaves.append(child_2)
        
        # Assign final cluster indices to leaves
        self.leaf_nodes_ = current_leaves
        for k, leaf in enumerate(current_leaves):
            self.tree_[leaf]['cluster_idx'] = k
    
    def get_clusters_at_level(self, level):
        """
        Get cluster assignments at a specific hierarchy level.
        
        Level 0 = all data in one cluster
        Level 1 = 2 clusters
        Level log2(K) = K clusters
        """
        # Find nodes at this level (or their descendants)
        nodes_at_level = []
        
        def traverse(node_id, current_level):
            if current_level == level or len(self.tree_[node_id]['children']) == 0:
                nodes_at_level.append(node_id)
            else:
                for child in self.tree_[node_id]['children']:
                    traverse(child, current_level + 1)
        
        traverse(0, 0)
        
        # Map final cluster labels to level-clusters
        level_labels = np.zeros(len(self.labels_), dtype=int)
        
        for level_cluster_id, node in enumerate(nodes_at_level):
            # Find all leaf descendants
            leaf_clusters = self._get_leaf_descendants(node)
            for leaf_cluster in leaf_clusters:
                mask = self.labels_ == leaf_cluster
                level_labels[mask] = level_cluster_id
        
        return level_labels, len(nodes_at_level)
    
    def _get_leaf_descendants(self, node_id):
        """Get all leaf cluster indices under a node."""
        if len(self.tree_[node_id]['children']) == 0:
            return [self.tree_[node_id]['cluster_idx']]
        
        descendants = []
        for child in self.tree_[node_id]['children']:
            descendants.extend(self._get_leaf_descendants(child))
        return descendants
    
    def print_hierarchy(self):
        """Print tree structure."""
        def print_node(node_id, indent=""):
            node = self.tree_[node_id]
            
            if len(node['children']) == 0:
                k = node['cluster_idx']
                size = np.sum(self.labels_ == k)
                print(f"{indent}Leaf Cluster {k} (n={size})")
            else:
                print(f"{indent}Internal Node (level={node['level']})")
                for child in node['children']:
                    print_node(child, indent + "  ")
        
        print("Cluster Hierarchy:")
        print_node(0)
 
 
# Demonstrate hierarchy
np.random.seed(42)
X = np.vstack([
    np.random.randn(100, 2) * 0.5 + [0, 0],
    np.random.randn(100, 2) * 0.5 + [3, 0],
    np.random.randn(100, 2) * 0.5 + [6, 0],
    np.random.randn(100, 2) * 0.5 + [9, 0],
    np.random.randn(100, 2) * 0.5 + [12, 0],
    np.random.randn(100, 2) * 0.5 + [15, 0],
    np.random.randn(100, 2) * 0.5 + [18, 0],
    np.random.randn(100, 2) * 0.5 + [21, 0],
])
 
bkm = BisectingKMeansWithHierarchy(n_clusters=8, selection='largest', random_state=42)
bkm.fit(X)
 
bkm.print_hierarchy()
 
print("\nMulti-resolution clustering:")
for level in range(4):
    labels_level, n_clusters = bkm.get_clusters_at_level(level)
    print(f"  Level {level}: {n_clusters} clusters, "
          f"sizes = {np.bincount(labels_level)}")

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78

# Scikit-learn includes BisectingKMeans as of version 1.1
from sklearn.cluster import BisectingKMeans as SklearnBisecting
from sklearn.preprocessing import StandardScaler
from sklearn.metrics import silhouette_score
import numpy as np
 
def production_bisecting_kmeans(X, K_max=10, random_state=42):
    """
    Production bisecting K-means with automated K selection.
    """
    # Standardize
    scaler = StandardScaler()
    X_scaled = scaler.fit_transform(X)
    
    print(f"Data: {X.shape[0]} samples × {X.shape[1]} features")
    print()
    
    # Find optimal K using silhouette
    results = []
    best_K = 2
    best_score = -1
    
    for K in range(2, K_max + 1):
        bkm = SklearnBisecting(
            n_clusters=K,
            init='k-means++',         # Initialization for bisections
            n_init=5,                  # Number of bisection trials
            bisecting_strategy='largest_cluster',  # or 'biggest_inertia'
            random_state=random_state
        )
        
        labels = bkm.fit_predict(X_scaled)
        
        inertia = bkm.inertia_
        silhouette = silhouette_score(X_scaled, labels)
        
        results.append({
            'K': K,
            'inertia': inertia,
            'silhouette': silhouette
        })
        
        print(f"K={K}: silhouette={silhouette:.4f}, inertia={inertia:.2e}")
        
        if silhouette > best_score:
            best_score = silhouette
            best_K = K
    
    print(f"\nOptimal K: {best_K} (silhouette={best_score:.4f})")
    
    # Train final model
    final_model = SklearnBisecting(
        n_clusters=best_K,
        init='k-means++',
        n_init=10,  # More trials for final
        random_state=random_state
    )
    final_model.fit(X_scaled)
    
    return {
        'model': final_model,
        'scaler': scaler,
        'optimal_K': best_K,
        'results': results,
        'cluster_centers_original': scaler.inverse_transform(final_model.cluster_centers_)
    }
 
 
# Example
np.random.seed(42)
X = np.vstack([
    np.random.randn(200, 5) + np.random.randn(5) * 3
    for _ in range(6)
])
 
result = production_bisecting_kmeans(X, K_max=10)
 
print(f"\nFinal cluster sizes: {np.bincount(result['model'].labels_)}")