1{
2  "op": "replace",
3  "filename": "causal-learn/bench/custom_algorithm.py",
4  "start_line": 3,
5  "end_line": 14,
6  "content": "from causallearn.graph.GraphNode import GraphNode\nfrom causallearn.graph.Endpoint import Endpoint\nfrom causallearn.graph.Edge import Edge\n\ndef run_causal_discovery(X: np.ndarray) -> GeneralGraph:\n    \"\"\"\n    Confidence-Calibrated Multi-Method Fusion (CCMF) for CPDAG Discovery.\n    \n    Novel algorithm that fuses constraint-based (PC at multiple significance\n    levels) and score-based (GES) methods using adaptive confidence weighting.\n    Multi-resolution PC provides edge stability estimates; GES brings score-based\n    optimality. Weighted consensus reduces both FP and FN vs any single method.\n    \"\"\"\n    n, p = X.shape\n    from causallearn.search.ConstraintBased.PC import pc\n    from causallearn.search.ScoreBased.GES import ges\n\n    # --- Phase 1: Multi-resolution estimation ---\n    graphs, wts = [], []\n\n    # Constraint-based at multiple significance levels for stability\n    for alpha in [0.01, 0.05, 0.1]:\n        try:\n            cg = pc(X, alpha=alpha, indep_test='fisherz', stable=True,\n                    show_progress=False)\n        except TypeError:\n            cg = pc(X, alpha=alpha, indep_test='fisherz', stable=True)\n        graphs.append(cg.G)\n        wts.append(1.0)\n\n    # Score-based (GES with BIC) - higher weight for Gaussian data\n    record = ges(X, score_func='local_score_BIC')\n    graphs.append(record['G'])\n    wts.append(2.5)\n\n    total_w = sum(wts)\n\n    # --- Phase 2: Confidence aggregation ---\n    adj_c = np.zeros((p, p))    # adjacency confidence\n    arr_c = np.zeros((p, p))    # arr_c[i,j] = weighted evidence for i->j\n\n    for g, w in zip(graphs, wts):\n        gn = g.get_nodes()\n        for i in range(p):\n            for j in range(i + 1, p):\n                e = g.get_edge(gn[i], gn[j])\n                if e is not None:\n                    adj_c[i, j] += w\n                    ep_i = e.get_proximal_endpoint(gn[i])\n                    ep_j = e.get_proximal_endpoint(gn[j])\n                    if ep_j == Endpoint.ARROW and ep_i == Endpoint.TAIL:\n                        arr_c[i, j] += w\n                    elif ep_i == Endpoint.ARROW and ep_j == Endpoint.TAIL:\n                        arr_c[j, i] += w\n\n    # --- Phase 3: Consensus CPDAG construction ---\n    nodes = [GraphNode(f\"X{i + 1}\") for i in range(p)]\n    G = GeneralGraph(nodes)\n\n    a_t = total_w * 0.45   # adjacency threshold (slightly below majority)\n    o_t = total_w * 0.25   # orientation threshold\n\n    for i in range(p):\n        for j in range(i + 1, p):\n            if adj_c[i, j] >= a_t:\n                fwd, bwd = arr_c[i, j], arr_c[j, i]\n                if fwd > bwd and fwd >= o_t:\n                    G.add_edge(Edge(nodes[i], nodes[j],\n                                    Endpoint.TAIL, Endpoint.ARROW))\n                elif bwd > fwd and bwd >= o_t:\n                    G.add_edge(Edge(nodes[i], nodes[j],\n                                    Endpoint.ARROW, Endpoint.TAIL))\n                else:\n                    G.add_edge(Edge(nodes[i], nodes[j],\n                                    Endpoint.TAIL, Endpoint.TAIL))\n\n    return G\n"
7}

1{
2  "op": "replace",
3  "filename": "causal-learn/bench/custom_algorithm.py",
4  "start_line": 3,
5  "end_line": 14,
6  "content": "from causallearn.graph.GraphNode import GraphNode\nfrom causallearn.graph.Endpoint import Endpoint\nfrom causallearn.graph.Edge import Edge\nfrom causallearn.utils.DAG2CPDAG import dag2cpdag\n\ndef run_causal_discovery(X: np.ndarray) -> GeneralGraph:\n    \"\"\"\n    Equal-Variance Iterative Sink Identification (EVIS).\n    Exploits equal error variance for full DAG identifiability in linear Gaussian SEMs.\n    Under equal variance, precision matrix diagonal ω_ii = σ^{-2}(1 + ||outgoing weights||^2),\n    so sink nodes (no children) minimize ω_ii. We iteratively peel sinks to recover the\n    causal ordering, then select parents via BIC, and convert the DAG to a CPDAG.\n    \"\"\"\n    n, p = X.shape\n    X_c = X - X.mean(axis=0)\n    S = X_c.T @ X_c / n\n    # Step 1: Causal ordering via iterative sink identification\n    remaining = list(range(p))\n    rev_order = []\n    for _ in range(p):\n        idx = np.array(remaining)\n        m = len(idx)\n        S_sub = S[np.ix_(idx, idx)]\n        Omega = np.linalg.inv(S_sub + 1e-8 * np.eye(m))\n        sink = int(np.argmin(np.diag(Omega)))\n        rev_order.append(remaining.pop(sink))\n    order = rev_order[::-1]\n    # Step 2: BIC forward-backward parent selection\n    logn = np.log(n)\n    parents = {}\n    for k in range(p):\n        node = order[k]\n        y = X_c[:, node]\n        cands = order[:k]\n        sel = []\n        best_bic = n * np.log(np.dot(y, y) / n + 1e-30)\n        # Forward\n        fwd = True\n        while fwd:\n            fwd = False\n            best_c, best_b = None, best_bic\n            for c in cands:\n                if c in sel:\n                    continue\n                trial = sel + [c]\n                Xp = X_c[:, trial]\n                coef = np.linalg.lstsq(Xp, y, rcond=None)[0]\n                resid = y - Xp @ coef\n                rv = np.dot(resid, resid) / n\n                bv = n * np.log(rv + 1e-30) + len(trial) * logn\n                if bv < best_b:\n                    best_b, best_c = bv, c\n            if best_c is not None:\n                sel.append(best_c)\n                best_bic = best_b\n                fwd = True\n        # Backward\n        bwd = True\n        while bwd:\n            bwd = False\n            for c in list(sel):\n                trial = [x for x in sel if x != c]\n                if trial:\n                    Xp = X_c[:, trial]\n                    coef = np.linalg.lstsq(Xp, y, rcond=None)[0]\n                    rv = np.dot(y - Xp @ coef, y - Xp @ coef) / n\n                else:\n                    rv = np.dot(y, y) / n\n                bv = n * np.log(rv + 1e-30) + len(trial) * logn\n                if bv < best_bic:\n                    sel.remove(c)\n                    best_bic = bv\n                    bwd = True\n                    break\n        parents[node] = sel\n    # Step 3: Build DAG and convert to CPDAG\n    nodes_list = [GraphNode(f\"X{i+1}\") for i in range(p)]\n    G = GeneralGraph(nodes_list)\n    for child, pa_list in parents.items():\n        for pa in pa_list:\n            G.add_edge(Edge(nodes_list[pa], nodes_list[child], Endpoint.TAIL, Endpoint.ARROW))\n    return dag2cpdag(G)\n"
7}

1{
2  "op": "replace",
3  "filename": "causal-learn/bench/custom_algorithm.py",
4  "start_line": 3,
5  "end_line": 84,
6  "content": "from causallearn.graph.GraphNode import GraphNode\nfrom causallearn.graph.Endpoint import Endpoint\nfrom causallearn.graph.Edge import Edge\nfrom causallearn.utils.DAG2CPDAG import dag2cpdag\nimport math\n\ndef _sel_pa(Xc, node, cands, logn, gamma):\n    \"\"\"EBIC-based forward-backward parent selection.\n    Uses Extended BIC (Chen & Chen 2008): BIC + 2*gamma*log(C(q,k))\n    where q=|candidates|, k=|selected|. This adaptively penalizes more\n    when selecting from larger candidate pools, reducing false positives.\"\"\"\n    n = Xc.shape[0]\n    y = Xc[:, node]\n    q = len(cands)\n    sel = []\n    def _sc(rv, k):\n        s = n * np.log(rv + 1e-30) + k * logn\n        if gamma > 0 and 0 < k <= q:\n            s += 2 * gamma * (math.lgamma(q + 1) - math.lgamma(k + 1) - math.lgamma(q - k + 1))\n        return s\n    best = _sc(np.dot(y, y) / n, 0)\n    # Forward selection\n    while True:\n        bc, bb = None, best\n        for c in cands:\n            if c in sel:\n                continue\n            t = sel + [c]\n            Xp = Xc[:, t]\n            cf = np.linalg.lstsq(Xp, y, rcond=None)[0]\n            rv = np.dot(y - Xp @ cf, y - Xp @ cf) / n\n            bv = _sc(rv, len(t))\n            if bv < bb:\n                bb, bc = bv, c\n        if bc is not None:\n            sel.append(bc)\n            best = bb\n        else:\n            break\n    # Backward elimination\n    while True:\n        rm = False\n        for c in list(sel):\n            t = [x for x in sel if x != c]\n            if t:\n                Xp = Xc[:, t]\n                cf = np.linalg.lstsq(Xp, y, rcond=None)[0]\n                rv = np.dot(y - Xp @ cf, y - Xp @ cf) / n\n            else:\n                rv = np.dot(y, y) / n\n            bv = _sc(rv, len(t))\n            if bv < best:\n                sel.remove(c)\n                best = bv\n                rm = True\n                break\n        if not rm:\n            break\n    return sel, best\n\ndef run_causal_discovery(X: np.ndarray) -> GeneralGraph:\n    \"\"\"\n    Equal-Variance Iterative Sink Identification with EBIC and Swap Refinement.\n    \n    Novel contributions:\n    1. Causal ordering via iterative precision-diagonal minimization exploiting\n       equal error variance: under equal Var(eps)=sigma^2, the precision diagonal\n       omega_ii = sigma^{-2}(1 + ||outgoing_weights_i||^2), so sinks minimize omega_ii.\n    2. EBIC parent selection adapts penalty to candidate pool size, reducing\n       false positives for later nodes in the ordering.\n    3. Adjacent-swap ordering refinement: greedily improves the ordering by\n       swapping adjacent pairs when total EBIC score decreases.\n    \"\"\"\n    n, p = X.shape\n    Xc = X - X.mean(axis=0)\n    S = Xc.T @ Xc / n\n    logn = np.log(n)\n    gamma = 0.5  # EBIC hyperparameter (standard choice from Chen & Chen 2008)\n    # Step 1: EVIS ordering via iterative sink identification\n    rem = list(range(p))\n    rord = []\n    for _ in range(p):\n        ix = np.array(rem)\n        m = len(ix)\n        Om = np.linalg.inv(S[np.ix_(ix, ix)] + 1e-8 * np.eye(m))\n        s = int(np.argmin(np.diag(Om)))\n        rord.append(rem.pop(s))\n    order = rord[::-1]  # source-first topological order\n    # Step 2: EBIC parent selection given ordering\n    pa, sc = {}, {}\n    for k in range(p):\n        pa[order[k]], sc[order[k]] = _sel_pa(Xc, order[k], order[:k], logn, gamma)\n    # Step 3: Ordering refinement via adjacent swaps (bubble-sort style)\n    for _ in range(10):\n        swapped = False\n        for k in range(p - 1):\n            a, b = order[k], order[k + 1]\n            old_sc = sc[a] + sc[b]\n            sel_b, sc_b = _sel_pa(Xc, b, order[:k], logn, gamma)\n            sel_a, sc_a = _sel_pa(Xc, a, order[:k] + [b], logn, gamma)\n            if sc_a + sc_b < old_sc - 1e-6:\n                order[k], order[k + 1] = b, a\n                pa[b], sc[b] = sel_b, sc_b\n                pa[a], sc[a] = sel_a, sc_a\n                swapped = True\n        if not swapped:\n            break\n    # Step 4: Build DAG and convert to CPDAG\n    nds = [GraphNode(f\"X{i + 1}\") for i in range(p)]\n    G = GeneralGraph(nds)\n    for ch in pa:\n        for pr in pa[ch]:\n            G.add_edge(Edge(nds[pr], nds[ch], Endpoint.TAIL, Endpoint.ARROW))\n    return dag2cpdag(G)\n"
7}

1{
2  "op": "replace",
3  "filename": "causal-learn/bench/custom_algorithm.py",
4  "start_line": 80,
5  "end_line": 80,
6  "content": "    gamma = 0.25  # Mild EBIC: balances false positive reduction (SF50) vs recall (Hard scenarios)"
7}