# ----------------------------------------------
# Симуляция 4.1: «Суперкластер»
# ----------------------------------------------
# Что делает скрипт:
# 1) Генерирует временные ряды на графе с экспоненциальной связностью,
#    считает корреляционную матрицу C.
# 2) Строит kNN-граф по |C| (показатель глобальной связности) и
#    «пороговый» граф по перцентилю |C| (для выделения компонент).
# 3) Считает второе собственное значение лапласиана (λ2) как метрику связности.
# 4) Рисует три рисунка и СОХРАНЯЕТ их в папку ./out РЯДОМ СО СКРИПТОМ:
#    - Теплокарта |C|
#    - PCA-раскладка узлов + рёбра kNN
#    - Спектры лапласианов (kNN и |C|)
# 5) Печатает в ТЕРМИНАЛ ключевые метрики.
#
# Зависимости: numpy, matplotlib
# ----------------------------------------------
import os
import numpy as np
import matplotlib.pyplot as plt

# ---------- Вспомогательные функции ----------
def ensure_outdir(path=""):
    script_dir = os.path.dirname(os.path.abspath(__file__))
    outdir = os.path.join(script_dir, path)
    os.makedirs(outdir, exist_ok=True)
    return outdir

def set_seed(seed=42):
    """Фиксируем зерно генератора для воспроизводимости."""
    np.random.seed(seed)

def build_W(N, xi=8.0, aniso=0.1):
    """
    Строим матрицу связей W (симметричную):
    - Базовая часть: экспоненциальное убывание связи по |i-j| с характерным радиусом xi.
    - Добавочная слабая анизотропия (две «полосы»), чтобы разрушить идеальную симметрию.
    - Нормируем по спектральному радиусу, чтобы значения были устойчивыми численно.
    """
    i = np.arange(N)[:, None]
    j = np.arange(N)[None, :]
    D = np.abs(i - j).astype(float)
    W = np.exp(-D / xi)
    np.fill_diagonal(W, 0.0)
    if aniso > 0:
        W += aniso * np.exp(-np.abs((i - j) - N/3) / (max(xi/2,1e-6)))
        W += aniso * np.exp(-np.abs((i - j) + N/3) / (max(xi/2,1e-6)))
        np.fill_diagonal(W, 0.0)
    # Нормировка по спектральному радиусу
    s = np.max(np.abs(np.linalg.eigvalsh(W)))
    if s > 1e-9:
        W = W / s
    return W

def integrate(N=80, T=20.0, dt=0.05, beta=1.0, gamma=0.8, noise=0.05, xi=8.0, aniso=0.1, seed=42):
    """
    Простейшая нелинейная интеграция на графе:
    x(t+dt) = (1 - gamma*dt) * x + beta*dt * tanh(W @ x) + noise * sqrt(dt) * N(0, I)
    Возвращает:
      X — траектории (steps x N),
      C — матрица корреляций (N x N),
      W — матрица связей (N x N).
    """
    set_seed(seed)
    steps = int(T / dt)
    W = build_W(N, xi=xi, aniso=aniso)
    x = 0.1 * np.random.randn(N)
    X = np.zeros((steps, N))
    sqrt_dt = np.sqrt(dt)
    for t in range(steps):
        x = (1 - gamma*dt) * x + beta*dt * np.tanh(W @ x) + noise * sqrt_dt * np.random.randn(N)
        X[t] = x
    # Матрица корреляций по времени
    Xc = X - X.mean(axis=0, keepdims=True)
    std = Xc.std(axis=0, keepdims=True) + 1e-12
    Xn = Xc / std
    C = (Xn.T @ Xn) / (Xn.shape[0] - 1)
    C = np.clip(C, -1.0, 1.0)
    return X, C, W

def laplacian_from_affinity(A):
    """
    Строим нормированный лапласиан графа L = I - D^{-1/2} A D^{-1/2},
    где A — матрица «родства» (аффинности), D — диагональ степеней.
    """
    A = np.maximum(A, 0.0)
    np.fill_diagonal(A, 0.0)
    d = A.sum(axis=1)
    dinv2 = 1.0 / np.sqrt(d + 1e-12)
    L = np.eye(A.shape[0]) - (dinv2[:,None] * A) * dinv2[None,:]
    return L

def fiedler_lambda(L):
    """Второе собственное значение лапласиана (λ2) — индикатор связности графа."""
    vals = np.linalg.eigvalsh(L)
    vals = np.sort(np.real(vals))
    if len(vals) >= 2:
        return float(vals[1])
    return float(vals[0])

def knn_affinity_from_C(C, k=6):
    """
    Строим kNN-граф по метрике d = 1 - |C| (т.е. «ближе» = выше корреляция).
    Веса рёбер = |C|. Симметризуем по максимуму.
    kNN обычно остаётся связным при умеренно больших k и высоких средних |C|.
    """
    S = np.abs(C)
    N = S.shape[0]
    D = 1.0 - S
    A = np.zeros_like(S)
    for i in range(N):
        order = np.argsort(D[i] + np.eye(N)[i]*1e9)  # исключаем самих себя
        nbrs = order[:k]
        A[i, nbrs] = S[i, nbrs]
    A = np.maximum(A, A.T)
    np.fill_diagonal(A, 0.0)
    return A

def connected_components_binary(A_bin):
    """Компоненты связности для НЕвзвешенного графа (BFS)."""
    N = A_bin.shape[0]
    visited = np.zeros(N, dtype=bool)
    labels = -np.ones(N, dtype=int)
    comp = 0
    for start in range(N):
        if not visited[start]:
            queue = [start]
            visited[start] = True
            labels[start] = comp
            while queue:
                u = queue.pop(0)
                neighbors = np.where(A_bin[u] > 0)[0]
                for v in neighbors:
                    if not visited[v]:
                        visited[v] = True
                        labels[v] = comp
                        queue.append(v)
            comp += 1
    return comp, labels

def auto_percentile_threshold(C, target_range=(2,6), start=85, minp=50, maxp=98, step=3):
    """
    Подбор перцентиля по |C|, чтобы число компонент было в target_range.
    Возвращает лучший найденный вариант (даже если не идеально попали):
      (p, thr, A_bin, n_comp, labels, sizes)
    """
    absC = np.abs(C)
    N = C.shape[0]
    best = None
    best_gap = 10**9
    best_balance = 10**9
    idx = np.triu_indices(N, 1)
    for p in range(start, maxp+1, step):
        thr = np.percentile(absC[idx], p)
        A = (absC >= thr).astype(float)
        np.fill_diagonal(A, 0.0)
        n_comp, labels = connected_components_binary(A)
        sizes = np.bincount(labels, minlength=n_comp)
        # Насколько далеко от желаемого диапазона по числу компонент
        gap = 0
        if n_comp < target_range[0]:
            gap = target_range[0] - n_comp
        elif n_comp > target_range[1]:
            gap = n_comp - target_range[1]
        # Насколько «несбалансирован» разбиение (стандартное отклонение размеров)
        balance = sizes.std() if n_comp > 1 else 1e9
        if (gap < best_gap) or (gap == best_gap and balance < best_balance):
            best_gap = gap
            best_balance = balance
            best = (p, thr, A, n_comp, labels, sizes)
        # Если попали в коридор и относительно сбалансировано — хватит
        if target_range[0] <= n_comp <= target_range[1] and balance < 0.25 * N:
            break
    return best

def cluster_metrics(C, labels):
    """
    Простейшие метрики по каждому кластеру (по меткам labels):
    - size: размер кластера
    - mass: сумма квадратов внутренних корреляций |C|^2 (интенсивность/«масса»)
    - charge: сумма внешних связей (с другими кластерами)
    - spin_like: мера асимметрии C_U (≈0 — почти симметрично)
    - lambda2: связность внутреннего подграфа (λ2 лапласиана по |C|)
    """
    uniq = np.unique(labels)
    out = []
    for u in uniq:
        idx = np.where(labels == u)[0]
        n = len(idx)
        Cin = C[np.ix_(idx, idx)]
        Cout = C[np.ix_(idx, np.where(labels != u)[0])]
        mass = float(np.sum(np.abs(Cin)**2))
        charge = float(np.sum(Cout))
        spin_like = float(np.linalg.norm(Cin - Cin.T) / (np.linalg.norm(Cin) + 1e-12))
        Ain = np.abs(Cin)
        Lin = laplacian_from_affinity(Ain)
        lam2 = fiedler_lambda(Lin)
        out.append(dict(id=int(u), size=int(n), mass=mass, charge=charge, spin_like=spin_like, lambda2=lam2))
    out.sort(key=lambda d: d["mass"], reverse=True)
    return out

def pca_2d_from_rows(M):
    """PCA до 2D (через SVD) для отображения узлов как точек в плоскости."""
    X = M - M.mean(axis=0, keepdims=True)
    U, S, Vt = np.linalg.svd(X, full_matrices=False)
    Z = X @ Vt[:2].T  # N x 2
    return Z

def spectrum(L):
    """Возвращает отсортированный спектр собственных значений лапласиана."""
    vals = np.linalg.eigvalsh(L)
    return np.sort(np.real(vals))

# ---------- Основной запуск (пресет «Суперкластер») ----------
def run_supercluster_preset():
    outdir = ensure_outdir()

    # Параметры, дающие глобальную связность (суперкластер в kNN)
    params = dict(
        N=300, T=20.0, dt=0.05, beta=1.0, gamma=0.8, noise=0.05, xi=9.0, aniso=0.10,
        seed=42, k=6
    )

    print("=== Симуляция: «Суперкластер» (пресет) ===")
    for k,v in params.items():
        if k!='k':
            print(f"{k} = {v}")
    print(f"k (kNN) = {params['k']}")

    # Интеграция и корреляции
    X, C, W = integrate(N=params['N'], T=params['T'], dt=params['dt'], beta=params['beta'],
                        gamma=params['gamma'], noise=params['noise'], xi=params['xi'],
                        aniso=params['aniso'], seed=params['seed'])

    # Граф kNN и λ2
    A_knn = knn_affinity_from_C(C, k=params['k'])
    L_knn = laplacian_from_affinity(A_knn)
    lam2_knn = fiedler_lambda(L_knn)

    # Пороговый граф по перцентилю |C| (для компонент)
    p_val, thr_val, A_bin, n_comp, labels, sizes = auto_percentile_threshold(
        C, target_range=(2,6), start=85, minp=50, maxp=98, step=3
    )

    # Лапласиан по |C| (как «глобальная жёсткость»)
    L_absC = laplacian_from_affinity(np.abs(C))
    lam2_absC = fiedler_lambda(L_absC)

    # --- Вывод метрик в терминал ---
    print("\n--- Глобальные метрики ---")
    print(f"Второе собственное λ2 (kNN-граф): {lam2_knn:.4f}  → чем больше, тем связнее граф")
    print(f"Второе собственное λ2 (аффинность |C|): {lam2_absC:.4f}")
    print(f"Перцентиль p={p_val}% (порог={thr_val:.3f}) → число компонент: {n_comp}")
    print(f"Размеры компонент: {sizes.tolist()}")

    # Метрики по кластерам (на основе компонент перцентильного графа)
    clusters = cluster_metrics(C, labels)
    print("\n--- Метрики кластеров (по пороговому графу) ---")
    for cl in clusters:
        print(f"Кластер {cl['id']:>2d}: "
              f"размер={cl['size']:>3d}, "
              f"масса={cl['mass']:.3f}, "
              f"связи-наружу={cl['charge']:.3f}, "
              f"анизо(spin_like)={cl['spin_like']:.5f}, "
              f"λ2-внутр={cl['lambda2']:.4f}")

    # --- Изображения ---

    # (1) Теплокарта |C|
    plt.figure()
    plt.imshow(np.abs(C), interpolation="nearest")
    plt.title("|C| — теплокарта (режим «Суперкластер»)")
    plt.xlabel("Узел")
    plt.ylabel("Узел")
    plt.colorbar(label="|C|")
    fig1_path = os.path.join(outdir, "supercluster_heatmap_absC.png")
    plt.savefig(fig1_path, dpi=180, bbox_inches="tight")
    plt.close()

    # (2) PCA + рёбра kNN
    coords = pca_2d_from_rows(C)
    plt.figure()
    rows, cols = np.where(A_knn > 0)
    for i, j in zip(rows, cols):
        if i < j:
            plt.plot([coords[i,0], coords[j,0]], [coords[i,1], coords[j,1]], alpha=0.3)
    for u in range(n_comp):
        idx = np.where(labels == u)[0]
        plt.scatter(coords[idx,0], coords[idx,1], label=f"комп. {u} (n={len(idx)})", s=20)
    plt.title("PCA узлов C + рёбра kNN")
    plt.xlabel("ПК1")
    plt.ylabel("ПК2")
    plt.legend(loc="best", fontsize=8)
    fig2_path = os.path.join(outdir, "supercluster_pca_knn.png")
    plt.savefig(fig2_path, dpi=180, bbox_inches="tight")
    plt.close()

    # (3) Спектры лапласианов (kNN и |C|)
    def spectrum(L):
        vals = np.linalg.eigvalsh(L)
        return np.sort(np.real(vals))

    spec_knn = spectrum(L_knn)
    spec_absC = spectrum(L_absC)

    plt.figure()
    plt.plot(spec_knn, marker='o', linestyle='None', label="спектр L (kNN)")
    plt.plot(spec_absC, marker='x', linestyle='None', label="спектр L (|C|)")
    plt.title("Спектры лапласианов")
    plt.xlabel("индекс собственного значения")
    plt.ylabel("значение")
    plt.legend()
    fig3_path = os.path.join(outdir, "supercluster_laplacian_spectra.png")
    plt.savefig(fig3_path, dpi=180, bbox_inches="tight")
    plt.close()

    print("\nСохранено:")
    print(f"- Теплокарта: {fig1_path}")
    print(f"- PCA + kNN: {fig2_path}")
    print(f"- Спектры:   {fig3_path}")

if __name__ == "__main__":
    run_supercluster_preset()

# --------------------------------------------------------------
# Симуляция 4.2c: «Сбалансированная многокластерность (спектральная)»
# --------------------------------------------------------------
# Честная часть: kNN по |C| + спектральная кластеризация Лапласиана.
# K выбирается АВТОМАТИЧЕСКИ по eigengap нормализованного лапласиана.
# Визуализация: кластерный layout (только для картинки) — группы узлов
# компактны и разведены друг от друга; можно управлять радиусом/джиттером.
# --------------------------------------------------------------

import os
import numpy as np
import matplotlib.pyplot as plt

# ---------- утилиты ----------
def ensure_outdir(path=""):
    script_dir = os.path.dirname(os.path.abspath(__file__))
    outdir = os.path.join(script_dir, path)
    os.makedirs(outdir, exist_ok=True)
    return outdir

def set_seed(seed=2024):
    np.random.seed(seed)

# ---------- генерация W с мягкой блочностью ----------
def build_W_grain(N, xi=3.0, aniso=0.15, n_groups=5, grain_strength=0.30, inter_atten=0.55):
    i = np.arange(N)[:, None]
    j = np.arange(N)[None, :]
    D = np.abs(i - j).astype(float)
    W = np.exp(-D / xi)
    np.fill_diagonal(W, 0.0)

    if aniso > 0:
        W += aniso * np.exp(-np.abs((i - j) - N/3) / (max(xi/2,1e-6)))
        W += aniso * np.exp(-np.abs((i - j) + N/3) / (max(xi/2,1e-6)))
        np.fill_diagonal(W, 0.0)

    sizes = [N // n_groups + (1 if r < (N % n_groups) else 0) for r in range(n_groups)]
    offs = np.cumsum([0] + sizes)
    group_id = np.zeros(N, dtype=int)
    for g in range(n_groups):
        lo, hi = offs[g], offs[g+1]
        group_id[lo:hi] = g
        W[lo:hi, lo:hi] *= (1.0 + grain_strength)

    G = (group_id[:,None] == group_id[None,:]).astype(float)
    W = W * (G + inter_atten*(1.0-G))
    np.fill_diagonal(W, 0.0)

    s = np.max(np.abs(np.linalg.eigvalsh(W)))
    if s > 1e-9:
        W = W / s
    return W

# ---------- динамика и корреляции ----------
def integrate(N=84, T=14.0, dt=0.05, beta=1.0, gamma=0.8, noise=0.10,
              xi=3.0, aniso=0.15, n_groups=5, grain_strength=0.30, inter_atten=0.55, seed=2024):
    set_seed(seed)
    steps = int(T / dt)
    W = build_W_grain(N, xi=xi, aniso=aniso, n_groups=n_groups,
                      grain_strength=grain_strength, inter_atten=inter_atten)
    x = 0.1 * np.random.randn(N)
    X = np.zeros((steps, N))
    sqrt_dt = np.sqrt(dt)
    for t in range(steps):
        x = (1 - gamma*dt) * x + beta*dt * np.tanh(W @ x) + noise * sqrt_dt * np.random.randn(N)
        X[t] = x

    Xc = X - X.mean(axis=0, keepdims=True)
    std = Xc.std(axis=0, keepdims=True) + 1e-12
    Xn = Xc / std
    C = (Xn.T @ Xn) / (Xn.shape[0] - 1)
    C = np.clip(C, -1.0, 1.0)
    return X, C, W

# ---------- граф и лапласиан ----------
def knn_affinity_from_C(C, k=3):
    S = np.abs(C)
    N = S.shape[0]
    D = 1.0 - S
    A = np.zeros_like(S)
    for i in range(N):
        order = np.argsort(D[i] + np.eye(N)[i]*1e9)
        nbrs = order[:k]
        A[i, nbrs] = S[i, nbrs]
    A = np.maximum(A, A.T)
    np.fill_diagonal(A, 0.0)
    return A

def laplacian_from_affinity(A):
    A = np.maximum(A, 0.0)
    np.fill_diagonal(A, 0.0)
    d = A.sum(axis=1)
    dinv2 = 1.0 / np.sqrt(d + 1e-12)
    L = np.eye(A.shape[0]) - (dinv2[:,None] * A) * dinv2[None,:]
    return L

def normalized_laplacian(A):
    d = A.sum(axis=1)
    dinv2 = 1.0/np.sqrt(d + 1e-12)
    return np.eye(A.shape[0]) - (dinv2[:,None]*A)*dinv2[None,:]

def fiedler_lambda(L):
    vals = np.linalg.eigvalsh(L)
    vals = np.sort(np.real(vals))
    return float(vals[1]) if len(vals) >= 2 else float(vals[0])

def spectrum(L):
    vals = np.linalg.eigvalsh(L)
    return np.sort(np.real(vals))

def select_K_by_eigengap(A, Kmin=2, Kmax=10):
    """
    Автовыбор числа кластеров по eigengap нормализованного лапласиана.
    Возвращает K в [Kmin, Kmax], где разрыв λ_k+1 - λ_k максимален.
    """
    Lsym = normalized_laplacian(A)
    vals = np.linalg.eigvalsh(Lsym)
    vals = np.sort(np.real(vals))
    gaps = []
    cand_K = []
    for K in range(Kmin, min(Kmax, len(vals)-1) + 1):
        gaps.append(vals[K] - vals[K-1])
        cand_K.append(K)
    if not cand_K:
        return max(2, min(5, len(vals)))
    return int(cand_K[int(np.argmax(gaps))])

# ---------- спектральная кластеризация ----------
def spectral_kmeans(A, K, iters=30, seed=123):
    d = A.sum(axis=1)
    dinv2 = 1.0/np.sqrt(d+1e-12)
    Lsym = np.eye(A.shape[0]) - (dinv2[:,None]*A)*dinv2[None,:]
    vals, vecs = np.linalg.eigh(Lsym)
    X = vecs[:, :K]
    X = X / (np.linalg.norm(X, axis=1, keepdims=True) + 1e-12)
    rng = np.random.default_rng(seed)
    N = X.shape[0]
    centers = X[rng.choice(N, K, replace=False)]
    labels = np.zeros(N, dtype=int)
    for _ in range(iters):
        dists = ((X[:,None,:] - centers[None,:,:])**2).sum(axis=2)
        labels = dists.argmin(axis=1)
        new_centers = np.stack([X[labels==k].mean(axis=0) if np.any(labels==k) else centers[k] for k in range(K)], axis=0)
        if np.allclose(new_centers, centers):
            break
        centers = new_centers
    return labels

def balance_score(sizes, N):
    sizes = np.array(sizes, dtype=float)
    std = sizes.std()
    max_frac = sizes.max()/N
    return - (std + 10.0*max_frac)

# ---------- метрики по кластерам ----------
def cluster_metrics(C, labels):
    uniq = np.unique(labels)
    out = []
    for u in uniq:
        idx = np.where(labels == u)[0]
        n = len(idx)
        Cin = C[np.ix_(idx, idx)]
        mass = float(np.sum(np.abs(Cin)**2))
        Ain = np.abs(Cin)
        Lin = laplacian_from_affinity(Ain)
        lam2 = fiedler_lambda(Lin)
        out.append(dict(id=int(u), size=int(n), mass=mass, lambda2=lam2))
    out.sort(key=lambda d: d["size"], reverse=True)
    return out

# ---------- базовое 2D-встраивание (опционально) ----------
def embed_2d_from_rows(M, seed=42):
    X = M - M.mean(axis=0, keepdims=True)
    try:
        import umap
        reducer = umap.UMAP(n_neighbors=25, min_dist=0.15, random_state=seed, metric="euclidean")
        Z = reducer.fit_transform(X)
    except Exception:
        try:
            from sklearn.manifold import TSNE
            Z = TSNE(n_components=2, perplexity=30, learning_rate="auto", init="pca",
                     random_state=seed).fit_transform(X)
        except Exception:
            U, S, Vt = np.linalg.svd(X, full_matrices=False)
            Z = X @ Vt[:2].T
    return Z

# ---------- ТОЛЬКО ДЛЯ ВИЗУАЛИЗАЦИИ: кластерный layout ----------
def spectral_embed_2d(A):
    d = A.sum(axis=1)
    dinv2 = 1.0/np.sqrt(d+1e-12)
    Lsym = np.eye(A.shape[0]) - (dinv2[:,None]*A)*dinv2[None,:]
    vals, vecs = np.linalg.eigh(Lsym)
    if vecs.shape[1] >= 3:
        X = vecs[:, 1:3]
    else:
        X = vecs[:, :2]
    X = X / (np.linalg.norm(X, axis=1, keepdims=True) + 1e-12)
    X = X - X.mean(axis=0, keepdims=True)
    return X

def clustered_layout(A, labels, gap=7.0, cluster_scale=3.0,
                     jitter_std=1.12, radial_noise=1.10,
                     anisotropy=0.20, seed=2024):
    """
    Делает координаты ТОЛЬКО для визуализации:
      • кластеры разводит по окружности (gap);
      • внутри кластера точки образуют облако (джиттер, радиальная «лохматость»,
        слабая анизотропия), масштаб контролируется cluster_scale.
    """
    rng = np.random.default_rng(seed)
    labels = np.asarray(labels)
    uniq = np.unique(labels)
    K = len(uniq)
    N = A.shape[0]

    # центры кластеров на окружности
    angles = np.linspace(0, 2*np.pi, K, endpoint=False)
    centers = np.c_[np.cos(angles), np.sin(angles)] * gap

    coords = np.zeros((N, 2))
    for ci, u in enumerate(uniq):
        idx = np.where(labels == u)[0]
        if len(idx) == 1:
            # одиночная точка — прямо в центр кластера
            coords[idx[0]] = centers[ci] + rng.normal(0, jitter_std, size=2)
            continue

        # компактная форма кластера из его внутрикластерного графа
        A_sub = A[np.ix_(idx, idx)]
        X = spectral_embed_2d(A_sub)                 # (m x 2)

        # нормируем на радиус 1
        rmax = np.max(np.linalg.norm(X, axis=1)) + 1e-9
        X = X / rmax

        # 1) радиальная «лохматость»: растягиваем/сжимаем каждую точку по её радиусу
        radii = np.linalg.norm(X, axis=1, keepdims=True)
        rad_scale = 1.0 + radial_noise * rng.normal(0, 1, size=(len(idx), 1))
        X = X * (rad_scale * (0.85 + 0.15*radii))    # чуть сильнее у внешних

        # 2) локальная случайная анизотропия и поворот
        theta = rng.uniform(0, 2*np.pi)
        c, s = np.cos(theta), np.sin(theta)
        R = np.array([[c, -s], [s, c]])             # поворот
        # лёгкое неравное масштабирование по осям
        sx = 1.0 + anisotropy * rng.normal(0, 0.6)
        sy = 1.0 + anisotropy * rng.normal(0, 0.6)
        S = np.diag([sx, sy])
        X = X @ (R @ S)

        # 3) мелкий изотропный джиттер
        X += rng.normal(0, jitter_std, size=X.shape)

        # держим кластер компактным и пропорциональным размеру
        rmax = np.max(np.linalg.norm(X, axis=1)) + 1e-9
        X = X / rmax
        X = X / (1.0 + 0.006*len(idx))
        X = X * cluster_scale

        # перенос в центр кластера на окружности
        coords[idx] = X + centers[ci]

    return coords

# ---------- основной пресет ----------
def run_balanced_multicluster_preset():
    outdir = ensure_outdir("")

    params = dict(
        N=300, T=14.0, dt=0.05,
        beta=1.0, gamma=0.8, noise=0.10,
        xi=3.0, aniso=0.15,
        n_groups=5, grain_strength=0.30, inter_atten=0.55,
        k=3,                # kNN (локальная аффинность)
        seed=2024
    )

    print("=== Симуляция: «Сбалансированная многокластерность» (спектральная) ===")
    for k,v in params.items():
        print(f"{k} = {v}")

    # (1) интеграция и корреляции
    X, C, W = integrate(N=params['N'], T=params['T'], dt=params['dt'], beta=params['beta'],
                        gamma=params['gamma'], noise=params['noise'], xi=params['xi'],
                        aniso=params['aniso'], n_groups=params['n_groups'],
                        grain_strength=params['grain_strength'], inter_atten=params['inter_atten'],
                        seed=params['seed'])

    # (2) kNN-граф и метрики связности
    A_knn = knn_affinity_from_C(C, k=params['k'])
    L_knn = laplacian_from_affinity(A_knn)
    lam2_knn = fiedler_lambda(L_knn)
    L_absC = laplacian_from_affinity(np.abs(C))
    lam2_absC = fiedler_lambda(L_absC)

    # (3) авт. выбор числа кластеров по eigengap
    K = select_K_by_eigengap(A_knn, Kmin=2, Kmax=10)
    labels = spectral_kmeans(A_knn, K, iters=40, seed=params['seed'])
    sizes = np.array(sorted([int(np.sum(labels==u)) for u in range(K)], reverse=True))
    n_comp = len(sizes)

    print("\n--- Глобальные метрики ---")
    print(f"λ2 (kNN): {lam2_knn:.4f}  → ближе к 0 — легче распадается")
    print(f"λ2 (|C|): {lam2_absC:.4f}")
    print(f"Автовыбранное число кластеров K = {K} (eigengap) → размеры: {sizes.tolist()}")

    # (4) кластерные метрики
    clusters = cluster_metrics(C, labels)
    print("\n--- Метрики кластеров (спектральное разбиение) ---")
    for cl in clusters:
        print(f"Кластер {cl['id']:>2d}: размер={cl['size']:>3d}, масса={cl['mass']:.1f}, λ2-внутр={cl['lambda2']:.4f}")

    # --- рисунки ---
    # (1) теплокарта |C|
    plt.figure()
    plt.imshow(np.abs(C), origin='lower', interpolation='nearest')
    plt.title("Теплокарта |C|")
    plt.colorbar(fraction=0.046, pad=0.04)
    fig1_path = os.path.join(outdir, "balanced_multicluster_heatmap.png")
    plt.savefig(fig1_path, dpi=180, bbox_inches="tight")
    plt.close()

    # (2) кластерная раскладка узлов + рёбра kNN (только визуализация)
    coords = clustered_layout(A_knn, labels, gap=7.0, cluster_scale=3.0)
    plt.figure()
    N = params['N']
    for i in range(N):
        for j in range(i+1, N):
            if A_knn[i,j] <= 0:
                continue
            same = (labels[i] == labels[j])
            a = float(max(0.03, min(0.6, A_knn[i,j])))
            if not same:
                a *= 0.12
            plt.plot([coords[i,0], coords[j,0]], [coords[i,1], coords[j,1]], alpha=a)
    for u in range(n_comp):
        idx = np.where(labels == u)[0]
        plt.scatter(coords[idx,0], coords[idx,1], label=f"кластер {u} (n={len(idx)})", s=28)
    plt.title("Кластерная раскладка узлов C + рёбра kNN")
    plt.xlabel("ось 1"); plt.ylabel("ось 2"); plt.legend(loc="best", fontsize=8)
    fig2_path = os.path.join(outdir, "balanced_multicluster_pca_knn.png")
    plt.savefig(fig2_path, dpi=180, bbox_inches="tight")
    plt.close()

    # (3) спектры лапласианов
    spec_knn = spectrum(L_knn)
    spec_absC = spectrum(L_absC)
    plt.figure()
    plt.plot(spec_knn, marker='o', linestyle='None', label="спектр L (kNN)")
    plt.plot(spec_absC, marker='x', linestyle='None', label="спектр L (|C|)")
    plt.title("Спектры лапласианов")
    plt.xlabel("индекс собственного значения"); plt.ylabel("значение"); plt.legend()
    fig3_path = os.path.join(outdir, "balanced_multicluster_laplacian_spectra.png")
    plt.savefig(fig3_path, dpi=180, bbox_inches="tight")
    plt.close()

    print("\nСохранено:")
    print(f"- Теплокарта: {fig1_path}")
    print(f"- Кластерная картинка: {fig2_path}")
    print(f"- Спектры:   {fig3_path}")

if __name__ == "__main__":
    run_balanced_multicluster_preset()

# --------------------------------------------------------------
# Фазовый переход: многокластер → интеграция
# Автоскан по межгрупповой связности (inter_atten = γ) и сохранение
# четырёх графиков: K(γ), Q(γ), λ2(γ), G(γ).
# Никаких CSV — только картинки в ./out рядом со скриптом.
# --------------------------------------------------------------

import os
import argparse
import numpy as np
import matplotlib.pyplot as plt

# ---------- утилиты ----------
def ensure_outdir(path=""):
    script_dir = os.path.dirname(os.path.abspath(__file__))
    outdir = os.path.join(script_dir, path)
    os.makedirs(outdir, exist_ok=True)
    return outdir

def set_seed(seed=2024):
    np.random.seed(seed)

# ---------- генерация W с мягкой блочностью ----------
def build_W_grain(N, xi=3.0, aniso=0.15, n_groups=5, grain_strength=0.30, inter_atten=0.55):
    i = np.arange(N)[:, None]
    j = np.arange(N)[None, :]
    D = np.abs(i - j).astype(float)
    W = np.exp(-D / xi)
    np.fill_diagonal(W, 0.0)

    if aniso > 0:
        W += aniso * np.exp(-np.abs((i - j) - N/3) / (max(xi/2,1e-6)))
        W += aniso * np.exp(-np.abs((i - j) + N/3) / (max(xi/2,1e-6)))
        np.fill_diagonal(W, 0.0)

    sizes = [N // n_groups + (1 if r < (N % n_groups) else 0) for r in range(n_groups)]
    offs = np.cumsum([0] + sizes)
    group_id = np.zeros(N, dtype=int)
    for g in range(n_groups):
        lo, hi = offs[g], offs[g+1]
        group_id[lo:hi] = g
        W[lo:hi, lo:hi] *= (1.0 + grain_strength)

    G = (group_id[:,None] == group_id[None,:]).astype(float)
    W = W * (G + inter_atten*(1.0-G))
    np.fill_diagonal(W, 0.0)

    s = np.max(np.abs(np.linalg.eigvalsh(W)))
    if s > 1e-9:
        W = W / s
    return W

# ---------- динамика и корреляции ----------
def integrate(N=300, T=14.0, dt=0.05, beta=1.0, gamma=0.8, noise=0.10,
              xi=3.0, aniso=0.15, n_groups=5, grain_strength=0.30, inter_atten=0.55, seed=2024):
    set_seed(seed)
    steps = int(T / dt)
    W = build_W_grain(N, xi=xi, aniso=aniso, n_groups=n_groups,
                      grain_strength=grain_strength, inter_atten=inter_atten)
    x = 0.1 * np.random.randn(N)
    X = np.zeros((steps, N))
    sqrt_dt = np.sqrt(dt)
    for t in range(steps):
        x = (1 - gamma*dt) * x + beta*dt * np.tanh(W @ x) + noise * sqrt_dt * np.random.randn(N)
        X[t] = x

    Xc = X - X.mean(axis=0, keepdims=True)
    std = Xc.std(axis=0, keepdims=True) + 1e-12
    Xn = Xc / std
    C = (Xn.T @ Xn) / (Xn.shape[0] - 1)
    C = np.clip(C, -1.0, 1.0)
    return X, C, W

# ---------- граф и лапласиан ----------
def knn_affinity_from_C(C, k=3):
    S = np.abs(C)
    N = S.shape[0]
    D = 1.0 - S
    A = np.zeros_like(S)
    for i in range(N):
        order = np.argsort(D[i] + np.eye(N)[i]*1e9)
        nbrs = order[:k]
        A[i, nbrs] = S[i, nbrs]
    A = np.maximum(A, A.T)
    np.fill_diagonal(A, 0.0)
    return A

def laplacian_from_affinity(A):
    A = np.maximum(A, 0.0)
    np.fill_diagonal(A, 0.0)
    d = A.sum(axis=1)
    dinv2 = 1.0 / np.sqrt(d + 1e-12)
    L = np.eye(A.shape[0]) - (dinv2[:,None] * A) * dinv2[None,:]
    return L

def normalized_laplacian(A):
    d = A.sum(axis=1)
    dinv2 = 1.0/np.sqrt(d + 1e-12)
    return np.eye(A.shape[0]) - (dinv2[:,None]*A)*dinv2[None,:]

def fiedler_lambda(L):
    vals = np.linalg.eigvalsh(L)
    vals = np.sort(np.real(vals))
    return float(vals[1]) if len(vals) >= 2 else float(vals[0])

# ---------- спектральная кластеризация ----------
def spectral_kmeans(A, K, iters=30, seed=123):
    d = A.sum(axis=1)
    dinv2 = 1.0/np.sqrt(d+1e-12)
    Lsym = np.eye(A.shape[0]) - (dinv2[:,None]*A)*dinv2[None,:]
    vals, vecs = np.linalg.eigh(Lsym)
    X = vecs[:, :K]
    X = X / (np.linalg.norm(X, axis=1, keepdims=True) + 1e-12)
    rng = np.random.default_rng(seed)
    N = X.shape[0]
    centers = X[rng.choice(N, K, replace=False)]
    labels = np.zeros(N, dtype=int)
    for _ in range(iters):
        dists = ((X[:,None,:] - centers[None,:,:])**2).sum(axis=2)
        labels = dists.argmin(axis=1)
        new_centers = np.stack([X[labels==k].mean(axis=0) if np.any(labels==k) else centers[k] for k in range(K)], axis=0)
        if np.allclose(new_centers, centers):
            break
        centers = new_centers
    return labels

def select_K_by_eigengap(A, Kmin=2, Kmax=10):
    Lsym = normalized_laplacian(A)
    vals = np.linalg.eigvalsh(Lsym)
    vals = np.sort(np.real(vals))
    gaps = []
    cand_K = []
    for K in range(Kmin, min(Kmax, len(vals)-1) + 1):
        gaps.append(vals[K] - vals[K-1])
        cand_K.append(K)
    if not cand_K:
        return max(2, min(5, len(vals)))
    return int(cand_K[int(np.argmax(gaps))])

# ---------- метрики ----------
def modularity_Q(A, labels):
    """Взвешенная модульность Ньюмана–Гирвана (векторизованно)."""
    A = np.maximum(A, 0.0)
    m2 = A.sum()  # = 2m
    if m2 <= 1e-12:
        return 0.0
    k = A.sum(axis=1)
    labels = np.asarray(labels)
    uniq = np.unique(labels)
    Q = 0.0
    for s in uniq:
        idx = (labels == s)
        e_s = float(A[np.ix_(idx, idx)].sum()) / m2
        a_s = float(k[idx].sum()) / m2
        Q += e_s - a_s*a_s
    return float(Q)

def giant_component_fraction_from_absC(C, percentile=88.0):
    """Доля узлов в крупнейшей компоненте порогового графа по |C| (вне диагонали)."""
    N = C.shape[0]
    A = np.zeros((N,N), dtype=int)
    S = np.abs(C).copy()
    triu = np.triu_indices(N, 1)
    vals = S[triu]
    thr = float(np.percentile(vals, percentile))
    A[(S >= thr) & (~np.eye(N, dtype=bool))] = 1
    A = np.maximum(A, A.T)
    seen = np.zeros(N, dtype=bool)
    best = 0
    for s in range(N):
        if seen[s]: continue
        q = [s]; seen[s]=True; cnt=1
        while q:
            u = q.pop()
            neigh = np.where(A[u]>0)[0]
            for v in neigh:
                if not seen[v]:
                    seen[v]=True
                    q.append(v)
                    cnt += 1
        if cnt > best: best = cnt
    return best/float(N)

# ---------- основной скан ----------
def run_phase_transition_sweep(
    inter_atten_values=None, percentile=88.0, T=14.0, N=300,
    verbose=False,
    xi=3.0, aniso=0.15, n_groups=5, grain_strength=0.30, k=3,
    beta=1.0, gamma=0.8, noise=0.10, seed=2024
):
    outdir = ensure_outdir("")
    print("=== Фазовый скан: многокластер → интеграция ===")
    print(f"N={N}, T={T}, k={k}, percentile={percentile}, seed={seed}")
    print(f"xi={xi}, aniso={aniso}, n_groups={n_groups}, grain_strength={grain_strength}")
    if inter_atten_values is None:
        inter_atten_values = np.linspace(0.35, 0.95, 12)

    Ks, L2s, L2abs, Qs, Gs = [], [], [], [], []
    xs = []

    for ia in inter_atten_values:
        if verbose:
            print(f"γ={ia:.3f} → интеграция и метрики...")
        X,C,W = integrate(N=N, T=T, dt=0.05, beta=beta, gamma=gamma, noise=noise,
                          xi=xi, aniso=aniso, n_groups=n_groups, grain_strength=grain_strength,
                          inter_atten=float(ia), seed=seed)
        A = knn_affinity_from_C(C, k=k)
        L_knn = laplacian_from_affinity(A)
        lam2_knn = fiedler_lambda(L_knn)
        L_absC = laplacian_from_affinity(np.abs(C))
        lam2_absC = fiedler_lambda(L_absC)

        K = select_K_by_eigengap(A, Kmin=2, Kmax=10)
        labels = spectral_kmeans(A, K, iters=30, seed=seed)
        Q = modularity_Q(A, labels)
        G = giant_component_fraction_from_absC(C, percentile=percentile)

        xs.append(float(ia)); Ks.append(K); L2s.append(lam2_knn); L2abs.append(lam2_absC); Qs.append(Q); Gs.append(G)
        if verbose:
            print(f"   K={K}, λ2(kNN)={lam2_knn:.4f}, Q={Q:.4f}, G={G:.3f}")

    xs = np.array(xs)

    # Плоты (каждый отдельно, без заданных цветов/стилей)
    plt.figure(); plt.plot(xs, Ks, marker='o'); plt.xlabel("inter_atten γ"); plt.ylabel("K (eigengap)")
    plt.title("Число кластеров K по eigengap vs γ")
    plt.savefig(os.path.join(outdir, "phase_K_vs_gamma.png"), dpi=180, bbox_inches="tight"); plt.close()

    plt.figure(); plt.plot(xs, Qs, marker='o'); plt.xlabel("inter_atten γ"); plt.ylabel("модульность Q")
    plt.title("Модульность Q vs γ")
    plt.savefig(os.path.join(outdir, "phase_Q_vs_gamma.png"), dpi=180, bbox_inches="tight"); plt.close()

    plt.figure(); plt.plot(xs, L2s, marker='o'); plt.xlabel("inter_atten γ"); plt.ylabel("λ2 (kNN)")
    plt.title("Второе собственное λ2(kNN) vs γ")
    plt.savefig(os.path.join(outdir, "phase_lambda2_vs_gamma.png"), dpi=180, bbox_inches="tight"); plt.close()

    plt.figure(); plt.plot(xs, Gs, marker='o'); plt.xlabel("inter_atten γ"); plt.ylabel("доля гигантской компоненты G")
    plt.title(f"Перколяция по |C| (порог p={percentile:.0f}%)")
    plt.savefig(os.path.join(outdir, "phase_G_vs_gamma.png"), dpi=180, bbox_inches="tight"); plt.close()

    print("Готово. Сохранены картинки в:", outdir)

# ---------- CLI ----------
def parse_args():
    ap = argparse.ArgumentParser(description="Фазовый переход: многокластер → интеграция")
    ap.add_argument("--quiet", action="store_true", help="не печатать прогресс")
    ap.add_argument("--gamma-start", type=float, default=0.35, help="начальное γ (inter_atten)")
    ap.add_argument("--gamma-end", type=float, default=0.95, help="конечное γ (inter_atten)")
    ap.add_argument("--points", type=int, default=12, help="число точек скана")
    ap.add_argument("--percentile", type=float, default=88.0, help="перцентиль для порога |C| в G")
    ap.add_argument("--T", type=float, default=14.0, help="время интеграции на точку")
    ap.add_argument("--N", type=int, default=300, help="число узлов")
    ap.add_argument("--seed", type=int, default=2024, help="сид")
    # продвинутые (оставим по умолчанию как в мультикластере)
    ap.add_argument("--xi", type=float, default=3.0)
    ap.add_argument("--aniso", type=float, default=0.15)
    ap.add_argument("--n-groups", type=int, default=5)
    ap.add_argument("--grain-strength", type=float, default=0.30)
    ap.add_argument("--k", type=int, default=3)
    ap.add_argument("--beta", type=float, default=1.0)
    ap.add_argument("--gamma", dest="gamma_decay", type=float, default=0.8)
    ap.add_argument("--noise", type=float, default=0.10)
    return ap.parse_args()

if __name__ == "__main__":
    args = parse_args()
    gammas = np.linspace(args.gamma_start, args.gamma_end, args.points)
    run_phase_transition_sweep(
        inter_atten_values=gammas,
        verbose=not args.quiet,
        percentile=args.percentile,
        T=args.T, N=args.N, xi=args.xi, aniso=args.aniso,
        n_groups=args.n_groups, grain_strength=args.grain_strength, k=args.k,
        beta=args.beta, gamma=args.gamma_decay, noise=args.noise, seed=args.seed
    )

# ---------------------------------------------------------------------------
# • Один запуск без аргументов
# • Строит метрики во времени (по скользящему окну |C|) и 3 снимка эмбеддинга
# • Стабильность:
#     - ограничение авто-K по eigengap: Kmax=4
#     - burn-in по времени (игнорируем самые ранние окна)
#     - медианное сглаживание метрик (окно=3)
# • (Опционально) сохраняет анимацию эмбеддинга в MP4 (FFmpeg) или GIF (fallback)
# • В терминал печатает ключевые метрики и пути к файлам
# ---------------------------------------------------------------------------

import os
import numpy as np
import matplotlib.pyplot as plt
import itertools

# ===== ПРЕСЕТЫ =====
N = 300
T = 14.0
dt = 0.05
beta = 1.0
gamma_decay = 0.8
noise = 0.10
seed = 2024

xi = 3.0
aniso = 0.15
n_groups = 5
grain_strength = 0.30
k = 3

# Выбранное состояние (около порога)
GAMMA_INTER = 0.62

# Скользящее окно и снимки
WIN_SEC = 3.0
STEP_FRAC = 0.01
BURNIN_SEC = 2.0            # игнорируем самые ранние окна
SNAP_FRACS = [0.25, 0.62, 0.95]

# Стабилизация кластеризации
KMAX_EIGENGAP = 4           # авто-K не выше 4
K_VISUAL = 3                # фиксированное K для раскраски снимков

# Рисование рёбер
DRAW_EDGES = True
INTRA_ONLY = True
EDGE_TOP_PERCENT = 85.0

# Анимация
SMOOTH_ALPHA = 0.99   # 0..1; ближе к 1 — плавнее, но «тяжелее»
MAKE_ANIMATION = True
ANIM_FPS = 25
ANIM_INTRA_ONLY = False          # показывать и межкластерные «мосты»
ANIM_EDGE_TOP_PERCENT = 95.0
ANIM_BASENAME = f"time_embed_anim_gamma_{GAMMA_INTER:.2f}"

# ===== Утилиты =====
import itertools

def procrustes_align(Z, Z_ref):
    """Ортогональное выравнивание текущих координат к предыдущим (снимает повороты/отражения)."""
    # центрируем
    Zc = Z - Z.mean(axis=0, keepdims=True)
    Rc = Z_ref - Z_ref.mean(axis=0, keepdims=True)
    # матрица кросс-ковариации и SVD
    M = Zc.T @ Rc
    U, _, Vt = np.linalg.svd(M)
    R = U @ Vt
    return Zc @ R  # выровненные координаты (как Z_ref по ориентации)

def relabel_to_match(prev_labels, new_labels, Z_ref, Z_curr, K=None):
    """Перенумеровывает кластеры текущего кадра так, чтобы соответствовать предыдущему."""
    if K is None:
        K = len(np.unique(prev_labels))
    # центры в прежнем и текущем кадре
    prev_cent = np.stack([Z_ref[prev_labels==k].mean(axis=0) for k in range(K)], axis=0)
    curr_cent = np.stack([Z_curr[new_labels==k].mean(axis=0) for k in range(K)], axis=0)
    best_perm, best_cost = None, np.inf
    for perm in itertools.permutations(range(K)):
        cost = 0.0
        for i, j in enumerate(perm):
            cost += np.sum((curr_cent[j] - prev_cent[i])**2)
        if cost < best_cost:
            best_cost, best_perm = cost, perm
    # применяем лучшую перестановку
    mapping = {j:i for i, j in enumerate(best_perm)}
    remapped = np.array([mapping[l] for l in new_labels], dtype=int)
    return remapped

def ensure_outdir(path=""):
    script_dir = os.path.dirname(os.path.abspath(__file__))
    outdir = os.path.join(script_dir, path)
    os.makedirs(outdir, exist_ok=True)
    return outdir

def set_seed(s):
    np.random.seed(s)

def rolling_median(x, w=3):
    x = np.asarray(x, dtype=float)
    if w <= 1 or x.size == 0:
        return x
    half = w // 2
    # паддинг краёв значениями краёв
    xp = np.pad(x, (half, half), mode='edge')
    out = np.empty_like(x)
    for i in range(x.size):
        out[i] = np.median(xp[i:i+w])
    return out

# ===== Генерация W =====
def build_W_grain(N, xi=xi, aniso=aniso, n_groups=n_groups, grain_strength=grain_strength, inter_atten=0.55):
    i = np.arange(N)[:, None]; j = np.arange(N)[None, :]
    D = np.abs(i - j).astype(float)
    W = np.exp(-D / xi); np.fill_diagonal(W, 0.0)
    if aniso > 0:
        W += aniso * np.exp(-np.abs((i - j) - N/3) / (max(xi/2,1e-6)))
        W += aniso * np.exp(-np.abs((i - j) + N/3) / (max(xi/2,1e-6)))
        np.fill_diagonal(W, 0.0)
    sizes = [N // n_groups + (1 if r < (N % n_groups) else 0) for r in range(n_groups)]
    offs = np.cumsum([0] + sizes)
    group_id = np.zeros(N, dtype=int)
    for g in range(n_groups):
        lo, hi = offs[g], offs[g+1]
        group_id[lo:hi] = g
        W[lo:hi, lo:hi] *= (1.0 + grain_strength)
    Gm = (group_id[:,None] == group_id[None,:]).astype(float)
    W = W * (Gm + inter_atten*(1.0-Gm))
    np.fill_diagonal(W, 0.0)
    s = np.max(np.abs(np.linalg.eigvalsh(W)))
    if s > 1e-9: W = W / s
    return W

# ===== Интеграция =====
def integrate_all(N=N, T=T, dt=dt, beta=beta, gamma=gamma_decay, noise=noise,
                  xi=xi, aniso=aniso, n_groups=n_groups, grain_strength=grain_strength,
                  inter_atten=GAMMA_INTER, seed=seed):
    set_seed(seed)
    steps = int(T/dt)
    W = build_W_grain(N, xi=xi, aniso=aniso, n_groups=n_groups,
                      grain_strength=grain_strength, inter_atten=inter_atten)
    x = 0.1*np.random.randn(N)
    X = np.zeros((steps, N))
    sqrt_dt = np.sqrt(dt)
    for t in range(steps):
        x = (1 - gamma*dt)*x + beta*dt*np.tanh(W @ x) + noise*sqrt_dt*np.random.randn(N)
        X[t] = x
    return X

def corr_from_window(Xseg):
    Xc = Xseg - Xseg.mean(axis=0, keepdims=True)
    std = Xc.std(axis=0, keepdims=True) + 1e-12
    Xn  = Xc / std
    C = (Xn.T @ Xn) / (Xn.shape[0] - 1)
    return np.clip(C, -1.0, 1.0)

# ===== Граф и метрики =====
def knn_affinity_from_C(C, k=k):
    S = np.abs(C); Nloc = S.shape[0]
    D = 1.0 - S; A = np.zeros_like(S)
    for i in range(Nloc):
        order = np.argsort(D[i] + np.eye(Nloc)[i]*1e9)
        nbrs = order[:k]; A[i, nbrs] = S[i, nbrs]
    A = np.maximum(A, A.T); np.fill_diagonal(A, 0.0)
    return A

def normalized_laplacian(A):
    d = A.sum(axis=1)
    dinv2 = 1.0/np.sqrt(d + 1e-12)
    return np.eye(A.shape[0]) - (dinv2[:,None]*A)*dinv2[None,:]

def laplacian_from_affinity(A):
    d = A.sum(axis=1)
    dinv2 = 1.0/np.sqrt(d + 1e-12)
    return np.eye(A.shape[0]) - (dinv2[:,None]*A)*dinv2[None,:]

def fiedler_lambda(L):
    vals = np.linalg.eigvalsh(L)
    vals = np.sort(np.real(vals))
    return float(vals[1]) if len(vals) >= 2 else float(vals[0])

def modularity_Q(A, labels):
    m = A.sum()/2.0 + 1e-12
    d = A.sum(axis=1)
    Q = 0.0
    for c in np.unique(labels):
        idx = np.where(labels==c)[0]
        As = A[np.ix_(idx, idx)].sum()
        ds = d[idx].sum()
        Q += (As/(2*m) - (ds/(2*m))**2)
    return float(Q)

def spectral_kmeans(A, K, iters=30, seed=seed):
    d = A.sum(axis=1)
    dinv2 = 1.0/np.sqrt(d+1e-12)
    Lsym = np.eye(A.shape[0]) - (dinv2[:,None]*A)*dinv2[None,:]
    vals, vecs = np.linalg.eigh(Lsym)
    Xemb = vecs[:, :K]
    Xemb = Xemb / (np.linalg.norm(Xemb, axis=1, keepdims=True) + 1e-12)
    rng = np.random.default_rng(seed)
    Nloc = Xemb.shape[0]
    centers = Xemb[rng.choice(Nloc, K, replace=False)]
    labels = np.zeros(Nloc, dtype=int)
    for _ in range(iters):
        dists = ((Xemb[:,None,:] - centers[None,:,:])**2).sum(axis=2)
        labels = dists.argmin(axis=1)
        new_centers = np.stack([Xemb[labels==k].mean(axis=0) if np.any(labels==k) else centers[k] for k in range(K)], axis=0)
        if np.allclose(new_centers, centers): break
        centers = new_centers
    return labels

def select_K_by_eigengap(A, Kmin=2, Kmax=KMAX_EIGENGAP):
    Lsym = normalized_laplacian(A)
    vals = np.linalg.eigvalsh(Lsym); vals = np.sort(np.real(vals))
    gaps = []; cand = []
    upper = int(min(Kmax, max(Kmin, len(vals)-1)))
    for K in range(Kmin, upper+1):
        gaps.append(vals[K]-vals[K-1]); cand.append(K)
    if not cand: return max(2, min(5, len(vals)-1))
    return int(cand[int(np.argmax(gaps))])

def giant_component_fraction_from_absC(C, percentile=88.0):
    Nloc = C.shape[0]; S = np.abs(C).copy()
    A = np.zeros((Nloc,Nloc), dtype=int)
    triu = np.triu_indices(Nloc, 1)
    thr = float(np.percentile(S[triu], percentile))
    A[(S >= thr) & (~np.eye(Nloc, dtype=bool))] = 1
    A = np.maximum(A, A.T)
    seen = np.zeros(Nloc, dtype=bool); best = 0
    for s in range(Nloc):
        if seen[s]: continue
        q = [s]; seen[s]=True; cnt=1
        while q:
            u = q.pop()
            for v in np.where(A[u]>0)[0]:
                if not seen[v]:
                    seen[v]=True; q.append(v); cnt+=1
        best = max(best, cnt)
    return best/float(Nloc)

# ===== Эмбеддинг (Diffusion Maps) =====
def embed_diffusion(A, dim=2, t=1.0):
    d = A.sum(axis=1) + 1e-12
    P = (A.T / d).T
    vals, vecs = np.linalg.eig(P)
    vals = np.real(vals); vecs = np.real(vecs)
    order = np.argsort(-np.abs(vals))
    vals = vals[order]; vecs = vecs[:, order]
    U = vecs[:, 1:dim+1]
    lam = vals[1:dim+1]
    U = U * (lam[np.newaxis, :]**t)
    U = (U - U.mean(axis=0)) / (U.std(axis=0) + 1e-9)
    return U

# ===== Рисование =====
def plot_embedding(A, Z, labels, title, path, draw_edges=True, intra_only=True, edge_top_percent=85.0,
                   xlim=None, ylim=None):
    plt.figure()
    ax = plt.gca()
    if draw_edges:
        wpos = A[A > 0.0]
        thr = np.percentile(wpos, edge_top_percent) if wpos.size>0 else np.inf
        wmax = wpos.max() if wpos.size>0 else 1.0
        Nloc = A.shape[0]
        for i in range(Nloc):
            for j in range(i+1, Nloc):
                wij = A[i,j]
                if wij <= thr: continue
                if intra_only and labels[i]!=labels[j]: continue
                alpha = max(0.05, min(0.6, float((wij-thr)/(wmax-thr+1e-9)+0.05)))
                plt.plot([Z[i,0], Z[j,0]], [Z[i,1], Z[j,1]], alpha=alpha)
    for u in np.unique(labels):
        idx = np.where(labels==u)[0]
        plt.scatter(Z[idx,0], Z[idx,1], s=22, alpha=0.9, edgecolors='none', label=f"кластер {int(u)} (n={len(idx)})")

    ax.set_aspect('equal', 'box')
    if xlim is not None: ax.set_xlim(*xlim)
    if ylim is not None: ax.set_ylim(*ylim)
    plt.title(title); plt.legend(loc="best", fontsize=8)
    plt.xlabel("ось 1"); plt.ylabel("ось 2")
    plt.savefig(path, dpi=180, bbox_inches="tight"); plt.close()

# ===== Анимация =====
def save_animation(frames, xlim, ylim, outdir, fps=8, intra_only=False, edge_top_percent=92.0, basename="anim"):
    import matplotlib.pyplot as plt
    import matplotlib.animation as animation

    fig, ax = plt.subplots()
    ax.set_aspect('equal', 'box')
    if xlim is not None: ax.set_xlim(*xlim)
    if ylim is not None: ax.set_ylim(*ylim)

    def draw_frame(idx):
        ax.clear()
        ax.set_aspect('equal', 'box')
        if xlim is not None: ax.set_xlim(*xlim)
        if ylim is not None: ax.set_ylim(*ylim)
        t_mid, A, Z, labels = frames[idx]

        # рёбра (только сильные)
        wpos = A[A > 0.0]
        thr = np.percentile(wpos, edge_top_percent) if wpos.size>0 else np.inf
        wmax = wpos.max() if wpos.size>0 else 1.0
        Nloc = A.shape[0]
        for i in range(Nloc):
            for j in range(i+1, Nloc):
                wij = A[i,j]
                if wij <= thr: continue
                if intra_only and labels[i]!=labels[j]: continue
                alpha = max(0.05, min(0.6, float((wij-thr)/(wmax-thr+1e-9)+0.05)))
                ax.plot([Z[i,0], Z[j,0]], [Z[i,1], Z[j,1]], alpha=alpha)

        # точки
        for u in np.unique(labels):
            idxs = np.where(labels==u)[0]
            ax.scatter(Z[idxs,0], Z[idxs,1], s=22, alpha=0.9, edgecolors='none', label=f"кл {int(u)} (n={len(idxs)})")

        ax.legend(loc="best", fontsize=8)
        ax.set_title(f"Эмерджентное пространство (t={t_mid:.2f}, γ={GAMMA_INTER:.2f})")
        ax.set_xlabel("ось 1"); ax.set_ylabel("ось 2")
        return []

    ani = animation.FuncAnimation(fig, draw_frame, frames=len(frames), interval=1000/fps, blit=False)

    # Сохраняем GIF через PillowWriter (без imageio/ffmpeg)
    try:
        from matplotlib.animation import PillowWriter
        gif_path = os.path.join(outdir, f"{basename}.gif")
        ani.save(gif_path, writer=PillowWriter(fps=fps))
        plt.close(fig)
        return gif_path
    except Exception:
        plt.close(fig)
        return None


# ===== Основной сценарий =====
def run():
    outdir = ensure_outdir("")
    print("=== Динамика во времени (стабилизировано) ===")
    print(f"N={N}, T={T}, k={k}, seed={seed}, γ={GAMMA_INTER}")
    print(f"xi={xi}, aniso={aniso}, n_groups={n_groups}, grain_strength={grain_strength}")
    print(f"окно={WIN_SEC}s, шаг={STEP_FRAC*100:.0f}%, burn-in={BURNIN_SEC}s")

    # Интеграция полного ряда
    X = integrate_all(inter_atten=GAMMA_INTER)

    steps = X.shape[0]
    WINS = int(WIN_SEC/dt)
    STRIDE = max(1, int(WINS*STEP_FRAC))
    BURN = int(BURNIN_SEC/dt)

    times = []
    K_series, Q_series, L2_series, G_series = [], [], [], []
    windows = []

    # Скользящие окна (с учётом burn-in)
    for t1 in range(BURN + WINS, steps+1, STRIDE):
        t0 = t1 - WINS
        Xseg = X[t0:t1]
        C = corr_from_window(Xseg)
        A = knn_affinity_from_C(C, k=k)

        Kauto = select_K_by_eigengap(A, Kmin=2, Kmax=KMAX_EIGENGAP)
        labels = spectral_kmeans(A, Kauto, iters=30, seed=seed)
        Q = modularity_Q(A, labels)
        L2 = fiedler_lambda(laplacian_from_affinity(A))
        G = giant_component_fraction_from_absC(C, percentile=88.0)

        t_mid = (t0 + t1)/2 * dt
        times.append(t_mid); K_series.append(Kauto); Q_series.append(Q); L2_series.append(L2); G_series.append(G)

        windows.append((t_mid, C, A))

    # Сглаживание метрик
    K_sm = rolling_median(K_series, 3)
    Q_sm = rolling_median(Q_series, 3)
    L2_sm = rolling_median(L2_series, 3)
    G_sm = rolling_median(G_series, 3)

    # --- Графики метрик (рисуем только сглаженные) ---
    def lineplot(t, y, ylabel, fname):
        plt.figure()
        plt.plot(t, y, marker='o', lw=1.5)
        plt.xlabel("время"); plt.ylabel(ylabel); plt.title(f"{ylabel} во времени (γ={GAMMA_INTER})")
        plt.grid(alpha=0.25); plt.savefig(os.path.join(outdir, fname), dpi=180, bbox_inches="tight"); plt.close()

    lineplot(times, K_sm, "K (eigengap, сглаж.)", f"time_K_vs_t_gamma_{GAMMA_INTER:.2f}.png")
    lineplot(times, Q_sm, "модульность Q (сглаж.)", f"time_Q_vs_t_gamma_{GAMMA_INTER:.2f}.png")
    lineplot(times, L2_sm, "λ₂(kNN) (сглаж.)", f"time_lambda2_vs_t_gamma_{GAMMA_INTER:.2f}.png")
    lineplot(times, G_sm, "доля G (p=88%, сглаж.)", f"time_G_vs_t_gamma_{GAMMA_INTER:.2f}.png")

    # --- Снимки эмбеддинга (3 шт.) ---
    mids = np.array([w[0] for w in windows])
    picks = []
    for frac in SNAP_FRACS:
        target = BURNIN_SEC + frac * (T - BURNIN_SEC)  # выбираем после burn-in
        if mids.size == 0:
            continue
        idx = int(np.argmin(np.abs(mids - target)))
        picks.append(idx)
    picks = sorted(set(picks))

    # Единые оси для снимков
    Zs = []
    EMB = []
    for idx in picks:
        _, C, A = windows[idx]
        Z = embed_diffusion(A, dim=2, t=1.0)
        labels_vis = spectral_kmeans(A, K_VISUAL, iters=40, seed=seed)
        EMB.append((idx, A, Z, labels_vis))
        Zs.append(Z)
    if Zs:
        allZ = np.vstack(Zs)
        xmin, xmax = float(allZ[:,0].min()), float(allZ[:,0].max())
        ymin, ymax = float(allZ[:,1].min()), float(allZ[:,1].max())
        dx, dy = xmax-xmin, ymax-ymin
        padx, pady = 0.08*dx, 0.08*dy
        xlim = (xmin-padx, xmax+padx); ylim = (ymin-pady, ymax+pady)
    else:
        xlim = ylim = None

    for idx, A, Z, labels_vis in EMB:
        t_mid, C, _ = windows[idx]
        # метрики по этому окну
        Kauto = select_K_by_eigengap(A, Kmin=2, Kmax=KMAX_EIGENGAP)
        L2 = fiedler_lambda(laplacian_from_affinity(A))
        G = giant_component_fraction_from_absC(C, percentile=88.0)
        sizes_vis = [int(np.sum(labels_vis==u)) for u in np.unique(labels_vis)]
        fname = os.path.join(outdir, f"time_embed_t_{t_mid:.2f}_gamma_{GAMMA_INTER:.2f}.png")
        title = (f"Эмерджентное пространство (DiffMaps), t={t_mid:.2f}, γ={GAMMA_INTER:.2f}; "
             f"K_auto={Kauto}, K_vis={K_VISUAL}, λ2={L2:.4f}, G={G:.3f}, "
             f"размеры(K_vis)={sorted(sizes_vis, reverse=True)}")
        plot_embedding(A, Z, labels_vis, title, fname,
                       draw_edges=DRAW_EDGES, intra_only=INTRA_ONLY, edge_top_percent=EDGE_TOP_PERCENT,
                       xlim=xlim, ylim=ylim)
        print(f"   [SNAP] t={t_mid:.2f}: K_auto={Kauto}, λ2={L2:.4f}, G={G:.3f} → {fname}")

    # --- Анимация ---
    # --- Анимация ---
    if MAKE_ANIMATION and len(windows) >= 2:
        # Собираем кадры с выравниванием и EMA-сглаживанием координат
        frames = []
        allZ_anim = []

        Z_ref = None
        labels_ref = None
        for (t_mid, C, A) in windows:
            # эмбеддинг текущего окна
            Z_raw = embed_diffusion(A, dim=2, t=1.0)
            # окраска (фикс. K_VISUAL) — затем стабилизируем перестановку меток
            labels_curr = spectral_kmeans(A, K_VISUAL, iters=30, seed=seed)

            if Z_ref is None:
                # первый кадр — берём как есть
                Z_smooth = Z_raw
            else:
                # 1) выравниваем ориентацию к предыдущему
                Z_aligned = procrustes_align(Z_raw, Z_ref)
                # 2) согласовываем метки по кластерам
                labels_curr = relabel_to_match(labels_ref, labels_curr, Z_ref, Z_aligned, K=K_VISUAL)
                # 3) EMA-сглаживание координат
                Z_smooth = SMOOTH_ALPHA * Z_ref + (1.0 - SMOOTH_ALPHA) * Z_aligned

            frames.append((t_mid, A, Z_smooth, labels_curr))
            allZ_anim.append(Z_smooth)
            Z_ref = Z_smooth
            labels_ref = labels_curr

        # Единые оси для всей анимации
        allZ_anim = np.vstack(allZ_anim)
        xmin, xmax = float(allZ_anim[:,0].min()), float(allZ_anim[:,0].max())
        ymin, ymax = float(allZ_anim[:,1].min()), float(allZ_anim[:,1].max())
        dx, dy = xmax - xmin, ymax - ymin
        padx, pady = 0.08*dx, 0.08*dy
        xlim_anim = (xmin - padx, xmax + padx)
        ylim_anim = (ymin - pady, ymax + pady)

        anim_path = save_animation(frames, xlim_anim, ylim_anim, outdir,
                                fps=ANIM_FPS, intra_only=ANIM_INTRA_ONLY,
                                edge_top_percent=ANIM_EDGE_TOP_PERCENT,
                                basename=ANIM_BASENAME)
        if anim_path:
            print("   [ANIM] Сохранена анимация:", anim_path)
        else:
            print("   [ANIM] Не удалось сохранить анимацию (нет подходящего writer'а)")


    print("Готово. Картинки в:", outdir)

if __name__ == "__main__":
    run()

# --------------------------------------------------------------
# • Три состояния γ: до / около / после порога (один запуск, без аргументов)
# • Терминал: K_auto, K_vis, размеры, λ2(kNN), G(p=88%), пути к PNG
# • Визуализация:
#   - авто-K ограничен MAX_K=4 (стабильнее)
#   - раскраска по фикс. K_VISUAL=3 (единый вид)
#   - рисуем только сильные ВНУТРИКЛАСТЕРНЫЕ рёбра
#   - ОДНИ И ТЕ ЖЕ ПРЕДЕЛЫ ОСЕЙ для всех кадров (LOCK_AXES=True)
#   - ДОП. КАДР для γ=0.90 с межкластерными рёбрами (видны «мосты»)
# --------------------------------------------------------------

import os
import numpy as np
import matplotlib.pyplot as plt

# ---------- базовые параметры ----------
N = 300
T = 14.0
dt = 0.05
beta = 1.0
gamma_decay = 0.8
noise = 0.10
seed = 2024

xi = 3.0
aniso = 0.15
n_groups = 5
grain_strength = 0.30
k = 3

MAX_K = 4          # авто-K по eigengap не выше 4
K_VISUAL = 3       # фиксированное K для окраски

# три характерные точки γ
gammas = [0.40, 0.62, 0.90]

# отрисовка рёбер (основные кадры)
DRAW_EDGES = True
INTRA_ONLY = True           # только внутрикластерные рёбра (по K_VISUAL)
EDGE_TOP_PERCENT = 85.0     # верхние 15% по силе

# общий масштаб осей для всех кадров
LOCK_AXES = True
AXIS_PADDING = 0.08         # доля от диапазона

# дополнительный кадр «перемычки» для наглядности интеграции
MAKE_ALT_VIEW = True
ALT_VIEW_GAMMA = 0.90
ALT_INTRA_ONLY = False
ALT_EDGE_TOP_PERCENT = 92.0
ALT_SUFFIX = "_bridges"

# ---------- утилиты ----------
def ensure_outdir(path=""):
    script_dir = os.path.dirname(os.path.abspath(__file__))
    outdir = os.path.join(script_dir, path)
    os.makedirs(outdir, exist_ok=True)
    return outdir

def set_seed(s=seed):
    np.random.seed(s)

# ---------- генерация W ----------
def build_W_grain(N, xi=xi, aniso=aniso, n_groups=n_groups, grain_strength=grain_strength, inter_atten=0.55):
    i = np.arange(N)[:, None]
    j = np.arange(N)[None, :]
    D = np.abs(i - j).astype(float)
    W = np.exp(-D / xi)
    np.fill_diagonal(W, 0.0)

    if aniso > 0:
        W += aniso * np.exp(-np.abs((i - j) - N/3) / (max(xi/2,1e-6)))
        W += aniso * np.exp(-np.abs((i - j) + N/3) / (max(xi/2,1e-6)))
        np.fill_diagonal(W, 0.0)

    sizes = [N // n_groups + (1 if r < (N % n_groups) else 0) for r in range(n_groups)]
    offs = np.cumsum([0] + sizes)
    group_id = np.zeros(N, dtype=int)
    for g in range(n_groups):
        lo, hi = offs[g], offs[g+1]
        group_id[lo:hi] = g
        W[lo:hi, lo:hi] *= (1.0 + grain_strength)

    Gm = (group_id[:,None] == group_id[None,:]).astype(float)
    W = W * (Gm + inter_atten*(1.0-Gm))
    np.fill_diagonal(W, 0.0)

    s = np.max(np.abs(np.linalg.eigvalsh(W)))
    if s > 1e-9:
        W = W / s
    return W

# ---------- динамика и корреляции ----------
def integrate(N=N, T=T, dt=dt, beta=beta, gamma=gamma_decay, noise=noise,
              xi=xi, aniso=aniso, n_groups=n_groups, grain_strength=grain_strength, inter_atten=0.55, seed=seed):
    set_seed(seed)
    steps = int(T / dt)
    W = build_W_grain(N, xi=xi, aniso=aniso, n_groups=n_groups,
                      grain_strength=grain_strength, inter_atten=inter_atten)
    x = 0.1 * np.random.randn(N)
    X = np.zeros((steps, N))
    sqrt_dt = np.sqrt(dt)
    for t in range(steps):
        x = (1 - gamma*dt) * x + beta*dt * np.tanh(W @ x) + noise * sqrt_dt * np.random.randn(N)
        X[t] = x

    Xc = X - X.mean(axis=0, keepdims=True)
    std = Xc.std(axis=0, keepdims=True) + 1e-12
    Xn = Xc / std
    C = (Xn.T @ Xn) / (Xn.shape[0] - 1)
    C = np.clip(C, -1.0, 1.0)
    return X, C, W

# ---------- граф, лапласиан ----------
def knn_affinity_from_C(C, k=k):
    S = np.abs(C)
    Nloc = S.shape[0]
    D = 1.0 - S
    A = np.zeros_like(S)
    for i in range(Nloc):
        order = np.argsort(D[i] + np.eye(Nloc)[i]*1e9)
        nbrs = order[:k]
        A[i, nbrs] = S[i, nbrs]
    A = np.maximum(A, A.T)
    np.fill_diagonal(A, 0.0)
    return A

def laplacian_from_affinity(A):
    d = A.sum(axis=1)
    dinv2 = 1.0 / np.sqrt(d + 1e-12)
    L = np.eye(A.shape[0]) - (dinv2[:,None] * A) * dinv2[None,:]
    return L

def fiedler_lambda(L):
    vals = np.linalg.eigvalsh(L)
    vals = np.sort(np.real(vals))
    return float(vals[1]) if len(vals) >= 2 else float(vals[0])

def normalized_laplacian(A):
    d = A.sum(axis=1)
    dinv2 = 1.0/np.sqrt(d + 1e-12)
    return np.eye(A.shape[0]) - (dinv2[:,None]*A)*dinv2[None,:]

# ---------- кластеризация и выбор K ----------
def spectral_kmeans(A, K, iters=30, seed=seed):
    d = A.sum(axis=1)
    dinv2 = 1.0/np.sqrt(d+1e-12)
    Lsym = np.eye(A.shape[0]) - (dinv2[:,None]*A)*dinv2[None,:]
    vals, vecs = np.linalg.eigh(Lsym)
    Xemb = vecs[:, :K]
    Xemb = Xemb / (np.linalg.norm(Xemb, axis=1, keepdims=True) + 1e-12)
    rng = np.random.default_rng(seed)
    Nloc = Xemb.shape[0]
    centers = Xemb[rng.choice(Nloc, K, replace=False)]
    labels = np.zeros(Nloc, dtype=int)
    for _ in range(iters):
        dists = ((Xemb[:,None,:] - centers[None,:,:])**2).sum(axis=2)
        labels = dists.argmin(axis=1)
        new_centers = np.stack([Xemb[labels==k].mean(axis=0) if np.any(labels==k) else centers[k] for k in range(K)], axis=0)
        if np.allclose(new_centers, centers):
            break
        centers = new_centers
    return labels

def select_K_by_eigengap(A, Kmin=2, Kmax=MAX_K):
    Lsym = normalized_laplacian(A)
    vals = np.linalg.eigvalsh(Lsym)
    vals = np.sort(np.real(vals))
    gaps = []
    cand_K = []
    upper = int(min(Kmax, max(Kmin, len(vals)-1)))
    for Kc in range(Kmin, upper + 1):
        gaps.append(vals[Kc] - vals[Kc-1])
        cand_K.append(Kc)
    if not cand_K:
        return max(2, min(5, len(vals)-1))
    return int(cand_K[int(np.argmax(gaps))])

# ---------- перколяция по |C| ----------
def giant_component_fraction_from_absC(C, percentile=88.0):
    Nloc = C.shape[0]
    A = np.zeros((Nloc,Nloc), dtype=int)
    S = np.abs(C).copy()
    triu = np.triu_indices(Nloc, 1)
    vals = S[triu]
    thr = float(np.percentile(vals, percentile))
    A[(S >= thr) & (~np.eye(Nloc, dtype=bool))] = 1
    A = np.maximum(A, A.T)
    seen = np.zeros(Nloc, dtype=bool)
    best = 0
    for s in range(Nloc):
        if seen[s]: continue
        q = [s]; seen[s]=True; cnt=1
        while q:
            u = q.pop()
            neigh = np.where(A[u]>0)[0]
            for v in neigh:
                if not seen[v]:
                    seen[v]=True
                    q.append(v)
                    cnt += 1
        if cnt > best: best = cnt
    return best/float(Nloc)

# ---------- диффузионные карты ----------
def embed_diffusion(A, dim=2, t=1.0):
    d = A.sum(axis=1) + 1e-12
    P = (A.T / d).T  # D^{-1}A
    vals, vecs = np.linalg.eig(P)
    vals = np.real(vals); vecs = np.real(vecs)
    order = np.argsort(-np.abs(vals))
    vals = vals[order]; vecs = vecs[:, order]
    U = vecs[:, 1:dim+1]           # пропускаем тривиальный λ≈1
    lam = vals[1:dim+1]
    U = U * (lam[np.newaxis, :]**t)
    # стандартизируем по осям (без «окружности»)
    U = (U - U.mean(axis=0)) / (U.std(axis=0) + 1e-9)
    return U

# ---------- рисование ----------
def plot_embedded(A, coords, labels, title, path,
                  draw_edges=True, intra_only=True, edge_top_percent=85.0,
                  xlim=None, ylim=None):
    plt.figure()
    ax = plt.gca()

    if draw_edges:
        wpos = A[A > 0.0]
        thr = np.percentile(wpos, edge_top_percent) if wpos.size > 0 else np.inf
        wmax = wpos.max() if wpos.size > 0 else 1.0
        Nloc = A.shape[0]
        for i in range(Nloc):
            for j in range(i+1, Nloc):
                wij = A[i, j]
                if wij <= thr:
                    continue
                if intra_only and labels[i] != labels[j]:
                    continue
                alpha = max(0.05, min(0.6, float((wij - thr) / (wmax - thr + 1e-9) + 0.05)))
                plt.plot([coords[i,0], coords[j,0]], [coords[i,1], coords[j,1]], alpha=alpha)

    uniq = np.unique(labels)
    for u in uniq:
        idx = np.where(labels==u)[0]
        plt.scatter(coords[idx,0], coords[idx,1], s=22, label=f"кластер {int(u)} (n={len(idx)})")

    ax.set_aspect('equal', 'box')
    if xlim is not None and ylim is not None:
        ax.set_xlim(*xlim); ax.set_ylim(*ylim)
    plt.title(title); plt.legend(loc="best", fontsize=8)
    plt.xlabel("ось 1"); plt.ylabel("ось 2")
    plt.savefig(path, dpi=180, bbox_inches="tight"); plt.close()

# ---------- основной сценарий ----------
def run():
    outdir = ensure_outdir("")
    print("=== Эмерджентное пространство (v3) ===")
    print(f"N={N}, T={T}, k={k}, seed={seed}")
    print(f"xi={xi}, aniso={aniso}, n_groups={n_groups}, grain_strength={grain_strength}")
    print(f"gammas={gammas}")

    # 1) Считаем всё и копим для общих осей
    results = []
    for g in gammas:
        print(f"γ={g:.3f} → интеграция, граф, встраивание...")
        X, C, W = integrate(inter_atten=float(g))
        A = knn_affinity_from_C(C, k=k)
        K_auto = select_K_by_eigengap(A, Kmin=2, Kmax=MAX_K)
        labels_vis = spectral_kmeans(A, K_VISUAL, iters=40, seed=seed)
        sizes_vis = sorted([int(np.sum(labels_vis==u)) for u in range(K_VISUAL)], reverse=True)
        L = laplacian_from_affinity(A)
        lam2 = fiedler_lambda(L)
        G = giant_component_fraction_from_absC(C, percentile=88.0)
        Z = embed_diffusion(A, dim=2, t=1.0)

        results.append(dict(
            gamma=g, A=A, Z=Z, labels=labels_vis,
            K_auto=K_auto, sizes=sizes_vis, lam2=lam2, G=G
        ))

    # 2) Общие пределы осей
    xlim = ylim = None
    if LOCK_AXES:
        allZ = np.vstack([r["Z"] for r in results])
        xmin, xmax = float(allZ[:,0].min()), float(allZ[:,0].max())
        ymin, ymax = float(allZ[:,1].min()), float(allZ[:,1].max())
        dx, dy = xmax - xmin, ymax - ymin
        xpad, ypad = AXIS_PADDING*dx, AXIS_PADDING*dy
        xlim = (xmin - xpad, xmax + xpad)
        ylim = (ymin - ypad, ymax + ypad)
        print(f"Оси зафиксированы: xlim={xlim}, ylim={ylim}")

    # 3) Рисуем кадры
    for r in results:
        g = r["gamma"]; A = r["A"]; Z = r["Z"]; labels = r["labels"]
        K_auto = r["K_auto"]; sizes = r["sizes"]; lam2 = r["lam2"]; G = r["G"]

        fname = os.path.join(outdir, f"emergent_space_gamma_{g:.2f}.png")
        title = (f"Эмерджентное пространство (DiffMaps), γ={g:.2f}; "
                 f"K_auto={K_auto}, K_vis={K_VISUAL}, размеры(K_vis)={sizes}")
        plot_embedded(A, Z, labels, title, fname,
                      draw_edges=DRAW_EDGES, intra_only=INTRA_ONLY, edge_top_percent=EDGE_TOP_PERCENT,
                      xlim=xlim, ylim=ylim)
        print(f"   K_auto={K_auto}, K_vis={K_VISUAL}, размеры={sizes}, λ2(kNN)={lam2:.4f}, G(p=88%)={G:.3f} → {fname}")

        # дополнительный кадр с «перемычками»
        if MAKE_ALT_VIEW and abs(g - ALT_VIEW_GAMMA) < 1e-9:
            fname2 = os.path.join(outdir, f"emergent_space_gamma_{g:.2f}{ALT_SUFFIX}.png")
            plot_embedded(A, Z, labels,
                          title + " (межкластерные рёбра)",
                          fname2,
                          draw_edges=True, intra_only=ALT_INTRA_ONLY, edge_top_percent=ALT_EDGE_TOP_PERCENT,
                          xlim=xlim, ylim=ylim)
            print(f"   → Доп. вид (межкластерные рёбра): {fname2}")

    print("Готово. Картинки в:", outdir)

if __name__ == "__main__":
    run()