機器學習叢集分析-從理論基礎到Python實作

def kmeans_clustering(data, k, max_iters=100):
    # 隨機選擇初始中心點
    centroids = data[np.random.choice(data.shape[0], k, replace=False)]
    
    for _ in range(max_iters):
        # 分配樣本到最近的中心點
        distances = np.sqrt(((data - centroids[:, np.newaxis])**2).sum(axis=2))
        labels = np.argmin(distances, axis=0)
        
        # 更新中心點
        new_centroids = np.array([data[labels == i].mean(axis=0) for i in range(k)])
        
        # 檢查收斂
        if np.all(centroids == new_centroids):
            break
            
        centroids = new_centroids
        
    return labels, centroids

def hierarchical_clustering(data, n_clusters):
    # 計算距離矩陣
    distances = np.zeros((len(data), len(data)))
    for i in range(len(data)):
        for j in range(len(data)):
            distances[i,j] = np.sqrt(np.sum((data[i] - data[j])**2))
    
    # 初始化叢集
    clusters = [[i] for i in range(len(data))]
    
    while len(clusters) > n_clusters:
        # 找出最近的兩個叢集
        min_dist = float('inf')
        merge_pairs = None
        
        for i in range(len(clusters)):
            for j in range(i+1, len(clusters)):
                dist = min([distances[x,y] for x in clusters[i] for y in clusters[j]])
                if dist < min_dist:
                    min_dist = dist
                    merge_pairs = (i,j)
        
        # 合併叢集
        clusters[merge_pairs[0]].extend(clusters[merge_pairs[1]])
        clusters.pop(merge_pairs[1])
    
    return clusters

def standardize_data(data):
    return (data - np.mean(data, axis=0)) / np.std(data, axis=0)

def elbow_method(data, max_k):
    distortions = []
    for k in range(1, max_k+1):
        labels, centroids = kmeans_clustering(data, k)
        distortion = sum(np.min(cdist(data, centroids, 'euclidean'), axis=1)) / data.shape[0]
        distortions.append(distortion)
    return distortions

def silhouette_score(data, labels):
    n_samples = len(data)
    scores = np.zeros(n_samples)
    
    for i in range(n_samples):
        # 計算樣本i與同叢集其他樣本的平均距離
        a = np.mean([np.sqrt(np.sum((data[i] - data[j])**2)) 
                    for j in range(n_samples) 
                    if labels[j] == labels[i] and i != j])
        
        # 計算樣本i與其他叢集的最小平均距離
        b = float('inf')
        for cluster in set(labels):
            if cluster != labels[i]:
                dist = np.mean([np.sqrt(np.sum((data[i] - data[j])**2)) 
                              for j in range(n_samples) 
                              if labels[j] == cluster])
                b = min(b, dist)
        
        scores[i] = (b - a) / max(a, b)
    
    return np.mean(scores)

from sklearn.cluster import DBSCAN
import numpy as np

# 建立 DBSCAN 模型
dbscan = DBSCAN(eps=0.5, min_samples=5)

# 假設 X 是我們的資料集
X = np.random.rand(100, 2)  # 示範用隨機資料
clusters = dbscan.fit_predict(X)

def simple_grid_clustering(data, grid_size):
    # 建立網格
    x_min, x_max = data[:, 0].min(), data[:, 0].max()
    y_min, y_max = data[:, 1].min(), data[:, 1].max()
    
    x_grid = np.linspace(x_min, x_max, grid_size)
    y_grid = np.linspace(y_min, y_max, grid_size)
    
    # 分配點到網格
    grid_counts = np.zeros((grid_size, grid_size))
    for point in data:
        x_idx = np.digitize(point[0], x_grid) - 1
        y_idx = np.digitize(point[1], y_grid) - 1
        grid_counts[x_idx, y_idx] += 1
        
    return grid_counts

from sklearn.mixture import GaussianMixture

# 建立 GMM 模型
gmm = GaussianMixture(n_components=3, random_state=42)

# 擬合模型並預測群集
clusters = gmm.fit_predict(X)

from sklearn.cluster import AgglomerativeClustering

# 建立階層式群集模型
hierarchical = AgglomerativeClustering(n_clusters=3)

# 進行群集分析
clusters = hierarchical.fit_predict(X)

import numpy as np
from scipy.spatial.distance import cdist
import matplotlib.pyplot as plt
from sklearn.datasets import make_blobs
from sklearn.metrics import adjusted_rand_score
from sklearn.cluster import KMeans

class KMeansClustering:
    def __init__(self, n_clusters=8, max_iter=300, tol=0.0001, random_state=0):
        self.n_clusters = n_clusters
        self.max_iter = max_iter
        self.tol = tol
        self.random_state = random_state

def _greedy_kmeans_plus_plus(self, X):
    np.random.seed(self.random_state)
    n_samples, n_features = X.shape
    n_local_trials = 2 + int(np.log(self.n_clusters))
    
    # 初始化中心點陣列
    centers = np.zeros((self.n_clusters, n_features))
    
    # 隨機選擇第一個中心點
    first_index = np.random.choice(np.arange(n_samples))
    centers[0] = X[first_index]
    
    # 計算到第一個中心點的距離
    sq_distances = cdist(X, centers[0].reshape(1, -1), 
                        metric='sqeuclidean').ravel()

# 讓我們深入解析這段機器學習聚類別演算法的程式碼，逐步瞭解其實作細節:

class KMeansClustering:
    def __init__(self, n_clusters=8, max_iter=300, tol=1e-4, random_state=None):
        # 初始化 K-Means 聚類別器引數
        self.n_clusters = n_clusters      # 群集數量
        self.max_iter = max_iter          # 最大迭代次數
        self.tol = tol                    # 收斂容忍度
        self.random_state = random_state  # 隨機種子
        
    def _greedy_kmeans_plus_plus(self, X):
        # K-Means++ 貪婪演算法初始化群集中心
        n_samples, n_features = X.shape
        centers = np.zeros((self.n_clusters, n_features))
        
        # 隨機選擇第一個中心點
        first_center = X[np.random.randint(n_samples)]
        centers[0] = first_center
        
        # 計算到第一個中心點的距離平方
        sq_distances = np.sum((X - first_center) ** 2, axis=1)
        
        # 選擇剩餘的中心點
        for i in range(1, self.n_clusters):
            min_cost = np.inf
            best_candidate_index = None
            min_new_sq_distances = None
            
            # 隨機選擇候選點
            n_candidates = 2 + int(np.log(i + 1))
            candidate_indices = np.random.randint(n_samples, size=n_candidates)
            
            # 評估每個候選點
            for candidate_index in candidate_indices:
                new_sq_distances = np.sum((X - X[candidate_index]) ** 2, axis=1)
                new_cost = np.sum(np.minimum(sq_distances, new_sq_distances))
                
                if new_cost < min_cost:
                    best_candidate_index = candidate_index
                    min_new_sq_distances = new_sq_distances
                    min_cost = new_cost
                    
            # 更新中心點和距離
            centers[i] = X[best_candidate_index]
            sq_distances = np.minimum(sq_distances, min_new_sq_distances)
            
        return centers

    def fit(self, X):
        # 訓練模型
        n_samples, n_features = X.shape
        self.inertia_ = np.inf  # 初始化群內平方和
        
        # 使用 K-Means++ 初始化中心點
        self.cluster_centers_ = self._greedy_kmeans_plus_plus(X)
        
        # Lloyd's algorithm 迭代最佳化
        for _ in range(self.max_iter):
            # 計算每個樣本到所有中心點的距離
            distances = cdist(X, self.cluster_centers_, metric='sqeuclidean')
            # 為每個樣本分配最近的中心點
            labels = np.argmin(distances, axis=1)
            # 計算新的群內平方和
            new_inertia = np.sum(np.min(distances, axis=1))
            
            # 更新群集中心
            new_centers = np.zeros((self.n_clusters, n_features))
            for k in range(self.n_clusters):
                new_centers[k] = np.mean(X[labels == k], axis=0)
                
            # 檢查收斂性
            if np.abs(new_inertia - self.inertia_) < self.tol:
                break
                
            self.inertia_ = new_inertia
            self.cluster_centers_ = new_centers
            
    def predict(self, X):
        # 預測新樣本的群集標籤
        distances = cdist(X, self.cluster_centers_, metric='sqeuclidean')
        predicted_labels = np.argmin(distances, axis=1)
        return predicted_labels

# 使用scikit-learn實作K-means聚類別
sk_kmeans = KMeans(n_clusters=8, n_init='auto', random_state=0)
sk_kmeans.fit(X1)
sk_kmeans_pred_res = sk_kmeans.predict(X1)
sk_kmeans_ari = adjusted_rand_score(y1, sk_kmeans_pred_res)
sk_kmeans_centroinds = sk_kmeans.cluster_centers_

print(f'Adjusted Rand Score for sk KMeans: {sk_kmeans_ari}', '', sep='\n')
print(sk_kmeans_centroinds, '', sep='\n')
print('prediction', sk_kmeans_pred_res, sep='\n')

plt.figure(figsize=(12, 5))
plt.subplot(1, 2, 1)
plt.subplots_adjust(wspace=0.2)
plt.scatter(X1[:, 0], X1[:, 1], c=kmeans_pred_res, cmap="rainbow")
plt.scatter(kmeans_centroinds[:, 0], kmeans_centroinds[:, 1], 
           marker="x", color="black", s=100)
plt.title("KMeans")
plt.xlabel("Feature 1")
plt.ylabel("Feature 2")

def single_linkage(cluster1, cluster2):
    return np.min([np.linalg.norm(a - b) for a in cluster1 for b in cluster2])

def complete_linkage(cluster1, cluster2):
    return np.max([np.linalg.norm(a - b) for a in cluster1 for b in cluster2])

def average_linkage(cluster1, cluster2):
    distances = [np.linalg.norm(a - b) for a in cluster1 for b in cluster2]
    return np.mean(distances)

def ward_linkage(cluster1, cluster2):
    c1_mean = np.mean(cluster1, axis=0)
    c2_mean = np.mean(cluster2, axis=0)
    n1, n2 = len(cluster1), len(cluster2)
    return (n1 * n2 / (n1 + n2)) * np.sum((c1_mean - c2_mean) ** 2)

class AgglomerativeClustering:
    def __init__(self, n_clusters=2, linkage='ward'):
        self.n_clusters = n_clusters
        self.linkage = linkage
        
    def fit(self, X):
        n_samples = X.shape[0]
        # 初始化叢集，每個樣本為一個叢集
        clusters = [[X[i]] for i in range(n_samples)]
        
        while len(clusters) > self.n_clusters:
            min_dist = float('inf')
            merge_i, merge_j = 0, 0
            
            # 尋找最近的兩個叢集
            for i in range(len(clusters)):
                for j in range(i + 1, len(clusters)):
                    dist = self._calculate_distance(clusters[i], clusters[j])
                    if dist < min_dist:
                        min_dist = dist
                        merge_i, merge_j = i, j
            
            # 合併叢集
            clusters[merge_i].extend(clusters[merge_j])
            clusters.pop(merge_j)
            
        return clusters

def initialize_distance_matrix(X):
    n_samples = X.shape[0]
    distances = np.zeros((n_samples, n_samples))
    for i in range(n_samples):
        for j in range(i + 1, n_samples):
            distances[i,j] = distances[j,i] = np.linalg.norm(X[i] - X[j])
    return distances

class HierarchicalAgglomerativeClustering:
    def __init__(self, n_clusters=2):
        self.n_clusters = n_clusters

    @staticmethod
    def _ward_distance(c1, c2):
        n1, n2 = len(c1), len(c2)
        c1_mean, c2_mean = np.mean(c1, axis=0), np.mean(c2, axis=0)
        sqeuclidean_dist = sqeuclidean(c1_mean, c2_mean)
        return (n1 * n2) / (n1 + n2) * sqeuclidean_dist

    @staticmethod
    def _update_labels(labels, min_cdist_idxs):
        labels[labels == min_cdist_idxs[1]] = min_cdist_idxs[0]
        labels[labels > min_cdist_idxs[1]] -= 1
        return labels

    def fit_predict(self, X):
        labels = np.arange(len(X))
        clusters = [[x] for x in X]
        
        while len(clusters) > self.n_clusters:
            min_cdist, min_cdist_idxs = np.inf, []
            for i in range(len(clusters) - 1):
                for j in range(i + 1, len(clusters)):
                    cdist = self._ward_distance(clusters[i], clusters[j])
                    if cdist < min_cdist:
                        min_cdist = cdist
                        min_cdist_idxs = (i, j)
            
            labels = self._update_labels(labels, min_cdist_idxs)
            clusters[min_cdist_idxs[0]].extend(
                clusters.pop(min_cdist_idxs[1]))
        
        return np.array(labels)

X2, y2 = make_blobs(n_samples=75, n_features=2, centers=5, random_state=0)

# 使用我們的實作
ac = HierarchicalAgglomerativeClustering(n_clusters=5)
ac_pred_res = ac.fit_predict(X2)
ac_ari = adjusted_rand_score(y2, ac_pred_res)

# 使用 scikit-learn 的實作
sk_ac = AgglomerativeClustering(n_clusters=5, linkage='ward')
sk_ac_pred_res = sk_ac.fit_predict(X2)
sk_ac_ari = adjusted_rand_score(y2, sk_ac_pred_res)

print(f'Adjusted Rand Score for sk AgglomerativeClustering: {sk_ac_ari}')
print('prediction', sk_ac_pred_res)

plt.figure(figsize=(12, 5))
plt.subplot(1, 2, 1)
plt.scatter(X2[:, 0], X2[:, 1], c=ac_pred_res, cmap='rainbow')
plt.title('AgglomerativeClustering')
plt.xlabel("Feature 1")
plt.ylabel("Feature 2")

linkage_matrix = linkage(X2, method='ward', metric='euclidean')
plt.figure(figsize=(10, 6))
dendrogram(linkage_matrix, color_threshold=10)

import numpy as np
import matplotlib.pyplot as plt
from scipy.sparse.linalg import eigsh
from scipy.sparse.csgraph import laplacian as csgraph_laplacian
from sklearn.datasets import make_moons
from sklearn.preprocessing import StandardScaler
from sklearn.cluster import KMeans
from sklearn.metrics.pairwise import rbf_kernel

class EnhancedSpectralClustering:
    def __init__(self, n_clusters=8, gamma=1.0, random_state=None):
        self.n_clusters = n_clusters
        self.gamma = gamma
        self.random_state = random_state

    @staticmethod
    def _normalize_eigenvectors(u):
        # 確保特徵向量的一致性
        max_indices = np.argmax(np.abs(u), axis=1)
        signs = np.sign(u[range(u.shape[0]), max_indices])
        return u * signs[:, np.newaxis]

    def _compute_spectral_embedding(self, affinity):
        # 計算正規化的拉普拉斯矩陣
        laplacian, degrees = csgraph_laplacian(
            affinity, 
            normed=True, 
            return_diag=True
        )
        laplacian *= -1

        # 初始化隨機向量
        initial_vector = np.random.RandomState(
            self.random_state
        ).uniform(-1, 1, laplacian.shape[0])

        # 計算特徵向量
        _, eigenvectors = eigsh(
            laplacian,
            k=self.n_clusters,
            sigma=1.0,
            which="LM",
            tol=0,
            v0=initial_vector
        )

        # 正規化處理
        normalized_eigenvectors = (
            eigenvectors.T[self.n_clusters::-1] / degrees
        )
        embedding = self._normalize_eigenvectors(normalized_eigenvectors)
        
        return embedding[:self.n_clusters].T

    def fit_predict(self, X):
        # 計算相似度矩陣
        affinity_matrix = rbf_kernel(X, gamma=self.gamma)
        
        # 取得譜嵌入
        embedding = self._compute_spectral_embedding(affinity_matrix)
        
        # 執行最終聚類別
        kmeans = KMeans(
            n_clusters=self.n_clusters,
            n_init='auto',
            random_state=self.random_state
        )
        
        return kmeans.fit_predict(embedding)

# 生成測試資料
X, y = make_moons(n_samples=200, noise=0.05, random_state=0)
X = StandardScaler().fit_transform(X)

# 執行聚類別
spectral = EnhancedSpectralClustering(n_clusters=2, gamma=1.0, random_state=0)
labels = spectral.fit_predict(X)

# 自定義的光譜分群實作
sc = ArpackSpectralClustering(n_clusters=2, gamma=10, random_state=0)
sc_pred_res = sc.fit_predict(X3)
sc_ari = adjusted_rand_score(y3, sc_pred_res)

# scikit-learn的光譜分群實作
sk_sc = SpectralClustering(n_clusters=2, gamma=10, random_state=0)
sk_sc_pred_res = sk_sc.fit_predict(X3)
sk_sc_ari = adjusted_rand_score(y3, sk_sc_pred_res)

from sklearn.cluster import DBSCAN
import numpy as np

# 初始化DBSCAN
dbscan = DBSCAN(eps=0.5, min_samples=5)

# 假設我們有一個資料集X
X = np.random.rand(100, 2)  # 生成範例資料

# 進行聚類別
clusters = dbscan.fit_predict(X)

# 取得聚類別標籤
labels = dbscan.labels_

class DBSCANClustering:
    def __init__(self, eps=0.5, min_samples=5, metric='euclidean', 
                 algorithm='auto', leaf_size=30):
        self.eps = eps  # 鄰域半徑
        self.min_samples = min_samples  # 最小樣本數
        self.metric = metric  # 距離度量方式
        self.algorithm = algorithm  # 鄰近點搜尋演算法
        self.leaf_size = leaf_size  # 葉子節點大小

    @staticmethod
    def _dbscan_inner(core_points, neighborhoods, labels):
        label, queue = 0, []
        for point in range(len(core_points)):
            if not core_points[point] or labels[point] != -1:
                continue
            labels[point] = label
            queue.append(point)

            while queue:
                current_point = queue.pop(0)
                if core_points[current_point]:
                    current_neighbors = neighborhoods[current_point]
                    for neighbor in current_neighbors:
                        if labels[neighbor] == -1:
                            labels[neighbor] = label
                            queue.append(neighbor)
            label += 1

def fit_predict(self, X):
    nn = NearestNeighbors(n_neighbors=self.min_samples, 
                         radius=self.eps, 
                         metric=self.metric,
                         algorithm=self.algorithm, 
                         leaf_size=self.leaf_size)
    
    neighborhoods = nn.fit(X).radius_neighbors(X, return_distance=False)
    labels = np.full(len(X), -1, dtype=np.intp)
    core_points = np.asarray([len(n) >= self.min_samples for n in neighborhoods])
    self._dbscan_inner(core_points, neighborhoods, labels)
    self.labels_ = labels
    return self.labels_

# DBSCAN 基礎實作
dbscan = DBSCANClustering(eps=0.3, min_samples=3)
dbscan_pred_res = dbscan.fit_predict(X4)
dbscan_ari = adjusted_rand_score(y4, dbscan_pred_res)

# scikit-learn 的 DBSCAN 實作
sk_dbscan = DBSCAN(eps=0.3, min_samples=3)
sk_dbscan_pred_res = sk_dbscan.fit_predict(X4)
sk_dbscan_ari = adjusted_rand_score(y4, sk_dbscan_pred_res)

import matplotlib.pyplot as plt
import numpy as np
from sklearn.cluster import DBSCAN

# 建立分群模型
dbscan = DBSCAN(eps=0.5, min_samples=5)

# 視覺化分群結果
def visualize_clusters(X, labels, title):
    plt.figure(figsize=(8, 6))
    plt.scatter(X[:, 0], X[:, 1], c=labels, cmap='rainbow')
    plt.title(title)
    plt.xlabel('特徵一')
    plt.ylabel('特徵二')
    plt.colorbar()
    plt.show()

# 進行分群預測
cluster_labels = dbscan.fit_predict(X)

from sklearn.neighbors import NearestNeighbors

def find_optimal_eps(X, n_neighbors=5):
    neigh = NearestNeighbors(n_neighbors=n_neighbors)
    neigh.fit(X)
    distances, _ = neigh.kneighbors(X)
    return np.sort(distances[:, -1])

# 計算最佳eps值
distances = find_optimal_eps(X)

from sklearn.preprocessing import StandardScaler
from joblib import parallel_backend

# 資料預處理
scaler = StandardScaler()
X_scaled = scaler.fit_transform(X)

# 平行運算最佳化
with parallel_backend('threading', n_jobs=-1):
    dbscan = DBSCAN(eps=0.5, min_samples=5)
    labels = dbscan.fit_predict(X_scaled)

from sklearn.metrics import silhouette_score

def evaluate_clustering(X, labels):
    # 計算分群品質
    if len(np.unique(labels)) > 1:
        score = silhouette_score(X, labels)
        return score
    return None

# 自適應密度估計
def adaptive_dbscan(X, eps_range, min_samples_range):
    best_score = -1
    best_params = None
    
    for eps in eps_range:
        for min_samples in min_samples_range:
            dbscan = DBSCAN(eps=eps, min_samples=min_samples)
            labels = dbscan.fit_predict(X)
            score = evaluate_clustering(X, labels)
            
            if score and score > best_score:
                best_score = score
                best_params = (eps, min_samples)
                
    return best_params

from sklearn.cluster import AffinityPropagation
import numpy as np

# 建立示範資料
X = np.array([[1, 2], [2, 3], [2, 2], [4, 6], [4, 5], [5, 5]])

# 初始化親和傳播模型
af = AffinityPropagation(damping=0.8, preference=-50)

# 執行聚類別
clusters = af.fit_predict(X)

from sklearn.preprocessing import StandardScaler

scaler = StandardScaler()
X_scaled = scaler.fit_transform(X)
af = AffinityPropagation(damping=0.8)
clusters = af.fit_predict(X_scaled)

# 責任值計算
def calculate_responsibility(similarity_matrix, availability_matrix):
    responsibility = np.zeros_like(similarity_matrix)
    for i in range(len(similarity_matrix)):
        for k in range(len(similarity_matrix)):
            max_value = float('-inf')
            for k_prime in range(len(similarity_matrix)):
                if k_prime != k:
                    value = similarity_matrix[i][k_prime] + availability_matrix[i][k_prime]
                    max_value = max(max_value, value)
            responsibility[i][k] = similarity_matrix[i][k] - max_value
    return responsibility

# 可達性計算
def calculate_accessibility(responsibility_matrix):
    accessibility = np.zeros_like(responsibility_matrix)
    for i in range(len(responsibility_matrix)):
        for k in range(len(responsibility_matrix)):
            sum_positive = 0
            for i_prime in range(len(responsibility_matrix)):
                if i_prime != i and i_prime != k:
                    sum_positive += max(0, responsibility_matrix[i_prime][k])
            accessibility[i][k] = min(0, responsibility_matrix[k][k] + sum_positive)
    return accessibility

# 群集中心識別
def identify_exemplars(responsibility_matrix, accessibility_matrix):
    exemplars = set()
    n = len(responsibility_matrix)
    for i in range(n):
        max_value = float('-inf')
        max_k = -1
        for k in range(n):
            current = responsibility_matrix[i][k] + accessibility_matrix[i][k]
            if current > max_value:
                max_value = current
                max_k = k
        if max_k == i:  # 如果節點選擇自己作為中心
            exemplars.add(i)
    return exemplars

def compute_similarity(X):
    # 計算負的歐氏距離平方作為相似度
    distances = euclidean_distances(X, squared=True)
    return -distances

def update_responsibility(S, A, damping):
    # 計算責任度矩陣
    AS = A + S
    row_max = np.max(AS, axis=1)
    row_max_indices = np.argmax(AS, axis=1)
    AS[np.arange(AS.shape[0]), row_max_indices] = -np.inf
    row_second_max = np.max(AS, axis=1)
    
    R = S - row_max.reshape(-1, 1)
    R[np.arange(len(R)), row_max_indices] = S[np.arange(len(R)), row_max_indices] - row_second_max
    
    return (1 - damping) * R + damping * R_old

def initialize_affinity_propagation(X, preference=None):
    # 計算相似度矩陣
    S = compute_similarity(X)
    
    if preference is None:
        # 使用相似度矩陣的中位數作為偏好值
        preference = np.median(S)
    
    # 設定對角線元素為偏好值
    np.fill_diagonal(S, preference)
    
    return S

def iterate_messages(S, max_iter=200, convergence_iter=15, damping=0.5):
    n_samples = S.shape[0]
    A = np.zeros((n_samples, n_samples))  # 可用度矩陣
    R = np.zeros((n_samples, n_samples))  # 責任度矩陣
    
    # 用於收斂檢查的變數
    convergence = np.zeros((convergence_iter, n_samples))
    
    for iter_count in range(max_iter):
        # 更新責任度
        R = update_responsibility(S, A, damping)
        
        # 更新可用度
        A = update_availability(R, damping)
        
        # 檢查收斂性
        exemplars = np.where(np.diag(A + R) > 0)[0]
        convergence[iter_count % convergence_iter] = exemplars
        
        if iter_count >= convergence_iter:
            if np.all(convergence[0] == convergence[1:]):
                break

def affinity_propagation_clustering(similarity_matrix, max_iter=200, convergence_iter=15, damping=0.5):
    n_samples = similarity_matrix.shape[0]
    samples_indexes = np.arange(n_samples)
    
    # 初始化矩陣
    responsibility_matrix = np.zeros((n_samples, n_samples))
    availability_matrix = np.zeros((n_samples, n_samples))
    exemplars_convergence_matrix = np.zeros((n_samples, convergence_iter))
    
    # 處理退化情況，加入隨機擾動
    rng = np.random.RandomState()
    noise = (np.finfo(similarity_matrix.dtype).eps * similarity_matrix + 
             np.finfo(similarity_matrix.dtype).tiny * 100) * rng.standard_normal(size=(n_samples, n_samples))
    similarity_matrix += noise
    
    # 主要迭代過程
    for iter in range(max_iter):
        # 計算責任矩陣 (Responsibility Matrix)
        temp_matrix = availability_matrix + similarity_matrix
        max_indexes = np.argmax(temp_matrix, axis=1)
        max_values = temp_matrix[samples_indexes, max_indexes]
        temp_matrix[samples_indexes, max_indexes] = -np.inf
        second_max_values = np.max(temp_matrix, axis=1)
        
        # 更新責任值
        np.subtract(similarity_matrix, max_values[:, None], temp_matrix)
        max_responsibility = similarity_matrix[samples_indexes, max_indexes] - second_max_values
        temp_matrix[samples_indexes, max_indexes] = max_responsibility
        
        # 應用衰減因子
        temp_matrix *= (1 - damping)
        responsibility_matrix = responsibility_matrix * damping + temp_matrix
        
        # 計算可用性矩陣 (Availability Matrix)
        np.maximum(responsibility_matrix, 0, temp_matrix)
        temp_matrix.flat[::n_samples + 1] = responsibility_matrix.flat[::n_samples + 1]
        temp_matrix -= np.sum(temp_matrix, axis=0)
        
        # 更新對角線元素
        diag_availability = np.diag(temp_matrix).copy()
        temp_matrix.clip(0, np.inf, temp_matrix)
        temp_matrix.flat[::n_samples + 1] = diag_availability
        
        # 應用衰減因子到可用性矩陣
        temp_matrix *= (1 - damping)
        availability_matrix = availability_matrix * damping - temp_matrix
        
        # 檢查收斂性
        exemplar = (np.diag(availability_matrix) + np.diag(responsibility_matrix)) > 0
        exemplars_convergence_matrix[:, iter % convergence_iter] = exemplar
        
        # 判斷是否收斂
        if iter >= convergence_iter:
            if check_convergence(exemplars_convergence_matrix, n_samples):
                break
    
    # 確定最終的叢集
    return finalize_clusters(similarity_matrix, exemplar)

# 人工智慧叢集（Clustering）實作與分析 - 親和傳播（Affinity Propagation）演算法

class AffinityPropagationClustering:
    def __init__(self, preference=None, convergence_iter=15, max_iter=200, damping=0.5, random_state=None):
        """
        初始化親和傳播演算法的引數
        
        引數說明：
        - preference: 選擇範例點作為叢集中心的偏好值
        - convergence_iter: 收斂迭代次數
        - max_iter: 最大迭代次數 
        - damping: 阻尼係數(0.5-1之間)
        - random_state: 隨機數種子
        """
        self.preference = preference
        self.convergence_iter = convergence_iter
        self.max_iter = max_iter
        self.damping = damping
        self.random_state = random_state

    def _affinity_propagation_inner(self, S, preference, convergence_iter, max_iter, damping, random_state):
        """
        親和傳播演算法的核心實作
        
        引數說明：
        - S: 相似度矩陣
        - preference: 偏好值 
        - convergence_iter: 收斂迭代次數
        - max_iter: 最大迭代次數
        - damping: 阻尼係數
        - random_state: 隨機數種子

        回傳:
        - cluster_centers_indices: 叢集中心的索引
        - labels: 每個樣本的叢集標籤
        """
        if preference is None:
            preference = np.median(S)

        n_samples = S.shape[0]

        # 初始化責任矩陣(R)和可用性矩陣(A) 
        A = np.zeros((n_samples, n_samples))
        R = np.zeros((n_samples, n_samples))
        
        # 主要迭代過程
        for it in range(max_iter):
            # 更新責任值
            Rold = R.copy()
            AS = A + S
            
            # 計算每個點對其他點的責任值
            Y = AS.max(axis=1)
            for i in range(n_samples):
                AS[i,i] = -np.inf
                Y2 = AS[i,:].max()
                R[i,:] = S[i,:] - Y[i]
                R[i,AS[i,:].argmax()] = S[i,AS[i,:].argmax()] - Y2
            
            R = (1-damping) * R + damping * Rold
            
            # 更新可用性值 
            Aold = A.copy()
            Rp = np.maximum(R, 0)
            for i in range(n_samples):
                Rp[i,i] = R[i,i]
            
            A = np.sum(Rp, axis=0).reshape(1,-1) - Rp
            dA = np.diag(A)
            A = np.minimum(A, 0)
            
            for i in range(n_samples):
                A[i,i] = dA[i]
                
            A = (1-damping) * A + damping * Aold

        # 確定叢集中心和標籤
        E = (A + R) > 0
        K = np.sum(E, axis=0)
        cluster_centers_indices = np.where(K>0)[0]
        labels = -1 * np.ones(n_samples)
        
        for i in range(n_samples):
            if E[i, i]:
                labels[i] = i
                continue
            
            exemplars = np.where(E[i,:])[0]
            if len(exemplars) == 0:
                continue
            
            labels[i] = exemplars[np.argmax(S[i,exemplars])]

        return cluster_centers_indices, labels

    def fit_predict(self, X):
        """
        對輸入資料進行擬合並回傳叢集標籤
        
        引數:
        X: 輸入資料矩陣
        
        回傳:
        labels: 每個樣本的叢集標籤
        """
        # 計算相似度矩陣
        self.affinity_matrix_ = -euclidean_distances(X, squared=True)

from sklearn.cluster import AffinityPropagation

def optimize_affinity_propagation(X, damping=0.5, max_iter=200):
    # 初始化親和性傳播模型
    ap = AffinityPropagation(
        damping=damping,
        max_iter=max_iter,
        random_state=42
    )
    
    # 執行叢集分析
    cluster_labels = ap.fit_predict(X)
    
    return ap, cluster_labels