1{
2  "op": "replace",
3  "filename": "causal-learn/bench/custom_algorithm.py",
4  "start_line": 3,
5  "end_line": 14,
6  "content": "# =====================================================================\n# Non-Gaussianity-Adaptive Causal Discovery (NACD)\n# Novel: adaptive nonlinear function + BIC-Lasso parent selection\n# =====================================================================\ndef run_causal_discovery(X: np.ndarray) -> np.ndarray:\n    \"\"\"\n    Input:  X of shape (n_samples, n_variables)\n    Output: adjacency matrix B of shape (n_variables, n_variables)\n            B[i, j] != 0  means j -> i  (follows causal-learn convention)\n    \"\"\"\n    import warnings\n    warnings.filterwarnings('ignore')\n    from sklearn.utils import check_array\n    from sklearn.linear_model import LassoLarsIC\n\n    X = check_array(X, copy=True).astype(np.float64)\n    n, p = X.shape\n    Xs = (X - X.mean(0)) / (X.std(0) + 1e-10)\n\n    # === Novel 1: Adaptive nonlinear function based on kurtosis ===\n    # Light-tailed noise (Uniform, kurt < 0): x*exp(-x^2/2) captures bounded support\n    # Heavy-tailed noise (Laplace/Exponential, kurt > 0): tanh is robust\n    kurt = np.median(np.mean(Xs ** 4, axis=0) - 3.0)\n    if kurt < -0.5:\n        def g(x):\n            return x * np.exp(-0.5 * x * x)\n    else:\n        def g(x):\n            return np.tanh(x)\n\n    # === DirectLiNGAM Causal Ordering (vectorized) ===\n    remaining = list(range(p))\n    order = []\n    R = Xs.copy()\n\n    for _ in range(p):\n        nc = len(remaining)\n        if nc <= 1:\n            if nc == 1:\n                order.append(remaining[0])\n            break\n\n        vR = np.sum(R * R, axis=0) / n  # variances\n        bc, bs = 0, np.inf\n\n        for c in range(nc):\n            if vR[c] < 1e-10:\n                continue\n            xj = R[:, c]\n            # Regression betas for all columns on candidate root xj\n            betas = (xj @ R) / (n * vR[c])\n            # Residuals after removing linear effect of xj\n            res = R - np.outer(xj, betas)\n            # Standardize for dependence computation\n            xjs = (xj - xj.mean()) / (xj.std() + 1e-10)\n            rn = (res - res.mean(0)) / (res.std(0) + 1e-10)\n            # Non-linear covariance dependence measure\n            d = np.abs(g(xjs) @ rn / n) + np.abs(xjs @ g(rn) / n)\n            d[c] = 0.0  # exclude self\n            s = d.sum()\n            if s < bs:\n                bs, bc = s, c\n\n        order.append(remaining[bc])\n        # Residualize: remove effect of chosen root variable\n        xm = R[:, bc]\n        bm = (xm @ R) / (n * (vR[bc] + 1e-10))\n        R = (R - np.outer(xm, bm))[:, [i for i in range(nc) if i != bc]]\n        remaining = [remaining[i] for i in range(nc) if i != bc]\n\n    # === Novel 2: Adjacency with BIC-regularized Lasso parent selection ===\n    # LassoLarsIC selects sparsity via BIC, adapting to each node's\n    # effective dimensionality without manual threshold tuning\n    B = np.zeros((p, p))\n\n    for k in range(1, p):\n        j = order[k]\n        anc = [order[i] for i in range(k)]\n        na = len(anc)\n        if na == 0:\n            continue\n\n        y = Xs[:, j]\n        Xa = Xs[:, anc]\n\n        try:\n            model = LassoLarsIC(criterion='bic', max_iter=10000)\n            model.fit(Xa, y)\n            coefs = model.coef_\n        except Exception:\n            # Fallback: OLS + BIC-consistent threshold\n            coefs = np.linalg.lstsq(Xa, y, rcond=None)[0]\n            res = y - Xa @ coefs\n            rv = np.var(res) + 1e-10\n            se = np.sqrt(rv / (np.var(Xa, axis=0) * n + 1e-10)) + 1e-10\n            coefs = np.where(np.abs(coefs) / se > np.sqrt(np.log(n)), coefs, 0.0)\n\n        for idx in range(na):\n            if abs(coefs[idx]) > 1e-8:\n                B[j, anc[idx]] = coefs[idx]\n\n    return B\n# ====================================================================="
7}

1{
2  "op": "replace",
3  "filename": "causal-learn/bench/custom_algorithm.py",
4  "start_line": 3,
5  "end_line": 14,
6  "content": "# =====================================================================\n# Two-Stage LiNGAM: Non-Gaussianity Ordering + BIC-Lasso Parent Selection\n# =====================================================================\ndef run_causal_discovery(X: np.ndarray) -> np.ndarray:\n    \"\"\"\n    Stage 1: DirectLiNGAM recovers causal ordering via non-Gaussianity.\n    Stage 2: BIC-penalized LARS-Lasso selects parents; OLS debiases coefficients.\n    B[i, j] != 0 means j -> i.\n    \"\"\"\n    import warnings; warnings.filterwarnings('ignore')\n    from sklearn.utils import check_array\n    from sklearn.linear_model import LassoLarsIC\n    from causallearn.search.FCMBased.lingam import DirectLiNGAM\n\n    X = check_array(X)\n    n_samples, n_features = X.shape\n\n    # Stage 1: Causal ordering via DirectLiNGAM\n    dlingam = DirectLiNGAM()\n    dlingam.fit(X)\n    order = dlingam.causal_order_\n\n    # Stage 2: BIC-Lasso parent selection + OLS debiasing\n    B = np.zeros((n_features, n_features))\n    for i in range(1, len(order)):\n        target = order[i]\n        parents = list(order[:i])\n        if not parents:\n            continue\n        y = X[:, target]\n        Xp = X[:, parents]\n        try:\n            lasso = LassoLarsIC(criterion='bic')\n            lasso.fit(Xp, y)\n            sel = np.abs(lasso.coef_) > 1e-10\n            if np.any(sel):\n                sel_idx = [parents[j] for j in range(len(parents)) if sel[j]]\n                Xs = X[:, sel_idx]\n                beta = np.linalg.lstsq(Xs, y, rcond=None)[0]\n                for j, s in enumerate(sel_idx):\n                    B[target, s] = beta[j]\n        except Exception:\n            for j, s in enumerate(parents):\n                v = dlingam.adjacency_matrix_[target, s]\n                if abs(v) > 0:\n                    B[target, s] = v\n    return B\n# =====================================================================\n"
7}

1{
2  "op": "replace",
3  "filename": "causal-learn/bench/custom_algorithm.py",
4  "start_line": 3,
5  "end_line": 50,
6  "content": "# =====================================================================\n# Robust DirectLiNGAM with Backward Stepwise Edge Refinement\n# Novel: combines non-Gaussianity-based DAG discovery (DirectLiNGAM)\n# with partial F-test backward elimination. Edges must pass both\n# LiNGAM's independence criterion AND statistical significance for\n# variance explanation, providing complementary validation layers.\n# =====================================================================\ndef run_causal_discovery(X: np.ndarray) -> np.ndarray:\n    \"\"\"\n    Stage 1: DirectLiNGAM for non-Gaussianity-based ordering and initial DAG.\n    Stage 2: Backward stepwise edge refinement using partial F-tests — each\n    parent must significantly reduce residual variance to be retained.\n    B[i, j] != 0 means j -> i.\n    \"\"\"\n    import warnings; warnings.filterwarnings('ignore')\n    from sklearn.utils import check_array\n    from scipy import stats\n    from causallearn.search.FCMBased.lingam import DirectLiNGAM\n\n    X = check_array(X)\n    n_samples, n_features = X.shape\n\n    # Stage 1: DirectLiNGAM (ordering + non-Gaussianity-based edge estimation)\n    model = DirectLiNGAM()\n    model.fit(X)\n    B = model.adjacency_matrix_.copy()\n\n    # Stage 2: Backward stepwise refinement via partial F-tests\n    # For each node with multiple parents, iteratively remove the least\n    # significant parent until all remaining parents are significant.\n    alpha = 0.01  # Conservative significance level\n\n    for i in range(n_features):\n        parents = list(np.where(B[i, :] != 0)[0])\n        if len(parents) <= 1:\n            continue\n\n        y = X[:, i]\n        current = list(parents)\n\n        changed = True\n        while changed and len(current) > 1:\n            changed = False\n            Xc = X[:, current]\n            beta = np.linalg.lstsq(Xc, y, rcond=None)[0]\n            rss_full = np.sum((y - Xc @ beta) ** 2)\n            df_full = max(n_samples - len(current), 1)\n            mse_full = rss_full / df_full\n\n            # Find the least significant parent via partial F-test\n            worst_pval = -1.0\n            worst_idx = -1\n            for j in range(len(current)):\n                mask = list(range(len(current)))\n                mask.remove(j)\n                Xr = Xc[:, mask]\n                beta_r = np.linalg.lstsq(Xr, y, rcond=None)[0]\n                rss_red = np.sum((y - Xr @ beta_r) ** 2)\n\n                f_stat = max((rss_red - rss_full), 0.0) / (mse_full + 1e-15)\n                pval = 1.0 - stats.f.cdf(f_stat, 1, df_full)\n\n                if pval > worst_pval:\n                    worst_pval = pval\n                    worst_idx = j\n\n            if worst_pval > alpha and worst_idx >= 0:\n                del current[worst_idx]\n                changed = True\n\n        # Zero out edges that were removed\n        removed = set(parents) - set(current)\n        for r in removed:\n            B[i, r] = 0.0\n\n    return B\n# =====================================================================\n"
7}

1{
2  "n": 2
3}