S07: From Atom-Mapped Reactions to DPO Rules

import logging
import os
import warnings

# SynKit reports progress with tqdm. Outside a live Jupyter session tqdm falls
# back to a text bar that writes one stderr line per update, which buries the
# actual results. Bars stay on when you run this notebook yourself.
os.environ.setdefault("TQDM_DISABLE", "1")
logging.disable(logging.INFO)
warnings.filterwarnings("ignore")

import rdkit
import pandas as pd
import networkx as nx
from pathlib import Path
from tqdm.auto import tqdm
from synedu.Utils import load_database

print("RDKit version:", rdkit.__version__)
print("NetworkX version:", nx.__version__)

DATA_DIR = Path("data")
DATA_PATH = DATA_DIR / "data_aam.json.gz"
data = pd.DataFrame(
    load_database(DATA_PATH)[:100]
)  # 100 reactions for manageable runtime
display(data.head())
print(data.shape)

from tqdm.auto import tqdm
from synkit.Chem.Reaction.standardize import Standardize
from synkit.Chem.Reaction.canon_rsmi import CanonRSMI

tqdm.pandas(desc="Canonicalizing")


def std_canon(aam: str):
    if aam is None:
        return None
    try:
        std = Standardize()
        canon = CanonRSMI(wl_iterations=4)
        aam_std = std.fit(aam, remove_aam=False)
        canon_out = canon.canonicalise(aam_std)
        return getattr(canon_out, "canonical_rsmi", None)
    except Exception:
        return None


tqdm.pandas(desc="Canonicalizing — rxn_mapper")
data['rxn_canon'] = data['rxn_mapper'].progress_apply(std_canon)

tqdm.pandas(desc="Canonicalizing — graphormer")
data['graph_canon'] = data['graphormer'].progress_apply(std_canon)

tqdm.pandas(desc="Canonicalizing — local_mapper")
data['local_canon'] = data['local_mapper'].progress_apply(std_canon)

import pandas as pd
from typing import Tuple, List, Union, Sequence


def extract_aam(
    df: pd.DataFrame,
    keys: Sequence[str] = ("rxn_mapper", "graphormer", "local_mapper"),
    strip: bool = True,
    require_not_none: bool = True,
    return_indices: bool = False,
    inplace: bool = False,
) -> Union[
    Tuple[pd.DataFrame, pd.DataFrame],
    Tuple[pd.DataFrame, List[int], pd.DataFrame, List[int]],
]:
    k0, k1, k2 = keys
    s0, s1, s2 = df[k0].copy(), df[k1].copy(), df[k2].copy()
    if strip:
        mask0 = s0.map(lambda x: isinstance(x, str))
        mask1 = s1.map(lambda x: isinstance(x, str))
        mask2 = s2.map(lambda x: isinstance(x, str))
        if mask0.any():
            s0.loc[mask0] = s0.loc[mask0].str.strip()
        if mask1.any():
            s1.loc[mask1] = s1.loc[mask1].str.strip()
        if mask2.any():
            s2.loc[mask2] = s2.loc[mask2].str.strip()
    if require_not_none:
        non_null = s0.notna() & s1.notna() & s2.notna()
    else:
        non_null = pd.Series(True, index=df.index)
    eq_mask = non_null & s0.eq(s1) & s1.eq(s2)
    out_df = df.loc[eq_mask].copy()
    non_out_df = df.loc[~eq_mask].copy()
    if not out_df.empty:
        out_df = out_df.assign(aam=s0.loc[eq_mask].values)
        out_df = out_df.drop(columns=[k0, k1, k2], errors="ignore")
    if inplace:
        df.drop(df.index, inplace=True)
        for col in out_df.columns:
            df[col] = out_df[col].values
        df.reset_index(drop=True, inplace=True)
    if return_indices:
        return out_df, out_df.index.to_list(), non_out_df, non_out_df.index.to_list()
    return out_df, non_out_df


out, rest = extract_aam(data, keys=("rxn_canon", "graph_canon", "local_canon"))
print("Equivalent maps:", len(out))
print("Non equivalent maps:", len(rest))

import pandas as pd
from typing import Sequence, Tuple, Optional
from synkit.Chem.Reaction.aam_validator import AAMValidator


def refinement(
    df: pd.DataFrame,
    aam_cols: Sequence[str] = ("rxn_canon", "graph_canon", "local_canon"),
    validator: Optional[AAMValidator] = None,
) -> Tuple[pd.DataFrame, pd.DataFrame]:
    if validator is None:
        validator = AAMValidator()
    cols = list(aam_cols)
    if len(cols) < 2:
        raise ValueError("aam_cols must contain at least two columns")
    missing = [c for c in cols if c not in df.columns]
    if missing:
        raise KeyError(f"missing columns: {missing}")
    df2 = df.copy()
    df2[cols] = df2[cols].map(lambda x: x.strip() if isinstance(x, str) else x)
    ref = cols[0]
    solved_idx, unsolved_idx = [], []
    for idx, row in df2[cols].iterrows():
        ref_val = row[ref]
        if ref_val is None or pd.isna(ref_val):
            unsolved_idx.append(idx)
            continue
        ok = True
        for c in cols[1:]:
            val = row[c]
            if val is None or pd.isna(val):
                ok = False
                break
            try:
                if not bool(validator.smiles_check(ref_val, val)):
                    ok = False
                    break
            except Exception:
                ok = False
                break
        (solved_idx if ok else unsolved_idx).append(idx)
    solved = df2.loc[solved_idx].copy()
    solved = solved.drop(columns=[c for c in cols if c != ref], errors="ignore")
    solved = solved.rename(columns={ref: "aam"})
    unsolved = df2.loc[unsolved_idx].copy()
    return solved, unsolved


solve, unsolve = refinement(rest)

df_aam = pd.concat([out, solve], axis=0)
df_aam.shape

df_aam.drop_duplicates(subset='aam', inplace=True)

from synedu.Utils import rsmi_to_its
from synedu.Utils.its_vis import visualize_its
import matplotlib.pyplot as plt


def _safe_rsmi_to_its(rsmi, core=False):
    """Wrap rsmi_to_its to silently return None on bad/unmapped SMILES."""
    try:
        if not isinstance(rsmi, str) or ">>" not in rsmi:
            return None
        return rsmi_to_its(rsmi, core=core)
    except Exception:
        return None


tqdm.pandas(desc="ITS construction")
df_aam["ITS"] = df_aam["aam"].progress_apply(lambda s: _safe_rsmi_to_its(s, core=False))

tqdm.pandas(desc="Reaction center construction")
df_aam["RC"] = df_aam["aam"].progress_apply(lambda s: _safe_rsmi_to_its(s, core=True))

# Drop rows where ITS/RC could not be built
_before = len(df_aam)
df_aam.dropna(subset=["ITS", "RC"], inplace=True)
df_aam.reset_index(drop=True, inplace=True)
print(f"Kept {len(df_aam)}/{_before} reactions with valid ITS/RC")

fig, ax = plt.subplots(1, 2, figsize=(20, 8))
_idx = df_aam.index[1]
visualize_its(df_aam["ITS"][_idx], show_node_labels=False, ax=ax[0], show_legend=False)
visualize_its(df_aam["RC"][_idx], show_node_labels=False, ax=ax[1], show_legend=False)
ax[0].set_title("ITS graph", fontsize=24, fontweight="bold")
ax[1].set_title("Reaction center", fontsize=24, fontweight="bold")
plt.tight_layout()
plt.show()

from synedu.Utils import cluster_its_graphs

_rc_graphs = df_aam["RC"].tolist()
_labels = cluster_its_graphs(
    _rc_graphs,
    node_attrs=["element", "formal_charge"],
    edge_attrs=["order"],
)
df_aam["class"] = _labels
result = df_aam.to_dict("records")
print(f"Total reactions: {len(result)},  unique rule classes: {len(set(_labels))}")
result[0]

from synedu.Utils import stratified_sample

rc_records = stratified_sample(result, key="class")  # one record per class
rc = [r["RC"] for r in rc_records]
print(f"Representative reaction centers: {len(rc)}")

import networkx as nx
import matplotlib.pyplot as plt
from synedu.Utils.reaction import build_its
from synedu.Utils.its_vis import visualize_its


def its_decompose(its: nx.Graph):
    """Decompose ITS back into reactant and product molecular graphs."""
    r, p = nx.Graph(), nx.Graph()
    for n, data in its.nodes(data=True):
        r_attrs, p_attrs = {}, {}
        for k, v in data.items():
            if isinstance(v, tuple) and len(v) == 2:
                r_attrs[k], p_attrs[k] = v[0], v[1]
            else:
                r_attrs[k] = p_attrs[k] = v
        r_attrs["atom_map"] = p_attrs["atom_map"] = n
        r.add_node(n, **r_attrs)
        p.add_node(n, **p_attrs)
    for u, v, data in its.edges(data=True):
        br, bp = data.get("order", (0.0, 0.0))
        if float(br) > 0:
            r.add_edge(u, v, order=float(br))
        if float(bp) > 0:
            p.add_edge(u, v, order=float(bp))
    return r, p


def h_expand_its(its: nx.Graph) -> nx.Graph:
    """
    Expand implicit H transfers in an ITS into explicit H nodes.

    Only *balanced* donor-acceptor pairs are made explicit.
    Explicit H atoms already present in the ITS (from mapped SMILES) are
    skipped — their bonds are already encoded as ITS edges, so treating them
    as donors/acceptors would create spurious H-H bonds.
    """
    r, p = its_decompose(its)

    donors, acceptors = [], []
    for n in sorted(set(r.nodes()) & set(p.nodes())):
        if r.nodes[n].get("element") == "H":
            continue  # already explicit — skip to avoid H-H bonds
        hr = int(r.nodes[n].get("hcount", 0) or 0)
        hp = int(p.nodes[n].get("hcount", 0) or 0)
        # GetTotalNumHs() counts explicit H neighbours — subtract them
        # so we only track *implicit* transfers (explicit H nodes are
        # already encoded as ITS edges and must not be double-counted).
        expl_r = sum(1 for m in r.neighbors(n) if r.nodes[m].get("element") == "H")
        expl_p = sum(1 for m in p.neighbors(n) if p.nodes[m].get("element") == "H")
        implicit_r = max(0, hr - expl_r)
        implicit_p = max(0, hp - expl_p)
        donors.extend([n] * max(0, implicit_r - implicit_p))
        acceptors.extend([n] * max(0, implicit_p - implicit_r))

    n_shared = min(len(donors), len(acceptors))
    if n_shared == 0:
        return its  # no balanced transfers to expand

    next_id = max(set(r.nodes()) | set(p.nodes()), default=0) + 1
    for i in range(n_shared):
        hid = next_id + i
        _h_attrs = dict(
            element="H", formal_charge=0, aromatic=False, hcount=0, atom_map=hid
        )
        r.add_node(hid, **_h_attrs)
        r.add_edge(donors[i], hid, order=1.0)
        r.nodes[donors[i]]["hcount"] = max(
            0, (r.nodes[donors[i]].get("hcount", 0) or 0) - 1
        )
        p.add_node(hid, **_h_attrs)
        p.add_edge(acceptors[i], hid, order=1.0)
        p.nodes[acceptors[i]]["hcount"] = max(
            0, (p.nodes[acceptors[i]].get("hcount", 0) or 0) - 1
        )

    return build_its(r, p)


rc_h = [h_expand_its(g) for g in rc]
print(f"H-expanded reaction centers: {len(rc_h)}")

# ── Show comparison: first RC with a shared H transfer ───────────────
_show_idx = next(
    (i for i, g in enumerate(rc_h) if g.number_of_nodes() > rc[i].number_of_nodes()), 0
)
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 5), facecolor="white")
visualize_its(
    rc[_show_idx],
    ax=ax1,
    title="RC — implicit H",
    show_legend=False,
    layout="kamada_kawai",
)
visualize_its(
    rc_h[_show_idx],
    ax=ax2,
    title="RC — explicit H transfer",
    show_legend=True,
    layout="kamada_kawai",
)
plt.suptitle(
    "Reaction center before / after hydrogen expansion", fontsize=11, fontweight="bold"
)
plt.tight_layout()
plt.show()

_with_h = sum(1 for a, b in zip(rc, rc_h) if b.number_of_nodes() > a.number_of_nodes())
print(f"{_with_h}/{len(rc)} rule classes have explicit H transfer")

from itertools import permutations as _perms
import networkx as nx
import matplotlib.pyplot as plt
from synedu.Utils import rsmi_to_its, wl_hash
from synedu.Utils.reaction import build_its
from synedu.Utils.its_vis import visualize_its


def h_expand_all(its: nx.Graph):
    """
    Return all structurally distinct H-expanded ITS graphs.

    Tries every permutation of donor-acceptor pairings;
    uses the WL hash only as a bucket key, then de-duplicates
    each bucket by exact labeled-graph isomorphism.
    """
    r, p = its_decompose(its)
    donors, acceptors = [], []
    for n in sorted(set(r.nodes()) & set(p.nodes())):
        if r.nodes[n].get("element") == "H":
            continue
        hr = int(r.nodes[n].get("hcount", 0) or 0)
        hp = int(p.nodes[n].get("hcount", 0) or 0)
        expl_r = sum(1 for m in r.neighbors(n) if r.nodes[m].get("element") == "H")
        expl_p = sum(1 for m in p.neighbors(n) if p.nodes[m].get("element") == "H")
        implicit_r = max(0, hr - expl_r)
        implicit_p = max(0, hp - expl_p)
        donors.extend([n] * max(0, implicit_r - implicit_p))
        acceptors.extend([n] * max(0, implicit_p - implicit_r))

    if not donors or not acceptors:
        return [its]

    n_shared = min(len(donors), len(acceptors))
    donors = donors[:n_shared]
    acceptors = acceptors[:n_shared]
    results, seen = [], {}

    for acc_perm in _perms(acceptors):
        r_c, p_c = r.copy(), p.copy()
        nxt = max(set(r.nodes()) | set(p.nodes()), default=0) + 1
        for i, (d, a) in enumerate(zip(donors, acc_perm)):
            hid = nxt + i
            _h = dict(
                element="H", formal_charge=0, aromatic=False, hcount=0, atom_map=hid
            )
            r_c.add_node(hid, **_h)
            r_c.add_edge(d, hid, order=1.0)
            r_c.nodes[d]["hcount"] = max(0, (r_c.nodes[d].get("hcount", 0) or 0) - 1)
            p_c.add_node(hid, **_h)
            p_c.add_edge(a, hid, order=1.0)
            p_c.nodes[a]["hcount"] = max(0, (p_c.nodes[a].get("hcount", 0) or 0) - 1)
        cand = build_its(r_c, p_c)
        h = wl_hash(cand, node_attrs=["element", "formal_charge"], edge_attrs=["order"])
        bucket = seen.setdefault(h, [])
        node_match = lambda a, b: (
            a.get("element") == b.get("element")
            and a.get("formal_charge") == b.get("formal_charge")
        )
        edge_match = lambda a, b: a.get("order") == b.get("order")
        if any(
            nx.is_isomorphic(
                cand,
                previous,
                node_match=node_match,
                edge_match=edge_match,
            )
            for previous in bucket
        ):
            continue
        bucket.append(cand)
        results.append(cand)
    return results


# ── Real reaction: trichloroacetonitrile + N-hydroxy compound ────────────
_rsmi_amb = (
    "[CH3:1][O:2][C:3](=[O:4])[N:5]([OH:6])[c:7]1[cH:8][cH:9][cH:10][cH:11][cH:12]1"
    ".[N:18]#[C:17][C:14]([Cl:13])([Cl:15])[Cl:16]"
    ">>"
    "[CH3:1][O:2][C:3](=[O:4])[NH:5][c:7]1[cH:8][cH:9][cH:10][cH:11][c:12]1"
    "[NH:18][C:17](=[O:6])[C:14]([Cl:13])([Cl:15])[Cl:16]"
)

_rc_amb = rsmi_to_its(_rsmi_amb, core=True)

print("RC summary:")
for n in sorted(_rc_amb.nodes()):
    d = _rc_amb.nodes[n]
    hc = d.get("hcount", (0, 0))
    role = (
        "donor"
        if isinstance(hc, tuple) and hc[0] > hc[1]
        else "acceptor" if isinstance(hc, tuple) and hc[1] > hc[0] else ""
    )
    if role:
        print(f"  atom {n:2d} ({d.get('element')})  hcount={hc}  → {role}")

candidates = h_expand_all(_rc_amb)
print(f"\n→ {len(candidates)} structurally distinct expansion(s)")

_wl_fn = lambda g: wl_hash(
    g, node_attrs=["element", "formal_charge"], edge_attrs=["order"]
)
for i, cand in enumerate(candidates):
    print(
        f"  Pairing {i+1}: {cand.number_of_nodes()} nodes, "
        f"{cand.number_of_edges()} edges  WL={_wl_fn(cand)[:12]}…"
    )

# ── Visualize: RC + both expansions ─────────────────────────────────────
_n = len(candidates)
_pair_titles = [
    "Pairing A\nO6→N5,  C12→N18\nfirst correspondence topology",
    "Pairing B\nO6→N18,  C12→N5\nalternative correspondence topology",
]

fig, axes = plt.subplots(1, _n + 1, figsize=(6.5 * (_n + 1), 6), facecolor="white")

visualize_its(
    _rc_amb,
    ax=axes[0],
    title="Ambiguous RC (implicit H)\n2 donors · 2 acceptors",
    show_legend=False,
    show_edge_labels=True,
    layout="kamada_kawai",
)
axes[0].text(
    0.5,
    -0.05,
    "Donors: O6 (N–OH) · C12 (Ar–H)\nAcceptors: N5 (amine) · N18 (nitrile N)",
    ha="center",
    va="top",
    transform=axes[0].transAxes,
    fontsize=8,
    color="#333",
    bbox=dict(boxstyle="round,pad=0.3", fc="lightyellow", alpha=0.95),
)

for i, cand in enumerate(candidates):
    ttl = _pair_titles[i] if i < len(_pair_titles) else f"Pairing {i+1}"
    visualize_its(
        cand,
        ax=axes[i + 1],
        title=ttl,
        show_legend=(i == _n - 1),
        show_edge_labels=True,
        layout="kamada_kawai",
    )
    axes[i + 1].text(
        0.5,
        -0.05,
        f"WL: {_wl_fn(cand)[:14]}…",
        ha="center",
        va="top",
        transform=axes[i + 1].transAxes,
        fontsize=7,
        color="#666",
    )

plt.suptitle(
    "Ambiguous H expansion — 2 non-isomorphic reaction centers\n"
    "(trichloroacetonitrile + N-hydroxy compound)",
    fontsize=11,
    fontweight="bold",
)
plt.tight_layout()
plt.show()

from synedu.Utils import wl_hash

tqdm.pandas(desc="WL hash calculation")
_wl = lambda g: wl_hash(
    g, node_attrs=["element", "formal_charge"], edge_attrs=["order"]
)
df_aam["wl_hash"] = df_aam["RC"].progress_apply(_wl)
df_aam["wl_hash_its"] = df_aam["ITS"].progress_apply(_wl)

df_aam['wl_hash'].value_counts()

import networkx as nx
import matplotlib.pyplot as plt
from synedu.Utils.its_vis import visualize_its

# ── Minimal demo RC: C-O broken, C-N formed, C-C preserved, + H transfer ──
_demo = nx.Graph()
_demo.add_node(
    1, element="C", formal_charge=(0, 0), hcount=(0, 0), aromatic=(False, False)
)
_demo.add_node(
    2, element="O", formal_charge=(0, 0), hcount=(1, 0), aromatic=(False, False)
)
_demo.add_node(
    3, element="N", formal_charge=(0, 0), hcount=(0, 1), aromatic=(False, False)
)
_demo.add_node(
    4, element="C", formal_charge=(0, 0), hcount=(0, 0), aromatic=(False, False)
)
_demo.add_node(
    5, element="H", formal_charge=(0, 0), hcount=(0, 0), aromatic=(False, False)
)

_demo.add_edge(1, 2, order=(1.0, 0.0))  # C-O  broken    → left
_demo.add_edge(1, 3, order=(0.0, 1.0))  # C-N  formed    → right
_demo.add_edge(1, 4, order=(1.0, 1.0))  # C-C  preserved → context
_demo.add_edge(2, 5, order=(1.0, 0.0))  # O-H  broken    → left
_demo.add_edge(3, 5, order=(0.0, 1.0))  # N-H  formed    → right

# ── Edge classification ────────────────────────────────────────────────────
_BOND_LBL = {1.0: "1", 1.5: "#", 2.0: "2", 3.0: "3"}

print(f"{'Edge':>7}  {'order':>12}  GML block")
print("-" * 40)
for u, v, d in _demo.edges(data=True):
    br, bp = float(d["order"][0]), float(d["order"][1])
    eu, ev = _demo.nodes[u]["element"], _demo.nodes[v]["element"]
    if br == bp and br > 0:
        block = "context"
    elif br > 0 and bp == 0:
        block = "left"
    elif br == 0 and bp > 0:
        block = "right"
    else:
        block = "left + right"  # bond order changes (e.g. single→double)
    print(f"{eu:>2}-{ev:<2}  {str(d['order']):>14}  {block}")

# ── Visualise ─────────────────────────────────────────────────────────────
fig, ax = plt.subplots(figsize=(5, 4), facecolor="white")
visualize_its(
    _demo,
    ax=ax,
    title="Demo RC — edge classification",
    show_legend=True,
    show_edge_labels=True,
    layout="kamada_kawai",
)
plt.tight_layout()
plt.show()

def _bond_label(order: float) -> str:
    """Map a scalar bond order to its GML label string."""
    return {1.0: "1", 1.5: "#", 2.0: "2", 3.0: "3"}.get(
        round(float(order), 2), str(order)
    )


def _node_label(attrs: dict) -> str:
    """Return element symbol with optional charge suffix (e.g. N+, O-)."""
    elem = attrs.get("element", "*")
    chg_raw = attrs.get("formal_charge", 0)
    # formal_charge can be a scalar or an ITS (r, p) tuple — take reactant side
    chg = int((chg_raw[0] if isinstance(chg_raw, tuple) else chg_raw) or 0)
    if chg > 0:
        return f"{elem}+"
    if chg < 0:
        return f"{elem}-"
    return str(elem)


# ── Quick tests ──────────────────────────────────────────────────────────
print("_bond_label tests:")
for b, expected in [(1.0, "1"), (2.0, "2"), (1.5, "#"), (3.0, "3")]:
    out = _bond_label(b)
    print(f"  _bond_label({b}) = {out!r}  {'✓' if out == expected else '✗'}")

print("\n_node_label tests:")
cases = [
    ({"element": "C", "formal_charge": 0}, "C"),
    ({"element": "N", "formal_charge": 1}, "N+"),
    ({"element": "O", "formal_charge": -1}, "O-"),
    ({"element": "N", "formal_charge": (1, 1)}, "N+"),  # ITS tuple form
]
for attrs, expected in cases:
    out = _node_label(attrs)
    print(f"  _node_label({attrs}) = {out!r}  {'✓' if out == expected else '✗'}")

import networkx as nx
from typing import Union


def its_to_gml(
    its: nx.Graph,
    rule_id: Union[int, str] = 0,
    include_context_edges: bool = True,
) -> str:
    """Convert an ITS / H-expanded RC graph to a MØD GML rule string.

    :param its: ITS graph (edge order = (b_r, b_p) tuples).
    :param rule_id: Identifier written to the ``ruleID`` field.
    :param include_context_edges: If True, preserved bonds are emitted in
        the ``context`` block (required by MØD for bond-order context).
    :returns: GML rule string.
    """
    left_edges, right_edges, ctx_edges = [], [], []

    for u, v, d in its.edges(data=True):
        order = d.get("order", (0.0, 0.0))
        br, bp = (
            (float(order[0]), float(order[1]))
            if isinstance(order, (tuple, list))
            else (float(order), float(order))
        )

        if abs(br - bp) < 1e-9:  # preserved bond
            if br > 0 and include_context_edges:
                ctx_edges.append((u, v, _bond_label(br)))
        else:
            if br > 0:
                left_edges.append((u, v, _bond_label(br)))  # broken
            if bp > 0:
                right_edges.append((u, v, _bond_label(bp)))  # formed

    lines = [f'rule [', f'  ruleID "{rule_id}"', '  left [']
    for u, v, lbl in left_edges:
        lines.append(f'    edge [ source {u} target {v} label "{lbl}" ]')
    lines.append('  ]')

    lines.append('  context [')
    for n in sorted(its.nodes()):
        lines.append(f'    node [ id {n} label "{_node_label(its.nodes[n])}" ]')
    for u, v, lbl in ctx_edges:
        lines.append(f'    edge [ source {u} target {v} label "{lbl}" ]')
    lines.append('  ]')

    lines.append('  right [')
    for u, v, lbl in right_edges:
        lines.append(f'    edge [ source {u} target {v} label "{lbl}" ]')
    lines.append('  ]')
    lines.append(']')
    return '\n'.join(lines)


# ── Apply to the demo RC and inspect the output ──────────────────────────
gml_demo = its_to_gml(_demo, rule_id="demo", include_context_edges=True)
print(gml_demo)

# Generate one GML rule string per unique H-expanded reaction center
rules = [
    its_to_gml(g, rule_id=i, include_context_edges=True) for i, g in enumerate(rc_h)
]
print(f"Generated {len(rules)} GML rules")
print("\nFirst rule:")
print(rules[0] if rules else "(empty)")

# Print a GML rule with H transfer (if one exists), else print the first rule
_ex = next((r for r in rules if "H" in r), rules[0] if rules else "(no rules)")
print(_ex)