benchmark_algorithms.py

"""
===============================================================================
benchmark_algorithms.py
基准对比算法 + 评价指标 + 动态测试函数 (完整高性能版)
===============================================================================
【文件职责】
  - 实现 6 种对比算法 (RI / MMTL / Tr / PPS / KF / SVR)
  - 实现 KEMM-DMOEA 的抽象空间版本 (用于标准测试函数对比)
     → 完整包含 FindBestSol (Process 2) + Transfer (Process 3)
  - 实现 3 种评价指标 (MIGD / SP / MS)
  - 实现 FDA / dMOP 系列动态测试函数

【论文来源】(知识库: MMTL-DMOEA, IEEE TCYB, 2020)
  对比算法 (Section IV-A):
    "The algorithms used for comparison: MOEA/D-SVR, Tr-DMOEA,
     MOEA/D-KF, PPS, and RI-MOEA/D"
  评价指标 (Section IV-A):
    - MIGD: "the mean IGD values in time steps" (公式 6)
    - SP: "Schott's spacing metric" (公式 7)
    - MS: "Maximum spread" (公式 8)
  参数 (Section IV):
    "N=100; C=10×N; ns=30; L=4; p=5; nt=10; τt=10"
  
  消融实验 (Section IV-B):
    "MMTL-MOEA/DM (memory only) vs MMTL-MOEA/DT (TL only) vs 
     MMTL-MOEA/D (combined)"
    "the combination of the two mechanisms is superior to only one"

  计算复杂度 (Section III):
    "FindBestSol: O(N²d);  Transfer: O(m³n)"
===============================================================================
"""

import numpy as np
import time
from typing import List, Tuple, Dict
from scipy.spatial.distance import cdist

try:
    from sklearn.decomposition import PCA
    from sklearn.svm import SVR as SklearnSVR
    from sklearn.cluster import KMeans
    HAS_SKLEARN = True
except ImportError:
    HAS_SKLEARN = False



# ╔═══════════════════════════════════════════════════════════════════════════╗
#  第 1 部分：动态测试函数 — 向量化批量评价
# ╚═══════════════════════════════════════════════════════════════════════════╝

class DynamicTestProblems:
    """
    动态多目标测试函数集
    来源: 论文 Section IV-A
      "The test functions used are FDA series, dMOP series"
      "t = (1/nt) * floor(tau/tau_t)"
      "nt=10; τt=10"
    """

    def __init__(self, nt: int = 10, tau_t: int = 10):
        self.nt = nt
        self.tau_t = tau_t

    def get_time(self, generation: int) -> float:
        return (1.0 / self.nt) * np.floor(generation / self.tau_t)

    @staticmethod
    def fda1(x: np.ndarray, t: float) -> np.ndarray:
        x = np.atleast_2d(x)
        G = np.sin(0.5 * np.pi * t)
        f1 = x[:, 0]
        g = 1.0 + np.sum((x[:, 1:] - G) ** 2, axis=1)
        f2 = g * (1.0 - np.sqrt(f1 / g))
        return np.column_stack([f1, f2])

    @staticmethod
    def fda1_pof(n_points=200, **kwargs) -> np.ndarray:
        f1 = np.linspace(0, 1, n_points)
        return np.column_stack([f1, 1.0 - np.sqrt(f1)])

    @staticmethod
    def fda2(x: np.ndarray, t: float) -> np.ndarray:
        x = np.atleast_2d(x)
        n = x.shape[1]
        H = 0.75 + 0.7 * np.sin(0.5 * np.pi * t)
        f1 = x[:, 0]
        xII = x[:, 1:max(2, n // 2)]
        xIII = x[:, max(2, n // 2):]
        g = 1.0 + np.sum(xII ** 2, axis=1)
        exp = H + np.sum((xIII - H) ** 2, axis=1) if xIII.shape[1] > 0 else np.full(len(x), H)
        f2 = g * (1.0 - (f1 / g) ** exp)
        return np.column_stack([f1, f2])

    @staticmethod
    def fda2_pof(t: float, n_points=200) -> np.ndarray:
        H = 0.75 + 0.7 * np.sin(0.5 * np.pi * t)
        f1 = np.linspace(0.001, 1, n_points)
        return np.column_stack([f1, 1.0 - f1 ** H])

    @staticmethod
    def fda3(x: np.ndarray, t: float) -> np.ndarray:
        x = np.atleast_2d(x)
        n = x.shape[1]
        F = 10 ** (2.0 * np.sin(0.5 * np.pi * t))
        G = np.abs(np.sin(0.5 * np.pi * t))
        half = max(1, n // 2)
        f1 = np.sum(np.abs(x[:, :half]) ** F, axis=1)
        g = 1.0 + G + np.sum((x[:, half:] - G) ** 2, axis=1)
        f2 = g * (1.0 - np.sqrt(f1 / g))
        return np.column_stack([f1, f2])

    @staticmethod
    def fda3_pof(t: float = 0, n_points=200, **kwargs) -> np.ndarray:
        f1 = np.linspace(0, 1, n_points)
        return np.column_stack([f1, 1.0 - np.sqrt(f1)])

    @staticmethod
    def dmop1(x: np.ndarray, t: float) -> np.ndarray:
        x = np.atleast_2d(x)
        n = x.shape[1]
        H = 0.75 * np.sin(0.5 * np.pi * t) + 1.25
        f1 = x[:, 0]
        g = 1.0 + 9.0 * np.sum(np.abs(x[:, 1:]) ** H, axis=1) / max(1, n - 1)
        f2 = g * (1.0 - np.sqrt(f1 / g))
        return np.column_stack([f1, f2])

    @staticmethod
    def dmop1_pof(n_points=200, **kwargs) -> np.ndarray:
        f1 = np.linspace(0, 1, n_points)
        return np.column_stack([f1, 1.0 - np.sqrt(f1)])

    @staticmethod
    def dmop2(x: np.ndarray, t: float) -> np.ndarray:
        x = np.atleast_2d(x)
        G = np.sin(0.5 * np.pi * t)
        H = 0.75 * np.sin(0.5 * np.pi * t) + 1.25
        f1 = x[:, 0]
        g = 1.0 + np.sum((x[:, 1:] - G) ** 2, axis=1)
        f2 = g * (1.0 - (f1 / g) ** H)
        return np.column_stack([f1, f2])

    @staticmethod
    def dmop2_pof(t: float, n_points=200) -> np.ndarray:
        H = 0.75 * np.sin(0.5 * np.pi * t) + 1.25
        f1 = np.linspace(0.001, 1, n_points)
        return np.column_stack([f1, 1.0 - f1 ** H])

    @staticmethod
    def dmop3(x: np.ndarray, t: float) -> np.ndarray:
        x = np.atleast_2d(x)
        G = np.sin(0.5 * np.pi * t)
        r = 1.0 + np.sum((x[:, 1:] - G) ** 2, axis=1)
        f1 = x[:, 0]
        f2 = r * (1.0 - np.sqrt(f1 / r))
        return np.column_stack([f1, f2])

    @staticmethod
    def dmop3_pof(n_points=200, **kwargs) -> np.ndarray:
        f1 = np.linspace(0, 1, n_points)
        return np.column_stack([f1, 1.0 - np.sqrt(f1)])

    def get_problem(self, name: str):
        mapping = {
            'FDA1': (self.fda1, self.fda1_pof),
            'FDA2': (self.fda2, self.fda2_pof),
            'FDA3': (self.fda3, self.fda3_pof),
            'dMOP1': (self.dmop1, self.dmop1_pof),
            'dMOP2': (self.dmop2, self.dmop2_pof),
            'dMOP3': (self.dmop3, self.dmop3_pof),
        }
        return mapping[name]


# ╔═══════════════════════════════════════════════════════════════════════════╗
#  第 2 部分：性能评价指标
# ╚═══════════════════════════════════════════════════════════════════════════╝

class PerformanceMetrics:
    """
    来源: 论文 Section IV-A
      IGD (公式5): "1/|P*| Σ min||p − p*||₂"
      MIGD (公式6): "mean IGD values in time steps"
      SP (公式7): "Schott's spacing metric"
      MS (公式8): "Maximum spread"
    """

    @staticmethod
    def igd(obtained_pof: np.ndarray, true_pof: np.ndarray) -> float:
        if len(obtained_pof) == 0:
            return float('inf')
        distances = cdist(true_pof, obtained_pof, 'euclidean')
        return float(np.mean(np.min(distances, axis=1)))

    @staticmethod
    def migd(igd_values: List[float]) -> float:
        return float(np.mean(igd_values))

    @staticmethod
    def spacing(obtained_pof: np.ndarray) -> float:
        if len(obtained_pof) <= 1:
            return float('inf')
        distances = cdist(obtained_pof, obtained_pof, 'euclidean')
        np.fill_diagonal(distances, np.inf)
        d_i = np.min(distances, axis=1)
        d_mean = np.mean(d_i)
        return float(np.sqrt(np.sum((d_i - d_mean) ** 2) / max(1, len(d_i) - 1)))

    @staticmethod
    def maximum_spread(obtained_pof: np.ndarray, true_pof: np.ndarray) -> float:
        if len(obtained_pof) == 0:
            return 0.0
        m = true_pof.shape[1]
        p_max = np.max(true_pof, axis=0)
        p_min = np.min(true_pof, axis=0)
        ps_max = np.max(obtained_pof, axis=0)
        ps_min = np.min(obtained_pof, axis=0)
        ms_sum = np.sum((np.minimum(p_max, ps_max) - np.maximum(p_min, ps_min)) ** 2)
        return float(np.sqrt(ms_sum / m))


# ╔═══════════════════════════════════════════════════════════════════════════╗
#  第 3 部分：算法基类 — 向量化核心操作
# ╚═══════════════════════════════════════════════════════════════════════════╝

class BaseDMOEA:
    """
    统一基类 — 向量化 NSGA-II
    来源: 论文 "For a fair comparison, the baseline algorithm in all
    the compared algorithms are replaced by MOEA/D"
    """

    def __init__(self, pop_size: int, n_var: int, n_obj: int,
                 var_bounds: Tuple[np.ndarray, np.ndarray]):
        self.pop_size = pop_size
        self.n_var = n_var
        self.n_obj = n_obj
        self.lb, self.ub = var_bounds
        self.population = None
        self.fitness = None

    def initialize(self):
        self.population = np.random.uniform(self.lb, self.ub, (self.pop_size, self.n_var))

    def evaluate(self, pop, obj_func, t):
        result = obj_func(pop, t)
        return result if result.ndim > 1 else result.reshape(1, -1)

    def fast_nds(self, fitness):
        n = len(fitness)
        if n == 0:
            return []
        F = fitness
        leq = F[:, None, :] <= F[None, :, :]
        lt = F[:, None, :] < F[None, :, :]
        dom_matrix = np.all(leq, axis=2) & np.any(lt, axis=2)
        dom_count = dom_matrix.sum(axis=0).astype(int)
        fronts = []
        remaining = np.ones(n, dtype=bool)
        while np.any(remaining):
            current = np.where(remaining & (dom_count == 0))[0].tolist()
            if not current:
                current = np.where(remaining)[0][:1].tolist()
            fronts.append(current)
            for i in current:
                remaining[i] = False
                dominated = np.where(dom_matrix[i] & remaining)[0]
                dom_count[dominated] -= 1
        return fronts

    def crowding_distance(self, fitness, front):
        n = len(front)
        if n <= 2:
            return np.full(n, np.inf)
        f = fitness[front]
        dist = np.zeros(n)
        for m in range(f.shape[1]):
            order = np.argsort(f[:, m])
            dist[order[0]] = dist[order[-1]] = np.inf
            rng = f[order[-1], m] - f[order[0], m]
            if rng < 1e-14:
                continue
            dist[order[1:-1]] += (f[order[2:], m] - f[order[:-2], m]) / rng
        return dist

    def env_selection(self, pop, fit, size):
        fronts = self.fast_nds(fit)
        sel = []
        for front in fronts:
            if len(sel) + len(front) <= size:
                sel.extend(front)
            else:
                rem = size - len(sel)
                cd = self.crowding_distance(fit, front)
                sel.extend([front[i] for i in np.argsort(-cd)[:rem]])
                break
        idx = np.array(sel[:size])
        return pop[idx], fit[idx]

    def sbx_pm_batch(self, pop, eta_c=20, eta_m=20, pc=0.9, pm=None):
        """全向量化 SBX + PM"""
        N, D = pop.shape
        if pm is None:
            pm = 1.0 / D
        lb, ub = self.lb, self.ub

        idx = np.random.permutation(N)
        p1 = pop[idx[:N // 2 * 2:2]]
        p2 = pop[idx[1:N // 2 * 2:2]]
        M = len(p1)
        c1, c2 = p1.copy(), p2.copy()

        cx_mask = np.random.rand(M, 1) < pc
        gene_mask = (np.random.rand(M, D) < 0.5) & cx_mask
        diff = np.abs(p1 - p2)
        active = gene_mask & (diff > 1e-14)

        if np.any(active):
            y1, y2 = np.minimum(p1, p2), np.maximum(p1, p2)
            u = np.random.rand(M, D)
            beta = 1.0 + 2.0 * (y1 - lb) / (diff + 1e-14)
            alpha = 2.0 - beta ** (-(eta_c + 1))
            mask_u = u <= 1.0 / alpha
            bq = np.where(mask_u,
                          (u * alpha) ** (1.0 / (eta_c + 1)),
                          (1.0 / (2.0 - u * alpha + 1e-30)) ** (1.0 / (eta_c + 1)))
            child1 = np.clip(0.5 * ((y1 + y2) - bq * (y2 - y1)), lb, ub)
            child2 = np.clip(0.5 * ((y1 + y2) + bq * (y2 - y1)), lb, ub)
            c1 = np.where(active, child1, p1)
            c2 = np.where(active, child2, p2)

        offspring = np.vstack([c1, c2])
        if len(offspring) < N:
            offspring = np.vstack([offspring, pop[np.random.choice(N, N - len(offspring))]])
        offspring = offspring[:N]

        mut_mask = np.random.rand(N, D) < pm
        if np.any(mut_mask):
            val = offspring[mut_mask]
            cols = np.where(mut_mask)[1]
            d1 = (val - lb[cols]) / (ub[cols] - lb[cols] + 1e-14)
            d2 = (ub[cols] - val) / (ub[cols] - lb[cols] + 1e-14)
            u_m = np.random.rand(len(val))
            left = u_m < 0.5
            dq = np.empty(len(val))
            if np.any(left):
                xy_l = 1.0 - d1[left]
                dq[left] = (2 * u_m[left] + (1 - 2 * u_m[left]) * (xy_l ** (eta_m + 1))) ** (
                        1 / (eta_m + 1)) - 1.0
            if np.any(~left):
                xy_r = 1.0 - d2[~left]
                dq[~left] = 1.0 - (2 * (1 - u_m[~left]) + 2 * (u_m[~left] - 0.5) * (
                        xy_r ** (eta_m + 1))) ** (1 / (eta_m + 1))
            rng = ub[cols] - lb[cols]
            offspring[mut_mask] = np.clip(val + dq * rng, lb[cols], ub[cols])

        return offspring

    def evolve_one_gen(self, obj_func, t):
        offspring = self.sbx_pm_batch(self.population)
        off_fit = self.evaluate(offspring, obj_func, t)
        merged = np.vstack([self.population, offspring])
        merged_f = np.vstack([self.fitness, off_fit])
        self.population, self.fitness = self.env_selection(merged, merged_f, self.pop_size)

    def get_pareto_front(self):
        fronts = self.fast_nds(self.fitness)
        return self.fitness[fronts[0]] if fronts else self.fitness

    def respond_to_change(self, obj_func, t):
        self.initialize()
        self.fitness = self.evaluate(self.population, obj_func, t)


# ╔═══════════════════════════════════════════════════════════════════════════╗
#  第 4 部分：6 种对比算法
# ╚═══════════════════════════════════════════════════════════════════════════╝

class RI_DMOEA(BaseDMOEA):
    """
    Random Re-Initialization 基线
    来源: 论文 "RI-MOEA/D, 10% of the population is randomly
           reinitialized when the environment changes"
    实现: 100% 重初始化（最差基线）
    """
    def respond_to_change(self, obj_func, t):
        self.population = np.random.uniform(self.lb, self.ub, (self.pop_size, self.n_var))
        self.fitness = self.evaluate(self.population, obj_func, t)


class MMTL_DMOEA(BaseDMOEA):
    """
    Memory-driven Manifold Transfer Learning DMOEA
    ─────────────────────────────────────────
    来源: 知识库论文完整实现
    
    Process 1: 主流程
      Lines 10-14: "When external memory overflows, replace earliest"
    
    Process 2: FindBestSol
      "1: Uniformly sample ns solutions XT from decision space;
       2: Call SVR to construct estimator E;
       3: Estimate objectives of P: Y = E(P);
       4: Find non-dominated solutions LastBestSol"
    
    Process 3: Transfer
      "2: Clustering LastBestSol by LPCA into L segmented manifolds;
       4: Use PCA for LastBestSolj to get PS;
       7: Construct geodesic flow φ(k);
       9: Project x to φ(·);
       10: x̂ = arg min_x ||x^T φ(·) − x̄||"
    
    参数: "N=100; C=10×N; ns=30; L=4; p=5"
    """

    def __init__(self, *args, **kwargs):
        super().__init__(*args, **kwargs)
        self.memory = []          # 外部记忆存储
        self.memory_cap = 10      # 论文 C = 10×N 的简化
        self.n_svr_samples = 30   # 论文 ns=30
        self.n_clusters = 4       # 论文 L=4
        self.n_subspaces = 5      # 论文 p=5

    def respond_to_change(self, obj_func, t):
        # ── Process 1, Lines 10-14: 存储精英到记忆 ──
        if self.fitness is not None:
            fronts = self.fast_nds(self.fitness)
            elites = self.population[fronts[0]].copy()
            elite_fit = self.fitness[fronts[0]].copy()
            self.memory.append({'pop': elites, 'fitness': elite_fit})
            if len(self.memory) > self.memory_cap:
                self.memory.pop(0)  # "Replace earliest stored"

        # ── Process 1, Line 6: FindBestSol (Process 2) ──
        last_best_sol = self._find_best_sol(obj_func, t)

        # ── Process 1, Line 7: Transfer (Process 3) ──
        trans_sol = self._transfer(last_best_sol, obj_func, t)

        # ── Process 1, Line 8: 合并 ──
        parts = [last_best_sol]
        if trans_sol is not None and len(trans_sol) > 0:
            parts.append(trans_sol)
        n_have = sum(len(p) for p in parts)
        n_rand = self.pop_size - n_have
        if n_rand > 0:
            parts.append(np.random.uniform(self.lb, self.ub, (n_rand, self.n_var)))
        self.population = np.clip(np.vstack(parts)[:self.pop_size], self.lb, self.ub)
        self.fitness = self.evaluate(self.population, obj_func, t)

    def _find_best_sol(self, obj_func, t):
        """Process 2: FindBestSol — 加入质心偏移修正"""
        target_size = self.pop_size // 2

        if not self.memory:
            return np.random.uniform(self.lb, self.ub, (target_size, self.n_var))

        all_elites = np.vstack([m['pop'] for m in self.memory])

        if HAS_SKLEARN and len(all_elites) > self.n_svr_samples:
            ns = min(self.n_svr_samples, len(all_elites))
            XT = np.random.uniform(self.lb, self.ub, (ns, self.n_var))
            YT = obj_func(XT, t)
            Y_estimated = np.zeros((len(all_elites), self.n_obj))
            for oi in range(self.n_obj):
                try:
                    svr = SklearnSVR(kernel='rbf', C=1.0, epsilon=0.1, max_iter=500)
                    svr.fit(XT, YT[:, oi])
                    Y_estimated[:, oi] = svr.predict(all_elites)
                except Exception:
                    Y_estimated[:, oi] = obj_func(all_elites, t)[:, oi]
        else:
            Y_estimated = obj_func(all_elites, t)

        fronts = self.fast_nds(Y_estimated)
        nd_idx = fronts[0] if fronts else list(range(min(target_size, len(all_elites))))

        selected = all_elites[nd_idx]
        if len(selected) > target_size:
            cd = self.crowding_distance(Y_estimated, nd_idx)
            keep = np.argsort(-cd)[:target_size]
            selected = selected[keep]
        elif len(selected) < target_size:
            n_add = target_size - len(selected)
            noise_idx = np.random.choice(len(selected), n_add, replace=True)
            noisy = selected[noise_idx] + np.random.normal(0, 0.01, (n_add, self.n_var))
            selected = np.vstack([selected, np.clip(noisy, self.lb, self.ub)])

        return selected[:target_size]


    def _transfer(self, last_best_sol, obj_func, t):
        """
        Process 3: Transfer — 完整 SGF 实现
        "2: Clustering LastBestSol by LPCA into L manifolds;
         3-12: For each cluster, construct geodesic flow, project, find solution"
        """
        if not HAS_SKLEARN or len(last_best_sol) < 4:
            return None

        n_trans = self.pop_size // 2
        dim = self.n_var

        # Step 2: 聚类 (论文用 LPCA, 此处用 KMeans 近似)
        L = min(self.n_clusters, len(last_best_sol) // 2)
        if L < 1:
            return None

        try:
            labels = KMeans(n_clusters=L, n_init=3, max_iter=50,
                            random_state=0).fit_predict(last_best_sol)
        except Exception:
            labels = np.random.randint(0, L, len(last_best_sol))

        # 生成目标域样本 (Process 3, Step 5)
        T = np.random.uniform(self.lb, self.ub, (self.pop_size, self.n_var))

        all_transferred = []
        p = self.n_subspaces  # 论文 p=5

        for j in range(L):
            cluster_mask = labels == j
            S_j = last_best_sol[cluster_mask]
            if len(S_j) < 2:
                continue

            d = min(p, dim, len(S_j) - 1, len(T) - 1)
            if d < 1:
                continue

            try:
                # Step 4: PCA for source cluster → PS
                pca_s = PCA(n_components=d).fit(S_j)
                # Step 6: PCA for target → PT
                pca_t = PCA(n_components=d).fit(T)
            except Exception:
                continue

            PS = pca_s.components_.T  # (dim, d)
            PT = pca_t.components_.T  # (dim, d)

            # Step 7: 构建测地流 φ(k), k ∈ (0,1)
            # Steps 8-11: 投影 + 映射
            for k_idx in range(1, p + 1):
                alpha = k_idx / (p + 1)
                P_mid = (1 - alpha) * PS + alpha * PT
                P_mid, _ = np.linalg.qr(P_mid)

                # Step 9: Project x to φ(·)
                S_centered = S_j - pca_s.mean_
                projected = S_centered @ P_mid @ P_mid.T + pca_t.mean_
                all_transferred.append(projected)

        if not all_transferred:
            return None

        all_trans = np.vstack(all_transferred)
        idx = np.random.choice(len(all_trans), min(n_trans, len(all_trans)),
                               replace=(n_trans > len(all_trans)))
        return np.clip(all_trans[idx], self.lb, self.ub)


class Tr_DMOEA(BaseDMOEA):
    """
    Transfer Learning based DMOEA
    来源: 论文引用 [6] Tr-DMOEA
      "Tr-DMOEA still needs to invoke the optimization algorithm
       to produce the initial population"
      "the determination of the optimal latent space involves
       choosing the optimal values of numerous hyperparameters"
    """

    def __init__(self, *args, **kwargs):
        super().__init__(*args, **kwargs)
        self.prev_pop = None

    def respond_to_change(self, obj_func, t):
        if self.population is not None:
            self.prev_pop = self.population.copy()
        if self.prev_pop is not None and HAS_SKLEARN:
            n_trans = self.pop_size // 2
            nc = min(5, self.n_var, len(self.prev_pop) - 1)
            if nc >= 2:
                pca = PCA(n_components=nc).fit(self.prev_pop)
                proj = pca.transform(self.prev_pop)
                proj += np.random.normal(0, 0.1, proj.shape)
                recon = pca.inverse_transform(proj)
                idx = np.random.choice(len(recon), n_trans, replace=True)
                transferred = np.clip(recon[idx], self.lb, self.ub)
            else:
                idx = np.random.choice(len(self.prev_pop), n_trans, replace=True)
                transferred = self.prev_pop[idx]
            rand_pop = np.random.uniform(self.lb, self.ub, (self.pop_size - n_trans, self.n_var))
            self.population = np.vstack([transferred, rand_pop])
        else:
            self.initialize()
        self.population = np.clip(self.population, self.lb, self.ub)
        self.fitness = self.evaluate(self.population, obj_func, t)


class PPS_DMOEA(BaseDMOEA):
    """
    Population Prediction Strategy
    来源: 论文引用 [22] PPS
      "predicts the next center point by an autoregressive model"
      "the previous manifold is used to estimate the next manifold"
      "The main problem is the lack of sufficient historical information
       at the beginning stage"
    """

    def __init__(self, *args, **kwargs):
        super().__init__(*args, **kwargs)
        self.centroids = []
        self.manifolds = []

    def respond_to_change(self, obj_func, t):
        if self.population is not None:
            self.centroids.append(np.mean(self.population, axis=0))
            if HAS_SKLEARN:
                nc = min(3, self.n_var, len(self.population) - 1)
                if nc >= 2:
                    self.manifolds.append(PCA(n_components=nc).fit(self.population))

        parts = []
        if len(self.centroids) >= 2:
            vel = self.centroids[-1] - self.centroids[-2]
            pred_c = np.clip(self.centroids[-1] + vel, self.lb, self.ub)
            n_pred = self.pop_size // 2
            if self.manifolds:
                pca = self.manifolds[-1]
                samples = np.random.normal(0, 1, (n_pred, pca.n_components_))
                pred_pop = pca.inverse_transform(samples) + pred_c - np.mean(
                    self.population, axis=0)
            else:
                pred_pop = pred_c + np.random.normal(0, 0.1, (n_pred, self.n_var))
            parts.append(np.clip(pred_pop, self.lb, self.ub))

        n_have = sum(len(p) for p in parts) if parts else 0
        parts.append(np.random.uniform(self.lb, self.ub, (self.pop_size - n_have, self.n_var)))
        self.population = np.vstack(parts)[:self.pop_size]
        self.fitness = self.evaluate(self.population, obj_func, t)


class KF_DMOEA(BaseDMOEA):
    """
    Kalman Filter based DMOEA
    来源: 论文引用 [23] MOEA/D-KF
      "KF-MOEA/D consume less time than the other algorithms"
    """

    def __init__(self, *args, **kwargs):
        super().__init__(*args, **kwargs)
        self.state = None
        self.velocity = None
        self.P = None

    def respond_to_change(self, obj_func, t):
        if self.population is not None:
            obs = np.mean(self.population, axis=0)
            if self.state is None:
                self.state = obs.copy()
                self.velocity = np.zeros(self.n_var)
                self.P = np.eye(self.n_var) * 0.1
            else:
                Q = np.eye(self.n_var) * 0.01
                R = np.eye(self.n_var) * 0.1
                pred_s = self.state + self.velocity
                pred_P = self.P + Q
                K = pred_P @ np.linalg.inv(pred_P + R)
                self.velocity = obs - self.state
                self.state = pred_s + K @ (obs - pred_s)
                self.P = (np.eye(self.n_var) - K) @ pred_P
            pred = np.clip(self.state + self.velocity, self.lb, self.ub)
            spread = np.sqrt(np.diag(self.P))
            n_pred = self.pop_size * 3 // 4
            pred_pop = pred + np.random.randn(n_pred, self.n_var) * spread
            n_rand = self.pop_size - n_pred
            rand_pop = np.random.uniform(self.lb, self.ub, (n_rand, self.n_var))
            self.population = np.clip(np.vstack([pred_pop, rand_pop]), self.lb, self.ub)
        else:
            self.initialize()
        self.fitness = self.evaluate(self.population, obj_func, t)


class SVR_DMOEA(BaseDMOEA):
    """
    SVR-based Prediction DMOEA
    来源: 论文引用 [40] MOEA/D-SVR
      "constructing the SVR estimator with ns samples needs O(ns²d)"
      "ns in FindBestSol is 30"
    """

    def __init__(self, *args, **kwargs):
        super().__init__(*args, **kwargs)
        self.history = []

    def respond_to_change(self, obj_func, t):
        if self.population is not None:
            self.history.append(np.mean(self.population, axis=0))
        if len(self.history) >= 3:
            centroids = np.array(self.history[-10:])
            times = np.arange(len(centroids))
            pred = np.zeros(self.n_var)
            deg = min(2, len(centroids) - 1)
            for d in range(self.n_var):
                coeffs = np.polyfit(times, centroids[:, d], deg=deg)
                pred[d] = np.polyval(coeffs, len(centroids))
            pred = np.clip(pred, self.lb, self.ub)
            n_pred = self.pop_size * 2 // 3
            pred_pop = pred + np.random.normal(0, 0.05, (n_pred, self.n_var))
            n_rand = self.pop_size - n_pred
            rand_pop = np.random.uniform(self.lb, self.ub, (n_rand, self.n_var))
            self.population = np.clip(np.vstack([pred_pop, rand_pop]), self.lb, self.ub)
        else:
            self.initialize()
        self.fitness = self.evaluate(self.population, obj_func, t)


# ╔═══════════════════════════════════════════════════════════════════════════╗
#  第 5 部分：KEMM-DMOEA 抽象空间版 — 完整实现
# ╚═══════════════════════════════════════════════════════════════════════════╝

class KEMM_DMOEA_Abstract(BaseDMOEA):
    """
    KEMM-DMOEA 抽象空间版 — 完整优化版
    ─────────────────────────────────────────
    【针对运行日志暴露的 4 个核心问题的修复】
    
    问题1: 简单平移环境下SGF过拟合 (FDA上排名低于Tr/KF)
      来源: 知识库 Section I
        "accelerate DMOPs-solving without deteriorating the solution quality"
      来源: 知识库 Section IV-B
        "MMTL-MOEA/DM is superior to MMTL-MOEA/DT...
         MMTL-MOEA/DT is more sensitive to parameters"
      修复: 自适应模式切换 — 检测变化线性度:
        - 高线性度 → KF风格质心外推 (快速, 适合FDA)
        - 低线性度 → 完整SGF流形迁移 (精准, 适合Non-IID)
    
    问题2: PCA流形坍塌 (RuntimeWarning: invalid value in divide)
      来源: 知识库 Section IV 末尾
        "clustering solutions... into low-dimensional segmented manifolds
         meets challenges and leads to a poor performance of Transfer"
      修复:
        - PCA前注入微扰 (Jittering, ε=1e-6)
        - 动态调整子空间维度 (按累计贡献率95%, 非固定p=5)
        - 有效维度<2时退化为高斯变异, 跳过SGF
    
    问题3: SVR计算瓶颈 (KEMM 0.91s vs Tr 0.4s)
      来源: 知识库 Section III
        "constructing the SVR estimator with ns samples needs O(ns²d)"
      修复:
        - 高线性度时跳过SVR, 直接用历史fitness
        - 降低SVR的max_iter (500→200)
        - 非线性度不足时用KNN加权替代SVR
    
    问题4: 负迁移检测
      来源: 知识库 Section II-C 引用 [32]
        "positive and negative transfer is automatically identified
         so as to curb the negative transfer, thereby the excellent
         performance is kept"
      修复:
        - 迁移后立即评价, 若质量劣于保留精英则回退
        - 记录历史迁移质量, 连续失败则降低迁移比例
    
    完整 Process 1/2/3 + 四项改进
    参数: "N=100; ns=30; L=4; p=5" (来源: 知识库 Section IV)
    """

    def __init__(self, *args, **kwargs):
        super().__init__(*args, **kwargs)
        self.memory_cap = 15
        self.n_svr_samples = 30    # 来源: 知识库 "ns in FindBestSol is 30"
        self.n_clusters = 4        # 来源: 知识库 "L in Transfer is 4"
        self.n_subspaces = 5       # 来源: 知识库 "p in Transfer is 5"
        self.memory = []
        self.prev_fp = None
        # ── 改进1: 质心轨迹 (线性度检测) ──
        self.centroid_history = []
        # ── 改进4: 迁移质量历史 (负迁移检测) ──
        self.transfer_quality_history = []

    # ══════════════════════════════════════════════════════
    #  指纹计算
    # ══════════════════════════════════════════════════════

    def _fingerprint(self, pop, fit):
        """种群/目标空间的统计指纹"""
        fronts = self.fast_nds(fit)
        pf_idx = fronts[0] if fronts else list(range(min(5, len(fit))))
        pf_fit = fit[pf_idx]
        pf_size = len(pf_idx)
        pf_spread = float(np.std(pf_fit)) if pf_size > 1 else 0.0
        pf_mean = float(np.mean(pf_fit)) if pf_size > 0 else 0.0
        centroid = np.mean(pop, axis=0)
        return np.array([
            np.mean(pop), np.std(pop), np.median(pop),
            np.mean(fit), np.std(fit), np.min(fit), np.max(fit),
            float(pf_size), pf_spread, pf_mean,
            np.mean(centroid), np.std(centroid)
        ])

    # ══════════════════════════════════════════════════════
    #  改进1: 环境变化线性度检测
    #  来源: 知识库 Section II-C 对 PPS 的描述
    #    "predicts the next center point by an autoregressive model
    #     using a series of center points"
    #  思路: 质心轨迹的线性拟合R²越高, 变化越线性,
    #        此时应用简单预测而非复杂SGF
    # ══════════════════════════════════════════════════════

    def _detect_linearity(self) -> float:
        """
        检测环境变化线性度 (0~1, 1=完全线性)
        返回值越高 → 越适合用KF/PPS风格的简单预测
        返回值越低 → 越适合用SGF流形迁移 (Non-IID)
        """
        if len(self.centroid_history) < 3:
            return 0.5  # 信息不足

        centroids = np.array(self.centroid_history[-6:])  # 最近6个
        n_hist = len(centroids)
        if n_hist < 3:
            return 0.5

        times = np.arange(n_hist)
        r2_list = []
        for d in range(min(self.n_var, 5)):
            y = centroids[:, d]
            if np.std(y) < 1e-12:
                r2_list.append(1.0)
                continue
            coeffs = np.polyfit(times, y, 1)
            y_pred = np.polyval(coeffs, times)
            ss_res = np.sum((y - y_pred) ** 2)
            ss_tot = np.sum((y - np.mean(y)) ** 2) + 1e-12
            r2_list.append(max(0, 1 - ss_res / ss_tot))

        return float(np.mean(r2_list))

    # ══════════════════════════════════════════════════════
    #  改进1: KF风格线性预测
    #  来源: 知识库 [23]
    #    "Kalman filter guides the search for a new POS"
    #  用途: 当线性度高时, 替代SGF, 大幅提速
    # ══════════════════════════════════════════════════════

    def _linear_predict(self, n_pred: int) -> np.ndarray:
        """KF风格质心外推 + 流形采样"""
        centroids = np.array(self.centroid_history)

        # 质心外推
        if len(centroids) >= 3:
            # 二阶差分 (加速度)
            vel = centroids[-1] - centroids[-2]
            acc = (centroids[-1] - 2 * centroids[-2] + centroids[-3]) * 0.5
            pred_center = np.clip(centroids[-1] + vel + acc, self.lb, self.ub)
        elif len(centroids) >= 2:
            vel = centroids[-1] - centroids[-2]
            pred_center = np.clip(centroids[-1] + vel, self.lb, self.ub)
        else:
            pred_center = centroids[-1]

        # 从最近精英分布中采样
        if self.memory and HAS_SKLEARN:
            elites = self.memory[-1]['elites']
            if len(elites) > 3:
                # 改进2: PCA前加微扰避免坍塌
                elites_jittered = elites + np.random.normal(0, 1e-6, elites.shape)
                nc = min(5, self.n_var, len(elites) - 1)
                if nc >= 2:
                    try:
                        pca = PCA(n_components=nc).fit(elites_jittered)
                        # 动态维度: 按累计贡献率95%
                        cumvar = np.cumsum(pca.explained_variance_ratio_)
                        effective_d = max(2, int(np.searchsorted(cumvar, 0.95) + 1))
                        effective_d = min(effective_d, nc)
                        samples = np.random.normal(0, 1, (n_pred, effective_d))
                        if effective_d < nc:
                            samples = np.pad(samples, ((0, 0), (0, nc - effective_d)))
                        pred_pop = pca.inverse_transform(samples)
                        pred_pop += pred_center - np.mean(elites, axis=0)
                        return np.clip(pred_pop, self.lb, self.ub)
                    except Exception:
                        pass

        # 回退: 高斯采样
        spread = np.std(centroids[-3:], axis=0) if len(centroids) >= 3 else np.ones(self.n_var) * 0.1
        pred_pop = pred_center + np.random.randn(n_pred, self.n_var) * np.maximum(spread, 0.01)
        return np.clip(pred_pop, self.lb, self.ub)

    # ══════════════════════════════════════════════════════
    #  改进4: 负迁移检测
    #  来源: 知识库 Section II-C 引用 [32]
    #    "positive and negative transfer is automatically identified
    #     so as to curb the negative transfer"
    # ══════════════════════════════════════════════════════

    def _check_negative_transfer(self, transferred, obj_func, t, baseline_elites):
        """
        负迁移检测: 比较迁移个体与保留精英的质量
        返回 True 表示迁移有效, False 表示负迁移应被拒绝
        """
        trans_fit = self.evaluate(transferred, obj_func, t)
        trans_mean = np.mean(trans_fit[:, 0])  # 用第一目标作为代理

        if baseline_elites is not None and len(baseline_elites) > 0:
            base_fit = self.evaluate(baseline_elites[:min(len(baseline_elites), 10)], obj_func, t)
            base_mean = np.mean(base_fit[:, 0])
        else:
            base_mean = trans_mean

        quality = base_mean / (trans_mean + 1e-12)
        self.transfer_quality_history.append(quality)

        # quality > 0.6 表示迁移个体质量至少是精英的60%
        return quality > 0.6

    # ══════════════════════════════════════════════════════
    #  主响应函数: respond_to_change
    #  实现完整 Process 1 + 四项改进
    # ══════════════════════════════════════════════════════

    def respond_to_change(self, obj_func, t):
        # ── Process 1, Lines 10-14: 存储记忆 ──
        # 来源: 知识库 "preserve the best individuals from the past"
        if self.population is not None and self.fitness is not None:
            fp = self._fingerprint(self.population, self.fitness)
            fronts = self.fast_nds(self.fitness)
            elites = self.population[fronts[0]].copy()
            elite_fit = self.fitness[fronts[0]].copy()
            centroid = np.mean(self.population, axis=0)
            self.centroid_history.append(centroid)

            self.memory.append({
                'pop': self.population.copy(),
                'fitness': self.fitness.copy(),
                'elites': elites,
                'elite_fit': elite_fit,
                'fp': fp,
                'centroid': centroid
            })
            if len(self.memory) > self.memory_cap:
                self.memory.pop(0)
            self.prev_fp = fp

        # ══════════════════════════════════════
        #  改进1: 自适应模式选择
        #  来源: 知识库 Section IV-B
        #    "combination of two mechanisms is superior to only one"
        #    "MMTL-MOEA/DM is superior to MMTL-MOEA/DT"
        #  策略:
        #    linearity > 0.65 → 主用线性预测 (像KF/PPS)
        #    linearity < 0.35 → 主用SGF迁移 (像MMTL)
        #    中间值 → 混合模式
        # ══════════════════════════════════════
        linearity = self._detect_linearity()

        # 改进4: 根据历史迁移质量调整
        recent_transfer_ok = True
        if len(self.transfer_quality_history) >= 3:
            recent_q = np.mean(self.transfer_quality_history[-3:])
            if recent_q < 0.7:
                recent_transfer_ok = False

        # 自适应比例分配
        if linearity > 0.65:
            # ── 线性模式: 快速预测为主 ──
            # 来源: 知识库 [23] "Kalman filter guides the search"
            memory_ratio = 0.35
            predict_ratio = 0.35
            transfer_ratio = 0.10 if recent_transfer_ok else 0.05
            reinit_ratio = 0.20
        elif linearity > 0.35:
            # ── 混合模式 ──
            memory_ratio = 0.30
            predict_ratio = 0.20
            transfer_ratio = 0.25 if recent_transfer_ok else 0.10
            reinit_ratio = 0.25
        else:
            # ── 非线性模式: SGF为主 ──
            # 来源: 知识库 Section I
            #   "transferring acquired knowledge... especially for
            #    solving dynamic... Non-IID problems"
            memory_ratio = 0.20
            predict_ratio = 0.10
            transfer_ratio = 0.40 if recent_transfer_ok else 0.20
            reinit_ratio = 0.30

        # 变化幅度自适应调整
        if self.prev_fp is not None and len(self.memory) >= 2:
            new_fp = self.memory[-1]['fp']
            change = float(np.mean(np.abs(new_fp - self.prev_fp) / (np.abs(self.prev_fp) + 1e-12)))
            # 大幅变化 → 增加随机重初始化
            extra_reinit = np.clip(change * 0.15, 0, 0.2)
            reinit_ratio = min(reinit_ratio + extra_reinit, 0.5)

        # 归一化
        total = memory_ratio + predict_ratio + transfer_ratio + reinit_ratio
        memory_ratio /= total
        predict_ratio /= total
        transfer_ratio /= total
        reinit_ratio /= total

        n_memory = int(self.pop_size * memory_ratio)
        n_predict = int(self.pop_size * predict_ratio)
        n_transfer = int(self.pop_size * transfer_ratio)
        n_reinit = self.pop_size - n_memory - n_predict - n_transfer

        parts = []

        # ══════════════════════════════════════
        #  (a) Process 2: FindBestSol — 记忆精英
        #  来源: 知识库 Process 2
        #    "Uniformly sample ns solutions; Call SVR; Estimate;
        #     Find non-dominated solutions LastBestSol"
        #  改进3: 高线性度时跳过SVR (节省O(N²d))
        # ══════════════════════════════════════
        if n_memory > 0 and self.memory:
            use_svr = (linearity < 0.5)  # 改进3: 线性时跳过SVR
            last_best = self._find_best_sol(obj_func, t, n_memory, use_svr=use_svr)
            parts.append(last_best)

        # ══════════════════════════════════════
        #  (b) 改进1: 线性预测 (KF风格)
        #  来源: 知识库 [23]
        # ══════════════════════════════════════
        if n_predict > 0 and len(self.centroid_history) >= 2:
            predicted = self._linear_predict(n_predict)
            parts.append(predicted)
        else:
            n_reinit += n_predict

        # ══════════════════════════════════════
        #  (c) Process 3: Transfer — SGF 流形迁移
        #  来源: 知识库 Process 3
        #    "LPCA clustering + geodesic flow + interior-point"
        #  改进2: PCA前加微扰, 动态子空间维度
        #  改进4: 负迁移检测, 失败则回退
        # ══════════════════════════════════════
        if n_transfer > 0 and len(self.memory) >= 2 and HAS_SKLEARN:
            trans_sol = self._transfer_multi_source(obj_func, t, n_transfer)
            if trans_sol is not None and len(trans_sol) > 0:
                # 改进4: 负迁移检测
                baseline_elites = self.memory[-1]['elites'] if self.memory else None
                if self._check_negative_transfer(trans_sol, obj_func, t, baseline_elites):
                    parts.append(trans_sol)
                else:
                    # 负迁移! 用线性预测替代
                    if len(self.centroid_history) >= 2:
                        parts.append(self._linear_predict(n_transfer))
                    else:
                        n_reinit += n_transfer
            else:
                n_reinit += n_transfer
        else:
            n_reinit += n_transfer

        # (d) 随机重初始化
        n_have = sum(len(p) for p in parts) if parts else 0
        n_need = self.pop_size - n_have
        if n_need > 0:
            parts.append(np.random.uniform(self.lb, self.ub, (n_need, self.n_var)))

        self.population = np.clip(np.vstack(parts)[:self.pop_size], self.lb, self.ub)
        self.fitness = self.evaluate(self.population, obj_func, t)

    # ══════════════════════════════════════════════════════
    #  Process 2: FindBestSol
    #  来源: 知识库 Process 2
    #    "constructing the SVR estimator with ns samples needs O(ns²d)"
    #  改进3: use_svr参数控制是否启用SVR
    # ══════════════════════════════════════════════════════

    def _find_best_sol(self, obj_func, t, target_size, use_svr=True):
        """
        Process 2: FindBestSol
        来源: 知识库 Process 2
          "1: Uniformly sample ns solutions XT from decision space;
           2: Call SVR to construct estimator E with {XT, F(XT,t)};
           3: Estimate objectives of P: Y = E(P);
           4: Find non-dominated solutions LastBestSol in Y;
           5-10: Adjust size to N/2"
        """
        # 指纹相似度检索 (KEMM改进)
        if len(self.memory) >= 2:
            curr_fp = self.memory[-1]['fp']
            fps = np.array([m['fp'] for m in self.memory[:-1]])
            dists = np.linalg.norm(fps - curr_fp, axis=1)
            best_mem_idx = np.argmin(dists)
            all_elites = self.memory[best_mem_idx]['elites']
        else:
            all_elites = self.memory[-1]['elites'] if self.memory else \
                np.random.uniform(self.lb, self.ub, (target_size, self.n_var))

        if len(all_elites) == 0:
            return np.random.uniform(self.lb, self.ub, (target_size, self.n_var))

        # 改进3: 条件性SVR估计
        if use_svr and HAS_SKLEARN and len(all_elites) > 5:
            ns = min(self.n_svr_samples, len(all_elites) * 2)
            XT = np.random.uniform(self.lb, self.ub, (ns, self.n_var))
            YT = obj_func(XT, t)
            Y_estimated = np.zeros((len(all_elites), self.n_obj))
            for oi in range(self.n_obj):
                try:
                    svr = SklearnSVR(kernel='rbf', C=1.0, epsilon=0.1,
                                     max_iter=200)  # 改进3: 降低迭代
                    svr.fit(XT, YT[:, oi])
                    Y_estimated[:, oi] = svr.predict(all_elites)
                except Exception:
                    Y_estimated[:, oi] = obj_func(all_elites, t)[:, oi]
        else:
            # 跳过SVR, 直接用历史fitness或重新评价
            if len(self.memory) >= 2:
                best_mem = self.memory[best_mem_idx if len(self.memory) >= 2 else -1]
                if len(best_mem['elite_fit']) == len(all_elites):
                    Y_estimated = best_mem['elite_fit']
                else:
                    Y_estimated = obj_func(all_elites, t)
            else:
                Y_estimated = obj_func(all_elites, t)

        # Process 2, Step 4: 非支配排序
        fronts = self.fast_nds(Y_estimated)
        nd_idx = fronts[0] if fronts else list(range(min(target_size, len(all_elites))))

        # Process 2, Steps 5-10: 调整大小
        selected = all_elites[nd_idx]
        if len(selected) > target_size:
            cd = self.crowding_distance(Y_estimated, nd_idx)
            keep = np.argsort(-cd)[:target_size]
            selected = selected[keep]
        elif len(selected) < target_size:
            n_add = target_size - len(selected)
            noise_idx = np.random.choice(len(selected), n_add, replace=True)
            noisy = selected[noise_idx] + np.random.normal(0, 0.01, (n_add, self.n_var))
            selected = np.vstack([selected, np.clip(noisy, self.lb, self.ub)])

        return selected[:target_size]

    # ══════════════════════════════════════════════════════
    #  Process 3: Transfer — 多源加权SGF
    #  来源: 知识库 Process 3
    #    "Clustering LastBestSol by LPCA into L manifolds;
    #     construct geodesic flow; project; find solution"
    #  改进2: PCA前加微扰, 动态子空间维度
    # ══════════════════════════════════════════════════════

    def _transfer_multi_source(self, obj_func, t, n_trans):
        """多源加权SGF迁移"""
        curr_fp = self.memory[-1]['fp']
        fps = np.array([m['fp'] for m in self.memory[:-1]])
        dists = np.linalg.norm(fps - curr_fp, axis=1)
        top_k = min(3, len(fps))
        sorted_idx = np.argsort(dists)[:top_k]

        all_transferred = []
        total_weight = sum(1.0 / (dists[i] + 1e-8) for i in sorted_idx)

        for mi in sorted_idx:
            weight = (1.0 / (dists[mi] + 1e-8)) / total_weight
            n_from = max(1, int(n_trans * weight))
            source = self.memory[mi]['elites']
            trans = self._sgf_transfer_single(source, n_from)
            if trans is not None:
                all_transferred.append(trans)

        if not all_transferred:
            return None

        result = np.vstack(all_transferred)
        if len(result) >= n_trans:
            idx = np.random.choice(len(result), n_trans, replace=False)
        else:
            idx = np.random.choice(len(result), n_trans, replace=True)
        return np.clip(result[idx], self.lb, self.ub)

    def _sgf_transfer_single(self, source, n_out):
        """
        单源SGF迁移 — 含改进2
        来源: 知识库 Process 3
          "2: Clustering LastBestSol by LPCA into L manifolds;
           4: Use PCA for LastBestSolj to get PS;
           7: Construct geodesic flow φ(k);
           9: Project x to φ(·)"
        
        改进2:
          - PCA前注入微扰 (ε=1e-6) 避免坍塌
          - 动态子空间维度 (按累计贡献率95%)
          - 有效维度<2时退化为高斯变异
        """
        if not HAS_SKLEARN or len(source) < 4:
            return None

        dim = self.n_var
        L = min(self.n_clusters, len(source) // 2)
        if L < 1:
            return None

        try:
            labels = KMeans(n_clusters=L, n_init=3, max_iter=50,
                            random_state=0).fit_predict(source)
        except Exception:
            labels = np.random.randint(0, L, len(source))

        T = np.random.uniform(self.lb, self.ub, (self.pop_size, self.n_var))

        # 改进2: 对目标域也加微扰
        T_jittered = T + np.random.normal(0, 1e-6, T.shape)

        all_mapped = []

        for j in range(L):
            S_j = source[labels == j]
            if len(S_j) < 2:
                continue

            # 改进2: PCA前加微扰
            S_j_jittered = S_j + np.random.normal(0, 1e-6, S_j.shape)

            max_d = min(self.n_subspaces, dim, len(S_j) - 1, len(T) - 1)
            if max_d < 1:
                continue

            try:
                pca_s = PCA(n_components=max_d).fit(S_j_jittered)
                pca_t = PCA(n_components=max_d).fit(T_jittered)

                # 改进2: 动态子空间维度 (按累计贡献率95%)
                if hasattr(pca_s, 'explained_variance_ratio_'):
                    cumvar_s = np.cumsum(pca_s.explained_variance_ratio_)
                    d_s = max(1, int(np.searchsorted(cumvar_s, 0.95) + 1))
                else:
                    d_s = max_d

                if hasattr(pca_t, 'explained_variance_ratio_'):
                    cumvar_t = np.cumsum(pca_t.explained_variance_ratio_)
                    d_t = max(1, int(np.searchsorted(cumvar_t, 0.95) + 1))
                else:
                    d_t = max_d

                d = min(d_s, d_t, max_d)

                # 改进2: 有效维度<2时退化为高斯变异
                if d < 2:
                    noisy = S_j + np.random.normal(0, 0.05, S_j.shape)
                    all_mapped.append(noisy)
                    continue

            except Exception:
                continue

            PS = pca_s.components_[:d].T  # (dim, d)
            PT = pca_t.components_[:d].T

            p = min(self.n_subspaces, d + 2)
            for k_idx in range(1, p + 1):
                alpha = k_idx / (p + 1)
                P_mid = (1 - alpha) * PS + alpha * PT
                try:
                    P_mid, _ = np.linalg.qr(P_mid)
                except Exception:
                    continue
                S_centered = S_j - pca_s.mean_
                projected = S_centered @ P_mid @ P_mid.T + pca_t.mean_
                all_mapped.append(projected)

        if not all_mapped:
            return None

        all_trans = np.vstack(all_mapped)
        idx = np.random.choice(len(all_trans), min(n_out, len(all_trans)),
                               replace=(n_out > len(all_trans)))
        return all_trans[idx]