t-SNE 与 UMAP：非线性降维可视化

# PCA 在瑞士卷数据上的失败示例
from sklearn.datasets import make_swiss_roll
from sklearn.decomposition import PCA
import matplotlib.pyplot as plt

X, _ = make_swiss_roll(n_samples=1500, noise=0.05)
pca = PCA(n_components=2)
X_pca = pca.fit_transform(X)

plt.scatter(X_pca[:, 0], X_pca[:, 1], c=_, cmap="viridis")
plt.title("PCA on Swiss Roll: structure destroyed")
plt.show()

# 非线性降维正确展开瑞士卷
from sklearn.manifold import TSNE

tsne = TSNE(n_components=2, perplexity=30, random_state=42)
X_tsne = tsne.fit_transform(X)

plt.scatter(X_tsne[:, 0], X_tsne[:, 1], c=_, cmap="viridis")
plt.title("t-SNE on Swiss Roll: structure preserved")
plt.show()

# 手动实现 t-SNE 的高斯条件概率
import numpy as np
from scipy.spatial.distance import pdist, squareform

def compute_high_dim_affinity(X, sigma):
    """计算高维空间中基于高斯的条件概率 p(j|i)"""
    n = X.shape[0]
    sq_dists = squareform(pdist(X, "sqeuclidean"))
    P = np.exp(-sq_dists / (2 * sigma**2))
    np.fill_diagonal(P, 0)
    P = P / P.sum(axis=1, keepdims=True)
    P = np.maximum(P, 1e-12)
    return P

X = np.random.randn(100, 50)
P = compute_high_dim_affinity(X, sigma=1.0)
print("Per-row probability sums:", P.sum(axis=1)[:5])

# t 分布低维相似度计算
def compute_low_dim_affinity(Y):
    """计算低维空间中基于 t 分布的联合概率 q(j|i)"""
    n = Y.shape[0]
    sq_dists = squareform(pdist(Y, "sqeuclidean"))
    Q_num = 1.0 / (1.0 + sq_dists)
    np.fill_diagonal(Q_num, 0)
    Q = Q_num / Q_num.sum()
    Q = np.maximum(Q, 1e-12)
    return Q

Y = np.random.randn(100, 2)
Q = compute_low_dim_affinity(Y)
print("Q sum (should be ~1.0):", Q.sum())

# perplexity 网格搜索对比
from sklearn.datasets import load_digits
from sklearn.manifold import TSNE
import matplotlib.pyplot as plt

digits = load_digits()
X, y = digits.data, digits.target

fig, axes = plt.subplots(2, 3, figsize=(15, 10))
for i, perp in enumerate([5, 15, 30, 50, 75, 100]):
    tsne = TSNE(n_components=2, perplexity=perp, 
                random_state=42, n_iter=1000)
    X_emb = tsne.fit_transform(X[:300])
    ax = axes[i // 3, i % 3]
    sc = ax.scatter(X_emb[:, 0], X_emb[:, 1], 
                    c=y[:300], cmap="tab10", s=10)
    ax.set_title(f"perplexity={perp}")
plt.tight_layout()
plt.show()

# learning_rate 与 early_exaggeration 调优
perplexity = 30
results = {}
for lr in [50, 200, 500, 1000]:
    for ee in [4, 12, 24]:
        tsne = TSNE(
            n_components=2, perplexity=perplexity,
            learning_rate=lr, early_exaggeration=ee,
            n_iter=2000, random_state=42
        )
        X_emb = tsne.fit_transform(X[:300])
        # 用 KL divergence 评估质量
        results[(lr, ee)] = tsne.kl_divergence_

best = min(results, key=results.get)
print(f"Best: lr={best[0]}, ee={best[1]}, "
      f"KL={results[best]:.2f}")

# UMAP 核心：模糊单形集构建
import numpy as np
from sklearn.neighbors import NearestNeighbors

def compute_fuzzy_simplicial(X, k=15):
    """构建高维模糊单形集（简化版）"""
    n = X.shape[0]
    nn = NearestNeighbors(n_neighbors=k+1).fit(X)
    distances, indices = nn.kneighbors(X)
    distances = distances[:, 1:]  # 去掉自身
    indices = indices[:, 1:]
    
    # 局部连通度 rho 和宽度 sigma
    rho = distances[:, 0]  # 到最近邻的距离
    sigma = np.ones(n)
    for i in range(n):
        lo, hi = 0, 100
        for _ in range(64):  # 二分搜索
            mid = (lo + hi) / 2
            val = np.sum(np.exp(-(distances[i] - rho[i]) / mid))
            if val > np.log2(k):
                hi = mid
            else:
                lo = mid
        sigma[i] = hi
    
    # 成员度计算
    memberships = np.exp(-(distances - rho.reshape(-1,1)) / sigma.reshape(-1,1))
    return memberships, indices

X = np.random.randn(500, 20)
memberships, indices = compute_fuzzy_simplicial(X)
print("Membership matrix shape:", memberships.shape)

# UMAP 实战：完整流程
import umap
from sklearn.datasets import make_moons

X, y = make_moons(n_samples=500, noise=0.1)
reducer = umap.UMAP(
    n_components=2,
    n_neighbors=15,
    min_dist=0.1,
    metric="euclidean",
    random_state=42
)
embedding = reducer.fit_transform(X)

import matplotlib.pyplot as plt
plt.scatter(embedding[:, 0], embedding[:, 1], 
            c=y, cmap="Spectral", s=20)
plt.title("UMAP on Moons Dataset")
plt.show()

# 三种方法在同一数据集上的对比
from sklearn.datasets import load_digits
from sklearn.decomposition import PCA
from sklearn.manifold import TSNE
import umap
import time

X, y = load_digits().data, load_digits().target

# PCA
t0 = time.time()
X_pca = PCA(n_components=2).fit_transform(X)
print(f"PCA: {time.time()-t0:.2f}s")

# t-SNE
t0 = time.time()
X_tsne = TSNE(n_components=2, perplexity=30, 
              random_state=42).fit_transform(X)
print(f"t-SNE: {time.time()-t0:.2f}s")

# UMAP
t0 = time.time()
X_umap = umap.UMAP(n_components=2, random_state=42
                   ).fit_transform(X)
print(f"UMAP: {time.time()-t0:.2f}s")

# 量化对比：信任度（trustworthiness）评估
from sklearn.metrics import trustworthiness

methods = {
    "PCA": X_pca,
    "t-SNE": X_tsne,
    "UMAP": X_umap
}

for name, emb in methods.items():
    tw = trustworthiness(X, emb, n_neighbors=15)
    print(f"{name}: trustworthiness = {tw:.4f}")

# UMAP 通常获得最高的 trustworthiness
# 因为它同时保留了局部和全局结构

# 完整的降维可视化 pipeline
from sklearn.preprocessing import StandardScaler
from sklearn.decomposition import PCA
import umap
import matplotlib.pyplot as plt

def visualize_pipeline(X, labels, feature_names=None):
    """完整的降维可视化流程"""
    # Step 1: 标准化
    scaler = StandardScaler()
    X_scaled = scaler.fit_transform(X)
    
    # Step 2: PCA 探索方差
    pca = PCA()
    pca.fit(X_scaled)
    cumvar = np.cumsum(pca.explained_variance_ratio_)
    n_comp = np.argmax(cumvar > 0.95) + 1
    print(f"95% variance explained by {n_comp} components")
    
    # Step 3: PCA 降维预处理
    X_reduced = PCA(n_components=min(n_comp, 50)
                    ).fit_transform(X_scaled)
    
    # Step 4: UMAP 最终可视化
    reducer = umap.UMAP(n_neighbors=15, min_dist=0.1,
                        random_state=42)
    embedding = reducer.fit_transform(X_reduced)
    
    # Step 5: 可视化
    plt.figure(figsize=(10, 8))
    sc = plt.scatter(embedding[:, 0], embedding[:, 1],
                     c=labels, cmap="tab20", s=15, alpha=0.7)
    plt.colorbar(sc, label="Cluster")
    plt.title("UMAP Visualization (PCA preprocessed)")
    plt.tight_layout()
    plt.show()
    
    return embedding

# 量化评估：trustworthiness + silhouette
from sklearn.metrics import (trustworthiness, 
                              silhouette_score)
from sklearn.neighbors import NearestNeighbors

def evaluate_embedding(X_high, X_low, labels):
    """综合评估降维质量"""
    # Trustworthiness: 局部结构保留程度
    tw = trustworthiness(X_high, X_low, n_neighbors=15)
    
    # Silhouette: 聚类清晰度
    sil = silhouette_score(X_low, labels)
    
    # Neighborhood preservation
    nn_high = NearestNeighbors(n_neighbors=15
                               ).fit(X_high).kneighbors()[1]
    nn_low = NearestNeighbors(n_neighbors=15
                              ).fit(X_low).kneighbors()[1]
    preserved = np.mean([
        len(set(nn_high[i]) & set(nn_low[i])) / 15
        for i in range(len(X_high))
    ])
    
    print(f"Trustworthiness: {tw:.4f}")
    print(f"Silhouette: {sil:.4f}")
    print(f"Neighbor preservation: {preserved:.4f}")
    return {"tw": tw, "sil": sil, "np": preserved}

# MNIST 可视化：t-SNE + UMAP 对比
from sklearn.datasets import fetch_openml
import umap
from sklearn.manifold import TSNE
import matplotlib.pyplot as plt

# 加载 MNIST
mnist = fetch_openml("mnist_784", version=1, as_frame=False)
X, y = mnist.data, mnist.target.astype(int)

# 随机采样 7000 张（t-SNE 计算量大）
rng = np.random.RandomState(42)
idx = rng.choice(len(X), 7000, replace=False)
X_sample, y_sample = X[idx], y[idx]

# t-SNE（较慢）
tsne = TSNE(n_components=2, perplexity=50, 
            n_iter=1500, random_state=42)
X_tsne = tsne.fit_transform(X_sample)

# UMAP（较快）
umap_model = umap.UMAP(n_neighbors=10, min_dist=0.1,
                       random_state=42)
X_umap = umap_model.fit_transform(X_sample)

# 对比可视化
fig, axes = plt.subplots(1, 2, figsize=(16, 7))
for ax, emb, title in [(axes[0], X_tsne, "t-SNE"),
                        (axes[1], X_umap, "UMAP")]:
    sc = ax.scatter(emb[:, 0], emb[:, 1], 
                    c=y_sample, cmap="tab10", 
                    s=5, alpha=0.6)
    ax.set_title(title, fontsize=14)
    ax.set_xticks([])
    ax.set_yticks([])
plt.tight_layout()
plt.show()

# 词嵌入可视化：GloVe + UMAP
import numpy as np
import umap
import matplotlib.pyplot as plt
from sklearn.metrics.pairwise import cosine_distances

# 加载预训练的 GloVe 词向量（示例简化版）
word_vectors = {
    "king": np.random.randn(50),
    "queen": np.random.randn(50),
    "prince": np.random.randn(50),
    "cat": np.random.randn(50),
    "dog": np.random.randn(50),
    "bird": np.random.randn(50),
    "red": np.random.randn(50),
    "blue": np.random.randn(50),
    "green": np.random.randn(50),
    "one": np.random.randn(50),
    "two": np.random.randn(50),
    "three": np.random.randn(50),
}

words = list(word_vectors.keys())
X = np.array(list(word_vectors.values()))
categories = ["royal"]*3 + ["animal"]*3 + ["color"]*3 + ["number"]*3

# UMAP 降维
reducer = umap.UMAP(n_components=2, n_neighbors=5,
                    min_dist=0.3, metric="cosine",
                    random_state=42)
embedding = reducer.fit_transform(X)

# 可视化 + 标注
fig, ax = plt.subplots(figsize=(12, 8))
cmap = {"royal": "gold", "animal": "forestgreen",
        "color": "crimson", "number": "royalblue"}
for word, (x, y), cat in zip(words, embedding, categories):
    ax.scatter(x, y, c=cmap[cat], s=100)
    ax.annotate(word, (x, y), fontsize=10,
                fontweight="bold")
ax.set_title("Word Embedding UMAP Visualization")
ax.grid(True, alpha=0.3)
plt.show()